cant get result Yalmip

[Chokri Sandy]

Hello friends, the only result I get is NaN, any help is appreciated

Thank you

 

 

 

 

clear

all;

clc

%% Define Matrix

A1 = [0 1; 26.3679 -0.2839]; B1 = [0 -0.2633]'; C1 = [1 0];

A2 = [0 1; 25.7826 -0.2833]; B2 = [0 -0.2469]'; C2 = [1 0];

A3 = [0 1; 24.1023 -0.2818]; B3 = [0 -0.2002]'; C3 = [1 0];

A4 = [0 1; 21.5219 -0.2802]; B4 = [0 -0.1299]'; C4 = [1 0];

Ap = [0 0.0050; 0.1076 0.0014]; Bp = [0 0.0039]'; Cp = [0 0];

Q11 = [0.01 0; 0 0.01]; Q22 = [0.01 0; 0 0.01];

%% Declaration of variables

P22 = sdpvar(2,2);

Z1 = sdpvar(2,1); Z2 = sdpvar(2,1); Z3 = sdpvar(2,1); Z4 = sdpvar(2,1);

rho = sdpvar(1,1);

%% Formulation of the LMI matrix

U11 = [P22*A1+A1'*P22-Z1*C1-C1'*Z1'+Ap'*Ap+Q22 P22;

P22 -(3+(1/rho))^-1*eye(2)];

U22 = [P22*A2+A2'*P22-Z2*C2-C2'*Z2'+Ap'*Ap+Q22 P22;

P22 -(3+(1/rho))^-1*eye(2)];

U33 = [P22*A3+A3'*P22-Z3*C3-C3'*Z3'+Ap'*Ap+Q22 P22;

P22 -(3+(1/rho))^-1*eye(2)];

U44 = [P22*A4+A4'*P22-Z4*C4-C4'*Z4'+Ap'*Ap+Q22 P22;

P22 -(3+(1/rho))^-1*eye(2)];

F = [P22>zeros(2)];

F = F+[U11<zeros(4)];

F = F+[U22<zeros(4)];

F = F+[U33<zeros(4)];

F = F+[U44<zeros(4)];

ops = sdpsettings('verbose',1,'warning',0,'showprogress','1','dimacs','1');

 

 

sol=solvesdp(F,rho,ops);

solver_time=sol.solvertime;

P22 = double(P22);rho=double(rho);

Z1=double(Z1);Z2=double(Z2);Z3=double(Z3);Z4=double(Z4);

L1 = inv(P22)*Z1;L2 = inv(P22)*Z2;L3 = inv(P22)*Z3;L4 = inv(P22)*Z3;

[Johan Löfberg]

1. Strict inequalities is not supported (typically ill-posed problem since there exist no minimizer only an infimizer)
http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Blog.Prepare-your-code

2. there are no public solvers which will solve your problem, since it nasty nonlinear nonconvex terms involving rho.

However, why the nasty parametrization? Introduce a new variable alpha to replace (3+(1/rho))^-1, and then minimize alpha instead (the nasty function in rho is monotonic). Once you have the optimal alpha, you can easily compute the corresponding rho

[Chokri Sandy]

thank Mr Johan, once I modified the cod like you recommended, I found a solution equal to rho= 0.0096 but using matlab lmi toolbox the solution is 0.6084 and using cvx is infeasable ?. One more question, is it possible to solve a set of linear matrix inequality without objective function . Thank you again for your help
best regards 

[Johan Löfberg]

The fact you have such different solutions indicates that either the problem is horribly ill-conditioned, or you have implemented something wrong. Perhaps you code post your code and data

Of course you can solve a problem without any objective. A solver will then just return some feasible solution. Which one depends on the algorithm

[Chokri Sandy]

Dear Johan, attached to this post three versions of the code for the same problem , using LMI matlab tool box, YALMIP and CVX  I got different result in YALMIP then in LMI toolbox and no feasible solution in CVX Please help

 

 

 

[Johan Löfberg]

Why on earth are you posting code in a word document, or a pdf document?? Code is pure text, so why don't you just attach the .m file, or, if you for some reason don't want to post a .m file, in a .txt file.

[Chokri Sandy]

Sorry I am just not thinking anymore, attached are the .m files.  The author claim that using LMI toolbox he got  rho^2 = 0.0684 and P22 = [1.9318  -0.2116;  -0.2116 0.2038] . Thank you again for all your help.

Best regards

[Johan Löfberg]

The cvx code fails since you explicitly have to tell CVX that you are defining an sdp by cvx_begin sdp

Note, in the YALMIP code, you don't have to write U11 <= zeros(4) etc, U11 <= 0 is sufficient. Why are you turning on dimacs error measures? and why are you turning off warnings?

