SOSTools - Getting example to work?

EDIT: Changed example to a more suitable one for discussion (smaller):

I am trying to solve example 4.2 from this paper: Analysis of Non-polynomial Systems using the Sum of Squares Decomposition. Here the system

$$\begin{split} \dot{x}_1 &= x_2 \\ \dot{x}_2 &= \sin(x_1)\cos(x_1) - 10\sin(x_1) \\ \end{split}$$

is investigated (i.e. a polynomial Lyapunov function is searched for). Note that the physical parameters are chosen as suggested in the paper (all parameters are set to $1$ except $g = 10$. The system is first rewritten like

$$\begin{split} \dot{x}_1 &= x_2 \\ \dot{x}_2 &= u_1u_2 - 10u_1 \\ \dot{u}_1 &= x_2u_2 \\ \dot{u}_2 &= -x_2u_1 \end{split}$$

with the equality constraint

$$u_1^2 + u_2^2 - 1 = 0 \,.$$

Then the following Lyapunov function candidate is used:

$$V = a_1x_2^2 + a_2u_1^2 + a_3u_2^2 + a_4u_2 + a_5$$

with $a_3 + a_4 + a_5 = 0$. Then the search is performed for $V$ such that

$$V - \epsilon_1(1 - u_2) - \epsilon_2x_2^2 \geq 0$$

with $\epsilon_1 = 0.1, \epsilon_2 = 0.1$. Then the following Lyapunov function is derived:

$$V = 0.33445x_2^2 + 1.4615u_1^2 + 1.7959u_2^2 - 6.689u_2 + 4.8931 \,.$$

Problems:

My code (see above) works, but produces a different Lyapunov function... why is that? I think I used everything like explained, is there something else wrong?

clear;
syms x1 x2 u1 u2;
vars = [x1; x2; u1; u2];

f = [x2; u1*u2 - 10*u1; x2*u2; -x2*u1];
prog = sosprogram(vars);

[prog, V] = sospolyvar(prog, [1; u2; u2^2; u1^2; x2^2], 'wscoeff');
[prog, s1] = sospolyvar(prog, [x1^2; x2^2; u1^2; u2^2], 'wscoeff');

g1 = u1^2 + u2^2 - 1;

prog = sosineq(prog, V - 0.1*(1 - u2) - 0.1*x2^2);

dV = -(diff(V, x1)*f(1) + diff(V, x2)*f(2) + diff(V, u1)*f(3) + diff(V, u2)*f(4)) + s1*g1;
prog = sosineq(prog, dV);

prog = soseq(prog, coeff_3 + coeff_2 + coeff_1);
prog = sossolve(prog);
SOLV = sosgetsol(prog, V)


Also I don't see how the mentioned $V$ qualifies as a Lyapunov function... the following code creates a surface plot of its orbital derivative:

syms x1 x2 real;

f(1) = x2;
f(2) = sin(x1)*cos(x1) - 10*sin(x1);

V = 0.33445*x2^2 + 1.4615*sin(x1)^2 + 1.7959*cos(x1)^2 - 6.689*cos(x1) + 4.8931;
V_dot = diff(V, x1)*f(1) + diff(V, x2)*f(2);

x1_vec = -0.5:0.1:0.5;
x2_vec = -0.5:0.1:0.5;

V_map = zeros(length(x1_vec), length(x2_vec));
V_dot_map = zeros(length(x1_vec), length(x2_vec));

for k = 1:length(x1_vec)
for h = 1:length(x2_vec)
V_map(k, h) = subs(V, [x1, x2], [x1_vec(k), x2_vec(h)]);
V_dot_map(k, h) = subs(V_dot, [x1, x2], [x1_vec(k), x2_vec(h)]);
end
end

close all;
figure; surf(x1_vec, x2_vec, V_map.', 'edgecolor', 'none'); title('V');
figure; surf(x1_vec, x2_vec, V_dot_map.', 'edgecolor', 'none'); title('V_{dot}');


From the plot one can see that $\dot{V}$ is clearly positive in some regions (symmetric to the negative part) and increasing, for example along the $x_1 = x_2$ line, which is not supposed to happen?

Any help is appreciated.

• Your code does not in any sense encode the (in-)equality constraints in the model. Those must be appended to the model through the use of the positivstellensatz with multipliers (or similiar strategies) – Johan Löfberg Aug 13 '17 at 18:54
• Hm ok, could you provide a simple example or a reference how to do this? Because the examples in the paper don't go into detail on this point... – SampleTime Aug 13 '17 at 19:51
• Some example in the bottom of this page yalmip.github.io/tutorial/sumofsquaresprogramming. Has some references also, and the code is very similiar to sostools. In the paper, the polynomials g, sigma and lambda are the multipliers – Johan Löfberg Aug 14 '17 at 5:15
• Hi again, I got it to work, thanks. However, the results are strange: In the paper, example 4.2 (the whirling pendulum), the Lyapunov function given is $V = 0.33445x_2^2 + 1.4615\sin(x_1)^2 + 1.7959\cos(x_1)^2 - 6.689\cos(x_1) + 4.8931$. However, the derivative along the trajectories of this function is not negative. For example along the line $x_1 = x_2$ it is $> 0$ and increasing... what is wrong here? – SampleTime Aug 15 '17 at 17:59
• Along that line the vector flow looks almost parallel to the level curves, so I guess the requirements on negative definiteness is not strong enough, or there are numerical problems in the solution. Defining how definite something is can be really tricky, as you have to figure out which growth you should have on the decay (quadratic, quartic, mixed, etc). It is a mess/art to do stability analysis using SOS on general problems – Johan Löfberg Aug 16 '17 at 6:57

There has to be a typo somewhere, either in the paper, or in your code.

I tried to reconstruct the Lyapunov function (using YALMIP) and it is not possible to get that Lyapunov function.

First try where I try to find coefficients in the Lyapunov function as close as possible to the claimed data. The closest solution is far from the claimed

sdpvar x1 x2 u1 u2;

f = [x2; u1*u2 - 10*u1; x2*u2; -x2*u1];

c1 = sdpvar(5,1);
c2 = sdpvar(4,1);
V = c1'*[1; u2; u2^2; u1^2; x2^2];
s1 = c2'*[x1^2; x2^2; u1^2; u2^2]
g1 = u1^2 + u2^2 - 1;

Model = [sos(V-0.1*(1-u2)-0.1*x2^2) + sos(-jacobian(V,[x1;x2;u1;u2])*f + s1*g1)];
Model = [Model, sum(c1(1:3)) == 0];
optimize(Model,norm(c1 - [4.8931;-6.689;1.7959;1.4615;0.33445]),[],[c1;c2])
value(c1)


The closest I can get is quickly found by solving a mixed-integer program (YALMIP has a built-in framework for MISDP) where I try to maximize the number of matching coefficients

fits = binvar(5,1);
Model = [sos(V-0.1*(1-u2)-0.1*x2^2) + sos(-jacobian(V,[x1;x2;u1;u2])*f + s1*g1)];
Model = [Model, sum(c1(1:3)) == 0];
Model = [Model, implies(fits, c1 == [4.8931;-6.689;1.7959;1.4615;0.33445]),-100 <= c1 <= 100];
optimize(Model,-sum(fits),[],[c1;c2])
value(c1)


2 of the 5 coefficients match

Hence, my summary would be that the reported coefficients are wrong, or reported with to few significant digits.

Indeed, the following is infeasible

Vclaimed = 0.33445*x2^2 + 1.4615*u1^2 + 1.7959*u2^2 - 6.689*u2 + 4.8931;
Model = [sos(Vclaimed -0.1*(1-u2)-0.1*x2^2) + sos(-jacobian(Vclaimed ,[x1;x2;u1;u2])*f + s1*g1)];
optimize(Model,[],[],c2)


The fact that you get a completely different solution sounds natural as you're not minimizing any objective, and the solution to the problem is not unique.

To compute a verifiable solution, we can try to find an integer or rational solution. An integer solution is not feasible, but we can search for a rational solution.

Model = [sos(V-0.1*(1-u2)-0.1*x2^2) + sos(-jacobian(V,[x1;x2;u1;u2])*f+s1*g1)];
aux = sdpvar(5,1);
Model = [Model, sum(c1(1:3)) == 0,integer(aux), c1 == aux/20]
optimize(Model,[],[],[c1;c2])
format long
value(c1)


The solution we are interested in is the rational solution which is given by aux/20 (c1 will have errors from solver tolerances in the equality). We can compute the value of $\dot{V}$ in some random points and see that it is 0 everywhere

Vdot = replace(-jacobian(V,[x1;x2;u1;u2])*f,c1,value(aux)/20)
for i = 1:100
x = randn(2,1);
assign([x1 x2 u1 u2], [x(1) x(2) sin(x(1)) cos(x(1))]);
value(Vdot)
end

• First of all thanks for your detailed answer. However, it is impossible that 4 of 5 coefficients are correct and only the constant term $a_5$ is wrong. Thats simply due to the constraint $a_3 + a_4 + a_5 = 0$. If $a_3 = 1.7959$ and $a_4 = −6.689$ are correct, then $a_5$ would have no other choice than being $4.8931$... using the code from my question I get a solution with a residual norm of 1.0569e-10 and a feasration of 0.9976, so I would expect the result to actually be a valid Lyapunov function... however, looking at its plot it is not valid. – SampleTime Aug 17 '17 at 17:53
• Missed that constraint. If you add sum(c1(1:3)==0) you will see that it at most can match 2 coefficients, i.e., even worse so the reported numbers are just not to be trusted. – Johan Löfberg Aug 17 '17 at 18:00
• This is strange... looking at the example, there is even a second Lyapunov function given, which is parametrizable: $V = 7.3601 + 0.33033x_2^2 - 6.6066u_2 + 7.0698\dot{\theta}_a^2 + 0.3304\dot{\theta_a}^2u_2^2$. Here $\dot{\theta}_a^2$ is a constant value that should be $< 10$. However, regardless of the value I plug in for $\dot{\theta}_a^2$, the time derivative $\dot{V}$ gets positive for some $x$... – SampleTime Aug 17 '17 at 18:27
• I added code to compute a rigorous certificate, i.e. a solution where $\dot{V}$ turns out to be identical to 0. – Johan Löfberg Aug 17 '17 at 19:02
• ...no wonder thought that $\dot{V}=0$ is the best one can get, as the system isn't asymptotically stable but enters a limit-cycle directly – Johan Löfberg Aug 17 '17 at 19:45