.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/plot_narx_msa.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code or to run this example in your browser via JupyterLite. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_plot_narx_msa.py: =========================== Multi-step-ahead NARX model =========================== .. currentmodule:: fastcan In this example, we will compare one-step-ahead NARX and multi-step-ahead NARX. .. GENERATED FROM PYTHON SOURCE LINES 10-14 .. code-block:: Python # Authors: The fastcan developers # SPDX-License-Identifier: MIT .. GENERATED FROM PYTHON SOURCE LINES 15-28 Nonlinear system ---------------- `Duffing equation `_ is used to generate simulated data. The mathematical model is given by .. math:: \ddot{y} + 0.1\dot{y} - y + 0.25y^3 = u where :math:`y` is the output signal and :math:`u` is the input signal, which is :math:`u(t) = 2.5\cos(2\pi t)`. The phase portraits and the vector field of the Duffing equation are shown below. .. GENERATED FROM PYTHON SOURCE LINES 28-103 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from scipy.integrate import odeint def duffing_equation(y, t): """Non-autonomous system""" # y1 is displacement and y2 is velocity y1, y2 = y # u is sinusoidal input u = 2.5 * np.cos(2 * np.pi * t) # dydt is derivative of y1 and y2 dydt = [y2, -0.1 * y2 + y1 - 0.25 * y1**3 + u] return dydt def auto_duffing_equation(y, t): """Autonomous system""" y1, y2 = y dydt = [y2, -0.1 * y2 + y1 - 0.25 * y1**3] return dydt dur = 10 n_samples = 1000 y0 = None if y0 is None: n_init = 10 x0 = np.linspace(0, 2, n_init) y0_y = np.cos(np.pi * x0) y0_x = np.sin(np.pi * x0) y0 = np.c_[y0_x, y0_y] else: n_init = len(y0) t = np.linspace(0, dur, n_samples) sol = np.zeros((n_init, n_samples, 2)) for i in range(n_init): sol[i] = odeint(auto_duffing_equation, y0[i], t) # Phase portraits for i in range(n_init): plt.plot(sol[i, :, 0], sol[i, :, 1], c="tab:blue") # Vector field y_min = np.min(sol[:, :, 0]) - 0.2 y_max = np.max(sol[:, :, 0]) + 0.2 dot_y_min = np.min(sol[:, :, 1]) - 0.2 dot_y_max = np.max(sol[:, :, 1]) + 0.2 y, dot_y = np.meshgrid( np.linspace(y_min, y_max, 30), np.linspace(dot_y_min, dot_y_max, 30) ) ddot_y = auto_duffing_equation([y, dot_y], 0)[1] plt.streamplot( y, dot_y, dot_y, ddot_y, color=(0.5, 0.5, 0.5, 0.3), density=1.5, minlength=0.02, maxlength=0.1, linewidth=0.5, arrowsize=0.5, ) plt.xlim(y_min, y_max) plt.ylim(dot_y_min, dot_y_max) plt.title("Phase portraits and vector field of Duffing equation") plt.xlabel("y(t)") plt.ylabel("dy/dt(t)") plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_narx_msa_001.png :alt: Phase portraits and vector field of Duffing equation :srcset: /auto_examples/images/sphx_glr_plot_narx_msa_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 104-109 Generate training-test data --------------------------- In the phase portraits, it is shown that the system has two stable equilibria. We use one to generate training data and the other to generate test data. .. GENERATED FROM PYTHON SOURCE LINES 109-132 .. code-block:: Python # 10 s duration with 0.01 Hz sampling time, # so 1000 samples in total for each measurement dur = 10 n_samples = 1000 t = np.linspace(0, dur, n_samples) # External excitation is the same for each measurement u = 2.5 * np.cos(2 * np.pi * t).reshape(-1, 1) # Small additional white noise rng = np.random.default_rng(12345) e_train_0 = rng.normal(0, 0.0004, n_samples) e_test = rng.normal(0, 0.0004, n_samples) # Solve differential equation to get displacement as y # Initial condition at displacement 0.6 and velocity 0.8 sol = odeint(duffing_equation, [0.6, 0.8], t) y_train_0 = sol[:, 0] + e_train_0 # Initial condition at displacement 0.6 and velocity -0.8 sol = odeint(duffing_equation, [0.6, -0.8], t) y_test = sol[:, 0] + e_test .. GENERATED FROM PYTHON SOURCE LINES 133-140 One-step-head VS. multi-step-ahead NARX --------------------------------------- First, we use :meth:`make_narx` to obtain the reduced NARX model. Then, the NARX model will be fitted with one-step-ahead predictor and multi-step-ahead predictor, respectively. Generally, the training of one-step-ahead (OSA) NARX is faster, while the multi-step-ahead (MSA) NARX is more accurate. .. GENERATED FROM PYTHON SOURCE LINES 140-185 .. code-block:: Python from sklearn.metrics import r2_score from fastcan.narx import make_narx max_delay = 3 narx_model = make_narx( X=u, y=y_train_0, n_terms_to_select=5, max_delay=max_delay, poly_degree=3, verbose=0, ) def plot_prediction(ax, t, y_true, y_pred, title): ax.plot(t, y_true, label="true") ax.plot(t, y_pred, label="predicted") ax.legend() ax.set_title(f"{title} (R2: {r2_score(y_true, y_pred):.5f})") ax.set_xlabel("t (s)") ax.set_ylabel("y(t)") # OSA NARX narx_model.fit(u, y_train_0) y_train_0_osa_pred = narx_model.predict(u, y_init=y_train_0[:max_delay]) y_test_osa_pred = narx_model.predict(u, y_init=y_test[:max_delay]) # MSA NARX narx_model.fit(u, y_train_0, coef_init="one_step_ahead") y_train_0_msa_pred = narx_model.predict(u, y_init=y_train_0[:max_delay]) y_test_msa_pred = narx_model.predict(u, y_init=y_test[:max_delay]) fig, ax = plt.subplots(2, 2, figsize=(8, 6)) plot_prediction(ax[0, 0], t, y_train_0, y_train_0_osa_pred, "OSA NARX on Train 0") plot_prediction(ax[0, 1], t, y_train_0, y_train_0_msa_pred, "MSA NARX on Train 0") plot_prediction(ax[1, 0], t, y_test, y_test_osa_pred, "OSA NARX on Test") plot_prediction(ax[1, 1], t, y_test, y_test_msa_pred, "MSA NARX on Test") fig.tight_layout() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_narx_msa_002.png :alt: OSA NARX on Train 0 (R2: 0.69558), MSA NARX on Train 0 (R2: 0.99997), OSA NARX on Test (R2: -8.81038), MSA NARX on Test (R2: 0.72033) :srcset: /auto_examples/images/sphx_glr_plot_narx_msa_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 186-196 Multiple measurement sessions ----------------------------- The plot above shows that the NARX model cannot capture the dynamics at the left equilibrium shown in the phase portraits. To improve the performance, let us append another measurement session to the training data to include the dynamics of both equilibria. Here, we need to insert (at least max_delay number of) `np.nan` to indicate the model that the original training data and the appended data are from different measurement sessions. The plot shows that the prediction performance of the NARX on test data has been largely improved. .. GENERATED FROM PYTHON SOURCE LINES 196-234 .. code-block:: Python e_train_1 = rng.normal(0, 0.0004, n_samples) # Solve differential equation to get displacement as y # Initial condition at displacement 0.5 and velocity -1 sol = odeint(duffing_equation, [0.5, -1], t) y_train_1 = sol[:, 0] + e_train_1 u_all = np.r_[u, [[np.nan]] * max_delay, u] y_all = np.r_[y_train_0, [np.nan] * max_delay, y_train_1] narx_model = make_narx( X=u_all, y=y_all, n_terms_to_select=5, max_delay=max_delay, poly_degree=3, verbose=0, ) narx_model.fit(u_all, y_all) y_train_0_osa_pred = narx_model.predict(u, y_init=y_train_0[:max_delay]) y_train_1_osa_pred = narx_model.predict(u, y_init=y_train_1[:max_delay]) y_test_osa_pred = narx_model.predict(u, y_init=y_test[:max_delay]) narx_model.fit(u_all, y_all, coef_init="one_step_ahead") y_train_0_msa_pred = narx_model.predict(u, y_init=y_train_0[:max_delay]) y_train_1_msa_pred = narx_model.predict(u, y_init=y_train_1[:max_delay]) y_test_msa_pred = narx_model.predict(u, y_init=y_test[:max_delay]) fig, ax = plt.subplots(3, 2, figsize=(8, 9)) plot_prediction(ax[0, 0], t, y_train_0, y_train_0_osa_pred, "OSA NARX on Train 0") plot_prediction(ax[0, 1], t, y_train_0, y_train_0_msa_pred, "MSA NARX on Train 0") plot_prediction(ax[1, 0], t, y_train_1, y_train_1_osa_pred, "OSA NARX on Train 1") plot_prediction(ax[1, 1], t, y_train_1, y_train_1_msa_pred, "MSA NARX on Train 1") plot_prediction(ax[2, 0], t, y_test, y_test_osa_pred, "OSA NARX on Test") plot_prediction(ax[2, 1], t, y_test, y_test_msa_pred, "MSA NARX on Test") fig.tight_layout() plt.show() .. image-sg:: /auto_examples/images/sphx_glr_plot_narx_msa_003.png :alt: OSA NARX on Train 0 (R2: 0.78968), MSA NARX on Train 0 (R2: 0.99951), OSA NARX on Train 1 (R2: 0.82392), MSA NARX on Train 1 (R2: 0.99944), OSA NARX on Test (R2: 0.77897), MSA NARX on Test (R2: 0.97492) :srcset: /auto_examples/images/sphx_glr_plot_narx_msa_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 3.808 seconds) .. _sphx_glr_download_auto_examples_plot_narx_msa.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: lite-badge .. image:: images/jupyterlite_badge_logo.svg :target: ../lite/lab/index.html?path=auto_examples/plot_narx_msa.ipynb :alt: Launch JupyterLite :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_narx_msa.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_narx_msa.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_narx_msa.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_