Verifying feability of LMIs using seduli (YALMIP interface)

[Shaik]

Hello,

        I would like to thank Johan Lofberg 

for the YALMIP which eased my LMI formulation, and of course I will refer in my results. 

       Coming to the point, I tried to solve a controller synthesis problem via the feasibility of LMIs in YALMIP interface. My program is given below and my issue is that the eig_max (second line from the bottom of the code), which is the maximum of eigen values of LMIs, is of the order 10^-12. Can i still consider the LMIs as feasible? Why I am asking is, when I use the controller of the resulting program (K1,K2), I obtain the stable dynamics. 

X =sdpvar(2,2,2,'symmetric'); Y =sdpvar(1,2,2);

ld_L=[0.2 0.8;0.8 0.2];ld_H=ld_L;

H1=[2 10;0 2]; H2=[2 0;10 2]; B1=-[1 0]'; B2=-[0 1]'; %system dynamics

ld_L=sqrt(ld_L);ld_H=sqrt(ld_H);

lmi1=blkvar;

lmi1(1,1)=-X(:,:,1);

lmi1(1,3)=(ld_L(1,1))*(H1*X(:,:,1)+B1*Y(:,:,1));

lmi1(2,2)= -X(:,:,2);

lmi1(2,3)=(ld_L(1,2))*(H1*X(:,:,1)+B1*Y(:,:,1)); 

lmi1(3,3)=- -X(:,:,1);  

lmi1=sdpvar(lmi1);

lmi2=blkvar;

lmi2(1,1)=-X(:,:,1);

lmi2(1,3)=(ld_H(1,1))*(H1*X(:,:,1)+B1*Y(:,:,1));

lmi2(2,2)=-X(:,:,2);

lmi2(2,3)=(ld_H(1,2))*(H1*X(:,:,1)+B1*Y(:,:,1)); 

lmi2(3,3)=-X(:,:,1);

lmi2=sdpvar(lmi2);  

lmi3=blkvar;

lmi3(1,1)=-X(:,:,1);

lmi3(1,3)= (ld_L(2,1))*(H2*X(:,:,2)+B2*Y(:,:,2));

lmi3(2,2)=-X(:,:,2);

lmi3(2,3)= (ld_L(2,2))*(H2*X(:,:,2)+B2*Y(:,:,2)); 

lmi3(3,3)=-X(:,:,2);

lmi3=sdpvar(lmi3);

lmi4=blkvar;

lmi4(1,1)=-X(:,:,1);

lmi4(1,3)= (ld_H(2,1))*(H2*X(:,:,2)+B2*Y(:,:,2));

lmi4(2,2)=-X(:,:,2);

lmi4(2,3)=  (ld_H(2,2))*(H2*X(:,:,2)+B2*Y(:,:,2)); 

lmi4(3,3)=-X(:,:,2);

lmi4=sdpvar(lmi4);

F=set([X(:,:,1)>0,X(:,:,2)>0,lmi1<0,lmi2<0,lmi3<0,lmi4<0]);

options = sdpsettings('solver','sedumi', 'verbose',1,'debug',1);

solution =solvesdp(F,[],options) %#ok<*NOPTS>

X=double(X);Y=double(Y);

lmi1=double(lmi1);lmi2=double(lmi2);lmi3=double(lmi3);lmi4=double(lmi4);

eig_max= max([eig(lmi1);eig(lmi2);eig(lmi3);eig(lmi4)])      

K1=Y(:,:,1)/(X(:,:,1));  K2=Y(:,:,2)/(X(:,:,2)); % stablilizing controller

[Johan Löfberg]

Ouch, pretty nasty code to read with all those blkvars and 3-dimensional matrices.

Some details first

1. Square matrices are symmetric by default, so that flag is redundant here

2. There is no notion of strict inequalities, in fact, if you use a recent version of YALMIP, it should warn you when you use strict inequalities. See recent post on the YALMIP Wiki http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Blog.Prepare-your-code (however, sedumi has a bug which turns off all warnings in MATLAB, so if you have run sedumi after starting matlab, these warnings will never be seen)

3. You use the SET operator which has been declared obsolete for 5 years or so.

Anyway, the problem formulation is not well-posed. As formulated now, X=0 and Y=0 is feasible. Indeed, SeDuMi essentially returns zero on most of your variables. Since you most likely want positive definite X, you can dehomogenize the problem by adding an explicit positive definiteness constraint. Since the problem is homogeneous, the scale does not matter, hence a reasonable model is

F=[X(:,:,1)>=eye(2),X(:,:,2)>=eye(2),lmi1<=0,lmi2<=0,lmi3<=0,lmi4<=0];

This leads to an infeasible problem though, indicating that there is no solution with positive definite X.

[Shaik]

Thanks Johan,

                    As per your comment, I changed the constraints as

F=[X(:,:,1)>=eye(2),X(:,:,2)>=eye(2),lmi1<=0,lmi2<=0,lmi3<=0,lmi4<=0];

This resulted in the eigen values of X as

 eig(X(:,:,1)) =  1.0e-011 *   [0.099401830024720    0.161626763973713]

 eig(X(:,:,2)) =   [0.000000000001032    2.703464226753521]

which are positive, but most of them are close to zero. Can't I conclude that the solution is feasible?? Also the controllers K1=[2 10], K2=[10 2] lead to stable dynamics.

[Johan Löfberg]

Well, obviously SeDuMi has failed to return a solution. Both X are far from feasible (they should be larger than I). SeDuMi simply returns the current iterate it has when it detects infeasibility (look at the error code you get from solvesdp, it is problable 1, i.e., infeasible)

The fact that K1 and K2 are stabilizing is mere luck. The computation of K is essentially noise/noise

[Shaik]

Thanks for your reply, which is  a lot helpful.

[Johan Löfberg]

BTW, your problem is infeasible no matter what you do. The pair (H1,B1) is not stabilizable (you're stuck with an eigenvalue at 2)