droste89 / sostab

A sum-of-squares toolbox for stability analysis

SOStab

A sum-of-squares Matlab toolbox for stability analysis.

The goal of the toolbox is to facilitate the use of SoS programming for calculating an approximation of the region of attraction (RoA) of a dynamical system. The only input needed are an equilibrium point, the admissible set and the dynamics of the system.

SOStab requires YALMIP, as well as a semidefinite solver. Mosek is used by default, but it can be replaced by any other solver, provided they are installed and interfaced through YALMIP.

This Readme illustrates how to use the class with several examples and then describe the properties and the method of the class in the last section.

Installation

Only the SOStab.m file is needed for the toolbox to work. You can put it in your working directory or add it to your Matlab path, to have it available everywhere:

addpath SOStab

Simple example

Initialize the class

Initialize with an equilibrium point $x_{eq}$ and a range $\Delta x$, that defines the admissible set: $[x_{eq}-\Delta x, x_{eq}+\Delta x]$. The two vectors must have the same size. You can change the default values for the solver used and the verbose parameter of the optimization calls.

toy = SOStab([0;0],[1;1]);
%toy.solver = 'mosek';
%toy.verbose = 2;

Define the dynamics of the system

Define $f$ such that $\dot{x} = f(x)$, the vector must have the same size as $x_{eq}$. The property x of the class is a YALMIP sdpvar, which is a symbolic representation of the state variable. It has the same size as the problem. The property t is also an sdpvar, representing time that can be useful here. For example, here we implement a simple radial system:

