Hi,
I am currently trying to find a Hinf controller via LMI. I have found a controller using matlab's control toolbox, but I would like to do this by solving the LMIs with yalmip. I feel like I am cheating when I am using the matlab functions and want to do it myself.
From matlab function I get an answer. But I do not get anything out from my yalmip solution, everything is NaN. Can anyone see where I have made mistakes?
Underneath is the code.

Thanks in advanced

Tore
--------------------
clear all
x=6;% states
y=3;% outputs
u=1;% inputs
A=[ 0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
-4.2192 7.1833e-006 8.4361e-006 -0.24951 -2.2901 -2.2901
-7.9894e-011 -6.3708e-011 172.6 7.351e-010 0.00012081 1.237
-0.00064632 4.1633e-006 -195 -0.028875 -0.37408 -1.7712];
B=[ 0
0
0
-9.9032
6.3169e-008
-1.4971];
C=[ 0 0 0 0 9.549 9.549
0 0 0 0 926.3 0
1 0 0 0 0 0];
D=[ 0
0
0];
%% Formulation the generalized plant
% Given an LTI plant P with state-space equations:
% dx/dt = A * x + B1 * w + B2 * u
% zinf = C1 * x + D11 * w + D12 * u
% y = C2 * x + D21 * w + D22 * u
B1=[1 0 0 0 0 1]';
B2=B;
C1=[1 1 0 2 1 1];% Hinf performance
D11=zeros(1,1);
D12=ones(1,u)*2;
C2=C;
D21=ones(y,1);
D22=D;
P=[A B1 B2;C1 D11 D12;C2 D21 D22];% Plant
%% LMI
gama0=10;% Upper bound
gama=sdpvar(1);
R=sdpvar(x);
S=sdpvar(x);
Ak=sdpvar(x);
Bk=sdpvar(x,y);
Ck=sdpvar(u,x);
Dk=sdpvar(u,y);
C_1=[A*R+R*A'+B2*Ck+Ck'*B2' Ak'+A+B2*Dk*C2 B1+B2*Dk*D21 (C1*R+D12*Ck)';
(Ak'+A+B2*Dk*C2)' A'*S+S*A+Bk*C2+C2'*Bk' S*B1+Bk*D21 (C1+D12*Dk*C2)';
(B1+B2*Dk*D21)' (S*B1+Bk*D21)' -1 (D11+D12*Dk*D21)';
C1*R+D12*Ck C1+D12*Dk*C2 D11+D12*Dk*D21 -gama^2]<0;
C_2=[gama^2<gama0^2];
F=[C_1,C_2];
solvesdp(F,gama^2);
R=double(R);
S=double(S);
Ak=double(Ak);
Bk=double(Bk);
Ck=double(Ck);
Dk=double(Dk);
Acl=[A+B2*Dk*C2 B2*Ck;Bk*C2 Ak];
eig(Acl)
------------------------------

The problem is that you do not have a linear matrix inequality, since you have a gama^2. This is easily fixed by performing the optimization problem in a variable which models gama^2,

I.e define
gamasquared=sdpvar(1);

and replace all occurances of gama^2 with this new linear variable.

However, there are more issues in your code. In YALMIP, square matrices are by default symmetric,. That means you currently force Ak to be symmetric. I don't think you want that. Hence, you should use

Ak = sdpvar(x,x,'full')

and similiar on all matrices which should be fully parameterized
users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Tutorials.Basics

BTW, your constraint C_2 has no meaning. Either the optimal solution satisfies this, or the problem is infeasible. It only makes sense if you don't have an objective.

Finally, the problem is still not solved. This is either due to some small thing I missed, or due to the fact that your model is insanely badly scaled. I heard a great quote recently, "Solvers use 10 digits of accuracy, five for actual reasonable modeling, and five extra to try to cover for bad modeling from users". Your A matrix ranges over 13 digits!

Hence, as I test, I balance the problem
S = balreal(ss(A,B,C,D))
A = S.A;
B = S.B;
C = S.C;
D = S.D;

Slight improvement in scaling. Still, I get crap solutions

Don't you want R and S to be positive definite?

I add that, and I can compute a solution which is strictly feasible. However, the eigenvalues of Ac1 are not all negative, so maybe you've modelled something wrong (I suspect you want Acl to be Hurwitz)

I think you should also have [R eye(x);eye(x) S] >0. Still, Acl is not Hurwitz despite a strictly feasible solution.

Hi again!

Thanks for good coments!
The reason for the funnylooking model is that it is the direct outcome form the wind turbine software FAST. FAST linearizes a nonlinear WT model, so often there is loooots of digits.
My plan was to make mixed h2/hinf control via LMIs, and now it works. In the code above is missing some vital parts regarding some change of variables tricks, but now it works!
For your info I give the code underneath!

Tore

---------------
clear all
clf
x=6;% states
y=3;% outputs
u=1;% inputs
A=[ 0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
-4.2192 7.1833e-006 8.4361e-006 -0.24951 -2.2901 -2.2901
-7.9894e-011 -6.3708e-011 172.6 7.351e-010 0.00012081 1.237
-0.00064632 4.1633e-006 -195 -0.028875 -0.37408 -1.7712];
B=[ 0
0
0
-9.9032
6.3169e-008
-1.4971];
C=[ 0 0 0 0 9.549 9.549
0 0 0 0 926.3 0
1 0 0 0 0 0];
D=[ 0
0
0];
%% Formulation the generalized plant
% Given an LTI plant P with state-space equations:
% dx/dt = A * x + B1 * w + B2 * u
% z2 = C12 * x + D1 * w + D2 * u
% zinf = C1i * x + D1i * w + D2i * u
% y = C2 * x + D21 * w + D22 * u
B1=[1 0 0 0 0 1]';
B2=B;
C12=[0 1 5 1 15 0];% H2 performance
D1=zeros(1,1);
D2=ones(1,u)*10;
% C1i=[10 0 0 0 0 5];% Hinf performance
C1i=[1 0 0 1e-1 1 1e-1];% Hinf performance
D1i=zeros(1,1);
D2i=ones(1,u)*2;
C2=C;
D21=ones(y,1);
D22=D;
%% LMI
alpha=10; % Hinf weight
beta=1; % H2 weight
gama0=16; % Hinf bound
eta0=40; % H2 bound
gama=sdpvar(1);
X=sdpvar(x);
Y=sdpvar(x);
W=sdpvar(1);
Ahat=sdpvar(x,x,'full');
Bhat=sdpvar(x,y);
Chat=sdpvar(u,x);
Dhat=sdpvar(u,y);
C_1=[A*X+X*A'+B2*Chat+(B2*Chat)' Ahat'+A+B2*Dhat*C2 B1+B2*Dhat*D21 (C1i*X+D2i*Chat)';
(Ahat'+A+B2*Dhat*C2)' A'*Y+Y*A+Bhat*C2+C2'*Bhat' Y*B1+Bhat*D21 (C1i+D2i*Dhat*C2)';
(B1+B2*Dhat*D21)' (Y*B1+Bhat*D21)' -gama (D1i+D2i*Dhat*D21)';
C1i*X+D2i*Chat C1i+D2i*Dhat*C2 D1i+D2i*Dhat*D21 -gama]<0;
C_2=[X eye(x);
eye(x) Y]>0;
% C_3=[lambdaij*[X eye(x);eye(x) Y]+myij*[A*X+B2*Chat*A+B2*Dhat*C2;Ahat Y*A+Bhat*C2]+myji*[X*A'+Chat'*B2' Ahat';(A+B2*Dhat*C2)' At*Y+C2'*Bhat']]<0;
C_3=[W C12*X+D2*Chat C12+D2*Dhat*C2;
(C12*X+D2*Chat)' X eye(x);
(C12+D2*Dhat*C2)' eye(x) Y]>0;
C_4=[trace(W)<eta0];
C_5=[gama<gama0];
F=[C_1,C_2,C_3,C_4,C_5];
obj=alpha*gama+beta*trace(W);
solvesdp(F,obj);
X=double(X);
Y=double(Y);
Ahat=double(Ahat);
Bhat=double(Bhat);
Chat=double(Chat);
Dhat=double(Dhat);
[u,sd,v]=svd(eye(x)-X*Y); % factorize I-XY
sd2=diag(sqrt(diag(sd)));
M=u*sd2;
Nt=sd2*v';
N=Nt';
Dk=Dhat;
Ck=(Chat-Dk*C2*X)*inv(M');
Bk=inv(N)*(Bhat-Y*B2*Dk);
Ak=inv(N)*(Ahat-Nt'*Bk*C2*X-Y*B2*Ck*M'-Y*(A+B*Dk*C2)*X)*inv(M');
% Constructing the closed loop
Bn=[B1 B];
Cn=[C12;C1i;C2];
Dn=[D1 D2;D1i D2i;D21 D22];
P=ltisys(A,Bn,Cn,Dn);
K=pck(Ak,Bk,Ck,Dk);
clsys=slft(P,K,1,1);
[Acl,Bcl,Ccl,Dcl]=unpck(clsys);
eig(Acl)
------------------------------

I don't get a stable Acl when I run your code

I also get one pole with location: 1.0274e-005, but I assume it is so close to zeros that assume it is ok. Next step for me is to implement pole placement constraints.