H-infinity state feedback controller

[Mohit]

Hi,

I am trying to design a full state feedback h-infinity controller. The code is given below. I am not getting the correct answer for the problem. I have verified that by calculating the infinity norm of the closed loop system for the controller obtained from solving LMI problem.

Also, I have found out the variables which satisfies the constraints and still gives me a value of the norm smaller than what is obtained from the LMI problem.

The code goes like this:

----------------------------------------------------------------------------------------

k=90; m=10; c=5;

A=[0 1 0 0 ;

    -k/(2*m) -c/(2*m) 0 0;

    0 0 0 1;

    0 0 -k/m -c/m];

% Control input matrix

Bu = [0;1/m;0;0];

% Disturbance input matrix

Bw = [0; 1/m; 0;1/m];

% Ouput performance matrix

C = [0 -1 0 1]; 

Du = [0];

Dw = [0];

% Decision variable

Q = sdpvar(4,4,'sym');         

Y = sdpvar(1,4,'full');     

g = sdpvar(1,1);

% Solver Definition

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

% LMI definition

LMIs = [Q >=0];

LMIs = LMIs+[[(Q*A'+A*Q+Y'*Bu'+Bu*Y) Bw (Q*C'+Y'*Du'); Bw' -g Dw'; (C*Q+Du*Y) Dw -eye(1)] <=0];

diag = solvesdp(LMIs,g,options);      

check_LMI = checkset(LMIs)            

if check_LMI > 0 

   Q = double(Q);                     

   Y = double(Y);

   g = double(g);

   P = inv(Q)

   K = Y*P

   H_infty_norm = sqrt(g)

else disp('non feasible')

end

Ac=A+Bu*K;

Bc=Bw;

Cc=C+Du*K;

Dc=Dw;

G=ss(Ac,Bc,Cc,Dc);

Hnorm=norm(G,inf)

----------------------------------------------------------------------------------------------------------------------------

The output obtained is:

K =

  1.0e+002 *

   0.408193426197636  -1.629973303194311  -0.766755258053088   1.610783630944813

H_infty_norm =

   0.007210126360415

Hnorm =

   0.003707721768135

-------------------------------------------------------------------------------------------------------------------------

I have found the Q,Y and g which also satisfies the LMI constraints. The vale of g is 7.0810e-5, which is lesser than the value obtained above.

Code for that is given below:

clc

clear all

k=90; m=10; c=5;

A=[0 1 0 0 ;

    -k/(2*m) -c/(2*m) 0 0;

    0 0 0 1;

    0 0 -k/m -c/m];

% Control input matrix

Bu = [0;1/m;0;0];

% Disturbance input matrix

Bw = [0; 1/m; 0;1/m];

% Ouput performance matrix

C = [0 -1 0 1]; 

Du = [0];

Dw = [0];

K=[-k/(2) -c/(2) 0 0];

Ac=A+Bu*K;

Bc=Bw;

Cc=C+Du*K;

Dc=Dw;

G=ss(Ac,Bc,Cc,Dc);

Hnorm=norm(G,inf)

P= [  22.806803394759189   0.495132236091371 -22.806801858346688  -0.495132217887068

   0.495132236091371   2.591913924121580  -0.495132217886930  -2.591913753578595

 -22.806801858346688  -0.495132217886930  22.806803394760919   0.495132236091393

  -0.495132217887068  -2.591913753578595   0.495132236091393   2.591913924121776];

Q=inv(P)

Y=K*Q;

g= 7.0810e-5;

T=[(Q*A'+A*Q+Y'*Bu'+Bu*Y) Bw (Q*C'+Y'*Du'); Bw' -g Dw'; (C*Q+Du*Y) Dw -eye(1)]

[v,d]=eig(T)

-------------------------------------------------------------------------------------------------------------------------------------

I am not able figure out why Yalmip not giving me the correct value of the infinity norm. It is a convex optimization problem. So, there cannot be a case of local minima.

Does it have anything to do with the precision of the solver? I am using sdpt3.

Any kind of help is highly appreciated.

Thanks

Mohit 

[Johan Löfberg]

Your problem is ill-conditioned. The optimal value is very close to zero, but that happens when Q and Y tends to infinity (P goes singular). SDPT3 simply gives up when it thinks your solutions are too large. Other solvers yield other results

diag = solvesdp(LMIs,g,sdpsettings('solver','sedumi','verbose',0));

double(g) 

[norm(double(Q)) norm(double(Y))]

double(Y)*inv(double(Q))

diag = solvesdp(LMIs,g,sdpsettings('solver','sdpt3','verbose',0));

double(g) 

[norm(double(Q)) norm(double(Y))]

double(Y)*inv(double(Q))

diag = solvesdp(LMIs,g,sdpsettings('solver','mosek','verbose',0));

double(g) 

[norm(double(Q)) norm(double(Y))]

double(Y)*inv(double(Q))

ans =

   1.4417e-06

ans =

   1.0e+05 *

    0.3528    2.3030

