Why yalmip gives 'Numerical Problems' for a simple inequality

[harun tugal]

Hi,

I cound not figure out why yalmip gives numerical problems when n=2, yet it works fine when n is different from 2 

 

clear('yalmip')

opt = sdpsettings;

opt = sdpsettings(opt,'verbose',0); 

s=tf('s');  n=2; 

for i=1:1:n

Fi(i,1)=1/(s+3)^(i-1);

end

K = sdpvar(n,n,'symmetric');

[A,B,C,D]=ssdata(Fi);

P=sdpvar(size(A,1),size(A,2));

H=[A'*P+P*A P*B; B'*P zeros(size(P*B,2),size(B'*P,1))];

T=[C'*K*C C'*K*D ; D'*K*C D'*K*D];

F=set((H+T) > 0);

solution=solvesdp(F,[],opt); 

report=solution.problem

Thanks,

[Johan Löfberg]

First,

1. Don't use SET, it is obsolete
http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Commands.set

2. Don't use strict inequalities. Didn't you get warnings?
http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Blog.Prepare-your-code

Concerning your question, impossible to answer since you haven't told us which solver you are using. YALMIP does not create any such warnings, the solver does.

[harun tugal]

Hi Johan,

Thanks for suggestions, the solver is Sedumi. 

[Johan Löfberg]

Works here

 solvesdp(H+T>=0)
 
SeDuMi had unexplained problems, maybe due to linear dependence?
YALMIP tweaks the problem (adds 1e6 magnitude bounds on all variables) and restarts...
 
SeDuMi 1.3 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, theta = 0.250, beta = 0.500
eqs m = 4, order n = 11, dim = 13, blocks = 2
nnz(A) = 13 + 0, nnz(ADA) = 16, nnz(L) = 10
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            7.36E+06 0.000
  1 :   0.00E+00 5.51E+05 0.000 0.0748 0.9900 0.9900   1.00  1  1  2.0E+00
  2 :   0.00E+00 5.72E+02 0.000 0.0010 0.9999 0.9999   1.06  1  1  9.4E-01
  3 :   0.00E+00 5.77E-05 0.000 0.0000 1.0000 1.0000   1.00  1  1  5.9E-02
  4 :   0.00E+00 5.35E-10 0.000 0.0000 1.0000 1.0000   1.00  1  1  5.5E-07
  5 :   0.00E+00 4.00E-10 0.122 0.7476 0.9000 0.9000   0.63  4  4  5.3E-07
  6 :   0.00E+00 1.94E-10 0.000 0.4838 0.9000 0.9000   1.00  4  4  1.9E-07
  7 :   0.00E+00 3.78E-11 0.000 0.1955 0.9000 0.9000   0.97 10 11  4.3E-08
  8 :   0.00E+00 1.29E-11 0.000 0.3402 0.9000 0.9000  -0.75 10 11  1.5E-08
  9 :   0.00E+00 5.56E-13 0.000 0.0432 0.9900 0.9900   1.00  9  9  6.7E-10

iter seconds digits       c*x               b*y
  9      0.2   0.0  2.9998304023e-10  0.0000000000e+00
|Ax-b| =   4.8e-16, [Ay-c]_+ =   0.0E+00, |x|=  2.3e-16, |y|=  1.7e+06

Detailed timing (sec)
   Pre          IPM          Post
0.000E+00    1.250E-01    0.000E+00    
Max-norms: ||b||=0, ||c|| = 1000000,
Cholesky |add|=0, |skip| = 2, ||L.L|| = 1.

ans =

    yalmiptime: 0.140000000000001
    solvertime: 0.125000000000000
          info: 'Successfully solved (SeDuMi-1.3)'
       problem: 0

note what it says in the beginning though. The problem is too simple (or more exactly, it has problems with linear dependence) so SeDuMi gets confused. (you get the same response if you try to solve the following trivial problem, which is similiar to what you have)

sdpvar x y
solvesdp(x+y>=0,[],sdpsettings('solver','sedumi'))

You have to supply more infomation

[harun tugal]

Dear Johan,

As you said with the code below finally I get problem: 0

sdpvar x y
solvesdp(x+y>=0,[],sdpsettings('solver','sedumi'))

However, with the initial program I always get report=4 (numerical problems)   

SeDuMi had unexplained problems, maybe due to linear dependence?

YALMIP tweaks the problem (adds 1e6 magnitude bounds on all variables) and restarts...

 

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

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

eqs m = 4, order n = 11, dim = 13, blocks = 2

nnz(A) = 13 + 0, nnz(ADA) = 16, nnz(L) = 10

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

  0 :            7.36E+06 0.000

  1 :   0.00E+00 5.51E+05 0.000 0.0748 0.9900 0.9900   1.00  1  1  2.0E+00

  2 :   0.00E+00 5.72E+02 0.000 0.0010 0.9999 0.9999   1.06  1  1  9.4E-01

  3 :   0.00E+00 5.77E-05 0.000 0.0000 1.0000 1.0000   1.00  1  1  5.9E-02

  4 :   0.00E+00 4.92E-10 0.000 0.0000 1.0000 1.0000   1.00  1  1  5.0E-07

  5 :   0.00E+00 3.58E-10 0.133 0.7280 0.9000 0.9000   0.70  4  4  4.6E-07

  6 :   0.00E+00 1.75E-10 0.000 0.4883 0.9000 0.9000   1.01  4  4  1.7E-07

  7 :   0.00E+00 4.95E-11 0.000 0.2828 0.9000 0.9000   0.91  7 10  5.7E-08

  8 :   0.00E+00 2.61E-11 0.000 0.5285 0.9000 0.9000   0.85 10 10  2.9E-08

  9 :   0.00E+00 6.38E-12 0.000 0.2441 0.9000 0.9000   0.67  9  9  9.4E-09

 10 :   0.00E+00 4.68E-14 0.125 0.0073 0.9990 0.9990   0.96  9  9  7.2E-11

iter seconds digits       c*x               b*y

 10      0.1   0.1  2.7301201619e-10  0.0000000000e+00

|Ax-b| =   5.5e-16, [Ay-c]_+ =   0.0E+00, |x|=  2.3e-16, |y|=  1.8e+06

Detailed timing (sec)

   Pre          IPM          Post

9.998E-03    8.601E-02    4.999E-03    

Max-norms: ||b||=0, ||c|| = 1000000,

Cholesky |add|=0, |skip| = 2, ||L.L|| = 2.3418.

report =

     4

It is initial part of the stability analysis of a system with multipliers. Basically the definition of the program: Let  lambda(w)>0 and it can be defined as lambda(w)=Fi^* K Fi such that K=K^T. And with Kalman-Yakubovic-Lemma frequency depended inequality (FDI)  is transformed to linear matrix inequality (LMI) in order to be able to solve it with SDP.  

[Johan Löfberg]

As I said, this particular instance is simply not well-conditioned, and on your particular setup, SeDuMi gives you a hint about this

BTW, P=K=0 is feasible, so you should really find a better dehomogenized model

And since this is a KYP problem, don't you want P to be positive semidefinite?

[harun tugal]

Dear Johan,

Thank you very much for your particular interest ! I got the point about being not well-conditioned. As you suggest defining K or P as >= 0 solves the problem. 

Thanks indeed. 

[Johan Löfberg]

I think you have implemented something wrong.

From standard stability approaches etc, A stable (which you have) means A'*P+P*A is negative definite if P>0. You on the other hand are trying to get A'P+PA positive semidefinite.

[harun tugal]

I don't want to occupy you  with my problem, yet if it is ok for you I can attach a document that basically describes what I am trying to do  ?

[Johan Löfberg]

sure

[harun tugal]

Thanks a lot ! 

[Johan Löfberg]

As I suspected, you have switched the sign on the inequality.

And since the problem is homogenous in (P,K), you can safely replace the (practically impossible) < with <= -eye()

[Nasim En]

Hi,

the feas number should be a positive number or a negative?

[Johan Löfberg]

I presume you are using SeDuMi, for which the feas indicator converges to 1 for well-posed and numerically stable problems

[Nasim En]

yes, I'm using SeDuMi.

Thank you very much

[Nasim En]

[Nasim En]

Hi,

I followed the attached paper, but  I couldn't get the same result as reported in Example 1. 
I used MATLAB LMI toolbox and YALMIP but both were different.

I'm wondering if you could kindly see my code and help me to modify the LMI.

[Johan Löfberg]

You're solving a feasibility problem, so why would you expect to get the same solution.

BTW, there is no command solvesdp if you read the tutorials

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

and your problem is not well-posed as all variables = 0 is feasible

[Nasim En]

thank you

but in tutorials, the examples are for optimization and there isn't a feasibility problem.
how would the problem be well-posed? 

[Johan Löfberg]

I don't understand the question/claim

[Nasim En]

sorry,
I was confused. 
I'll study the tutorial again.

[Nasim En]

Hello Dear Johan

would you please help me?

why this message is displayed?

 "SeDuMi had unexplained problems, maybe due to linear dependence?

YALMIP tweaks the problem (adds 1e6 magnitude bounds on all variables) and restarts...

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

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

eqs m = 3716, order n = 9815, dim = 16967, blocks = 15

nnz(A) = 23383 + 0, nnz(ADA) = 2655948, nnz(L) = 3832046"

[Johan Löfberg]

Ill-posed problem, badly scaled problem, memory problems. sedumi simply crashed for some unknown reason

If you want to investigate further, you will have to put a debug break on line 41 in callsedumi and see what happens when you run that line in the command prompt

[Nasim En]

thank you so much

[Nasim En]

Hello Dear Johan

would you please help me?

why this message is displayed?

"Warning: Solver not applicable (sedumi)"

[Johan Löfberg]

You've formulated a problem which SeDuMi cannot solve (i.e., you do not have a linear SOCP/SDP-representable problem)

[Nasim En]

thank you so much

[Nasim En]

Dear Johan

I have tried many ways but the problem couldn't be feasible.  I want to make sure if something is wrong with my code or the LMI is ill conditioned. I appreciate your help.
This is the message.

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 = 395, order n = 1644, dim = 1670, blocks = 3

nnz(A) = 3569 + 0, nnz(ADA) = 3927, nnz(L) = 2438

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

  0 :            7.48E+04 0.000

  1 :   0.00E+00 2.49E+04 0.000 0.3335 0.9000 0.9000  -3.86  1  1  1.7E+04

  2 :   0.00E+00 9.73E+03 0.000 0.3902 0.9000 0.9000  -1.60  1  1  1.9E+04

  3 :   0.00E+00 2.62E+03 0.000 0.2692 0.9000 0.9000  -1.22  1  1  1.9E+04

  4 :   0.00E+00 1.58E+03 0.000 0.6048 0.9000 0.9000  -1.06  1  1  2.2E+04

  5 :   0.00E+00 5.00E+02 0.000 0.3154 0.9000 0.9000  -1.04  1  1  3.9E+04

  6 :   0.00E+00 2.55E+02 0.000 0.5105 0.9000 0.9000  -1.01  1  1  4.3E+04

  7 :   0.00E+00 1.12E+02 0.000 0.4409 0.9000 0.9000  -1.01  1  1  4.4E+04

  8 :   0.00E+00 4.78E+01 0.000 0.4250 0.9000 0.9000  -1.00  1  1  4.5E+04

  9 :   0.00E+00 1.98E+01 0.000 0.4141 0.9000 0.9000  -1.00  1  2  4.6E+04

 10 :   0.00E+00 7.98E+00 0.000 0.4032 0.9000 0.9000  -1.00  1  4  4.6E+04

 11 :   0.00E+00 3.14E+00 0.000 0.3930 0.9000 0.9000  -1.00  1  6  4.7E+04

 12 :   0.00E+00 1.19E+00 0.000 0.3811 0.9000 0.9000  -1.00  1  6  4.8E+04

 13 :   0.00E+00 4.44E-01 0.000 0.3713 0.9000 0.9000  -1.00  2 13  4.8E+04

 14 :   0.00E+00 1.61E-01 0.000 0.3622 0.9000 0.9000  -1.00  1 16  4.9E+04

 15 :   0.00E+00 5.68E-02 0.000 0.3537 0.9000 0.9000  -1.00  1 17  4.9E+04

 16 :   0.00E+00 1.97E-02 0.000 0.3458 0.9000 0.9000  -1.00  1 18  5.0E+04

 17 :   0.00E+00 6.65E-03 0.000 0.3382 0.9000 0.9000  -1.00  1 19  5.0E+04

 18 :   0.00E+00 2.22E-03 0.000 0.3338 0.9000 0.9000  -1.00  1 20  5.0E+04

 19 :   0.00E+00 7.13E-04 0.000 0.3215 0.9000 0.9000  -1.00  1 21  5.1E+04

 20 :   0.00E+00 2.28E-04 0.000 0.3193 0.9000 0.9000  -1.00  1 23  5.1E+04

 21 :   0.00E+00 7.14E-05 0.000 0.3137 0.9000 0.9000  -1.00  1 24  5.1E+04

 22 :   0.00E+00 2.20E-05 0.000 0.3084 0.9000 0.9000  -1.00  1 25  5.1E+04

 23 :   0.00E+00 6.69E-06 0.000 0.3035 0.9000 0.9000  -1.00  1 26  5.2E+04

 24 :   0.00E+00 1.98E-06 0.000 0.2958 0.9000 0.9000  -1.00  1 48  5.2E+04

 25 :   0.00E+00 5.93E-07 0.000 0.2997 0.9000 0.9000  -1.00  1 14  5.2E+04

 26 :   0.00E+00 1.70E-07 0.000 0.2860 0.9000 0.9000  -1.00  1 24  5.3E+04

 27 :   0.00E+00 4.77E-08 0.000 0.2815 0.9000 0.9000  -1.00  1 18  5.3E+04

 28 :   0.00E+00 3.43E-08 0.000 0.7182 0.9000 0.9000  -1.00  4 18  5.5E+04

 29 :   0.00E+00 3.00E-08 0.000 0.8744 0.9000 0.9000  -1.00  8 30  5.5E+04

Run into numerical problems.

Dual infeasible, primal improving direction found.

iter seconds  |Ax|    [Ay]_+     |x|       |y|

 29      0.7   1.4e-09   4.9e-13   4.4e+00   5.2e+01

Detailed timing (sec)

   Pre          IPM          Post

0.000E+00    2.960E-01    0.000E+00    

Max-norms: ||b||=0, ||c|| = 1,

Cholesky |add|=20, |skip| = 0, ||L.L|| = 2291.67.

ans = 

  struct with fields:

    yalmipversion: '20181012'

       yalmiptime: 0.0714

       solvertime: 0.3026

             info: 'Infeasible problem (SeDuMi-1.3)'

          problem: 1

[Johan Löfberg]

Without reproducible code all one can say that the problem appears infeasible

[Nasim En]

I will send you my code.

Thank you so much.

[Nasim En]

Hello Dear Johan

would you please help me?

in some cases, LMI is successfully solved and it's feasible but unfortunately placing the specified gains in the simulated model can not stabilize the system. what is the reason of this problem?

[Johan Löfberg]

Impossible to answer without the model.