local stability radius well defined. More examples to show

dynamicslab · Apr 19, 2023 · c68b873 · c68b873
1 parent 902bf57
commit c68b873
Showing 1 changed file with 36 additions and 16 deletions.
diff --git a/examples/17_trapping_inequality.ipynb b/examples/17_trapping_inequality.ipynb
@@ -18,7 +18,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 19,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -30,15 +30,16 @@
     "import matplotlib.gridspec as gridspec\n",
     "import scipy.io as sio\n",
     "from pysindy.utils import lorenz\n",
-    "\n",
+    "import sympy as sp\n",
+    "from sympy import limit, Symbol\n",
     "# ignore warnings\n",
     "import warnings\n",
     "warnings.filterwarnings(\"ignore\")"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 47,
+   "execution_count": 20,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -110,6 +111,21 @@
     "    print('Estimate of trapping region size, Rm = ', Rm)\n",
     "    print('Normalized trapping region size, Reff = ', Reff)\n",
     "\n",
+    "def get_trapping_radius(max_eigval, eps_Q, r, d):\n",
+    "    x = Symbol('x')\n",
+    "    delta = max_eigval ** 2 - 4 * np.sqrt(r ** 3) * eps_Q * np.linalg.norm(d, 2) / 3\n",
+    "    delta_func = max_eigval ** 2 - 4 * np.sqrt(r ** 3) * x * np.linalg.norm(d, 2) / 3\n",
+    "    if delta < 0:\n",
+    "        rad_trap = 0\n",
+    "        rad_stab = 0\n",
+    "    else:\n",
+    "        y_trap = -(3 / (2 * np.sqrt(r ** 3) * x)) * (max_eigval + sp.sqrt(delta_func))\n",
+    "        y_stab = (3 / (2 * np.sqrt(r ** 3) * x)) * (-max_eigval + sp.sqrt(delta_func))\n",
+    "        rad_trap = limit(y_trap, x, eps_Q, dir='+')\n",
+    "        rad_stab = limit(y_stab, x, eps_Q, dir='+')\n",
+    "    return rad_trap, rad_stab\n",
+    "\n",
+    "\n",
     "def check_stability_new(r, Xi, mod_matrix, sindy_opt, mean_val):\n",
     "    opt_m = sindy_opt.m_history_[-1]\n",
     "    PC_tensor = sindy_opt.PC_\n",
@@ -132,16 +148,13 @@
     "    d = C + np.dot(L, opt_m) + np.dot(np.tensordot(Q, opt_m, axes=([2], [0])), opt_m)\n",
     "    d = mod_matrix @ d\n",
     "    eps_Q = np.max(np.abs(Q_sum))\n",
-    "    delta = max_eigval ** 2 - 4 * np.sqrt(r ** 3) * eps_Q * np.linalg.norm(d, 2) / 3\n",
-    "    if delta < 0:\n",
-    "        Rm = 0\n",
-    "    else:\n",
-    "        Rm = (3 / (2 * np.sqrt(r ** 2) * eps_Q)) * (-max_eigval + np.sqrt(delta))\n",
+    "    Rm, DA = get_trapping_radius(max_eigval, eps_Q, r, d)\n",
     "    Reff = Rm / mean_val\n",
     "    print('Estimate of trapping region size, Rm = ', Rm)\n",
     "    print('Normalized trapping region size, Reff = ', Reff)\n",
+    "    print('Local stability size, DA = ', DA)\n",
     "\n",
-    "    return Rm\n",
+    "    return Rm, DA\n",
     "\n",
     "\n",
     "# use optimal m, calculate and plot the stability radius when the third-order\n",
