Error: Sedumi 1.3 on YALMIP: Dual infeasible, primal improving direction found

[Allan André do Nascimento]

Dear all,

I am new to both sedumi1.3 and Yalmip. I have been trying to run the following code:

clc
clear all


m1 = 290;
m2 = 59;
k1 = 16812;
k2 = 190000;
b1 = 1000;
alfa = 4.515e+13;
beta = 1;
gama = 1.545e+9;
tau = (1/30);
Ps=10342500;
A = 3.35e-4;


Aa1 = [0         1          0            0        0      0;        
   -k1/m1  -(b1)/m1    (k1/m1)      (b1/m1)   A/m1       0;        
    0         0          0            1        0         0;        
    k1/m2    b1/m2    -(k1+k2)/m2  -b1/m2    -A/m2       0;        
    0      -alfa*A       0          alfa*A   -beta   gama*sqrt(1.5*Ps); %0.36 e 0.8
    0         0          0            0        0        -1/tau ];


Aa2 = [0         1          0            0        0         0;        
      -k1/m1  -(b1)/m1    (k1/m1)      (b1/m1)   A/m1       0;        
       0         0          0            1        0         0;        
       k1/m2    b1/m2    -(k1+k2)/m2  -b1/m2    -A/m2       0;        
       0      -alfa*A       0          alfa*A   -beta  gama*sqrt(0.5*Ps);
       0         0          0            0        0        -1/tau ];
   




 Ba =[0                
    0                
    0                
    0              
    0                
    1/tau];


