ILMI algorithm

[Wajdi Belhaj]

I have written this code :
***********************************
Q = eye(2);
A=[0 1;1 0]; 
B=[1 0]';
C=[1 15];
delta = 1e-3;

a = sdpvar (1);
F = sdpvar(1,1);
P = sdpvar(2);
%X = sdpvar(2);

X = care(A,B*B',Q);

%i = 1;
while i <= 10 
RHS = [P zeros(2,1); zeros(1,2) 0];
LHS = [A'*P+P*A-X*B*B'*P-P*B*B'*X+X*B*B'*X (B'*P+F*C)' ; (B'*P+F*C) -eye(1)];

% Find optimal ai
ai = solvelocalgevp(LHS,RHS,[],a); 

% Now solve with that ai to recover the last feasible solution in the
% gevp code (could have been returned from there but it is easier to
% solve again...) 
% ops = sdpsettings('debug',1) ; 
Constraints = [LHS <= ai*RHS,P >= 0];
solvesdp(Constraints,trace(P))

if ai <= 0
P = double(P);
F = double(F);
break
end

Pi = double(P);
if norm(X-Pi) <= delta
P = [];
F = [];
break
end
X = Pi;
end
X
Pi= double(P) % extraire la matrice solution P
norm(X-Pi)
Fi= double(F) % extraire la matrice solution F
eig(A)
eig(A+B*Fi*C)
***********************************
with solvelocalgevp
***********************************
function t_upper = solvelocalgevp(A,B,C,t)
% Résoudre min t sous A(x)<=t*B(x), B>=0, C>=0

% Trouver des limites intiaux inf et sup 
%ops = sdpsettings('verbose',0); % En mettant verbose à 0, le solveur va exécuter un affichage minimal ops : variable options
ops = sdpsettings('debug',1);
sol = solvesdp([C>=0, A <= 0*B, B>=0],[],ops); % solvesdp(contraintes ,Objectif ,options)
if sol.problem == 0
% Ok, 0 est une limite sup , i.e. feasible 
t_upper = 0;
t_lower = -1;
while sol.problem == 0
t_lower = (t_lower)*2;
sol = solvesdp([C>=0, A <= t_lower*B, B>=0],[],ops);
end
else sol.problem = 1;
% Ok 0 est une limite inf, i.e. infeasible 
t_lower = 0;
t_upper = 1;
while sol.problem == 1
t_upper = t_upper*2;
sol = solvesdp([C>=0, A <= t_upper*B, B>=0],[],ops);
end
end
%else
% error('Hmm, you''d better robustify the code!')
%end

% OK, on a les limites , on commence la bissection 
while t_upper-t_lower >= 1e-3
t_mid = (t_upper+t_lower)/2;
sol = solvesdp([C>=0, A <= t_mid*B, B>=0],[],ops);
if sol.problem == 0
t_upper = t_mid; 
else
t_lower = t_mid;
end
%disp([t_lower t_upper]);
end
*************************************************************
I have obtained the following results :
iter seconds digits c*x b*y
21 0.4 Inf -7.8725074213e+001 -7.8725071780e+001
|Ax-b| = 4.7e-010, [Ay-c]_+ = 8.5E-010, |x|= 1.9e+004, |y|= 7.0e+001

Detailed timing (sec)
Pre IPM Post
3.120E-002 3.900E-001 0.000E+000 
Max-norms: ||b||=1, ||c|| = 7.989655e+001,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 2.44718.
ans = 
yalmiptime: 0.225999999999999
solvertime: 0.326999999999998
info: 'Successfully solved (SeDuMi-1.1)'
problem: 0
dimacs: [NaN NaN NaN NaN NaN NaN]
X =
8.938487227006524 3.738119536210800
3.738119536210800 46.878546847633686
Pi =
9.480842159422313 5.051249067025251
5.051249067025251 69.244229620159260
ans =
22.444410959621575
Fi =
-0.573452926049265
ans =
-1
1
ans =
-0.286726463024632 + 2.742185592942308i
-0.286726463024632 - 2.742185592942308i

but the exact results (in the article SOF stabilization : an ILMI approach - Cao et al. , 1998) are :
ai =
-0.0377 
Pi =
9.2557 6.48
6.4888 92.9376
Fi =
-0.7369

How is wrong with my code , please ?

[Johan Löfberg]

You exit immediately when ai turns negative, so there is no reason you should get anything optimal.

In additional to that error, you should really check the error code and see that you actually solve the problem successfully. This would be a better logic

while i <= 100 % Or what ever logic you want to control #iterations...

    RHS = [P zeros(2,1); zeros(1,2) 0];

    LHS = [A'*P+P*A-X*B*B'*P-P*B*B'*X+X*B*B'*X (B'*P+F*C)' ; (B'*P+F*C) -eye(1)];

    

    ai = solvelocalgevp(LHS,RHS,[],a)

       

    Constraints = [LHS <= ai*RHS,P >= 0];

    sol = solvesdp(Constraints,trace(P),sdpsettings('verbose',0));

    

    if sol.problem ~=0

        % For some reason, it didn't even find the old solution and failed,

        % either completely, or some numerical issues arose. In any case,

        % we give up, break and thus return latest successful solution
        % (or do something to mitigate these issues, for instance scale ai with 0.9999
        %  and try again, etc)

        break

    end

    
    % These were correct w.r.t ai, so we keep them

    Pi = double(P);

    Fi = double(Fi);

    

    if norm(X-Pi) <= delta       

        break

    end

    X = Pi;

end

X

norm(X-Pi)

eig(A)

eig(A+B*Fi*C)

You will note that you get a better result than what is reported in the paper. Here, I get -0.28

[Wajdi Belhaj]

it worked first of all . But when I repeated , it showed wrong results . I don't know if I should reset Sedumi and Yalmip in "Set Path" for Matlab or there is a problem in the code , thanks .

[Johan Löfberg]

What do you mean with wrong result?

Gives the same result repeatedly when rerunning the file here (with the changes I proposed). Note though that you get different results depending on solvers, due to different numerics leading to different paths in your local heuristic. SeDuMi quits at -0.28 while SDPT3 ends up at -6.05

[Wajdi Belhaj]

%end

% OK, on a les limites , on commence la bissection 

while t_upper-t_lower >= 1e-3

    t_mid = (t_upper+t_lower)/2;

    sol = solvesdp([C>=0, A <= t_mid*B, B>=0],[],ops);

    if sol.problem == 0

        t_upper = t_mid;       

    else

        t_lower = t_mid;

    end

    %disp([t_lower t_upper]);

end

************************************************************************************

Q = eye(2);

A=[0 1;1 0]; 

B=[1 0]';

C=[1 15];

delta = 1e-3;

a = sdpvar (1);

F = sdpvar(1);

P = sdpvar(2);

%X = sdpvar(2);

X = care(A,B*B',Q);

while i <= 100 % Nbre iterations...

    RHS = [P zeros(2,1); zeros(1,2) 0];

    LHS = [A'*P+P*A-X*B*B'*P-P*B*B'*X+X*B*B'*X (B'*P+F*C)' ; (B'*P+F*C) -eye(1)];

    

    ai = solvelocalgevp(LHS,RHS,[],a);

    

    

    %ops=sdpsettings('verbose',0,'warning',0);

    Constraints = [LHS <= ai*RHS,P >= 0];

    sol = solvesdp(Constraints,trace(P),sdpsettings('debug',1));

    

    if sol.problem ~=0

        

        break;

    end

    

    % 

    Pi = double(P);

    F = double(F);

    

    if norm(X-Pi) <= delta       

        break;

    end

    X = Pi;

end

ai

X

Fi=double(F)

Pi=double(P)

norm(X-Pi)

eig(A)

eig(A+B*Fi*C)

******************************************************************************************************

I have obtained the following results :

ai =

   0.037109375000000

X =

   9.894680859789526   1.932223148752961

   1.932223148752961  16.471470951978475

Fi =

  -0.178009124110935

Pi =

  10.073885697079936   1.948260616736751

   1.948260616736751  16.765692017211499

but the exact results with the same approach (in the article SOF stabilization : an ILMI approach - Cao et al. , 1998) are :

ai =
-0.0377 
Pi =
9.2557 6.48
6.4888 92.9376
Fi =
-0.7369

What is wrong mith my code or what I have to change to the code to be close to the final result , thanks .

[Johan Löfberg]

Why do you keep saying that the result in the referenced paper is "exact". What is "exact" about it? Clearly, you can compute solution with better (lower ai) values, as I just have showed you. Additionally, what you have is a local heuristic, so depending on solver, computer architecture, stopping criteria and general logic, position of the moon etc, you can converge to vastly different solutions.

BTW, there is a bug in both my code, and the code you copied from me (although you changed the code but introduced another bug). It should be

    Pi = double(P);

    Fi = double(F);

Still get the value -.28 using SeDuMi and -6 with SDPT3

[Wajdi Belhaj]

Given & (with realization ( At, Bt, Ct)) stabilizable via state feed-back.

Step 1. Select Q>0, and solve P from the following algebraic
Riccati equation
A'*P+P*A - P*B B'*P+Q=0.

Set i=1 and X1=P.

Step 2. Solve the following optimization problem for Pi , Ft and ai.
OP1: Minimize ai subject to the following LMI constraints
[A'*Pi+ Pi*A - Xi*B*B'*Pi - Pi*B*B'*Xi+Xi*B*B'*Xi - ai*Pi (B'*Pi+Ft*C)' ; (B'*Pi+Ft*C) -I] < 0 (1)
Pi = Pi' > 0 (2)

(eye(2)+(C*B*F3t)' + (C*B*F3t))>0 (3)
Denote a*i as the minimized value of ai.

Step 3. If a*i <= 0, F is a stabilizing static output feedback gain. Stop.

Step 4. Solve the following optimization proble for Pi and Ft.

OP2: Minimize trace(Pi) subject to the above LMI constraints (1)  (2) and (3) with ai = a*i . Denote P*i as the Pi that mini-mized trace (Pi).

Step 5. If || Xi*B - P*i*B || < delta , a prescribed tolerance, go to Step 6,
else set i = i+1 and Xi = P*i-1 , then go to Step 2.

Step 6. The system may not be stabilizable via static output feedback. Stop

********************************************************************************

%end

% OK, on a les limites , on commence la bissection 

while t_upper-t_lower >= 1e-4

    t_mid = (t_upper+t_lower)/2;

    sol = solvesdp([C>=0, A <= t_mid*B, B>=0],[],ops);

    if sol.problem == 0

        t_upper = t_mid;       

    else

        t_lower = t_mid;

    end

    %disp([t_lower t_upper]);

end

*********************************************************************************************************

A=[-0.0266 -36.6170 -18.8970 -32.0900 3.2509 -0.7626; 0.0001 -1.8997 0.9831 -0.0007 -0.1708 -0.0050; 0.0123 11.7200 -2.6316 0.0009 -31.6040 22.3960 ; 0 0 1 0 0 0; 0 0 0 0 -30 0; 0 0 0 0 0 -30];

B=[0 0; 0 0; 0 0; 0 0; 30 0; 0 30];

C=[0 1 0 0 0 0; 0 0 0 1 0 0];

[n,n]=size(A);

[n,m]=size(B);

[r,n]=size(C);

%Ft=[F1t F2t F3t];

C1t=[C zeros(2,2)];

C2t=[zeros(2,6) eye(2,2)];

C3t=[C*A zeros(2,2)];

At= [A zeros(6,2); C zeros(2,2)];

Bt = [B ;zeros(2,2)];

Ct = [C1t' C2t' C3t']';

Q = eye(8);

% A=[0 1;1 0]; 

% B=[1 0]';

% C=[1 15];

delta = 1e-5;

a = sdpvar (1);

Ft = sdpvar(2,6);

P = sdpvar(8,8);

%X = sdpvar(2);

eps = 0.01;

X = care(At,Bt*Bt',Q);

%i = 1;

%Y=eye(2);

%F3t=0;

%while (Y>zeros(2))  

%while i <= 100    

    RHS = [P zeros(8,2); zeros(2,8) zeros(2,2)];

    LHS = [At'*P+P*At-X*Bt*Bt'*P-P*Bt*Bt'*X+X*Bt*Bt'*X (Bt'*P+Ft*Ct)' ; (Bt'*P+Ft*Ct) -eye(2)];

    

    ops = sdpsettings('debug',1) ; 

    Constraints1 = [LHS <= a*RHS,P >= 0];

    solvesdp(Constraints1,a);

    

%     

% %     % Find optimal ai

%     ai = solvelocalgevp(LHS,RHS,[],a);  

     F1t=Fi(:,1:2);

     F2t=Fi(:,3:4);

     F3t=Fi(:,5:6);

     F3 = F3t*inv(eye(2)+ C*B*F3t);

     F2 = (eye(2)-F3*C*B)*F2t ;

     F1 = inv(eye(2)-F3*C*B)*F1t  ;

     Y =  eye(2)+(C*B*F3t)' + (C*B*F3t);    

    

    if ai <= 0

        Pi = double(P);

        Fi = double(Ft);

        break;%     end

%     

%     

%     if ai <= 0

%         Pi = double(P);

%         Fi = double(Ft);

%         break;

%     else

    else

    % Avec ce alpha i , résoudre le pb d'optimisation suivant : min trace(P) 

    %ai=double(a); 

    % ops = sdpsettings('verbose',0,'warning',0) ; 

     ops = sdpsettings('debug',1) ; 

    Constraints = [LHS <= ai*RHS,P >= 0,(eye(2)+(C*B*F3t)' + (C*B*F3t))>0];

    %Constraints = [LHS <= ai*RHS,P >= 0, eye(2)+(C*B*(Ft(:,5:6)))'+ (C*B*(Ft(:,5:6)))>0];

    sol = solvesdp(Constraints,trace(P),ops);

    

    

    Pi = double(P);

    Fi = double(Ft);

    

    Pi = double(P);

    Y = (norm(X-Pi) - delta );

    

   % if Y > 0

   %    disp ('La matrice I+C*B*F3t est singulière ') 

   % end 

    

    if (Y<=zeros(2))

    % norm(X-Pi) <= delta

       disp ('On ne peut pas décider par cet algo si le problème SOF est solvable')

       break;

    end

    X = Pi;

    %X=X+eps;

    %ai=ai-eps;

       

    

%      if sol.problem ~=0

%          

%          break;

%      end

    

    

    end

%     if sol.problem ~=0

%         

%         break;

%     end

    

    % 

    

%end

ai

%X

%Pi= double(P) % extraire la matrice solution P

norm(X-Pi)

Fi= double(Ft) % extraire la matrice solution F

F1t=Fi(:,1:2);

F2t=Fi(:,3:4);

F3t=Fi(:,5:6);

F3 = F3t*inv(eye(2)+ C*B*F3t)

F2 = (eye(2)-F3*C*B)*F2t 

F1 = inv(eye(2)-F3*C*B)*F1t  

eig(At)

%double(F1)

%double(F2)

%double(F3)

eig(At+Bt*Fi*Ct)

*********************************************************************************************

I have obtained the following results

ai =

   0.722351074218750

ans =

     0

Fi =

  1.0e+002 *

   2.347327881525589  -0.290145665756742   0.994771000946797   0.245590961617259   1.011059947220535  -0.972712801524232

  -0.010049072140165  -0.169981834683681   0.226684934190822   0.054409669962516  -0.081218856321405   0.069705942898811

F3 =

  1.0e+002 *

   1.011059947220535  -0.972712801524232

  -0.081218856321405   0.069705942898811

F2 =

  99.477100094679713  24.559096161725872

  22.668493419082207   5.440966996251622

F1 =

  1.0e+002 *

   2.347327881525589  -0.290145665756742

  -0.010049072140165  -0.169981834683681

ans =

                  0                     

                  0                     

 -5.675661980211077                     

  0.688638462613330 + 0.245532253157104i

  0.688638462613330 - 0.245532253157104i

 -0.259514945015590                     

-30.000000000000000                     

-30.000000000000000                     

ans =

  1.0e+002 *

 -5.425949134839419                     

 -0.169344631658648 + 0.205728294750513i

 -0.169344631658648 - 0.205728294750513i

 -0.024345896829333 + 0.015745136534263i

 -0.024345896829333 - 0.015745136534263i

 -0.000242356664840 + 0.001277369241297i

 -0.000242356664840 - 0.001277369241297i

 -0.000252435964747  

I wrote this code but I have a problem of constraint ( a warning message appears that One of the constraints evaluates to a FALSE LOGICAL variable + run into numerical problems ) + alpha i minimum is positive and it has to be negative .

Any help please , thanks.

[Johan Löfberg]

I guess you mean

 ai = solvelocalgevp(LHS,RHS,[],a);

and not

solvesdp(Constraints1,a)

since that program is bilinear....

Moving on, your code is completely messed up. On the following line, you do

F1t=Fi(:,1:2);

but Fi has not yet been defined so the code crashes (you probably run the code with old garbage in your work-space, thus not crashing)

Please clean up your code first.

[Wajdi Belhaj]

%     ops = sdpsettings('debug',1) ; 

%     Constraints1 = [LHS <= a*RHS,P >= 0];

%     solvesdp(Constraints1,a);

%     

%     

% %     % Find optimal ai

     ai = solvelocalgevp(LHS,RHS,[],a);  

     

     Fi=double(Ft);

     

     F1t=Fi(:,1:2);

     F2t=Fi(:,3:4);

     F3t=Fi(:,5:6);

     F3 = F3t*inv(eye(2)+ C*B*F3t);

     F2 = (eye(2)-F3*C*B)*F2t ;

     F1 = inv(eye(2)-F3*C*B)*F1t  ;

     Y =  eye(2)+(C*B*F3t)' + (C*B*F3t);    

    

    if ai <= 0

        Pi = double(P);

        Fi = double(Ft);

        break;%     end

%     

%     

%     if ai <= 0

%         Pi = double(P);

%         Fi = double(Ft);

%         break;

%     else

    else

    % Avec ce alpha i , résoudre le pb d'optimisation suivant : min trace(P) 

    %ai=double(a); 

    % ops = sdpsettings('verbose',0,'warning',0) ; 

    ops = sdpsettings('debug',1) ; 

    Constraints = [LHS <= ai*RHS,P >= 0,(eye(2)+(C*B*F3t)' + (C*B*F3t))>0];

    %Constraints = [LHS <= ai*RHS,P >= 0, eye(2)+(C*B*(Ft(:,5:6)))'+ (C*B*(Ft(:,5:6)))>0];

    sol = solvesdp(Constraints,trace(P),ops);

    

       %if sol.problem ~=0

           

       %    break;

       %end

       

    Pi = double(P);

    Fi = double(Ft);

    

    Pi = double(P);

    Y = (norm(X-Pi) - delta );

    

   % if Y > 0

   %    disp ('La matrice I+C*B*F3t est singulière ') 

   % end 

    

    if (Y<=zeros(2))

    % norm(X-Pi) <= delta

       disp ('On ne peut pas décider par cet algo si le problème SOF est solvable')

       break;

    end

    X = Pi;

    %X=X+eps;

    %ai=ai-eps;

       

    

%      if sol.problem ~=0

%          

%          break;

%      end

    

    

    end

    

%end

ai

%X

%Pi= double(P) % extraire la matrice solution P

norm(X-Pi)

Fi % extraire la matrice solution F

F1 

F2 

F3 

eig(At)

%double(F1)

%double(F2)

%double(F3)

eig(At+Bt*Fi*Ct)

I wrote this code but I have the same problem of constraint ( a warning message appears that One of the constraints evaluates to a FALSE LOGICAL variable + run into numerical problems ) + alpha i minimum is positive and it has to be negative .

What's wrong with my code please ? Thanks.

[Johan Löfberg]

This makes no sense

Constraints = [LHS <= ai*RHS,P >= 0,(eye(2)+(C*B*F3t)' + (C*B*F3t))>0];

The matrix (eye(2)+(C*B*F3t)' + (C*B*F3t)) is a fixed matrix, there are no decision variables in it. It evaluates to the identity matrix. This identity matrix is then compared with 0 (since you perform >) and the result is a logical matrix of size 2x2. Thus, you the constraint reads

Constraints = [LHS <= ai*RHS,P >= 0,[true false;false true]];

YALMIP gets confused and throws the error you see

BTW, once again, strict inequalities should not be used.