@@ -340,7 +353,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 33,
+   "execution_count": 21,
    "metadata": {},
    "outputs": [
     {
@@ -413,13 +426,12 @@
     "PQ_tensor = sindy_opt.PQ_\n",
     "Lenergy = np.tensordot(PL_tensor, Xi, axes=([3, 2], [0, 1]))\n",
     "Q = np.tensordot(PQ_tensor, Xi, axes=([4, 3], [0, 1]))\n",
-    "# Q = np.tensordot(sindy_opt.PQ_, Xi, axes=([4, 3], [0, 1]))\n",
     "print(np.max(np.abs((Q + np.transpose(Q, [1, 2, 0]) + np.transpose(Q, [2, 0, 1])))))"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 34,
+   "execution_count": 22,
    "metadata": {},
    "outputs": [
     {
@@ -453,7 +465,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 37,
+   "execution_count": 23,
    "metadata": {
     "scrolled": false
    },
@@ -487,6 +499,11 @@
       "Estimate of trapping region size, Rm =  105.19082697129966\n",
       "Normalized trapping region size, Reff =  4.4367936589400925\n",
       "Frobenius error =  0.7269189253531968\n",
+      "optimal m:  [-1.1598919  -0.12244458 36.74564581]\n",
+      "As eigvals:  [-10.0141668   -2.66230669  -0.93838016]\n",
+      "Estimate of trapping region size, Rm =  105.206801934060\n",
+      "Normalized trapping region size, Reff =  4.43746745926581\n",
+      "Local stability size, DA =  25812465313688.6\n",
       "Frobenius coefficient error =  1.4989984052599905\n",
       "Time-averaged derivative error =  1.0832957179907218e-05\n"
      ]
@@ -527,6 +544,7 @@
     "check_stability(r, Xi, np.eye(r), sindy_opt, mean_val)\n",
     "E_pred = np.linalg.norm(x_test - x_test_pred) / np.linalg.norm(x_test)\n",
     "print('Frobenius error = ', E_pred)\n",
+    "check_stability_new(r, Xi, np.eye(r), sindy_opt, mean_val)\n",
     "\n",
     "# compute relative Frobenius error in the model coefficients\n",
     "sigma = 10\n",
@@ -610,7 +628,9 @@
     "This problem is solved with a convex relaxation. Below, we relax the hard constraint slightly and instead solve \n",
     "$$ argmin_{\\mathbf{\\xi},\\mathbf m}\\|\\dot{\\mathbf a} - \\mathbf \\Theta(\\mathbf a) \\mathbf{\\xi}\\|^2 + \\gamma R(\\mathbf \\xi)  + \\eta \\lambda_1(\\mathbf A)  \\quad s.t. \\quad -\\epsilon_Q \\leq Q_{ijk} + Q_{jik} + Q_{kji} \\leq \\epsilon_Q.$$ \n",
     "This allows us to build locally Lyapunov stable models, and adjust the size of the local stability radius by varying $\\epsilon_Q$. A conservative estimate of the local stability is:\n",
-    "$$\\rho_+ = \\frac{3}{2r^{\\frac{3}{2}}\\epsilon_Q} \\left( \\sqrt{\\lambda^2_{\\text{max}}(\\textbf{A}_S) - \\frac{4r^{\\frac{3}{2}}\\epsilon_Q}{3}\\|\\mathbf{d}\\|_2} - \\lambda_{\\text{max}}(\\textbf{A}_S) \\right).$$"
+    "$$\\rho_+ = \\frac{3}{2r^{\\frac{3}{2}}\\epsilon_Q} \\left( \\sqrt{\\lambda^2_{\\text{max}}(\\textbf{A}_S) - \\frac{4r^{\\frac{3}{2}}\\epsilon_Q}{3}\\|\\mathbf{d}\\|_2} - \\lambda_{\\text{max}}(\\textbf{A}_S) \\right).$$\n",
+    "And the radius of the trapping region is given by:\n",
+    "$$\\rho_- = -\\frac{3}{2r^{\\frac{3}{2}}\\epsilon_Q} \\left( \\sqrt{\\lambda^2_{\\text{max}}(\\textbf{A}_S) - \\frac{4r^{\\frac{3}{2}}\\epsilon_Q}{3}\\|\\mathbf{d}\\|_2} + \\lambda_{\\text{max}}(\\textbf{A}_S) \\right).$$"
    ]
   },
   {
@@ -1413,7 +1433,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "sindy",
    "language": "python",
    "name": "python3"
   },
@@ -1444,7 +1464,7 @@
   },
   "vscode": {
    "interpreter": {
-    "hash": "f62b9ab411974d02b25bf2f89f4d54b9e7a8fbbc53d424ffba0504d4f9fff406"
+    "hash": "988acc18e7f5d05533957e37a841110745c8a9742da3b198a1f620e63121f7f6"
    }
   }
  },