how to code iterative lmi algorithms using yalmip

[Rahul Kant]

I tried to use first yalmip  but i don't know what went wrong and then again i started with lmi toolbox and this time my solution is not converging.
I would like to use the yalmip toolbox properly and get the converged solution.
I am to find a robust pid controller using the iterative lmi algorithm with three unknowns-P,F,alpha(scalar); two constraints.
first i have to minimize alpha and obtain the corresponding F value i.e. my controller. in case alpha doesn't minimizes i have to go for second minimization i.e. to minimize trace(P).

[Johan Löfberg]

You would have to specify your question better.

Here is an iterative nonsense problem

solvesdp(P >= eye(2));
for i =1:100
 solvesdp(P*double(P)+double(P)*P >=eye(2));
end

[Rahul Kant]

According to the algorithm,
1. A,B,C the state matrices are realized first.
2. Q>0 is selected and a ARE is solved for the unknown P
       A'P+PA-PBB'P+Q=0
3.Iteration starts i=1 and X(i)=P(obtained from above step); Three unknowns are there: P(i), F, alpha(i)( a scalar)
4.Solve the optimization for P(i),F,alpha(i):
     minimize alpha(i) subjected to LMI constraints:
    4.1      [ A'P(i)+P(i)A - X(i)BB'P(i) - P(i)BB'X(i) + X(i)BB'X(i) - alpha(i)P(i)       (B'P(i)+FC)';
                   B'P(i)+FC                                                                                  -I            ] < 0
    4.2  P(i)>0
5. alpha(i) is obtained and checked if alpha(i)<=0
     if yes then F is obtained as the desired controller
     if no then go for second minimization
6.Solve optimization for P(i) and F with new value of alpha(i):
      minimize trace(P(i)) subject to same above mentioned constraints 4.1 & 4.2.
7. P(i) is obtained and checked
       || X(i)B - P(i)B||< e(tolerance predefined)
       if yes then terminate
       if no then go for second iteration

[Johan Löfberg]

So where are you stuck?

Q = magiccomputation;
P = are(...)

goon = 1;
while goon
 X = P;
 [alpha,P,F] = solvegevp(X,A,...) % Implement bisection to optimize alpha. See decay-rate example on Wiki
 if alpha <= 0
  ?
 else
   [P,F] = solvelmi(X,alpha,A,...) % Just setup and solve 4.1
 end
 goon = norm(X-P*B) <= tolerance
end

[Rahul Kant]

I understood the bisection algorithm u suggested me to see, but I am not able to implement it to my work. I am showing here what i did:
rt=1.5e-3;
lt=30e-6;
r=76;
l=111.9e-3;
c1=62.855e-6;
f=60;
w0=2*pi*f;
A=[-rt/lt w0 0 -1/lt;w0 -rt/lt -2*w0 ((rt*c1*w0/l)-(w0/r));0 w0 -rt/l ((1/l)-(w0^2*c1));1/c1 0 -1/c1 -1/(r*c1)];
B=[1/lt 0 0 0]';
C=[0 0 0 1];
A1=[A zeros(4,1);C 0];
B1=[B;0];
C1=[C 0;zeros(1,4) eye(1);C*A 0];
Q=eye(5);
[P,L,G] = care(A1,B1,Q);
X1=P;
options = [1e-5,0,0,0,0];
M=sdpvar(5,5);
F=sdpvar(1,3);
a=sdpvar(1,1);
S=[M>0];
K=S+[A1'*M+M*A1-X1*(B1*B1')*M-M*(B1*B1')*X1+X1*(B1*B1')*X1+(B1'*M+F*C1)'*(B1'*M+F*C1)<a*M];
Sol=solvesdp(K,a)..

It was showing:
Sol =

    solvertime: 0
          info: 'No suitable solver'
       problem: -2
    yalmiptime: 0.2230
 no suitable solver. what shud be the suitable solver.

[Johan Löfberg]

You cannot solve the problem by simply plugging it in as it is stated, since it is bilinear in the variables, and thus you don't have any solver for it. You have not implemented the pseudo-code I gave you. Where is your while-loop!? You have to write those sub-function to setup and solve the bisection problem and the linear SDP problem

[Rahul Kant]

you mean to say that i need to write codes for the functions u mentioned i.e. solvegevp and solvelmi and implement the code u provided.

[Johan Löfberg]

Yes

[Rahul Kant]

Sir I am not able to code the sun functions as you mentioned. I am not getting it. here I have coded:

rt=1.5e-3;
lt=30e-6;
r=76;
l=111.9e-3;
c1=62.855e-6;
f=60;
w0=2*pi*f;
A=[-rt/lt w0 0 -1/lt;w0 -rt/lt -2*w0 ((rt*c1*w0/l)-(w0/r));0 w0 -rt/l ((1/l)-(w0^2*c1));1/c1 0 -1/c1 -1/(r*c1)];
B=[1/lt 0 0 0]';
C=[0 0 0 1];
A1=[A zeros(4,1);C 0];
B1=[B;0];
C1=[C 0;zeros(1,4) eye(1);C*A 0];
Q=eye(5);
[P,L,G] = care(A1,B1,Q);
X1=P;

M=sdpvar(5,5);
F=sdpvar(1,3);
a=sdpvar(1,1);
K=[M>0,A1'*M+M*A1-X1*(B1*B1')*M-M*(B1*B1')*X1+X1*(B1*B1')*X1+(B1'*M+F*C1)'*(B1'*M+F*C1)<a*M];
solvesdp(K,trace(M));
P0 = double(M);
a_lower = -max(eig(inv(P0)*A1'*M+M*A1-X1*(B1*B1')*M-M*(B1*B1')*X1+X1*(B1*B1')*X1+(B1'*M+F*C1)'*(B1'*M+F*C1)))/2
a_upper = a_lower*2;
K=[M>0,A1'*M+M*A1-X1*(B1*B1')*M-M*(B1*B1')*X1+X1*(B1*B1')*X1+(B1'*M+F*C1)'*(B1'*M+F*C1)<a_upper*M];
ops = sdpsettings('verbose',0,'warning',0);
sol = solvesdp(K,[],ops);
while (sol.problem==1)
    a_upper = a_upper*2;
    K=[M>0,A1'*M+M*A1-X1*(B1*B1')*M-M*(B1*B1')*X1+X1*(B1*B1')*X1+(B1'*M+F*C1)'*(B1'*M+F*C1)<a_upper*M];
    sol = solvesdp(K,[],ops);
end
tol = 0.01;
a_works = a_lower
while (a_upper-a_lower)>tol
  a_test = (a_upper+a_lower)/2;
  disp([a_lower a_upper a_test])
 K=[M>0,A1'*M+M*A1-X1*(B1*B1')*M-M*(B1*B1')*X1+X1*(B1*B1')*X1+(B1'*M+F*C1)'*(B1'*M+F*C1)<a_upper*M];
  sol = solvesdp(K,[],ops);
  if sol.problem==1
    a_upper = a_test;
  else
    a_lower = a_test;
    a_works = a_test;
    Mworks = double(M);
 end
end

[Rahul Kant]

 I am not very good at coding. Please help me code it .

[Johan Löfberg]

1. Don't you get a lot off warnings about the fact that you should not use strict inequalities?

2. It fails already in the first call to solvesdp (which you fail to see as you don't look at the diagnostics), since you have created a quadratic SDP inequality. Why have you done that? In your original form, it was linear since it had been Schur complemented compared to the one you use now.

3. Clean up your code into logical steps by writing separate functions which solves the gevp etc, as I recommended to begin with.

I am not going to code it for you as you will learn nothing from that (unless you pay consulting fees)