Lyapunov Function using sostool and Yalmip

[Muhammad Ali Murtaza]

Respected Professor Johan,

I was trying to reproduce an example given in documentation of sostool. Its about computing Lyapunov function. I was trying to compare solution using sostool for sosdemo2 example and using yalmip to solve the same problem.

System dynamics are given as follow

f = [-x1^3-x1*x3^2;
    -x2-x1^2*x2;
    -x3+3*x1^2*x3-3*x3/(x3^2+1)];

V is supposed to have only quadaratic terms. The candidate Lyapunov function must satisfy the following condition

V- (x1^2+x2^2+x3^2) > 0,

-(diff(V,x1)*f(1) + diff(V,x2)*f(2) + diff(V,x3)*f(3) )*(x3^2+1) > 0

The solution when solved using sostools give the following

7.02*x1^2 + 4.767*x2^2 + 1.984*x3^2

(They did reported in manual 5.5489x1^2 + 4.1068x2^2 + 1.7945x3^2, which may be o since solution looks closer). Code is provided in demo of sostool  

I code the following equivalent problem in yalmip

clear;
close all;

sdpvar x1 x2 x3;

% f = [-x1^3 - x1*x3^2;-x2 - x1^2*x2; -x3 + 3*x1^2*x3 - (3*x3)/(x3^2+1) ];
f = [-x1^3-x1*x3^2;
    -x2-x1^2*x2;
    -x3+3*x1^2*x3-3*x3/(x3^2+1)];


% v1 = monolist([x1 x2 x3],2);
% c1 = sdpvar(length(v1),1);   

% V = c1'*v1;
c1 = sdpvar(3,1);
V = c1'*[x1^2;x2^2;x3^2];

options = sdpsettings('solver','sedumi');
Model = [sos(V-1*(x1^2+ x2^2+ x3^2)), sos(-jacobian(V,[x1])*(x3^2+1)*f(1) + ...
    -jacobian(V,[x2])*(x3^2+1)*f(2) + -jacobian(V,[x3])*(x3^2+1)*f(3))];    
[sol,v,Q] = solvesos(Model,[],options,c1)
valc1 = value(c1)

It gives me PRIMAL_INFEASIBLE error

    yalmipversion: '20181012'
       yalmiptime: 0.5126
       solvertime: 0.1354
             info: 'Unbounded objective function (SeDuMi-1.3)'
          problem: 2

Can you let me know where I am doing wrong here. I tried with mosek and it gave me the same error. I did with sedumi since sostools uses sedumi.

[Johan Löfberg]

nonlinear fractions involving expressions are handled in a way that won't work when explicit representations are needed as in sos models, hence you will have to do the elimination of the fractional term manually, instead of simply multiplying the final expression with x3^2+1

or you can do 

sdpvar x1 x2 x3 t;

f = [-x1^3-x1*x3^2;

    -x2-x1^2*x2;

    -x3+3*x1^2*x3-3*x3/t];

c1 = sdpvar(3,1);

V = c1'*[x1^2;x2^2;x3^2];

options = sdpsettings('solver','sedumi');

x = [x1;x2;x3];

q = -jacobian(V,x)*f*t;

q = replace(q,t,x3^2+1)

Model = [sos(V-1*(x1^2+ x2^2+ x3^2)), sos(q)];    

[Muhammad Ali Murtaza]

Thank you Professor Johan for the guidance and correction. I am really obliged to you for taking the time out to correct our problem but is there a way to recognize the problem from yalmip error that this may be the cause of the problem.

[Johan Löfberg]

no, at the moment there is no check for this

[Muhammad Ali Murtaza]

Thank you Professor Johan.

Also the technique you mentioned above( declaring an additional variable and then using replace command before passing it to solvesos) always address this problem or is specific to this problem.

Finally one last question if we are using solvesos and we have few sdpvar declared but we want to fix some of the sdpvar variable and wan wantt to do optimization with rest of variables only. Mentioning only in decision variable of solvesos address that problem?  

I do have an sdpvar in my objective but I want to treat it as constant so if I dont specify it as my decision variable, will yalmip consider it as constant?

[Johan Löfberg]

It's specific to this problem.

The problem is that if you create something like g = y/f(x) where f(x) sin't a monomial, YALMIP will internally introduce a new variable z, return y/z and internally add a constraint [z == f(x)] to normalize the denominator. When you then create something by multiplying g with f(x)  to get rid of the numerator, the result will be g*f(x)  = (y/z)*f(x), i.e., the term is not eliminated, and the expression is impossibly a sos.

If you want a variable to be a constant, why even declare it as a decision variable? But sure, if you add x == thisvalue, it will automatically be detected as a parameter

[Muhammad Ali Murtaza]

Thank you Professor Johan for the detailed reply.