PySINDy_for_ControlledSystemIdentification.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/GirolamoOddo/AppliedMath_Notebooks/blob/main/PySINDy_for_ControlledSystemIdentification.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Vp45YnxEVNtd"
      },
      "source": [
        "#Overview PySINDy Colab   \n",
        "Application of pysindy in colab environment, here you can find:\n",
        "\n",
        "\n",
        "*   Train and Test SINDy Model w/o control inputs\n",
        "*   Train and Test SINDy Model w/ control inputs and\n",
        "    custom library\n",
        "\n",
        "This notebook is self-contained and can work with either synthetic or imported data.  \n",
        "Please note the following core dependencies:  \n",
        ">numpy: 1.23.5  \n",
        ">pysindy: 1.7.5  \n",
        ">pandas: 1.5.3\n",
        "\n",
        "For each case, there is also a plot section that highlights both the system and, if applicable, the control inputs.\n",
        "\n",
        "\n",
        "\n",
        "USEFUL RESOURCES:\n",
        "\n",
        "* OFFICIAL GIT: https://github.com/dynamicslab/pysindy  \n",
        "\n",
        "\n",
        "\n",
        "* OFFICIAL DOC: https://pysindy.readthedocs.io/en/latest/index.html  \n",
        "\n",
        "* UNOFFICIAL NOTES: https://uvadlc-notebooks.readthedocs.io/en/latest/tutorial_notebooks/DL2/Dynamical_Neural_Networks/Complete_DNN_2_1.html?highlight=PySINDy\n",
        "\n",
        "REFERENCES:\n",
        "\n",
        "* Brian M. de Silva, Kathleen Champion, Markus Quade, Jean-Christophe Loiseau, J. Nathan Kutz, and Steven L. Brunton., (2020). PySINDy: A Python package for the sparse identification of nonlinear dynamical systems from data. Journal of Open Source Software, 5(49), 2104, https://doi.org/10.21105/joss.02104\n",
        "\n",
        "* Kaptanoglu et al., (2022). PySINDy: A comprehensive Python package for robust sparse system identification. Journal of Open Source Software, 7(69), 3994, https://doi.org/10.21105/joss.03994\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0Zdlkedx9xI_"
      },
      "source": [
        "\n",
        "\n",
        "\n",
        "  \n",
        "  \n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CLdBVlvk06nz"
      },
      "source": [
        "# SINDy Download and Data Acqusition"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "TGM3okGScIpb",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "# @title Import SINDy (mandatory)\n",
        "\n",
        "%%capture\n",
        "\n",
        "!rm -rf pysindy\n",
        "!rm -rf sample_data\n",
        "\n",
        "!git clone https://github.com/dynamicslab/pysindy.git\n",
        "!cd pysindy/pysindy/.\n",
        "!ls\n",
        "!pip install pysindy\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "51Ywr_Jfb1Nr"
      },
      "outputs": [],
      "source": [
        "# @title Import Utilities (mandatory)\n",
        "import warnings\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import pysindy as ps\n",
        "import plotly.graph_objs as go\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import odeint\n",
        "from scipy.integrate import solve_ivp\n",
        "from sklearn.linear_model import Lasso\n",
        "from mpl_toolkits.mplot3d import Axes3D\n",
        "from pysindy.feature_library import PolynomialLibrary, FourierLibrary, IdentityLibrary\n",
        "\n",
        "\n",
        "#_______________________________________________________________________________\n",
        "# CHECK PACKAGES VERSION:\n",
        "#_______________________________________________________________________________\n",
        "#import os\n",
        "#packages = [\n",
        "#    'numpy',\n",
        "#    'pysindy',\n",
        "#    'pandas',\n",
        "#    'plotly'\n",
        "#    'scipy',\n",
        "#    'sklearn'\n",
        "#]\n",
        "#\n",
        "#for package in packages:\n",
        "#    try:\n",
        "#        version = os.popen(f\"pip show {package} | grep Version\").read().strip().split(': ')[1]\n",
        "#        print(f\"{package}: {version}\")\n",
        "#    except Exception as e:\n",
        "#        print(f\"Error getting version for {package}: {e}\")\n",
        "#\n",
        "#_______________________________________________________________________________\n",
        "# THE FOLLOWING CODE WORKS PROPERLY WITH THESE PACKAGES VERSION:\n",
        "#_______________________________________________________________________________\n",
        "#numpy: 1.23.5\n",
        "#pysindy: 1.7.5\n",
        "#pandas: 1.5.3\n",
        "#_______________________________________________________________________________"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7Wnxm6Tx5k1N"
      },
      "source": [
        "Now choose whether to: Load your own data or Generate a synthetic dataset."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "pMoAbl-VpyvC"
      },
      "outputs": [],
      "source": [
        "# @title Load your own data (if you have it)\n",
        "# Import custom data\n",
        "\n",
        "#from google.colab import drive\n",
        "#drive.mount('/content/drive',force_remount=True)\n",
        "\n",
        "# Correct the path for your Drive\n",
        "#train_data_path = '/content/drive/MyDrive/colab_notebooks/colab_set/colab_train'\n",
        "#test_data_path = '/content/drive/MyDrive/colab_notebooks/colab_set/colab_test'\n",
        "\n",
        "#x_train = pd.read_csv(train_data_path+\"train_sindy.csv\")\n",
        "#x_test = pd.read_csv(test_data_path+\"test_sindy.csv\")\n",
        "#dt = 0.01\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "Ihn03iuKus4m"
      },
      "outputs": [],
      "source": [
        "# @title Generate sintetic data (recommended option)\n",
        "# Simulate data\n",
        "dt = 0.00005\n",
        "time = np.arange(0, 20, dt)\n",
        "time_test = np.arange(0, 4, dt)\n",
        "coeff = np.array([1, 1, 1])\n",
        "\n",
        "# All 'Train data' are editable\n",
        "# Train data: Variables\n",
        "x = np.sin(time*4) + coeff[0]\n",
        "x = x.reshape(-1, 1)\n",
        "\n",
        "y = np.cos(time*4) + coeff[1]\n",
        "y = y.reshape(-1, 1)\n",
        "\n",
        "z = np.exp(-time) + coeff[2]\n",
        "z = z.reshape(-1, 1)\n",
        "\n",
        "x0 = np.array([x[0], y[0], z[0]])\n",
        "\n",
        "x_train = np.column_stack((x, y, z))\n",
        "x_train_norm = x_train / x_train.max(axis=0)\n",
        "\n",
        "# Train data: Controls\n",
        "u0 = time**(1/3) - np.sin(time)\n",
        "u0 = u0.reshape(-1, 1)\n",
        "\n",
        "u1 = time\n",
        "u1 = u1.reshape(-1, 1)\n",
        "\n",
        "u_train_control = np.column_stack((u0, u1))\n",
        "u_train_norm_control = u_train_control / u_train_control.max(axis=0)\n",
        "\n",
        "\n",
        "# Test data: Controls, Same Control Strategy\n",
        "u0_test = time_test**(1/3) - np.sin(time_test)\n",
        "u0_test = u0_test.reshape(-1, 1)\n",
        "\n",
        "u1_test = time_test\n",
        "u1_test = u1_test.reshape(-1, 1)\n",
        "\n",
        "u_test_control = np.column_stack((u0_test, u1_test))\n",
        "u_test_norm_control = u_test_control / u_test_control.max(axis=0)\n",
        "\n",
        "# Test data: Controls, Diff. Control Strategy\n",
        "u0_test_2 = np.exp(time_test)*np.sin(time_test**4)\n",
        "u0_test_2 = u0_test_2.reshape(-1, 1)\n",
        "\n",
        "u1_test_2 = np.exp(time_test**3)\n",
        "u1_test_2 = u1_test_2.reshape(-1, 1)\n",
        "\n",
        "u_test_control_2 = np.column_stack((u0_test_2, u1_test_2))\n",
        "u_test_norm_control_2 = u_test_control_2 / u_test_control_2.max(axis=0)\n",
        "\n",
        "#_____________________________________________________________________\n",
        "\n",
        "\n",
        "#print(x_train_norm)\n",
        "#print(u_train_control)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pSxyHGPg7Soi"
      },
      "source": [
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8kYkxucOZ8HY"
      },
      "source": [
        "# Train and Test SINDy w/o Control Actions"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Kij0v9Fd7dCn"
      },
      "source": [
        "In this section, the dynamic system will be modeled using SINDy, without any control actions. Therefore, the system evolves freely, and SINDy will identify the ODEs governing the system."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 542
        },
        "id": "X6-RXQkb-C_b",
        "outputId": "b2ae6308-9f4a-4909-bc30-34c6132cff92"
      },
      "outputs": [],
      "source": [
        "# @title View System Trajectory\n",
        "trace = go.Scatter3d(\n",
        "    x=x.flatten(), y=y.flatten(), z=z.flatten(),\n",
        "    mode='lines', name='Trajectory', line=dict(color='blue')\n",
        ")\n",
        "\n",
        "layout = go.Layout(\n",
        "    scene=dict(\n",
        "        bgcolor='black',\n",
        "         xaxis=dict(range=[-5, 5]),\n",
        "         yaxis=dict(range=[-5, 5]),\n",
        "         zaxis=dict(range=[0, 5])),\n",
        "    margin=dict(l=0, r=0, b=0, t=0))\n",
        "\n",
        "fig = go.Figure(data=[trace], layout=layout)\n",
        "fig.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZUgjVsfr8KFk"
      },
      "source": [
        "  ## Train and Test the SINDy model  \n",
        "Pay attention to the `threshold` value as it affects the complexity of the resulting system; the higher this threshold value is, the sparser the result will be."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "4nFERKYen_IP",
        "outputId": "33cb6528-dd3d-433a-8692-450d356a7a4c"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "(x)' = 4.001 y + -2.842 z + 0.479 sin(1 z) + -2.889 cos(1 z)\n",
            "(y)' = -4.000 x + 2.838 z + -0.473 sin(1 z) + 2.886 cos(1 z)\n",
            "(z)' = -0.364 z + 0.672 cos(1 z)\n"
          ]
        }
      ],
      "source": [
        "# @title Train the SINDy model\n",
        "# Train the model\n",
        "\n",
        "\n",
        "identity_lib = IdentityLibrary()\n",
        "poly_lib = PolynomialLibrary(degree=2)\n",
        "fourier_lib = FourierLibrary( include_sin=True, include_cos=True)\n",
        "combined_lib = identity_lib+fourier_lib\n",
        "\n",
        "feature_names = ['x','y','z']\n",
        "opt = ps.STLSQ(threshold=0.2) #0.2\n",
        "\n",
        "model = ps.SINDy(differentiation_method=ps.FiniteDifference(order=2), feature_names=feature_names, feature_library=combined_lib, optimizer=opt)\n",
        "model.fit(x_train, t=dt)\n",
        "model.print()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 542
        },
        "id": "4Z1szVEu669f",
        "outputId": "d80ffcc1-8c7d-4982-a1f3-9b70c022c667"
      },
      "outputs": [],
      "source": [
        "# @title Test the SINDy model\n",
        "# Simulate the model\n",
        "x0_new = x0.flatten()\n",
        "x_simulated = model.simulate(x0_new, t=time_test, integrator='odeint')\n",
        "\n",
        "x_sim = x_simulated[:,0]\n",
        "y_sim = x_simulated[:,1]\n",
        "z_sim = x_simulated[:,2]\n",
        "\n",
        "\n",
        "x0_new_2 = np.array([2, 3, 3])\n",
        "x_simulated_2 = model.simulate(x0_new_2, t=time_test, integrator='odeint')\n",
        "\n",
        "x_sim_2 = x_simulated_2[:,0]\n",
        "y_sim_2 = x_simulated_2[:,1]\n",
        "z_sim_2 = x_simulated_2[:,2]\n",
        "\n",
        "\n",
        "trace = go.Scatter3d(\n",
        "    x=x_sim.flatten(), y=y_sim.flatten(), z=z_sim.flatten(),\n",
        "    mode='lines', name='Sim Trajectory Same I.C.', line=dict(color='red'))\n",
        "\n",
        "layout = go.Layout(\n",
        "    scene=dict(\n",
        "        bgcolor='black',\n",
        "         xaxis=dict(range=[-5, 5]),\n",
        "         yaxis=dict(range=[-5, 5]),\n",
        "         zaxis=dict(range=[0, 5])),\n",
        "    margin=dict(l=0, r=0, b=0, t=0))\n",
        "\n",
        "fig = go.Figure(data=[trace], layout=layout)\n",
        "\n",
        "trace2 = go.Scatter3d(\n",
        "    x=x_sim_2.flatten(), y=y_sim_2.flatten(), z=z_sim_2.flatten(),\n",
        "    mode='lines', name='Sim Trajectory Diff I.C.', line=dict(color='green'))\n",
        "\n",
        "fig = go.Figure(data=[trace], layout=layout)\n",
        "\n",
        "trace1 = go.Scatter3d(\n",
        "    x=x.flatten(), y=y.flatten(), z=z.flatten(),\n",
        "    mode='lines', name='Training Trajectory', line=dict(color='blue'))\n",
        "\n",
        "fig.add_trace(trace2)\n",
        "fig.add_trace(trace1)\n",
        "fig.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "n2LW8uJILPHt"
      },
      "source": [
        "\n",
        "\n",
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8PHdQy8caTuX"
      },
      "source": [
        "# Train and Test SINDy w/ Control Actions and Custom Libraries"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QXMYzR73_3C-"
      },
      "source": [
        "In this section, we will do the same as in the previous section, but now the considered system will have control actions and will not evolve freely. In the SINDy model we want to obtain, it will be necessary to consider these inputs, which will be identified as 'u'."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "id": "i76p76SuPv9K",
        "outputId": "950d256d-21ed-4c63-c614-8db48484b570"
      },
      "outputs": [],
      "source": [
        "# @title View System Trajectory w/ Control\n",
        "trace = go.Scatter3d(\n",
        "    x=x.flatten(), y=y.flatten(), z=z.flatten(),\n",
        "    mode='lines', name='Trajectory', line=dict(color='blue'))\n",
        "\n",
        "layout = go.Layout(\n",
        "    scene=dict(\n",
        "         bgcolor='black',\n",
        "         xaxis=dict(range=[-5, 5]),\n",
        "         yaxis=dict(range=[-5, 5]),\n",
        "         zaxis=dict(range=[ 0, 5])),\n",
        "    margin=dict(l=0, r=0, b=0, t=0))\n",
        "\n",
        "fig = go.Figure(data=[trace], layout=layout)\n",
        "fig.show()\n",
        "\n",
        "plt.style.use('dark_background')\n",
        "\n",
        "fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(10, 6))\n",
        "\n",
        "ax1.plot(time, u0.flatten(), label='Control Signal 0', color='green')\n",
        "ax1.set_ylabel('Control Signal 0')\n",
        "ax1.grid(True)\n",
        "ax1.legend()\n",
        "\n",
        "ax2.plot(time, u1.flatten(), label='Control Signal 1', color='red')\n",
        "ax2.set_xlabel('Time')\n",
        "ax2.set_ylabel('Control Signal 1')\n",
        "ax2.grid(True)\n",
        "ax2.legend()\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "D_snxq1HYU93"
      },
      "source": [
        "## Build a Custom Candidate Library\n",
        "\n",
        "Here, we create a library by composing default SINDy libraries.  \n",
        "\n",
        "This part is not mandatory for the development of the SINDy model with control but is included for completeness. Furthermore, it is possible to generalize more and build a custom library from scratch.  \n",
        "\n",
        "Please note that using very large candidate function libraries tends to make the candidate matrix ill-conditioned or even singular. Therefore, it is advisable not to excessively expand the library unless in specific operational cases.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "id": "DSwW0NtsuVWG"
      },
      "outputs": [],
      "source": [
        "# @title Build a Custom Candidate Library Code\n",
        "# Build a custom library\n",
        "\n",
        "# Call default SINDy libraries\n",
        "poly_lib = PolynomialLibrary(degree=2)\n",
        "fourier_lib = ps.FourierLibrary(include_sin=True, include_cos=True)\n",
        "identity_lib = ps.IdentityLibrary()\n",
        "\n",
        "# Custom lib w/ custom function\n",
        "library_functions = [lambda x: np.tanh(x)]\n",
        "library_function_names = [\n",
        "    lambda x: \"tanh(\" + x + \")\"\n",
        "]\n",
        "tanh_lib = ps.CustomLibrary(\n",
        "    library_functions=library_functions,\n",
        "    function_names=library_function_names\n",
        ")\n",
        "\n",
        "# Sum libraries\n",
        "concat_lib = identity_lib+fourier_lib+poly_lib\n",
        "# Tensor libraries\n",
        "tensor_lib = identity_lib*fourier_lib*poly_lib\n",
        "# Custom library\n",
        "# NOTE: Using huge custom_lib make candidate matrix ill conditioned or singular\n",
        "custom_lib = concat_lib+tensor_lib\n",
        "combined_lib = identity_lib*fourier_lib+identity_lib"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3kkb_0_qPvbq"
      },
      "source": [
        " ## Train and Test the SINDy model w/ Control Actions\n",
        "The difference, compared to the previous code, lies in the availability of vectors with control actions.  \n",
        " `model.fit(x_train, u=u_train_control, t=dt)`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "kDsBxWkv7RKp",
        "outputId": "b2a21b44-f04a-49ce-993c-05bf66d445bf"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "(x)' = 0.756 y sin(1 z) + 0.366 y cos(1 z) + -2.445 z cos(1 z) + -1.511 u0 sin(1 z) + -1.084 u1 cos(1 z) + 3.165 y + -2.683 z + 1.272 u0 + 0.586 u1\n",
            "(y)' = -1.752 z sin(1 z) + 1.965 z cos(1 z) + 0.880 u0 sin(1 z) + 0.733 u1 cos(1 z) + -4.000 x + 4.415 z + -0.741 u0 + -0.396 u1\n",
            "(z)' = -0.352 z sin(1 z) + 0.538 z cos(1 z)\n"
          ]
        }
      ],
      "source": [
        "# @title Train the SINDy model w/ Control\n",
        "# Train the model w/ control inputs\n",
        "opt_contr = ps.STLSQ(threshold=0.32)\n",
        "\n",
        "model = ps.SINDy(feature_library=combined_lib,\n",
        "                 feature_names=['x', 'y', 'z', 'u0', 'u1'],\n",
        "                 optimizer=opt_contr)\n",
        "model.fit(x=x_train, u=u_train_control, t=dt)\n",
        "model.get_feature_names()\n",
        "model.print()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hvI8bOSt55jG"
      },
      "source": [
        "Note that since the control actions are totally invented, the SINDy model finds the low significance and removes them by increasing the sparsity treshold."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "cellView": "form",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "id": "mbA5x_oK7yNy",
        "outputId": "3ac1d66a-a483-4ec6-b807-97c9e58277d9"
      },
      "outputs": [],
      "source": [
        "# @title Test the SINDy model w/ Control\n",
        "# Simulate the model w/ control inputs\n",
        "\n",
        "# Same IC Same Control Strategy\n",
        "x0_new_cont = x0.flatten()\n",
        "x_simulated_cont = model.simulate(x0_new_cont, t=time_test, u=u_test_control, integrator='odeint')\n",
        "\n",
        "x_sim_cont = x_simulated_cont[:,0]\n",
        "y_sim_cont = x_simulated_cont[:,1]\n",
        "z_sim_cont = x_simulated_cont[:,2]\n",
        "\n",
        "# Diff IC Same Control Strategy\n",
        "x0_new_cont_2 = np.array([2, 3, 3])\n",
        "x_simulated_cont_2 = model.simulate(x0_new_cont_2, t=time_test,  u=u_test_control, integrator='odeint')\n",
        "\n",
        "x_sim_cont_2 = x_simulated_cont_2[:,0]\n",
        "y_sim_cont_2 = x_simulated_cont_2[:,1]\n",
        "z_sim_cont_2 = x_simulated_cont_2[:,2]\n",
        "\n",
        "# Diff IC Diff Control Strategy\n",
        "x0_new_cont_3 = np.array([2, 3, 3])\n",
        "x_simulated_cont_3 = model.simulate(x0_new_cont_3, t=time_test,  u=u_test_control_2, integrator='odeint')\n",
        "\n",
        "x_sim_cont_3 = x_simulated_cont_3[:,0]\n",
        "y_sim_cont_3 = x_simulated_cont_3[:,1]\n",
        "z_sim_cont_3 = x_simulated_cont_3[:,2]\n",
        "\n",
        "trace = go.Scatter3d(\n",
        "    x=x_sim_cont.flatten(), y=y_sim_cont.flatten(), z=z_sim_cont.flatten(),\n",
        "    mode='lines', name='Sim Trajectory Same I.C. Same Control', line=dict(color='red')\n",
        ")\n",
        "\n",
        "layout = go.Layout(\n",
        "    scene=dict(\n",
        "        bgcolor='black',\n",
        "         xaxis=dict(range=[-5, 5]),\n",
        "         yaxis=dict(range=[-5, 5]),\n",
        "         zaxis=dict(range=[0, 5])),\n",
        "    margin=dict(l=0, r=0, b=0, t=0))\n",
        "\n",
        "fig1 = go.Figure(data=[trace], layout=layout)\n",
        "\n",
        "trace1 = go.Scatter3d(\n",
        "    x=x.flatten(), y=y.flatten(), z=z.flatten(),\n",
        "    mode='lines', name='Training Trajectory', line=dict(color='blue'))\n",
        "\n",
        "             #(trace) #Sim Trajectory Same I.C. Same Control\n",
        "fig1.add_trace(trace1) #Training Traj\n",
        "fig1.show()\n",
        "\n",
        "\n",
        "# Second 3d plot w/ control actions\n",
        "trace2 = go.Scatter3d(\n",
        "    x=x_sim_cont_2.flatten(), y=y_sim_cont_2.flatten(), z=z_sim_cont_2.flatten(),\n",
        "    mode='lines', name='Sim Trajectory Diff I.C. Same Control', line=dict(color='green'))\n",
        "\n",
        "layout = go.Layout(\n",
        "    scene=dict(\n",
        "        bgcolor='black',\n",
        "         xaxis=dict(range=[-5, 5]),\n",
        "         yaxis=dict(range=[-5, 5]),\n",
        "         zaxis=dict(range=[0, 5])),\n",
        "    margin=dict(l=0, r=0, b=0, t=0))\n",
        "\n",
        "fig2 = go.Figure(data=[trace], layout=layout)\n",
        "\n",
        "trace3 = go.Scatter3d(\n",
        "    x=x_sim_cont_3.flatten(), y=y_sim_cont_3.flatten(), z=z_sim_cont_3.flatten(),\n",
        "    mode='lines', name='Sim Trajectory Diff I.C. Diff Control', line=dict(color='cyan'))\n",
        "\n",
        "fig2 = go.Figure(data=[trace2], layout=layout)\n",
        "\n",
        "             #(trace2) #Sim Trajectory Diff I.C. Same Control\n",
        "fig2.add_trace(trace3) #Sim Trajectory Diff I.C. Diff Control\n",
        "fig2.show()\n",
        "\n",
        "\n",
        "\n",
        "# Plot New Control Strategy\n",
        "plt.style.use('dark_background')\n",
        "\n",
        "fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(10, 6))\n",
        "\n",
        "ax1.plot(time_test, u0_test_2.flatten(), label=' New Control Signal 0', color='green')\n",
        "ax1.set_ylabel('Control Signal 0')\n",
        "ax1.grid(True)\n",
        "ax1.legend()\n",
        "\n",
        "ax2.plot(time_test, u1_test_2.flatten(), label='New Control Signal 1', color='red')\n",
        "ax2.set_xlabel('Time')\n",
        "ax2.set_ylabel('Control Signal 1')\n",
        "ax2.grid(True)\n",
        "ax2.legend()\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "xNlTykIk4r7K"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Notes on Noisy Data"
      ],
      "metadata": {
        "id": "hJXeyertoZvZ"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "There are mainly two ways in PySINDY to try to handle noisy signals:  \n",
        "\n",
        "* The first is to use 'Smooth' differentiation methods `ps.SmoothedFiniteDifference()`, as used follows:\n",
        "\n",
        "```\n",
        "identity_lib = IdentityLibrary()\n",
        "poly_lib = PolynomialLibrary(degree=2)\n",
        "fourier_lib = FourierLibrary( include_sin=True, include_cos=True)\n",
        "combined_lib = identity_lib+fourier_lib\n",
        "\n",
        "feature_names = ['x','y','z']\n",
        "opt = ps.STLSQ(threshold=0.2) #0.2\n",
        "\n",
        "model = ps.SINDy(differentiation_method=ps.SmoothedFiniteDifference(), feature_names=feature_names, feature_library=combined_lib, optimizer=opt)\n",
        "model.fit(x_train, t=dt)\n",
        "model.print()\n",
        "```\n",
        "\n",
        "\n",
        "* The second is to operate with multiple trajectories, i.e., assuming the signal is not cleanable, one proceeds by providing several trajectories to better understand the evolution of the dynamic system despite the noise. as follow:\n",
        "\n",
        "\n",
        "FROM OFFICIAL DOC:\n",
        "\n",
        " https://pysindy.readthedocs.io/en/latest/examples/1_feature_overview/example.html#Multiple-trajectories\n",
        "```\n",
        "if __name__ != \"testing\":\n",
        "    n_trajectories = 20\n",
        "    sample_range = (500, 1500)\n",
        "else:\n",
        "    n_trajectories = 2\n",
        "    sample_range = (10, 15)\n",
        "x0s = np.array([36, 48, 41]) * (np.random.rand(n_trajectories, 3) - 0.5) + np.array(\n",
        "    [0, 0, 25]\n",
        ")\n",
        "x_train_multi = []\n",
        "for i in range(n_trajectories):\n",
        "    x_train_multi.append(\n",
        "        solve_ivp(\n",
        "            lorenz, t_train_span, x0s[i], t_eval=t_train, **integrator_keywords\n",
        "        ).y.T\n",
        "    )\n",
        "\n",
        "model = ps.SINDy()\n",
        "model.fit(x_train_multi, t=dt)\n",
        "model.print()\n",
        "```\n",
        "\n"
      ],
      "metadata": {
        "id": "_anYLU6HqRiB"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "\n",
        "\n",
        "---\n",
        "\n"
      ],
      "metadata": {
        "id": "MhkLLM7F4vib"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Last Revision: G. Oddo 18/01/2024"
      ],
      "metadata": {
        "id": "3G98ViVa4StF"
      }
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "toc_visible": true,
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}