Loading code/1_NGRC_Lorenz_Basic_Model.ipynb +227 −60 Original line number Diff line number Diff line Loading @@ -92,7 +92,7 @@ "### 1.1 Discrete-Time Parameters\n", "\n", "Time is discretized using a fixed step size $dt$. \n", "Each continuous-time interval is converted into a number of discrete samples by dividing by $dt$.\n", "Each continuous time interval is converted into a number of discrete samples by dividing by $dt$.\n", "\n", "| Name | Symbol | Definition | Value |\n", "|------|--------|------------|--------|\n", Loading Loading @@ -293,17 +293,16 @@ }, { "cell_type": "markdown", "id": "3810a2f4", "id": "6657a8ad", "metadata": {}, "source": [ "### 2.1 Lorenz Time-Series Generation\n", "\n", "\n", "The Lorenz63 system is integrated over $[0, T_{\\text{total}}]$ \n", "and sampled at fixed time intervals $dt$.\n", "The Lorenz63 system is integrated over the time interval $[0, T_{\\text{total}}]$\n", "and sampled at uniform time steps $dt$.\n", "\n", "$$\n", "t_n = n\\,dt, \\quad n = 0, 1, \\dots, N_{\\text{total}}.\n", "t_n = n\\,dt, \\quad n = 0, 1, \\dots, N_{\\text{total}}\n", "$$\n", "\n", "**Initial condition**\n", Loading @@ -312,7 +311,7 @@ "(x(0), y(0), z(0)) = (17.6772,\\; 12.9314,\\; 43.9140)\n", "$$\n", "\n", "The discrete state vector is\n", "The discrete state vector is defined as\n", "\n", "$$\n", "\\mathbf{x}(t_n) =\n", Loading @@ -320,10 +319,11 @@ "x(t_n) \\\\\n", "y(t_n) \\\\\n", "z(t_n)\n", "\\end{bmatrix}.\n", "\\end{bmatrix}\n", "$$\n", "\n", "This produces the trajectories $x(t_n)$, $y(t_n)$, and $z(t_n)$.\n" "This produces the time series $x(t_n)$, $y(t_n)$, and $z(t_n)$.\n", "\n" ] }, { Loading Loading @@ -375,7 +375,7 @@ "+\n", "\\mathrm{Var}\\big(y(t_n)\\big)\n", "+\n", "\\mathrm{Var}\\big(z(t_n)\\big).\n", "\\mathrm{Var}\\big(z(t_n)\\big)\n", "$$\n", "\n", "This is used to normalize the prediction error.\n", Loading Loading @@ -422,7 +422,7 @@ "y(t_0) & \\cdots & y(t_{N_{\\text{total}}}) \\\\\n", "z(t_0) & \\cdots & z(t_{N_{\\text{total}}})\n", "\\end{bmatrix}\n", "\\in \\mathbb{R}^{3 \\times (N_{\\text{total}}+1)}.\n", "\\in \\mathbb{R}^{3 \\times (N_{\\text{total}}+1)}\n", "$$\n", "\n", "The matrix is partitioned sequentially along the time axis:\n", Loading @@ -436,7 +436,7 @@ "\\mathbf{X}_{\\text{train}}\n", "\\;\\big|\\;\n", "\\mathbf{X}_{\\text{test}}\n", "\\big],\n", "\\big]\n", "$$\n", "\n", "where\n", Loading @@ -446,7 +446,7 @@ "\\quad\n", "\\mathbf{X}_{\\text{train}} \\in \\mathbb{R}^{3 \\times N_{\\text{train}}},\n", "\\quad\n", "\\mathbf{X}_{\\text{test}} \\in \\mathbb{R}^{3 \\times N_{\\text{test}}}.\n", "\\mathbf{X}_{\\text{test}} \\in \\mathbb{R}^{3 \\times N_{\\text{test}}}\n", "$$\n" ] }, Loading @@ -459,18 +459,20 @@ "\n", "The NVAR feature vector at time $t$ is built as\n", "\n", "\n", "$$\n", "\\boldsymbol{\\phi}(\\mathbf{x}(t)) =\n", "\\mathbf{O}_{\\text{tot}}(t) =\n", "\\begin{bmatrix}\n", "\\text{Bias} \\\\\n", "\\text{Linear Delay Terms} \\\\\n", "\\text{Nonlinear Terms}\n", "\\end{bmatrix}.\n", "\\end{bmatrix}\n", "$$\n", "\n", "\n", "- The bias term is a constant equal to 1.\n", "- The linear features $\\boldsymbol{\\phi}_{\\text{lin}}(t)$ contain the current state and the delayed states.\n", "- The nonlinear features $\\boldsymbol{\\phi}_{\\text{nonlin}}(t)$ consist of products of the linear features.\n", "- The **linear features** $\\mathbf{O}_{\\text{lin}}(t)$ contain the current state variables and their delayed versions.\n", "- The **nonlinear features** $\\mathbf{O}_{\\text{nonlin}}(t)$ consist of polynomial combinations (products) of the linear features.\n", "\n", "The complete feature vector is formed by combining the bias, linear, and nonlinear terms, and this combined vector is used for training.\n", "\n" Loading @@ -485,8 +487,10 @@ "\n", "The linear feature vector at time $t$ is\n", "\n", "\n", "\n", "$$\n", "\\boldsymbol{\\phi}_{\\text{lin}}(t)\n", "\\mathbf{O}_{\\text{lin}}(t)\n", "=\n", "\\begin{bmatrix}\n", "x(t) \\\\\n", Loading @@ -495,10 +499,11 @@ "x(t-\\tau) \\\\\n", "y(t-\\tau) \\\\\n", "z(t-\\tau)\n", "\\end{bmatrix}.\n", "\\end{bmatrix}\n", "$$\n", "\n", "\n", "\n", "With $d=3$ variables and $k=2$ taps (current + delay),\n", "\n", "$$\n", Loading Loading @@ -545,17 +550,18 @@ "After removing the warmup samples, the NVAR training feature matrix is constructed as\n", "\n", "$$\n", "\\mathbf{\\Phi}_{\\text{train}} =\n", "\\mathbf{O}_{\\text{tot}}(t) =\n", "\\begin{bmatrix}\n", "1 \\\\\n", "\\boldsymbol{\\phi}_{\\text{lin}}(t) \\\\\n", "\\text{quadratic terms}\n", "\\end{bmatrix}.\n", "\\mathbf{O}_{\\text{lin}}(t) \\\\\n", "\\text{Nonlinear Terms}\n", "\\end{bmatrix}\n", "$$\n", "\n", "- The first row is the bias term.\n", "- The linear and quadratic features follow.\n", "- Only training time steps are included.\n" "- The linear and nonlinear features follow.\n", "- Only training time steps are included.\n", "\n" ] }, { Loading Loading @@ -604,7 +610,7 @@ "These form the nonlinear feature vector\n", "\n", "$$\n", "\\boldsymbol{\\phi}_{\\text{nonlin}}(t)\n", "\\mathbf{O}_{\\text{nonlin}}(t)\n", "=\n", "\\begin{bmatrix}\n", "x(t)^2 \\\\\n", Loading @@ -613,7 +619,7 @@ "\\vdots \\\\\n", "x(t)x(t-\\tau) \\\\\n", "y(t)z(t-\\tau)\n", "\\end{bmatrix}.\n", "\\end{bmatrix}\n", "$$\n", "\n", "\n", Loading Loading @@ -667,65 +673,67 @@ "- the linear delay features,\n", "- and the nonlinear (quadratic) features.\n", "\n", "\n", "\n", "$$\n", "\\mathbf{\\Phi}_{\\text{train}} =\n", "\\mathbf{O}_{\\text{tot}}(t) =\n", "\\begin{bmatrix}\n", "1 \\\\\n", "\\boldsymbol{\\phi}_{\\text{lin}}(t) \\\\\n", "\\boldsymbol{\\phi}_{\\text{nonlin}}(t)\n", "\\end{bmatrix}.\n", "\\mathbf{O}_{\\text{lin}}(t) \\\\\n", "\\mathbf{O}_{\\text{nonlin}}(t)\n", "\\end{bmatrix},\n", "\\qquad\n", "\\mathbf{O}_{\\text{tot}} \\in \\mathbb{R}^{28 \\times N_{\\text{train}}},\n", "$$\n", "\n", "Each column corresponds to one training time step.\n", "\n", "\n", "$$\n", "\\mathbf{\\Phi}_{\\text{train}} \\in \\mathbb{R}^{28 \\times N_{\\text{train}}}.\n", "$$\n" "where each column of $\\mathbf{O}_{\\text{tot}}$ corresponds to one training time step.\n" ] }, { "cell_type": "markdown", "id": "091d2571", "id": "8ddb7c0e", "metadata": {}, "source": [ "## 4. Training the Output Weights\n", "\n", "The output weight matrix $W_{\\text{out}}$ is learned using ridge regression. \n", "It maps the NVAR feature vector to the discrete state increment.\n", "It maps the NVAR observable vector to the discrete state increment.\n", "\n", "The model learns\n", "During training, the model learns the state increment\n", "\n", "$$\n", "\\Delta \\mathbf{x}(t)\n", "=\n", "\\mathbf{x}(t+1) - \\mathbf{x}(t).\n", "\\mathbf{x}(t+1) - \\mathbf{x}(t)\n", "$$\n", "\n", "Using the training feature matrix,\n", "Stacking the increments and observables over all training time steps gives\n", "\n", "$$\n", "\\Delta \\mathbf{X}_{\\text{train}}\n", "\\Delta \\mathbf{X}\n", "=\n", "W_{\\text{out}} \\, \\mathbf{\\Phi}_{\\text{train}}.\n", "W_{\\text{out}} \\, \\mathbf{O}_{\\text{tot}}\n", "$$\n", "\n", "where $\\mathbf{O}_{\\text{tot}}$ is the matrix whose columns are\n", "$\\mathbf{O}_{\\text{tot}}(t)$\n", "\n", "The ridge regression solution is\n", "\n", "$$\n", "W_{\\text{out}}\n", "=\n", "\\Delta \\mathbf{X}_{\\text{train}}\n", "\\, \\mathbf{\\Phi}_{\\text{train}}^{T}\n", "\\Delta \\mathbf{X}\n", "\\, \\mathbf{O}_{\\text{tot}}^{T}\n", "\\left(\n", "\\mathbf{\\Phi}_{\\text{train}} \\mathbf{\\Phi}_{\\text{train}}^{T}\n", "\\mathbf{O}_{\\text{tot}} \\mathbf{O}_{\\text{tot}}^{T}\n", "+ \\alpha I\n", "\\right)^{-1},\n", "\\right)^{-1}\n", "$$\n", "\n", "where $\\alpha$ is the regularization parameter.\n", "\n", "Only $W_{\\text{out}}$ is trained; the feature construction remains fixed.\n" "Only $W_{\\text{out}}$ is trained; the observable construction remains fixed.\n" ] }, { Loading Loading @@ -770,20 +778,20 @@ }, { "cell_type": "markdown", "id": "c02a357f", "id": "e937ec8e", "metadata": {}, "source": [ "## 5. Training Prediction\n", "## 5. Autonomous Prediction\n", "\n", "After learning $W_{\\text{out}}$, the model reconstructs the Lorenz trajectory\n", "by predicting the state increment.\n", "\n", "The model computes\n", "The model computes the predicted increment\n", "\n", "$$\n", "\\Delta \\mathbf{x}(t)\n", "\\Delta \\hat{\\mathbf{x}}(t)\n", "=\n", "W_{\\text{out}} \\, \\boldsymbol{\\phi}(t),\n", "W_{\\text{out}} \\, \\mathbf{O}_{\\text{tot}}(t),\n", "$$\n", "\n", "and the next state is obtained as\n", Loading @@ -791,14 +799,32 @@ "$$\n", "\\hat{\\mathbf{x}}(t+1)\n", "=\n", "\\mathbf{x}(t)\n", "\\hat{\\mathbf{x}}(t)\n", "+\n", "W_{\\text{out}} \\, \\boldsymbol{\\phi}(t).\n", "W_{\\text{out}} \\, \\mathbf{O}_{\\text{tot}}(t).\n", "$$\n", "\n", "Since the model learns the increment\n", "$\\Delta \\mathbf{x}(t) = \\mathbf{x}(t+1) - \\mathbf{x}(t)$,\n", "prediction is performed by adding this increment to the current state." "$\\Delta \\mathbf{x}(t) = \\mathbf{x}(t+1) - \\mathbf{x}(t)$\n", "during training, prediction is performed by adding the learned increment\n", "to the current predicted state.\n", "\n", "In autonomous mode, the model feeds back its own predicted state\n", "$\\hat{\\mathbf{x}}(t)$ at each time step.\n" ] }, { "cell_type": "markdown", "id": "0923d3f5", "metadata": {}, "source": [ "$$\n", "\\hat{\\mathbf{x}}(t+1)\n", "=\n", "\\hat{\\mathbf{x}}(t)\n", "+\n", "W_{\\text{out}} \\, \\mathbf{O}_{\\text{tot}}(t).\n", "$$\n" ] }, { Loading Loading @@ -938,14 +964,17 @@ "- Delay taps are shifted forward.\n", "- The next state is predicted as\n", "\n", "\n", "\n", "$$\n", "\\hat{\\mathbf{x}}(t+1)\n", "=\n", "\\hat{\\mathbf{x}}(t)\n", "+\n", "W_{\\text{out}} \\, \\boldsymbol{\\phi}(t).\n", "W_{\\text{out}} \\, \\mathbf{O}_{\\text{tot}}(t)\n", "$$\n", "\n", "\n", "The model uses its own previous prediction to generate the next state.\n", "\n", "#### Test Error (NRMSE)\n", Loading Loading @@ -1514,6 +1543,144 @@ "''' " ] }, { "cell_type": "code", "execution_count": 27, "id": "cb74bd11", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'\\n%matplotlib notebook\\n\\nimport matplotlib as mpl\\nmpl.rcParams[\\'animation.embed_limit\\'] = 200\\n\\nfrom matplotlib.animation import FuncAnimation\\nfrom IPython.display import HTML\\n\\n# Time settings\\nframes = plottime_pts\\ntrail = 100\\n\\n# Physical time window\\nt_error = t_eval[warmtrain_pts-1 : warmtrain_pts+plottime_pts-1]\\n\\n# Error computation\\ndist = np.sqrt(\\n np.sum(\\n (x[0:3, warmtrain_pts-1:warmtrain_pts+plottime_pts-1]\\n - x_test[0:3, :plottime_pts])**2,\\n axis=0\\n )\\n)\\n\\n# Lyapunov time marker\\ndt = t_eval[1] - t_eval[0]\\nt_boundary = t_eval[warmtrain_pts-1]\\nt_lyap = t_boundary + lyaptime_pts * dt\\n\\n# Create figure\\nfig = plt.figure(figsize=(13,6))\\nfig.suptitle(\"Actual vs Predicted Lorenz Trajectory\", y=0.98)\\n\\n# 2D Lorenz phase portrait (x vs z)\\nax_phase = fig.add_subplot(121)\\n\\nax_phase.set_xlim(np.min(x[0]), np.max(x[0]))\\nax_phase.set_ylim(np.min(x[2]), np.max(x[2]))\\n\\nax_phase.set_xlabel(\"x\")\\nax_phase.set_ylabel(\"z\")\\nax_phase.set_title(\"Phase Portrait (x vs z)\")\\nax_phase.grid(alpha=0.3)\\n\\ntrue_line, = ax_phase.plot([], [], color=\\'tab:blue\\', lw=2)\\ntrue_point, = ax_phase.plot([], [], \\'o\\', color=\\'tab:blue\\')\\n\\npred_line, = ax_phase.plot([], [], color=\\'tab:orange\\', lw=2, linestyle=\\'--\\')\\npred_point, = ax_phase.plot([], [], \\'o\\', color=\\'tab:orange\\')\\n\\nax_phase.legend([\"Ground Truth\", \"Prediction\"], loc=\"upper left\")\\n\\n# Prediction error plot\\nax_err = fig.add_subplot(122)\\n\\nax_err.set_xlim(t_error[0], t_error[-1])\\nax_err.set_ylim(0, np.max(dist)*1.1)\\n\\nax_err.set_title(\"Prediction Error\")\\nax_err.set_xlabel(\"Time\")\\nax_err.set_ylabel(\"Euclidean Distance\")\\nax_err.grid(alpha=0.3)\\n\\nax_err.axvline(t_lyap, linestyle=\\':\\', color=\\'tab:green\\', linewidth=1.2)\\nerr_line, = ax_err.plot([], [], color=\\'tab:red\\', lw=2, alpha=0.8)\\n\\n# Initialize animation\\ndef init():\\n true_line.set_data([], [])\\n pred_line.set_data([], [])\\n err_line.set_data([], [])\\n return true_line, pred_line, true_point, pred_point, err_line\\n\\n# Update animation\\ndef update(frame):\\n\\n start = max(0, frame - trail)\\n\\n # Actual trajectory\\n true_line.set_data(\\n x[0, warmtrain_pts-1+start:warmtrain_pts-1+frame],\\n x[2, warmtrain_pts-1+start:warmtrain_pts-1+frame]\\n )\\n\\n true_point.set_data(\\n [x[0, warmtrain_pts-1+frame]],\\n [x[2, warmtrain_pts-1+frame]]\\n )\\n\\n # Predicted trajectory\\n pred_line.set_data(\\n x_test[0, start:frame],\\n x_test[2, start:frame]\\n )\\n\\n pred_point.set_data(\\n [x_test[0, frame]],\\n [x_test[2, frame]]\\n )\\n\\n # Error evolution\\n err_line.set_data(t_error[:frame], dist[:frame])\\n\\n return true_line, pred_line, true_point, pred_point, err_line\\n\\n# Create animation\\nani = FuncAnimation(\\n fig,\\n update,\\n frames=frames,\\n init_func=init,\\n interval=30,\\n blit=False\\n)\\n\\nHTML(ani.to_jshtml())\\n'" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "'''\n", "%matplotlib notebook\n", "\n", "import matplotlib as mpl\n", "mpl.rcParams['animation.embed_limit'] = 200\n", "\n", "from matplotlib.animation import FuncAnimation\n", "from IPython.display import HTML\n", "\n", "# Time settings\n", "frames = plottime_pts\n", "trail = 100\n", "\n", "# Physical time window\n", "t_error = t_eval[warmtrain_pts-1 : warmtrain_pts+plottime_pts-1]\n", "\n", "# Error computation\n", "dist = np.sqrt(\n", " np.sum(\n", " (x[0:3, warmtrain_pts-1:warmtrain_pts+plottime_pts-1]\n", " - x_test[0:3, :plottime_pts])**2,\n", " axis=0\n", " )\n", ")\n", "\n", "# Lyapunov time marker\n", "dt = t_eval[1] - t_eval[0]\n", "t_boundary = t_eval[warmtrain_pts-1]\n", "t_lyap = t_boundary + lyaptime_pts * dt\n", "\n", "# Create figure\n", "fig = plt.figure(figsize=(13,6))\n", "fig.suptitle(\"Actual vs Predicted Lorenz Trajectory\", y=0.98)\n", "\n", "# 2D Lorenz phase portrait (x vs z)\n", "ax_phase = fig.add_subplot(121)\n", "\n", "ax_phase.set_xlim(np.min(x[0]), np.max(x[0]))\n", "ax_phase.set_ylim(np.min(x[2]), np.max(x[2]))\n", "\n", "ax_phase.set_xlabel(\"x\")\n", "ax_phase.set_ylabel(\"z\")\n", "ax_phase.set_title(\"Phase Portrait (x vs z)\")\n", "ax_phase.grid(alpha=0.3)\n", "\n", "true_line, = ax_phase.plot([], [], color='tab:blue', lw=2)\n", "true_point, = ax_phase.plot([], [], 'o', color='tab:blue')\n", "\n", "pred_line, = ax_phase.plot([], [], color='tab:orange', lw=2, linestyle='--')\n", "pred_point, = ax_phase.plot([], [], 'o', color='tab:orange')\n", "\n", "ax_phase.legend([\"Ground Truth\", \"Prediction\"], loc=\"upper left\")\n", "\n", "# Prediction error plot\n", "ax_err = fig.add_subplot(122)\n", "\n", "ax_err.set_xlim(t_error[0], t_error[-1])\n", "ax_err.set_ylim(0, np.max(dist)*1.1)\n", "\n", "ax_err.set_title(\"Prediction Error\")\n", "ax_err.set_xlabel(\"Time\")\n", "ax_err.set_ylabel(\"Euclidean Distance\")\n", "ax_err.grid(alpha=0.3)\n", "\n", "ax_err.axvline(t_lyap, linestyle=':', color='tab:green', linewidth=1.2)\n", "err_line, = ax_err.plot([], [], color='tab:red', lw=2, alpha=0.8)\n", "\n", "# Initialize animation\n", "def init():\n", " true_line.set_data([], [])\n", " pred_line.set_data([], [])\n", " err_line.set_data([], [])\n", " return true_line, pred_line, true_point, pred_point, err_line\n", "\n", "# Update animation\n", "def update(frame):\n", "\n", " start = max(0, frame - trail)\n", "\n", " # Actual trajectory\n", " true_line.set_data(\n", " x[0, warmtrain_pts-1+start:warmtrain_pts-1+frame],\n", " x[2, warmtrain_pts-1+start:warmtrain_pts-1+frame]\n", " )\n", "\n", " true_point.set_data(\n", " [x[0, warmtrain_pts-1+frame]],\n", " [x[2, warmtrain_pts-1+frame]]\n", " )\n", "\n", " # Predicted trajectory\n", " pred_line.set_data(\n", " x_test[0, start:frame],\n", " x_test[2, start:frame]\n", " )\n", "\n", " pred_point.set_data(\n", " [x_test[0, frame]],\n", " [x_test[2, frame]]\n", " )\n", "\n", " # Error evolution\n", " err_line.set_data(t_error[:frame], dist[:frame])\n", "\n", " return true_line, pred_line, true_point, pred_point, err_line\n", "\n", "# Create animation\n", "ani = FuncAnimation(\n", " fig,\n", " update,\n", " frames=frames,\n", " init_func=init,\n", " interval=30,\n", " blit=False\n", ")\n", "\n", "HTML(ani.to_jshtml())\n", "'''" ] }, { "cell_type": "code", "execution_count": 24, Loading
code/1_NGRC_Lorenz_Basic_Model.ipynb +227 −60 Original line number Diff line number Diff line Loading @@ -92,7 +92,7 @@ "### 1.1 Discrete-Time Parameters\n", "\n", "Time is discretized using a fixed step size $dt$. \n", "Each continuous-time interval is converted into a number of discrete samples by dividing by $dt$.\n", "Each continuous time interval is converted into a number of discrete samples by dividing by $dt$.\n", "\n", "| Name | Symbol | Definition | Value |\n", "|------|--------|------------|--------|\n", Loading Loading @@ -293,17 +293,16 @@ }, { "cell_type": "markdown", "id": "3810a2f4", "id": "6657a8ad", "metadata": {}, "source": [ "### 2.1 Lorenz Time-Series Generation\n", "\n", "\n", "The Lorenz63 system is integrated over $[0, T_{\\text{total}}]$ \n", "and sampled at fixed time intervals $dt$.\n", "The Lorenz63 system is integrated over the time interval $[0, T_{\\text{total}}]$\n", "and sampled at uniform time steps $dt$.\n", "\n", "$$\n", "t_n = n\\,dt, \\quad n = 0, 1, \\dots, N_{\\text{total}}.\n", "t_n = n\\,dt, \\quad n = 0, 1, \\dots, N_{\\text{total}}\n", "$$\n", "\n", "**Initial condition**\n", Loading @@ -312,7 +311,7 @@ "(x(0), y(0), z(0)) = (17.6772,\\; 12.9314,\\; 43.9140)\n", "$$\n", "\n", "The discrete state vector is\n", "The discrete state vector is defined as\n", "\n", "$$\n", "\\mathbf{x}(t_n) =\n", Loading @@ -320,10 +319,11 @@ "x(t_n) \\\\\n", "y(t_n) \\\\\n", "z(t_n)\n", "\\end{bmatrix}.\n", "\\end{bmatrix}\n", "$$\n", "\n", "This produces the trajectories $x(t_n)$, $y(t_n)$, and $z(t_n)$.\n" "This produces the time series $x(t_n)$, $y(t_n)$, and $z(t_n)$.\n", "\n" ] }, { Loading Loading @@ -375,7 +375,7 @@ "+\n", "\\mathrm{Var}\\big(y(t_n)\\big)\n", "+\n", "\\mathrm{Var}\\big(z(t_n)\\big).\n", "\\mathrm{Var}\\big(z(t_n)\\big)\n", "$$\n", "\n", "This is used to normalize the prediction error.\n", Loading Loading @@ -422,7 +422,7 @@ "y(t_0) & \\cdots & y(t_{N_{\\text{total}}}) \\\\\n", "z(t_0) & \\cdots & z(t_{N_{\\text{total}}})\n", "\\end{bmatrix}\n", "\\in \\mathbb{R}^{3 \\times (N_{\\text{total}}+1)}.\n", "\\in \\mathbb{R}^{3 \\times (N_{\\text{total}}+1)}\n", "$$\n", "\n", "The matrix is partitioned sequentially along the time axis:\n", Loading @@ -436,7 +436,7 @@ "\\mathbf{X}_{\\text{train}}\n", "\\;\\big|\\;\n", "\\mathbf{X}_{\\text{test}}\n", "\\big],\n", "\\big]\n", "$$\n", "\n", "where\n", Loading @@ -446,7 +446,7 @@ "\\quad\n", "\\mathbf{X}_{\\text{train}} \\in \\mathbb{R}^{3 \\times N_{\\text{train}}},\n", "\\quad\n", "\\mathbf{X}_{\\text{test}} \\in \\mathbb{R}^{3 \\times N_{\\text{test}}}.\n", "\\mathbf{X}_{\\text{test}} \\in \\mathbb{R}^{3 \\times N_{\\text{test}}}\n", "$$\n" ] }, Loading @@ -459,18 +459,20 @@ "\n", "The NVAR feature vector at time $t$ is built as\n", "\n", "\n", "$$\n", "\\boldsymbol{\\phi}(\\mathbf{x}(t)) =\n", "\\mathbf{O}_{\\text{tot}}(t) =\n", "\\begin{bmatrix}\n", "\\text{Bias} \\\\\n", "\\text{Linear Delay Terms} \\\\\n", "\\text{Nonlinear Terms}\n", "\\end{bmatrix}.\n", "\\end{bmatrix}\n", "$$\n", "\n", "\n", "- The bias term is a constant equal to 1.\n", "- The linear features $\\boldsymbol{\\phi}_{\\text{lin}}(t)$ contain the current state and the delayed states.\n", "- The nonlinear features $\\boldsymbol{\\phi}_{\\text{nonlin}}(t)$ consist of products of the linear features.\n", "- The **linear features** $\\mathbf{O}_{\\text{lin}}(t)$ contain the current state variables and their delayed versions.\n", "- The **nonlinear features** $\\mathbf{O}_{\\text{nonlin}}(t)$ consist of polynomial combinations (products) of the linear features.\n", "\n", "The complete feature vector is formed by combining the bias, linear, and nonlinear terms, and this combined vector is used for training.\n", "\n" Loading @@ -485,8 +487,10 @@ "\n", "The linear feature vector at time $t$ is\n", "\n", "\n", "\n", "$$\n", "\\boldsymbol{\\phi}_{\\text{lin}}(t)\n", "\\mathbf{O}_{\\text{lin}}(t)\n", "=\n", "\\begin{bmatrix}\n", "x(t) \\\\\n", Loading @@ -495,10 +499,11 @@ "x(t-\\tau) \\\\\n", "y(t-\\tau) \\\\\n", "z(t-\\tau)\n", "\\end{bmatrix}.\n", "\\end{bmatrix}\n", "$$\n", "\n", "\n", "\n", "With $d=3$ variables and $k=2$ taps (current + delay),\n", "\n", "$$\n", Loading Loading @@ -545,17 +550,18 @@ "After removing the warmup samples, the NVAR training feature matrix is constructed as\n", "\n", "$$\n", "\\mathbf{\\Phi}_{\\text{train}} =\n", "\\mathbf{O}_{\\text{tot}}(t) =\n", "\\begin{bmatrix}\n", "1 \\\\\n", "\\boldsymbol{\\phi}_{\\text{lin}}(t) \\\\\n", "\\text{quadratic terms}\n", "\\end{bmatrix}.\n", "\\mathbf{O}_{\\text{lin}}(t) \\\\\n", "\\text{Nonlinear Terms}\n", "\\end{bmatrix}\n", "$$\n", "\n", "- The first row is the bias term.\n", "- The linear and quadratic features follow.\n", "- Only training time steps are included.\n" "- The linear and nonlinear features follow.\n", "- Only training time steps are included.\n", "\n" ] }, { Loading Loading @@ -604,7 +610,7 @@ "These form the nonlinear feature vector\n", "\n", "$$\n", "\\boldsymbol{\\phi}_{\\text{nonlin}}(t)\n", "\\mathbf{O}_{\\text{nonlin}}(t)\n", "=\n", "\\begin{bmatrix}\n", "x(t)^2 \\\\\n", Loading @@ -613,7 +619,7 @@ "\\vdots \\\\\n", "x(t)x(t-\\tau) \\\\\n", "y(t)z(t-\\tau)\n", "\\end{bmatrix}.\n", "\\end{bmatrix}\n", "$$\n", "\n", "\n", Loading Loading @@ -667,65 +673,67 @@ "- the linear delay features,\n", "- and the nonlinear (quadratic) features.\n", "\n", "\n", "\n", "$$\n", "\\mathbf{\\Phi}_{\\text{train}} =\n", "\\mathbf{O}_{\\text{tot}}(t) =\n", "\\begin{bmatrix}\n", "1 \\\\\n", "\\boldsymbol{\\phi}_{\\text{lin}}(t) \\\\\n", "\\boldsymbol{\\phi}_{\\text{nonlin}}(t)\n", "\\end{bmatrix}.\n", "\\mathbf{O}_{\\text{lin}}(t) \\\\\n", "\\mathbf{O}_{\\text{nonlin}}(t)\n", "\\end{bmatrix},\n", "\\qquad\n", "\\mathbf{O}_{\\text{tot}} \\in \\mathbb{R}^{28 \\times N_{\\text{train}}},\n", "$$\n", "\n", "Each column corresponds to one training time step.\n", "\n", "\n", "$$\n", "\\mathbf{\\Phi}_{\\text{train}} \\in \\mathbb{R}^{28 \\times N_{\\text{train}}}.\n", "$$\n" "where each column of $\\mathbf{O}_{\\text{tot}}$ corresponds to one training time step.\n" ] }, { "cell_type": "markdown", "id": "091d2571", "id": "8ddb7c0e", "metadata": {}, "source": [ "## 4. Training the Output Weights\n", "\n", "The output weight matrix $W_{\\text{out}}$ is learned using ridge regression. \n", "It maps the NVAR feature vector to the discrete state increment.\n", "It maps the NVAR observable vector to the discrete state increment.\n", "\n", "The model learns\n", "During training, the model learns the state increment\n", "\n", "$$\n", "\\Delta \\mathbf{x}(t)\n", "=\n", "\\mathbf{x}(t+1) - \\mathbf{x}(t).\n", "\\mathbf{x}(t+1) - \\mathbf{x}(t)\n", "$$\n", "\n", "Using the training feature matrix,\n", "Stacking the increments and observables over all training time steps gives\n", "\n", "$$\n", "\\Delta \\mathbf{X}_{\\text{train}}\n", "\\Delta \\mathbf{X}\n", "=\n", "W_{\\text{out}} \\, \\mathbf{\\Phi}_{\\text{train}}.\n", "W_{\\text{out}} \\, \\mathbf{O}_{\\text{tot}}\n", "$$\n", "\n", "where $\\mathbf{O}_{\\text{tot}}$ is the matrix whose columns are\n", "$\\mathbf{O}_{\\text{tot}}(t)$\n", "\n", "The ridge regression solution is\n", "\n", "$$\n", "W_{\\text{out}}\n", "=\n", "\\Delta \\mathbf{X}_{\\text{train}}\n", "\\, \\mathbf{\\Phi}_{\\text{train}}^{T}\n", "\\Delta \\mathbf{X}\n", "\\, \\mathbf{O}_{\\text{tot}}^{T}\n", "\\left(\n", "\\mathbf{\\Phi}_{\\text{train}} \\mathbf{\\Phi}_{\\text{train}}^{T}\n", "\\mathbf{O}_{\\text{tot}} \\mathbf{O}_{\\text{tot}}^{T}\n", "+ \\alpha I\n", "\\right)^{-1},\n", "\\right)^{-1}\n", "$$\n", "\n", "where $\\alpha$ is the regularization parameter.\n", "\n", "Only $W_{\\text{out}}$ is trained; the feature construction remains fixed.\n" "Only $W_{\\text{out}}$ is trained; the observable construction remains fixed.\n" ] }, { Loading Loading @@ -770,20 +778,20 @@ }, { "cell_type": "markdown", "id": "c02a357f", "id": "e937ec8e", "metadata": {}, "source": [ "## 5. Training Prediction\n", "## 5. Autonomous Prediction\n", "\n", "After learning $W_{\\text{out}}$, the model reconstructs the Lorenz trajectory\n", "by predicting the state increment.\n", "\n", "The model computes\n", "The model computes the predicted increment\n", "\n", "$$\n", "\\Delta \\mathbf{x}(t)\n", "\\Delta \\hat{\\mathbf{x}}(t)\n", "=\n", "W_{\\text{out}} \\, \\boldsymbol{\\phi}(t),\n", "W_{\\text{out}} \\, \\mathbf{O}_{\\text{tot}}(t),\n", "$$\n", "\n", "and the next state is obtained as\n", Loading @@ -791,14 +799,32 @@ "$$\n", "\\hat{\\mathbf{x}}(t+1)\n", "=\n", "\\mathbf{x}(t)\n", "\\hat{\\mathbf{x}}(t)\n", "+\n", "W_{\\text{out}} \\, \\boldsymbol{\\phi}(t).\n", "W_{\\text{out}} \\, \\mathbf{O}_{\\text{tot}}(t).\n", "$$\n", "\n", "Since the model learns the increment\n", "$\\Delta \\mathbf{x}(t) = \\mathbf{x}(t+1) - \\mathbf{x}(t)$,\n", "prediction is performed by adding this increment to the current state." "$\\Delta \\mathbf{x}(t) = \\mathbf{x}(t+1) - \\mathbf{x}(t)$\n", "during training, prediction is performed by adding the learned increment\n", "to the current predicted state.\n", "\n", "In autonomous mode, the model feeds back its own predicted state\n", "$\\hat{\\mathbf{x}}(t)$ at each time step.\n" ] }, { "cell_type": "markdown", "id": "0923d3f5", "metadata": {}, "source": [ "$$\n", "\\hat{\\mathbf{x}}(t+1)\n", "=\n", "\\hat{\\mathbf{x}}(t)\n", "+\n", "W_{\\text{out}} \\, \\mathbf{O}_{\\text{tot}}(t).\n", "$$\n" ] }, { Loading Loading @@ -938,14 +964,17 @@ "- Delay taps are shifted forward.\n", "- The next state is predicted as\n", "\n", "\n", "\n", "$$\n", "\\hat{\\mathbf{x}}(t+1)\n", "=\n", "\\hat{\\mathbf{x}}(t)\n", "+\n", "W_{\\text{out}} \\, \\boldsymbol{\\phi}(t).\n", "W_{\\text{out}} \\, \\mathbf{O}_{\\text{tot}}(t)\n", "$$\n", "\n", "\n", "The model uses its own previous prediction to generate the next state.\n", "\n", "#### Test Error (NRMSE)\n", Loading Loading @@ -1514,6 +1543,144 @@ "''' " ] }, { "cell_type": "code", "execution_count": 27, "id": "cb74bd11", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'\\n%matplotlib notebook\\n\\nimport matplotlib as mpl\\nmpl.rcParams[\\'animation.embed_limit\\'] = 200\\n\\nfrom matplotlib.animation import FuncAnimation\\nfrom IPython.display import HTML\\n\\n# Time settings\\nframes = plottime_pts\\ntrail = 100\\n\\n# Physical time window\\nt_error = t_eval[warmtrain_pts-1 : warmtrain_pts+plottime_pts-1]\\n\\n# Error computation\\ndist = np.sqrt(\\n np.sum(\\n (x[0:3, warmtrain_pts-1:warmtrain_pts+plottime_pts-1]\\n - x_test[0:3, :plottime_pts])**2,\\n axis=0\\n )\\n)\\n\\n# Lyapunov time marker\\ndt = t_eval[1] - t_eval[0]\\nt_boundary = t_eval[warmtrain_pts-1]\\nt_lyap = t_boundary + lyaptime_pts * dt\\n\\n# Create figure\\nfig = plt.figure(figsize=(13,6))\\nfig.suptitle(\"Actual vs Predicted Lorenz Trajectory\", y=0.98)\\n\\n# 2D Lorenz phase portrait (x vs z)\\nax_phase = fig.add_subplot(121)\\n\\nax_phase.set_xlim(np.min(x[0]), np.max(x[0]))\\nax_phase.set_ylim(np.min(x[2]), np.max(x[2]))\\n\\nax_phase.set_xlabel(\"x\")\\nax_phase.set_ylabel(\"z\")\\nax_phase.set_title(\"Phase Portrait (x vs z)\")\\nax_phase.grid(alpha=0.3)\\n\\ntrue_line, = ax_phase.plot([], [], color=\\'tab:blue\\', lw=2)\\ntrue_point, = ax_phase.plot([], [], \\'o\\', color=\\'tab:blue\\')\\n\\npred_line, = ax_phase.plot([], [], color=\\'tab:orange\\', lw=2, linestyle=\\'--\\')\\npred_point, = ax_phase.plot([], [], \\'o\\', color=\\'tab:orange\\')\\n\\nax_phase.legend([\"Ground Truth\", \"Prediction\"], loc=\"upper left\")\\n\\n# Prediction error plot\\nax_err = fig.add_subplot(122)\\n\\nax_err.set_xlim(t_error[0], t_error[-1])\\nax_err.set_ylim(0, np.max(dist)*1.1)\\n\\nax_err.set_title(\"Prediction Error\")\\nax_err.set_xlabel(\"Time\")\\nax_err.set_ylabel(\"Euclidean Distance\")\\nax_err.grid(alpha=0.3)\\n\\nax_err.axvline(t_lyap, linestyle=\\':\\', color=\\'tab:green\\', linewidth=1.2)\\nerr_line, = ax_err.plot([], [], color=\\'tab:red\\', lw=2, alpha=0.8)\\n\\n# Initialize animation\\ndef init():\\n true_line.set_data([], [])\\n pred_line.set_data([], [])\\n err_line.set_data([], [])\\n return true_line, pred_line, true_point, pred_point, err_line\\n\\n# Update animation\\ndef update(frame):\\n\\n start = max(0, frame - trail)\\n\\n # Actual trajectory\\n true_line.set_data(\\n x[0, warmtrain_pts-1+start:warmtrain_pts-1+frame],\\n x[2, warmtrain_pts-1+start:warmtrain_pts-1+frame]\\n )\\n\\n true_point.set_data(\\n [x[0, warmtrain_pts-1+frame]],\\n [x[2, warmtrain_pts-1+frame]]\\n )\\n\\n # Predicted trajectory\\n pred_line.set_data(\\n x_test[0, start:frame],\\n x_test[2, start:frame]\\n )\\n\\n pred_point.set_data(\\n [x_test[0, frame]],\\n [x_test[2, frame]]\\n )\\n\\n # Error evolution\\n err_line.set_data(t_error[:frame], dist[:frame])\\n\\n return true_line, pred_line, true_point, pred_point, err_line\\n\\n# Create animation\\nani = FuncAnimation(\\n fig,\\n update,\\n frames=frames,\\n init_func=init,\\n interval=30,\\n blit=False\\n)\\n\\nHTML(ani.to_jshtml())\\n'" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "'''\n", "%matplotlib notebook\n", "\n", "import matplotlib as mpl\n", "mpl.rcParams['animation.embed_limit'] = 200\n", "\n", "from matplotlib.animation import FuncAnimation\n", "from IPython.display import HTML\n", "\n", "# Time settings\n", "frames = plottime_pts\n", "trail = 100\n", "\n", "# Physical time window\n", "t_error = t_eval[warmtrain_pts-1 : warmtrain_pts+plottime_pts-1]\n", "\n", "# Error computation\n", "dist = np.sqrt(\n", " np.sum(\n", " (x[0:3, warmtrain_pts-1:warmtrain_pts+plottime_pts-1]\n", " - x_test[0:3, :plottime_pts])**2,\n", " axis=0\n", " )\n", ")\n", "\n", "# Lyapunov time marker\n", "dt = t_eval[1] - t_eval[0]\n", "t_boundary = t_eval[warmtrain_pts-1]\n", "t_lyap = t_boundary + lyaptime_pts * dt\n", "\n", "# Create figure\n", "fig = plt.figure(figsize=(13,6))\n", "fig.suptitle(\"Actual vs Predicted Lorenz Trajectory\", y=0.98)\n", "\n", "# 2D Lorenz phase portrait (x vs z)\n", "ax_phase = fig.add_subplot(121)\n", "\n", "ax_phase.set_xlim(np.min(x[0]), np.max(x[0]))\n", "ax_phase.set_ylim(np.min(x[2]), np.max(x[2]))\n", "\n", "ax_phase.set_xlabel(\"x\")\n", "ax_phase.set_ylabel(\"z\")\n", "ax_phase.set_title(\"Phase Portrait (x vs z)\")\n", "ax_phase.grid(alpha=0.3)\n", "\n", "true_line, = ax_phase.plot([], [], color='tab:blue', lw=2)\n", "true_point, = ax_phase.plot([], [], 'o', color='tab:blue')\n", "\n", "pred_line, = ax_phase.plot([], [], color='tab:orange', lw=2, linestyle='--')\n", "pred_point, = ax_phase.plot([], [], 'o', color='tab:orange')\n", "\n", "ax_phase.legend([\"Ground Truth\", \"Prediction\"], loc=\"upper left\")\n", "\n", "# Prediction error plot\n", "ax_err = fig.add_subplot(122)\n", "\n", "ax_err.set_xlim(t_error[0], t_error[-1])\n", "ax_err.set_ylim(0, np.max(dist)*1.1)\n", "\n", "ax_err.set_title(\"Prediction Error\")\n", "ax_err.set_xlabel(\"Time\")\n", "ax_err.set_ylabel(\"Euclidean Distance\")\n", "ax_err.grid(alpha=0.3)\n", "\n", "ax_err.axvline(t_lyap, linestyle=':', color='tab:green', linewidth=1.2)\n", "err_line, = ax_err.plot([], [], color='tab:red', lw=2, alpha=0.8)\n", "\n", "# Initialize animation\n", "def init():\n", " true_line.set_data([], [])\n", " pred_line.set_data([], [])\n", " err_line.set_data([], [])\n", " return true_line, pred_line, true_point, pred_point, err_line\n", "\n", "# Update animation\n", "def update(frame):\n", "\n", " start = max(0, frame - trail)\n", "\n", " # Actual trajectory\n", " true_line.set_data(\n", " x[0, warmtrain_pts-1+start:warmtrain_pts-1+frame],\n", " x[2, warmtrain_pts-1+start:warmtrain_pts-1+frame]\n", " )\n", "\n", " true_point.set_data(\n", " [x[0, warmtrain_pts-1+frame]],\n", " [x[2, warmtrain_pts-1+frame]]\n", " )\n", "\n", " # Predicted trajectory\n", " pred_line.set_data(\n", " x_test[0, start:frame],\n", " x_test[2, start:frame]\n", " )\n", "\n", " pred_point.set_data(\n", " [x_test[0, frame]],\n", " [x_test[2, frame]]\n", " )\n", "\n", " # Error evolution\n", " err_line.set_data(t_error[:frame], dist[:frame])\n", "\n", " return true_line, pred_line, true_point, pred_point, err_line\n", "\n", "# Create animation\n", "ani = FuncAnimation(\n", " fig,\n", " update,\n", " frames=frames,\n", " init_func=init,\n", " interval=30,\n", " blit=False\n", ")\n", "\n", "HTML(ani.to_jshtml())\n", "'''" ] }, { "cell_type": "code", "execution_count": 24,
