{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Lyapunov Function Search"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Adapted from**: SOSTOOLS' SOSDEMO2 (See Section 4.2 of [SOSTOOLS User's Manual](http://sysos.eng.ox.ac.uk/sostools/sostools.pdf))"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(DynamicPolynomials.PolyVar{true}[x₁, x₂, x₃],)"
     },
     "metadata": {},
     "execution_count": 1
    }
   ],
   "cell_type": "code",
   "source": [
    "using DynamicPolynomials\n",
    "@polyvar x[1:3]"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "We define below the vector field $\\text{d}x/\\text{d}t = f$"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "3-element Vector{MultivariatePolynomials.RationalPoly{DynamicPolynomials.Polynomial{true, Int64}, DynamicPolynomials.Polynomial{true, Int64}}}:\n (-x₁³ - x₁x₃²) / (1)\n (-x₁²x₂ - x₂) / (1)\n (3x₁²x₃³ + 3x₁²x₃ - x₃³ - 4x₃) / (x₃² + 1)"
     },
     "metadata": {},
     "execution_count": 2
    }
   ],
   "cell_type": "code",
   "source": [
    "f = [-x[1]^3 - x[1] * x[3]^2,\n",
    "     -x[2] - x[1]^2 * x[2],\n",
    "     -x[3] - 3x[3] / (x[3]^2 + 1) + 3x[1]^2 * x[3]]"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "We need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v0.21.6/installation/#Supported-solvers) for a list of the available choices.\n",
    "We use `SOSModel` instead of `Model` to be able to use the `>=` syntax for Sum-of-Squares constraints."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using SumOfSquares\n",
    "using CSDP\n",
    "solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)\n",
    "model = SOSModel(solver);"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "We are searching for a Lyapunov function $V(x)$ with monomials $x_1^2$, $x_2^2$ and $x_3^2$.\n",
    "We first define the monomials to be used for the Lyapunov function:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "3-element Vector{DynamicPolynomials.Monomial{true}}:\n x₁²\n x₂²\n x₃²"
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "cell_type": "code",
   "source": [
    "monos = x.^2"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "We now define the Lyapunov function as a polynomial decision variable with these monomials:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(_[1])x₁² + (_[2])x₂² + (_[3])x₃²",
      "text/latex": "$$ (_[1])x_{1}^{2} + (_[2])x_{2}^{2} + (_[3])x_{3}^{2} $$"
     },
     "metadata": {},
     "execution_count": 5
    }
   ],
   "cell_type": "code",
   "source": [
    "@variable(model, V, Poly(monos))"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "We need to make sure that the Lyapunov function is strictly positive.\n",
    "We can do this with a constraint $V(x) \\ge \\epsilon (x_1^2 + x_2^2 + x_3^2)$,\n",
    "let's pick $\\epsilon = 1$:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(_[1] - 1)x₁² + (_[2] - 1)x₂² + (_[3] - 1)x₃² is SOS",
      "text/latex": "$$ (_[1] - 1)x_{1}^{2} + (_[2] - 1)x_{2}^{2} + (_[3] - 1)x_{3}^{2} \\text{ is SOS} $$"
     },
     "metadata": {},
     "execution_count": 6
    }
   ],
   "cell_type": "code",
   "source": [
    "@constraint(model, V >= sum(x.^2))"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "cell_type": "markdown",
   "source": [
    "We now compute $\\text{d}V/\\text{d}x \\cdot f$."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "((-2 _[1])x₁⁴x₃² + (-2 _[2])x₁²x₂²x₃² + (-2 _[1] + 6 _[3])x₁²x₃⁴ + (-2 _[1])x₁⁴ + (-2 _[2])x₁²x₂² + (-2 _[1] + 6 _[3])x₁²x₃² + (-2 _[2])x₂²x₃² + (-2 _[3])x₃⁴ + (-2 _[2])x₂² + (-8 _[3])x₃²) / (x₃² + 1)",
      "text/latex": "$$ \\frac{(-2 _[1])x_{1}^{4}x_{3}^{2} + (-2 _[2])x_{1}^{2}x_{2}^{2}x_{3}^{2} + (-2 _[1] + 6 _[3])x_{1}^{2}x_{3}^{4} + (-2 _[1])x_{1}^{4} + (-2 _[2])x_{1}^{2}x_{2}^{2} + (-2 _[1] + 6 _[3])x_{1}^{2}x_{3}^{2} + (-2 _[2])x_{2}^{2}x_{3}^{2} + (-2 _[3])x_{3}^{4} + (-2 _[2])x_{2}^{2} + (-8 _[3])x_{3}^{2}}{x_{3}^{2} + 1} $$"
     },
     "metadata": {},
     "execution_count": 7
    }
   ],
   "cell_type": "code",
   "source": [
    "using LinearAlgebra # Needed for `dot`\n",
    "dVdt = dot(differentiate(V, x), f)"
   ],
   "metadata": {},
   "execution_count": 7
  },
  {
   "cell_type": "markdown",
   "source": [
    "The denominator is $x[3]^2 + 1$ is strictly positive so the sign of `dVdt` is the\n",
    "same as the sign of its numerator."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(-2 _[1])x₁⁴x₃² + (-2 _[2])x₁²x₂²x₃² + (-2 _[1] + 6 _[3])x₁²x₃⁴ + (-2 _[1])x₁⁴ + (-2 _[2])x₁²x₂² + (-2 _[1] + 6 _[3])x₁²x₃² + (-2 _[2])x₂²x₃² + (-2 _[3])x₃⁴ + (-2 _[2])x₂² + (-8 _[3])x₃²",
      "text/latex": "$$ (-2 _[1])x_{1}^{4}x_{3}^{2} + (-2 _[2])x_{1}^{2}x_{2}^{2}x_{3}^{2} + (-2 _[1] + 6 _[3])x_{1}^{2}x_{3}^{4} + (-2 _[1])x_{1}^{4} + (-2 _[2])x_{1}^{2}x_{2}^{2} + (-2 _[1] + 6 _[3])x_{1}^{2}x_{3}^{2} + (-2 _[2])x_{2}^{2}x_{3}^{2} + (-2 _[3])x_{3}^{4} + (-2 _[2])x_{2}^{2} + (-8 _[3])x_{3}^{2} $$"
     },
     "metadata": {},
     "execution_count": 8
    }
   ],
   "cell_type": "code",
   "source": [
    "P = dVdt.num"
   ],
   "metadata": {},
   "execution_count": 8
  },
  {
   "cell_type": "markdown",
   "source": [
    "Hence, we constrain this numerator to be nonnegative:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(2 _[1])x₁⁴x₃² + (2 _[2])x₁²x₂²x₃² + (2 _[1] - 6 _[3])x₁²x₃⁴ + (2 _[1])x₁⁴ + (2 _[2])x₁²x₂² + (2 _[1] - 6 _[3])x₁²x₃² + (2 _[2])x₂²x₃² + (2 _[3])x₃⁴ + (2 _[2])x₂² + (8 _[3])x₃² is SOS",
      "text/latex": "$$ (2 _[1])x_{1}^{4}x_{3}^{2} + (2 _[2])x_{1}^{2}x_{2}^{2}x_{3}^{2} + (2 _[1] - 6 _[3])x_{1}^{2}x_{3}^{4} + (2 _[1])x_{1}^{4} + (2 _[2])x_{1}^{2}x_{2}^{2} + (2 _[1] - 6 _[3])x_{1}^{2}x_{3}^{2} + (2 _[2])x_{2}^{2}x_{3}^{2} + (2 _[3])x_{3}^{4} + (2 _[2])x_{2}^{2} + (8 _[3])x_{3}^{2} \\text{ is SOS} $$"
     },
     "metadata": {},
     "execution_count": 9
    }
   ],
   "cell_type": "code",
   "source": [
    "@constraint(model, P <= 0)"
   ],
   "metadata": {},
   "execution_count": 9
  },
  {
   "cell_type": "markdown",
   "source": [
    "The model is ready to be optimized by the solver:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "JuMP.optimize!(model)"
   ],
   "metadata": {},
   "execution_count": 10
  },
  {
   "cell_type": "markdown",
   "source": [
    "We verify that the solver has found a feasible solution:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "FEASIBLE_POINT::ResultStatusCode = 1"
     },
     "metadata": {},
     "execution_count": 11
    }
   ],
   "cell_type": "code",
   "source": [
    "JuMP.primal_status(model)"
   ],
   "metadata": {},
   "execution_count": 11
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can now obtain this feasible solution with:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "15.718362431619164x₁² + 12.28610732110516x₂² + 2.995845325502817x₃²",
      "text/latex": "$$ 15.718362431619164x_{1}^{2} + 12.28610732110516x_{2}^{2} + 2.995845325502817x_{3}^{2} $$"
     },
     "metadata": {},
     "execution_count": 12
    }
   ],
   "cell_type": "code",
   "source": [
    "value(V)"
   ],
   "metadata": {},
   "execution_count": 12
  },
  {
   "cell_type": "markdown",
   "source": [
    "---\n",
    "\n",
    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.2"
  },
  "kernelspec": {
   "name": "julia-1.7",
   "display_name": "Julia 1.7.2",
   "language": "julia"
  }
 },
 "nbformat": 4
}