toy.dynamics = (toy.x'*toy.x - 0.25)*toy.x;

Solve the outer approximation problem

Solve the outer approximation of the finite horizon RoA for a target set $\lbrace x \in \mathbb{R}^n : ||x-x_{eq}||< \varepsilon \rbrace$, with $\varepsilon=0.1$ and a time horizon $T=20$. The degree of the polynomials is defined by $d=6$. $d$ above 20 is useless and for high dimensions ($n > 6$), a standard computer RAM can't cope with $d$ above 10.

toy.SoS_out(6,20,0.1);

Solve the inner approximation problem

Solve the inner approximation of the RoA. Note that the results of this approximation are very sensitive to the relative size of the sets and the degree of the polynomials. You can check the relevance of the results by verifying the coefficients of $w$. If they are null for most of them, then this means that the solver was not able to find the optimum of the problem. You can try to increase the size of the target set or the degree of the polyomials.

toy.SoS_in(6,20,0.1);
witoy = value(toy.wcoef_inner);
witoy(1) = witoy(1)-1;
if witoy'*witoy > 0.00001
    toy.plot_roa(1,2,'inner',0,"x","y");
else
    disp("w not relevant")
end

Plot the RoA

Plot the projection of the RoA in two dimensions (for the last calculated approximation).

toy.plot_roa(1,2,'outer');
axis('equal');

The first two arguments indicate the indices of the variables to plot. The fourth (optionnal) argument indicates to plot the target set.

Plot the surface of $w$

Plot the surface of $w$ for the last calculated inner or outer approximation.

toy.plot_w(1,2, 'o');

Various other examples

Van der Pol

The Van der Pol oscillator is a simple 2 dimensional system:

$$ \begin{pmatrix} \dot{x}_1 \\ \dot{x}_2 \end{pmatrix} = \begin{pmatrix} -2x_2 \\ 0.8x_1 + 10(x_1^2-0.21)x_2 \end{pmatrix} $$

vdp = SOStab([0;0],[1.1;1.1]);
vdp.dynamics = [-2*vdp.x(2); 0.8*vdp.x(1) + 10*(vdp.x(1)^2-0.21)*vdp.x(2)];
vdp.SoS_out(12,1,0.5);
vovdp = value(vdp.vcoef_outer);
wovdp = value(vdp.wcoef_outer);
vdp.plot_roa(1,2,'outer',1,"x_1","x_2");
vdp.SoS_in(12,1,0.5);
vivdp = value(vdp.vcoef_inner);
wivdp = value(vdp.wcoef_inner);
wivdp(1) = wivdp(1)-1;
if wivdp'*wivdp>0.00001
    vdp.plot_roa(1,2,'inner',0,"x_1","x_2");
else
    disp("w not relevant")
end

Scaled pendulum

The pendulum is a simple example of "polynomialization" of a system.

$$\dot{\theta} = \sin(2\theta)=2\sin\theta\cos\theta $$

The system used as input of the toolbox will be $x= (\sin\theta, \cos\theta, \omega)$ (here $\omega=0$ is used for plotting in 2 dimensions), then we have:

$$\dot{x}= \begin{pmatrix}-2\sin\theta\cos^2\theta \\ 2\sin^2\theta\cos\theta \\ 0 \end{pmatrix}$$

Hence, the dynamic is indeed polynomial in the variables.

pen = SOStab([0;1;0],[1;1;1],[1,2]);
pen.dynamics = [-2*pen.x(1)*pen.x(2)^2;2*pen.x(1)^2*pen.x(2);0];
pen.SoS_out(6,40,0.5);
vopen = value(pen.vcoef_outer);
wopen = value(pen.wcoef_outer);
pen.plot_roa([1,2,1],3,'outer',1,"\theta","\omega");
pen.SoS_in(6,40,0.5);
vipen = value(pen.vcoef_inner);
wipen = value(pen.wcoef_inner);
wipen(1) = wipen(1)-1;
if wipen'*wipen > 0.00001
    pen.plot_roa([1,2],3,'inner',0,"\theta","\omega");
else
    disp("w not relevant")
end

Power system

Consider an electrical power network made of synchronous machines connected in a cycle, along with its second reduced order model (see M. Anghel, F. Milano, and A. Papachristodoulou, "Algorithmic construction of Lyapunov functions for power system stability analysis"):

$$ \begin{align*} \dot \theta_1 &= \omega_1, \qquad \dot \theta_2 = \omega_2 \\ \dot \omega_1 &= -\sin\theta_1 - 0.5\sin(\theta_1-\theta_2) - 0.4 , \omega_1 \\ \dot \omega_2 &= -0.5\sin\theta_2 - 0.5\sin(\theta_2-\theta_1) - 0.5 , \omega_2 + 0.05 \end{align*} $$

Similarly to the pendulum, it can be polynomialized we the change of variable $x= (\sin\theta_1, \cos\theta_1, \sin\theta_2, \cos\theta_2,\omega_1, \omega_2)$. The RoA of the system is calculated by:

eq = [sin(0.02);cos(0.02);sin(0.06);cos(0.06);0;0];
dev = [1;1;1;1;pi;pi];
ang_ind = [1,2;3,4];
pow_sys = SOStab(eq,dev,ang_ind);
pow_sys.dynamics = [pow_sys.x(5)*pow_sys.x(2);-pow_sys.x(5)*pow_sys.x(1);...
    pow_sys.x(6)*pow_sys.x(4);-pow_sys.x(6)*pow_sys.x(3);
    -pow_sys.x(1)-0.5*(pow_sys.x(1)*pow_sys.x(4)-pow_sys.x(2)*pow_sys.x(3))-0.4*pow_sys.x(5);
    -0.5*pow_sys.x(3)+0.5*(pow_sys.x(1)*pow_sys.x(4)-pow_sys.x(2)*pow_sys.x(3))-0.5*pow_sys.x(6)+0.05];
pow_sys.SoS_out(6,8,0.1);
vopow = value(pow_sys.vcoef_outer);
wopow = value(pow_sys.wcoef_outer);
wopow(1) = wopow(1)-1;
if wopow'*wopow > 0.00001
    pow_sys.plot_roa([1,2],[3,4],'o','1',"\theta_1","\theta_2");
else
    disp("w not relevant")

Properties and methods of the class

Properties of the class are the following:

• dimension: dimension of the problem (number of variables)
• x_eq: first argument of the call, the equilibrium state of the system
• delta_x: the range around the equilibrium defining the feasible set
• angle_eq: the vector $\theta_{eq}$ - the equilibrium angles - recalculated from the values of their sine and cosine (empty if no phases are involved)
• angle_ind: the indices of sine and cosine functions in the recast variable $x$, given as last input of SOStab call (empty if no phases are involved)
• x: a YALMIP sdpvar polynomial object, of the dimension of the problem. It represents the variable $x$ and is used by the user to define the dynamics of the system
• t: sdpvar polynomial of size 1, representing the time variable, which can be needed to define the dynamics of the system (if non-autonomous)
• D, the matrix $D$ of the variable change used in the toolbox to normalize the system
• invD, the inverse $D^{-1}$ of the matrix $D$
• solver, the choice of the solver used in the optimization, defined as Mosek by default
• verbose, the value of the verbose parameters of the YALMIP optimization calls, defined at 2 by default
• dynamics, a YALMIP polynomial defining the polynomial dynamics $f$ of the system

The class has seven methods:

• The class initialization SOStab takes as input the equilibrium point $x_{eq}$, the distance $\Delta x$ as two vectors, that must have the same length -- the dimension of the problem. It initializes an object of the class SOStab.

• moments is a method called inside the other methods, which calculates the integration of all the monomials on the normalized feasible set $[-1,1]^n$. The function has only one argument, the maximum degree of the monomials. It returns the vector $y$ such that the integration of a polynomial $p$ on the admissible set $X$ would be $\int_X p = y^\top v_p$ where $v_p$ is the vector of coefficients of $p$ in the basis of monomials (with Yalmip ordering convention if $n\geq 2$).

• SoS_out takes as input an optional matrix $A$ (by default $A = I_n$) and a real positive value $\varepsilon$ - defining the target $\mathbb{M} = \lbrace x \in \mathbb{R}^n : | A(x-x_{eq})|^2 \leq \varepsilon^2 \rbrace$ -, an even degree $d=2k$ for the polynomials and a time horizon $T$. It solves the strengthening of the optimization problem giving the RoA in bounded degree $d$ and returns its optimal value (an over approximation of the volume of the calculated RoA) and the coefficients of the polynomials $v^\star_k$ and $w^\star_k$ defining this approximation.

• SoS_in takes the same inputs and solves the counterpart of the previous problem, which is the inner approximation optimization problem. It returns the same outputs as the previous method. Currently, the inner approximation has some limitations on some of the instances tested, due to the algebraic representation of the boundary of the admissible state.

• plot_roa takes as inputs the two indices $i,j$ of the variables on which to project the ROA, a string which states which approximation to plot (outer or inner) and four optional arguments: an int indicating whether or not to plot the target set (1 or 0) - 0 by default -, two strings indicating the name of the axis (respectively for the first and second indices) - $x_i$ and $x_j$ by default - and a vector giving the size of the plotting mesh - $(40,40)$ by default. It plots the expected slice of the ROA, with all other variables at equilibrium.

  If both inner and outer approximations are called sequentially for the same variables, the two plots will appear on the same figure.

  If an angle $\theta$ is one of the plotted variable, the index argument should be a vector of 2 values: the index of $\sin\theta$ in the variables and the one of $\cos\theta$.

• plot_v and plot_w take the same inputs. They respectively plot the graph of $v^\star_k$ and $w^\star_k$ in 3D (with non-selected variables at equilibrium). \end{itemize}

sostab's Issues

Moments

Au lieu de définir y par une suite de 1 et réallouer la valeur 0 à tous les moments avec une puissance impaire, j'aurais plutôt définit y par une suite de 0, rien fait dès qu'il y a un moment avec une puissance impaire, et modifié la valeur de y seulement aux puissances paires, non ?

Code improvement

• Transform phi into x (unclear)
• Define a null dynamic inside the class creation
• separate properties between modifiable and non modifiable.
• Maybe add the angles to the initial call.
• Check angle indices are different

changement de taille des tableaux

Pour les contraintes SoS, Matlab me dit dans l'oreillette que changer la taille d'un tableau en cours de calcul est très cher. À la place, il me propose de créer dès le départ un tableau de la bonne taille (avec des 1 ou des 0) et d'en modifier les entrées dans les boucles. Ça peut être bien à faire pour gagner en efficacité

Separation sdpvar declaration & resolution

both methods SoS_out and SoS_in define the same objects (e.g. f, [-1,1]^n, the target set...). Very often both methods are called one after the other. It would be better to separate the definition of the variables that are common to the 2 methods, into a new method (e.g. obj.SoS_var), so that we do not repeat all variable declarations at each call of SoS_in and SoS_out. This way we will save computational time.

Issue in Initialization

Invalid default value for property 't' in class 'SOStab':
Undefined function 'sdpvar' for input arguments of type 'double'.

cas dimension 1

Que fait la toolbox si le système initial est de dimension 1?

run demo.mlx on R2019b version, got ERROR

run demo with R2019b version, got:

ERROR: variable_replace (line 31)
Nonlinear replacement not supported yet

ERROR:  sdpvar/replace (line 38)
    Z = variable_replace(X,Y,W);

ERROR:  SOStab/SoS_out (line 188)
				f = (obj.D) * replace(obj.dynamics, [obj.t;obj.x], [s*T;(obj.invD)*z+obj.x_eq]); % f defined thanks to variable change

ROA of SMIB system

We have tested this on a SMIB system, but something seems to go wrong. One with polynomialization returns ROA with volume = 0 and one with taylor expansion returns volume = 1.7 but both of them cannot plot the ROA. We are wondering why this will happen and codes are attached. Thank you for your patience and time if you could help us with this.
`clear
syms w theta
syms z1 z2
H = 0.5;
D = 0.1;
X = 0.1;
E = 1;
V = 1;
Pm = 0.9;
f = [
w;
(Pm - E * V * sin(theta) / X - D * w) / 2 / H;
];
x = [w; theta];
x0 = [0; asin(X * Pm)];
% sinθ = z1, cosθ = z2
f_tayl = taylor(f, x, x0, 'Order', 4);
% x_recast = [w; z1; z2];
% x0_recast = [0; sin(x0(2)); cos(x0(2))];

eq = x0;
dev = [1; 1];
% ang_ind = [2, 3];
SMIB = SOStab(eq, dev);
sysDyn = sym2SOStab('SMIB', f_tayl, x);
SMIB.dynamics = eval(sysDyn);
[vol, vc, wc] = SMIB.SoS_out(4,5,0.1);
vopow = value(SMIB.vcoef_outer);
wopow = value(SMIB.wcoef_outer);
SMIB.plot_roa(1, 2, 'outer');`

`clear
syms w theta
syms z1 z2
H = 0.5;
D = 0.1;
X = 0.1;
E = 1;
V = 1;
Pm = 0.9;
f = [
w;
(Pm - E * V * sin(theta) / X - D * w) / 2 / H;
];
x = [w; theta];
x0 = [0; asin(X * Pm)];
% sinθ = z1, cosθ = z2
f_recast = [
w * z2;
w * (-z1);
(Pm - E * V * z1 - D * w) / 2 / H;
];
x_recast = [w; z1; z2];
x0_recast = [0; sin(x0(2)); cos(x0(2))];

eq = x0_recast;
dev = [1; 1; 1;];
ang_ind = [2, 3];
SMIB = SOStab(eq, dev, ang_ind);
sysDyn = sym2SOStab('SMIB', f_recast, x_recast);
SMIB.dynamics = eval(sysDyn);
[vol, vc, wc] = SMIB.SoS_out(10,5,0.1);
vopow = value(SMIB.vcoef_outer);
wopow = value(SMIB.wcoef_outer);
SMIB.plot_roa(1, [2, 3], 'outer');`

function sysDyn = sym2SOStab(sysName, f, x) if ~isa(sysName, 'string') sysName = string(sysName); end f_str = string(f); x_str = string(x); [a, b] = size(x); nx = numel(x); dots = repmat(".", a, b); sysPre = join([repmat(sysName, a, b), dots], ""); % xPre = "x"; xNew_str = join([repmat("x(", a, b), join([string((1 : nx).'), repmat(")", a, b)]) ], ""); f_str_new = replace(f_str, x_str, join([sysPre, xNew_str], "")); sysDyn = f_str_new; sysDyn = join([sysDyn, repmat(";", nx, 1)]); sysDyn = join(sysDyn.'); sysDyn = join(["[", sysDyn, "]"]); end

Arguments à passer en option

La plupart du temps on définit A=eye(n). Il faudrait prendre ça par défaut et rendre l'input optionnel

Angular input minimization

• Remove the first column of the parsing matrix, and simply reconstruct the angular equilibria using the formula \theta = \sgn(\sin\theta)·\arccos(\cos\theta).
• Remove the last index in the plot_ methods: giving only the indices of the sin and cos allows to retrieve (using columns 2 and 3 of the parsing matrix) the index of the corresponding angle, and access to its equilibrium value using the above constructed angle_eq list.

BONUS: currently the default legend is "x_i" for the plot. Maybe we could add an if loop so that if the variable is an angle (which is easily identified by counting the index size in the input) so that the default legend for an angular value is "\theta_i" with i determined in the above instructions. But this is only cosmetics and can wait after paper submission.

Polir les arguments de SOStab

Je pense que c'est mieux si les chaînes de caractères nommant les variables n'apparaissent pas dans le premier appel de SOStab. L'idéal serait de ne les mettre que lorsqu'elles sont utiles (je te laisse voir où exactement, on peut s'appeler pour en discuter si tu veux)

Optional arguments

Clarify the use of optional arguments (in the plotting methods) with an input parser for example.

demo.mlx

Change notations of the notebook to match the paper. Include all 4 systems of the demo.m file (radial system, Van der Pol, pendulum, power system)

Create a proper README

• Define all class parameters and methods.
• Give an example.

Issue in simple Swing Equation - need help and some suggestions

First, really thanks for sharing this. However I find something that seems need improvements.
So I tried on Swing Equation, which should give me a similar result as using this method (Phase Portrait Plotter on 2D phase plane
[Version 1.0.3]
(https://www.mathworks.com/matlabcentral/fileexchange/110785-phase-portrait-plotter-on-2d-phase-plane#version_history_tab) (370 KB) by [Yu Zhang]

The codes I tried using your method are attached below (note some \times "*" operators are missing):

clf;
eq = [sin(pi/6);cos(pi/6);0];
dev = [1;1;4];
ang_ind = [1,2];
pow_sys = SOStab(eq,dev,ang_ind);
pow_sys.dynamics = [pow_sys.x(3)*pow_sys.x(2);-pow_sys.x(3)pow_sys.x(1);
0.5 - pow_sys.x(1) - 0.2pow_sys.x(3)
];

pow_sys.SoS_out(6,14,0.1);
vopow_sys = value(pow_sys.vcoef_outer);
wopow_sys = value(pow_sys.wcoef_outer);
pow_sys.plot_roa([1,2],3,'outer',0,"\theta","\omega");

The result looks like this, partially correct:

My comments are:

1. It seems your method is sensitive to the settings of D, T, and error (varepsilon). I tried difeerent settings but none of them give expected results and some of them seems not correct. Although maybe I'm wrong, could someone help on this?
2. You should remove the requirement on the initial acknowledgement of equilibrium point values. Refer to the method above I referred to.
3. The good thing of your method is it does NOT limited to the degree of 2 of ODE fcn.