Why lmilab fails: You made an error somewhere. Where? Well, that is exactly why we don't want to use lmilab. Messy to write and debug

[Johan Löfberg]

The optimal solution is  rho2^2

ans =

   9.1466e-05

Either you have implemented the model wrong (wrong data or something like that) or the author has implemented something wrong, or lmilab simply sucks and returns a far from optimal solution. Your YALMIP code is correct and implements the model you state.

Something is wrong somewhere though, with the claimed solution

F = [F, P22 == [1.9318  -0.2116;  -0.2116 0.2038]];

The optimal rho^2 is0.1034, which is different from the claimed value

[Chokri Sandy]

Dear Mr Johan

I agree LMI toolbox sucks I used it to check if I can get the same results as the author did, can you tell me please how fix CVX to make it run. Also how did you get rho^2 =0.1034 and the optimal objective value alpha =9.1466e-05 using my code and I am getting rho^2 =0.0096 and the optimal objective value alpha =9.2970e-03

I am confused

Thank you for your help

best regards

[Mark L. Stone]

In CVX, you can either invoke the CVX session with cvx_begn sdp as Yohan stated, and then all inequalities will be interpreted as LMI where eligible, or if invoked as cvx_begin, for each LMI instance, you write the LHS which is >= 0 as LMI in CVX as LHS == semidefinite(n) , where n is the dimension of the matrix  Your success with CVX will improve if you read the CVX User's Guide at http://cvxr.com/cvx/doc/CVX.pdf .

[Johan Löfberg]

I guess I got confused by your notation in the variable. The optimal value of the variable you call rho2 (which I gues means rho^2)is 0.0096 as you say.


The value 0.1034 was obtained when I forced YALMIP to find the P22 the "author" claimed was part of the optimal solution. This was just to show that the claimed solution does not give the value 0.0684, indicating there is either an error in your YALMIP code, or in the numbers reported in the paper by the authors. What paper are you trying to reproduce?

[Chokri Sandy]

Dear Mr Johan, Thank you for your quick reply and your help, attached is the paper I am trying to reproduce, the LMI system I am solving are equation [73], I need to be sure that the result are the same that the author claim so I can use them in my research ( Fuzzy control of flexible spacecraft). the initial data the author used are on page (582), also attached is the .m file for the whole problem.

 

Thank you for your support and your help

 

note : the author claim that     rho^2=0.6084   not as I mentioned in my last email (0.0684)

 

 

[Johan Löfberg]

The solution has to satisfy both 72 and 70. You have only implemented the constraints related to 72. Hence, your solution is most likely infeasible for the constraints in 70. Since the problem is nonconvex in all the variables, an iterative scheme such as the one proposed in the paper has to be used (basically fix rho, solve 72 to find P and Z, use that in 70 to find W. Decrease rho untill you cannot solve the two problems)

[Chokri Sandy]

regardless of what value of rho^2 ( 5....0.1), the first set have a feasible solution but the second set infeasible. there must be a way for Yalmip or sedumi to automatically find rho such that there is a feasible solution, on the website you mentioned that relaxation can be a way for non convex optimization , is it possible to relax the constrains in this case ? 

best regards

[Johan Löfberg]

If you remove the constraint on V44, the problem is feasible. This might indicate that there is something wrong with V44. Looking around, you definition of L4 is wrong. Still infeasible after that fix, but it shows that there are simple flaws in your code.

Regarding relaxation based methods etc. I would say that most likely is completely intractable here. A simple local iterative scheme as proposed is probably as good as it gets, unless you come up with a linearizing change of variables etc.

[Johan Löfberg]

Unfortunately, the proposed algorithm is pretty poor, since the success of the second step depends on the solution computed in step 1.

As you do it know where you search for any feasible P22, the obtained P22 tends to lead to an infeasible problem in the second step when SeDuMi is used. When the authors used LMILAB, it happened to work with the solution which LMILAB returned down to rho^2=0.6084.

However, if we tell the solver to return not just any P22, but the P22 closest to identity

sol=solvesdp(F,norm(P22-eye(2)),ops);

then the solution fits better for the second problem, and rho^2 can be selected below 0.6084 (0.4 is feasible)

Why norm(P22-I)? Because it worked...

[Chokri Sandy]

Dera Mr Löfberg

is it possible to impose (either in the constrain F = [ ...<=0...] or in the solver  sol=solvesdp(F,alpha,ops);) that the solution P22 have only positive eigvalues , if so how to do that ?

[Johan Löfberg]

P22 >= epsilon*eye(n), where you pick epsilon suitably?