examples/LPJ17e43.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Code to reproduce example 4.3 of\n",
    "\n",
    "B. Legat, R. M. Jungers, and P. A. Parrilo\n",
    "[**Certifying unstability of Switched Systems using Sum of Squares Programming**](https://arxiv.org/abs/1710.01814)\n",
    "\n",
    "This example is taken from example 2.8 of [PJ08] which is inspired from the construction of worst case for the CQLF by [AS98]. The JSR is 1.\n",
    "\n",
    "[PJ08] P. Parrilo and A. Jadbabaie\n",
    "*Approximation of the joint spectral radius using sum of squares.*\n",
    "Linear Algebra and its Applications, Elsevier, **2008**, 428, 2385-2402\n",
    "\n",
    "[AS98] Ando, T. and Shih, M.-h.\n",
    "*Simultaneous Contractibility*.\n",
    "SIAM Journal on Matrix Analysis & Applications, **1998**, 19, 487"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Hybrid System with automaton OneStateAutomaton(2)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using HybridSystems\n",
    "A1 = [1  0\n",
    "      1  0]\n",
    "A2 = [0  1\n",
    "      0 -1]\n",
    "s = discreteswitchedsystem([A1, A2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Choose one SDP solver (two examples are given below, you don't need to run both boxes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Mosek.MosekSolver(Any[(:LOG, 0)])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Mosek\n",
    "solver = MosekSolver(LOG=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using CSDP\n",
    "solver = CSDPSolver(printlevel=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using SwitchOnSafety"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.9999840270350961, 1.4143636144605014)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lb, ub = soslyapb(s, 1, solver=solver, tol=2e-4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see with `hassmp`, the atom extraction has found nothing using feasible dual solution obtained uring the binary search"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "false"
      ],
      "text/plain": [
       "false"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hassmp(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can try it again with the last feasible dual solution found using `sosextractcycle` but as expected it gives the same result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "true"
      ],
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psw = sosextractcycle(s, 1);\n",
    "isnull(psw)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.8408964152537146, 1.0001941594715107)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lb, ub = soslyapb(s, 2, solver=solver, tol=2e-4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see with `hassmp`, the atom extraction has found a proof that the JSR is above 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PSW(1.0, [1])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "getsmp(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can see below that this is obtained by the last dual solution found (e.g. with $γ$ closest to 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PSW(1.0, [1])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psw = sosextractcycle(s, 2);\n",
    "get(psw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that $\\mu_1 = 0.38237 \\delta_{(1, 1)} + 0.38299 \\delta_{(1, -1)}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Atomic measure on the variables DynamicPolynomials.PolyVar{true}[x1, x2] with 2 atoms:\n",
       " at [1.0, 1.0] with weight 0.38237048581016886\n",
       " at [1.0, -1.0] with weight 0.3829874730497317\n"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get(SwitchOnSafety.extractatomic(s, 2, 1, 1e-6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that $\\mu_2 = 0.38299 \\delta_{(1, 1)} + 0.38315 \\delta_{(1, -1)}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Atomic measure on the variables DynamicPolynomials.PolyVar{true}[x1, x2] with 2 atoms:\n",
       " at [1.0, 1.0] with weight 0.38298735054176786\n",
       " at [1.0, -1.0] with weight 0.3831451456900721\n"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "get(SwitchOnSafety.extractatomic(s, 2, 2, 1e-6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.1.0",
   "language": "julia",
   "name": "julia-1.1"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}