v=[-k1/m1  -(b1)/m1    (k1/m1)      (b1/m1)   (A/m1)       0];
Qi=[1 0 0 0 0 0;0 10 0 0 0 0;0 0 1 0 0 0;0 0 0 10 0 0;0 0 0 0 1e-11 0;0 0 0 0 0 10];%1e-11 0.5/1.5 %1e-12
Qb=((75)*(v'*v))+Qi;   %75 para 0.5/1.5      %Linha5 Coluna5colocar 1e-12 e coef qb=70 com 1.36/0.8 nas matrizes A1 e A2  
[V,D]=eig(Qb);
DD=sqrt(D);
%Vi=inv(V);
SQ=V*DD/V;%eigenvalue decomposition
R=(1e-2);
SR=sqrt(R);
Aa3=Aa1+Aa2;
zz=zeros(6,1);
yalmip('clear');
X = sdpvar(6);
Z = sdpvar(6);
M1=sdpvar(1,6);
M2=sdpvar(1,6);
F = [[((Aa1*X)+(X*Aa1')-(Ba*M1)-(M1'*Ba')),X*SQ,M1'*SR,M2'*SR;SQ*X,-eye(6),zz,zz;SR*M1,zz',-1,0;SR*M2,zz',0,-1]<0,...
     [((Aa2*X)+(X*Aa2')-(Ba*M1)-(M1'*Ba')),X*SQ,M1'*SR,M2'*SR;SQ*X,-eye(6),zz,zz;SR*M1,zz',-1,0;SR*M2,zz',0,-1]<0,...
     [((Aa3*X)+(X*Aa3')-(Ba*M1)-(Ba*M2)-(M1'*Ba')-(M2'*Ba')),X*SQ,M1'*SR,M2'*SR;SQ*X,-eye(6),zz,zz;SR*M1,zz',-1,0;SR*M2,zz',0,-1]<0,...
     [Z,eye(6);eye(6),X]>0,...
     X>0];
 des=trace(Z);
 sol=solvesdp(F,des);
Xf=double(X);
Pf= inv(Xf);
M11=double(M1);
M22=double(M2);
K1=M11*Pf;
K2=M22*Pf;


eigenproof=eig(Pf)

My problem is with the message I am receiving - "Run into numerical problems / Dual infeasible, primal improving direction found . I have tried to scale the states with a transformation matrix but I still run into numerical problems.

I decided then to test the feasibility of the closed loop system (with the K values found by this algorithm) with the Matlab LMI toolbox, which turned out to be feasible (t approximately -1e-14 and all eigenvalues of Pf positive).

First question: Does anyone have any clue on how to solve it, i mean make the problem more treatable numerically (if so, could you please explain in a very clear manner, as I am new to everything).

Second question: Can I trust this matlab toolbox answer (LMIs are feasible)? I would say I can because, it seems that the solution is almost unfeasible but matlab uses somehow more "relaxed constraints" whereas sedumi uses some slightly stringent constraints to my problem ( I have also tried to change the defaul settings but it is was unsucessful). If not, why?

Thank you all very much.

Best regards,

Allan    

[Johan Löfberg]

It could very well be that the solution is feasible. SeDuMi doesn't like the numbers though, as the optimal value is so enormously large, it thinks it is about to diverge to infinity, i.e., unbounded, and gives up and claims unbounded/infeasible dual

By naively scaling the objective with a small number, you can trick sedumi to return other diagnostics, and also this time return a feasible solution. A typical sign of a numerically poor model

>> sol=solvesdp(F,des*1e-6,sdpsettings('solver','sedumi'))

SeDuMi 1.32 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.

Alg = 2: xz-corrector, theta = 0.250, beta = 0.500

eqs m = 54, order n = 61, dim = 769, blocks = 6

nnz(A) = 879 + 0, nnz(ADA) = 2412, nnz(L) = 1233

 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec

  0 :            3.22E+25 0.000

  1 :  -1.14E-05 8.55E+24 0.000 0.2652 0.9000 0.9000   2.60  5  1  1.3E+13

  2 :  -2.74E-05 2.61E+24 0.000 0.3050 0.9000 0.9000   0.62  9  8  5.3E+12

  3 :  -5.32E-05 8.88E+23 0.000 0.3409 0.9000 0.9000   0.30  5  8  2.7E+12

  4 :  -9.31E-05 3.24E+23 0.000 0.3650 0.9000 0.9000   0.10  8  5  1.6E+12

 84 :  -2.85E+00 9.06E-06 0.000 0.9587 0.9000 0.9000  -0.63  4  1  1.6E-02

 85 :  -2.96E+00 8.74E-06 0.184 0.9648 0.9000 0.9000  -0.67  4  1  1.6E-02

 86 :  -3.18E+00 8.58E-06 0.435 0.9808 0.9000 0.9000  -1.79  4  1  2.0E-02

Run into numerical problems.

iter seconds digits       c*x               b*y

 86      2.4  -1.4  7.5223927317e+01 -3.1836328531e+00

|Ax-b| =   1.7e-02, [Ay-c]_+ =   1.3E+00, |x|=  3.7e+05, |y|=  9.3e+10

No sensible solution found.

Detailed timing (sec)

   Pre          IPM          Post

8.006E-03    1.754E+00    1.006E-03    

Max-norms: ||b||=1.000000e-06, ||c|| = 2,

Cholesky |add|=2, |skip| = 1, ||L.L|| = 6.51893e+09.

sol = 

    yalmiptime: 0.0779

    solvertime: 1.7641

          info: 'Numerical problems (SeDuMi-1.3)'

       problem: 4

>> check(F)

 

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|   ID|          Constraint|   Primal residual|   Dual residual|

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

|   #1|   Matrix inequality|                 0|     -1.3194e-16|

|   #2|   Matrix inequality|                 0|     -1.0959e-16|

|   #3|   Matrix inequality|                 0|     -6.7492e-17|

|   #4|   Matrix inequality|        1.9146e-06|      3.5461e-12|

|   #5|   Matrix inequality|        3.0078e-06|      5.5001e-12|

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

[Allan André do Nascimento]

Thanks Mr. Lofberg, always very helpful!

A good feasible solution is all I need. Is the matlab proof enough?

[Johan Löfberg]

Look at the solution. Does it make sense. Can you see something from it. For instance, if you look at X, it tells you a story. Why are almost all elements zero, except one single diagonal element. Maybe that tells you how you should rescale the system, or that you have uncontrollable states, or...

Is the MATLAB proof enough? Well, if it satisfies your needs, then you should be satisfied. However, you'll have to judge it from it's practical use. Look at the solution. It is not uncommon that people just compute solutions and report them, but never actually look at them from a practical point of view. Does a control law with numbers in the range 10^15 make sense, if you actually intend to use it, and similar stuff.

[Allan André do Nascimento]

Excellent, will do it Mr. Lofberg! Thank you very much once again.