ans =

   1.0e+03 *

    0.0447   -1.6203   -0.0888    1.6225

ans =

   5.1986e-05

ans =

   1.0e+03 *

    0.4760    2.2886

ans =

   40.8193 -162.9973  -76.6755  161.0784

ans =

   5.7191e-07

ans =

   1.0e+05 *

    0.5680    4.5166

ans =

   1.0e+03 *

    0.0449   -1.2813   -0.0894    1.2803

[Mohit]

Hi Johan

Yes, you are absolutely right. Indeed the optimal value of the problem is zero. I see that the singularity of P is causing the problem.

The actual problem was:

Minimize: g 

Subject to: P>0 and

               [(A+Bu*K)' *P+P*(A+Bu*K)    P*Bw      (C+Du*K)'  ;

                                 *                         -g*I            Dw'       ;            < 0

                                 *                            *              -I         ]

* - denotes the symmetric values and I is the identity matrix of the appropriate size. The variables of the optimization are P (sym), K (full) and g.

Since the second constraint is not linear ( as it has Bu*K*P terms in it) it was pre and post multiplied with diag{K^-1,I, I}

and later on the following substitutions were made: Q= inv(P), Y=K*Q, to obtain the following LMI in Q, Y and g

 

Minimize: g 

Subject to: Q>0 and

               [(QA'+AQ+BuY+Y'Bu'    Bw      QC'+Y'Du'  ;

                                 *                -g*I            Dw'     ;            < 0

                                 *                  *              -I         ]

Is there anyway to shift the optimal of the optimization by some transformation so that the singularity of P can be avoided? 

Or if there is any other way in which I can pose this problem so that the singularity of the P will not come into picture.

 

[Johan Löfberg]

The only thing that comes to my mind is a paper by a colleague of mine on (Anders Helmersson, Employing Kronecker Canonical Form for

LMI-Based Synthesis Problems)

http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=06104363

[Mohit]

I will have a look at that paper.

Thank you very much for your prompt reply. Really helped me a lot.

[Ilhan Polat]

I am quickly jumping to conclusion but as far as I can see your plant is not controllable as the second mass-spring-damper system is decoupled from the first. Hence your system is not minimal. But you need a generalized plant for Hinf theory. Hence the gamma level you find is not relevant.

[Mohit]

Hi IIhan

I am trying to find a controller which matches the response of the first mass-spring-damper-system with that of the the second mass-spring-damper-system based on Hinf theory. The second system is totally independent. It is like matching trajectory of the first system with that of the second. I want to use LMI for that. Will it work if I first find the minimal realization of the system by using minreal in matlab and then find the Hinf norm using LMI?

Do you have any suggestions on that.

[SKS]

Hi,

I am trying to design a full state feedback h-infinity controller. The code is given below. I am not getting the controller values. please find the mistake and correct me, i am hopeful for your optimiscic response.

clc
clear all
close all
yalmip 'clear'
X=sdpvar(4);
Y=sdpvar(4);
Ka=sdpvar(4,4,'full');
Kb=sdpvar(4,4,'full');
Kc=sdpvar(2,4,'full');
Kd=sdpvar(2,4,'full');
W=1.1;
ls=0.0286e-5; lm=0.0032e-2; lr=0.0286e-5;
a=[ls 0 lm 0; 0 ls 0 lm; lm 0 lr 0; 0 lm 0 lr];
b=inv(a);
lsi=3.1252; lmi=0.0279; w0=2*pi*50; delw=2*pi*.04*50; rs=0.01; rr=0.01;
A=[-rs/ls w0 rs*lm/ls 0; -w0 -rs/ls 0 rs*lm/ls; -lsi*rs/ls (w0*lsi-lmi*delw*lm/ls) (lsi*rs*lm/ls-rr) (lmi*delw*lm^2/ls-lmi*delw*lm);  (lmi*w0*lm/ls-w0*lsi) -rs*lsi/ls (lmi*delw*lr-lmi*delw*lm^2/ls) (lsi*rs*lm/ls-lmi*rr)]
B1=[1 0 0 0; 0 1 0 0; lsi 0 0 0; 0 lsi 0 0]
B2=[0 0; 0 0; -lmi 0; 0 -lmi]
C1=[0 0 -1 0; 0 0 0 -1]
C2=eye(4)
D11=[0 0 1 0; 0 0 0 1]
D12=[0 0; 0 0]
D21=zeros(4)
D22=zeros(2)
Ny=null([B2' D12'], 'r')
ny12=zeros(4,6);ny21=zeros(4);ny22=eye(6);
n1=cat(1,Ny,ny21);n2=cat(1,ny12,ny22);
Ny1=cat(2,n1,n2);
% Nx=null([C2 D21], 'r')
Nx=eye(4)
F=[[Nx zeros(4,6);zeros(6,4) eye(6)]'*[X*A+A'*X X*B1 C1';B1'*X -W*eye(4) D11';C1 D11 -W*eye(2)]*[Nx zeros(4,6);zeros(6,4) eye(6)]<0,Ny1'*[A*Y+Y*A' Y*C1' B1;C1*Y -W*eye(2) D11;B1' D11' -W*eye(4)]*Ny1<0,[Y eye(4);eye(4) X]>0,X>0,Y>0]
sol=solvesdp(F)
X=double(X)
Y=double(Y)
W=double(W)
% X=vpa(X,8)
% % % E=eig(vpa(P))
% Y=vpa(Y,8)
X =[  211711.94,   310.62937, -67712.046,   -75.11958;
 310.62937,   107.59296, -99.381282, -0.33541042;
 -67712.046,  -99.381282,  21666.755,   24.037673;
 -75.11958, -0.33541042,  24.037673, 0.097628741]
 
 
Y =[ 108.77908,  -0.11658389,    1.9857134, 0.026509699; -0.11658389,    106.85555, -0.093638681,   1.8558274;
  1.9857134, -0.093638681,    2498468.2,   2.3378861;
0.026509699,    1.8558274,    2.3378861,   2498468.3]
W=vpa(W)
Rx=chol(X)
 Ry=chol(Y)
%   Rx=vpa(Rx)
%   Ry=vpa(Ry)
 Z=Ry*Rx';
 [U,S,V]=svd(Z)
 sigma2=S*S'
 row2=sigma2-eye(4)
 row=sqrt(row2)
 T=((S)^(-1/2))*V'*Rx
 X3=S
 X2=T'*row

% p=[X*A+A'*X A'*X2 X*B1 C1'; X2'*A zeros(4) X2'*B1 zeros(4,2); B1'*X B1'*X2 -W*eye(4) D11'; C1 zeros(2,4) D11 -W*eye(2)]
% q=[X*B2 X2; X2'*B2 X3; zeros(4,2) zeros(4); D12 zeros(2,4)]*[Kd Kc;Kb Ka]*[C2 zeros(4) D21 zeros(4,2); zeros(4) eye(4) zeros(4) zeros(4,2)]
% r=[C2 zeros(4) D21 zeros(4,2); zeros(4) eye(4) zeros(4) zeros(4,2)]'*[Kd Kc;Kb Ka]'*[X*B2 X2; X2'*B2 X3; zeros(4,2) zeros(4); D12 zeros(2,4)]'
% Fp=(plus(p,q))

Q=[X*A+A'*X A'*X2 X*B1 C1';
    X2'*A zeros(4) X2'*B1 zeros(4,2);
    B1'*X B1'*X2 -W*eye(4) D11';
    C1 zeros(2,4) D11 -W*eye(2)];
R=[X*B2*Kd*C2+X2*Kb*C2 X*B2*Kc+X2*Ka X*B2*Kd*D21+X2*Kb*D21 zeros(4,2);
    X2'*B2*Kd*C2+X2*Kb*C2 X2'*B2*Kc+X3*Ka X2'*B2*Kd*D21+X3*Kb*D21 zeros(4,2);
    zeros(4) zeros(4) zeros(4) zeros(4,2);
    D12*Kd*C2 D12*Kc D12*Kd*D21 zeros(2)];
S=[C2'*Kd'*(X*B2)'+C2'*Kb'*X2' C2'*Kd'*(X2'*B2)'+C2'*Kb'*X3' zeros(4) C2'*Kd'*D12';
    Kc'*(X*B2)'+Ka'*X2' Kc'*(X2'*B2)'+Ka'*X3' zeros(4) Kc'*D12';
    D21'*Kd'*(X*B2)'+D21'*Kb'*X2' D21'*Kd'*(X2'*B2)'+D21'*Kb'*X3' zeros(4) D21'*Kd'*D12';
    zeros(2,4) zeros(2,4) zeros(2,4) zeros(2)];
Fp=[plus(Q,R,S)<0,Ka>0,Kb>0,Kc>0,Kd>0];
sol=solvesdp(r<0,Ka>0,Kb>0,Kc>0,Kd>0)
Ka=double(Ka)
Kb=double(Kb)
Kc=double(Kc)
Kd=double(Kd)

 

Any kind of help is highly appreciated.

Thanks

satendra kumar

PHD Scholar
 

[Johan Löfberg]

Well, the very first problem you solve is infeasible, and you are not checking that before the code proceeds

solvesdp is obsolete it is called optimize now, and strict inequalities are not supported (and yalmip is screaming at you about this)

https://yalmip.github.io/tutorial/basics/

already you first constraint is infeasible....

optimize(F(1))

ans = 

  struct with fields:

    yalmiptime: 0.4154

    solvertime: 0.3876

          info: 'Infeasible problem (SDPT3-4)'

       problem: 1

in addition to that (or perhaps the cause for all the issues), you data is numerically horrible

A =

   1.0e+05 *

   -0.3497    0.0031    0.0000         0

   -0.0031   -0.3497         0    0.0000

   -1.0927    0.0094    0.0000    0.0000

   -0.0000   -1.0927   -0.0000    0.0000