SOSTools - Getting example to work?

Question

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.

Answer 1

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