Run into numerical problems

Anh Lam Do

unread,

Feb 8, 2016, 10:33:39 AM2/8/16

to YALMIP

Dear Yohan Löfberg,

I implemented the below code to optimize the Hinf norm of a closed-loop systems (under some constraints) (with Yalmip+Sedumi). And what I obtained is the warning "Run into numerical problem". Cound you please say if the model is well-posed? And how to avoid the "Run into numerical problem"?

Indeed, when I check a gain the conditions with the obtained cost value "alpha" (minimal), I see that the constraint "X>0" is not satisfied.

I would prefer have a bigger minimum 'alpha' (not optimal) than seeing this kind of problem.

Thank you very much in advance.

Kind regards,

Lam

Code:

----------------------------------------------------------

A =[-7.1379   -0.8860         0         0
    22.3116   -7.6465         0         0
         0    1.0000          0         0
   19.4444   10.0000    19.4444         0];

B= [3.4687
   39.6043
         0
         0];

Bw=[-0.5045
         0
         0
         0];

Cz =[ 0 0 0 1];

H11 =[ 0 0 0 1];

h011 = 0.3000;

sizeN=size(A,1);

sizeU=size(B,2);

Q = sdpvar(sizeN,sizeN,'symmetric');

R = sdpvar(sizeU,sizeN,'full');

gamma=sdpvar(1,1,'full');

% constraint 1

M=[Q*A'+A*Q+R'*B'+B*R Q*Cz' Bw;

Cz*Q -1 0;

Bw' 0 -gamma]; %M<0

% Constraint 2

X11 = Q;

X12 = (Q*A'+R'*B')*H11';

X21 = X12';

X22 = h011*h011;

X = [X11 X12;X21 X22]; %X>0

% Constraint 3

alpha1=0.1;

alpha2=10;

W1=A*Q+Q*A'+B*R+R'*B'+2*alpha1*Q; % W1<0

W2=A*Q+Q*A'+B*R+R'*B'+2*alpha2*Q; % W2>0

F = [Q>0,M<0,X>0,W1<0,W2>0]

%%% Find feasible solution minimizing gamma

ops = sdpsettings('warning',1,'verbose',1,'solver','sedumi');

solution = solvesdp(F,gamma,ops);

------------------------End of code-----------------------------

Results: Run into numerical problems!! (see below)

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 = 15, order n = 24, dim = 110, blocks = 6
nnz(A) = 160 + 0, nnz(ADA) = 225, nnz(L) = 120
it :     b*y       gap    delta rate   t/tP* t/tD*   feas cg cg prec
0 :            2.88E+002 0.000
1 : -2.15E-001 7.69E+001 0.000 0.2673 0.9000 0.9000   1.51 1 1 3.5E+001
2 : -4.18E-001 2.49E+001 0.000 0.3231 0.9000 0.9000   0.95 1 1 1.3E+001
3 : -7.61E-001 6.77E+000 0.000 0.2724 0.9000 0.9000   0.37 1 1 5.5E+000
4 : -1.21E+000 2.23E+000 0.000 0.3297 0.9000 0.9000   0.07 1 1 3.2E+000
5 : -1.85E+000 8.04E-001 0.000 0.3602 0.9000 0.9000   0.12 1 1 1.9E+000
6 : -2.52E+000 3.19E-001 0.000 0.3962 0.9000 0.9000   0.06 1 1 1.3E+000
7 : -3.53E+000 1.28E-001 0.000 0.4026 0.9000 0.9000   0.20 1 1 8.0E-001
8 : -4.41E+000 4.90E-002 0.000 0.3821 0.9000 0.9000 -0.08 1 1 6.5E-001
9 : -6.28E+000 1.70E-002 0.000 0.3470 0.9000 0.9000   0.14 1 1 3.5E-001
10 : -8.14E+000 5.53E-003 0.000 0.3248 0.9000 0.9000 -0.02 1 1 2.4E-001
11 : -1.06E+001 1.89E-003 0.000 0.3415 0.9000 0.9000   0.05 1 2 1.4E-001
12 : -1.30E+001 6.89E-004 0.000 0.3654 0.9000 0.9000 -0.07 1 2 1.0E-001
13 : -1.69E+001 2.39E-004 0.000 0.3470 0.9000 0.9000   0.08 1 2 5.8E-002
14 : -2.03E+001 8.79E-005 0.000 0.3673 0.9000 0.9000 -0.03 1 2 4.1E-002
15 : -2.52E+001 3.05E-005 0.000 0.3466 0.9000 0.9000   0.08 1 5 2.3E-002
16 : -2.94E+001 1.17E-005 0.000 0.3829 0.9000 0.9000   0.00 2 6 1.6E-002
17 : -3.54E+001 3.90E-006 0.000 0.3346 0.9000 0.9000   0.14 1 7 8.4E-003
18 : -4.00E+001 1.48E-006 0.000 0.3791 0.9000 0.9000   0.10 1 7 5.3E-003
19 : -4.57E+001 4.52E-007 0.000 0.3057 0.9000 0.9000   0.27 1 8 2.4E-003
20 : -4.93E+001 1.59E-007 0.000 0.3526 0.9000 0.9000   0.34 1 8 1.2E-003
21 : -5.24E+001 3.55E-008 0.000 0.2223 0.9000 0.9000   0.58 1 8 3.4E-004
22 : -5.34E+001 2.62E-009 0.000 0.0739 0.9900 0.9900   0.79 1 11 2.8E-005
23 : -5.34E+001 6.37E-010 0.000 0.2430 0.9000 0.9000   0.97 16 19 7.1E-006
Run into numerical problems.

iter seconds digits c*x b*y
23 0.3 4.1 -5.3447794086e+001 -5.3443066471e+001
|Ax-b| = 1.1e-005, [Ay-c]_+ = 8.8E-008, |x|= 5.1e+005, |y|= 1.1e+004

Detailed timing (sec)
Pre IPM Post
5.801E-002 2.340E-001 2.200E-002
Max-norms: ||b||=1, ||c|| = 1.009029e+000,
Cholesky |add|=0, |skip| = 3, ||L.L|| = 10.7835.

Johan Löfberg

unread,

Feb 8, 2016, 11:03:11 AM2/8/16

to YALMIP

If you look at the solution, you'll see that the optimal solution tends to a solution where a row/column of Q is zero, and thus optimal tends to singular.

Standard trick is to solve the problem tro get the optimal value, and then solve a feasibility problem where you relax required optimality, and use the fact that interior-point solvers typically return well behaved interior solutions (analytic centers or similar) for pure feasibility problems

solution = solvesdp(F,gamma,ops);
solution = solvesdp([F,gamma <= 1.05*value(gamma)],[],ops);

BTW, you cannot have strict inequalities in your model. YALMIP is telling you that

Anh Lam Do

unread,

Feb 12, 2016, 10:21:55 AM2/12/16

to YALMIP

Thank you very much for your solution. It really helps me.

Kind regards,

Lam

Viet Long

unread,

Oct 23, 2018, 3:40:35 PM10/23/18

to YALMIP

Dear Prof. Johan Löfberg,

I have a basic question about your answer, we have the LMI conditions follows as:

#=================================================================#

F = [Q>0,M<0,X>0,W1<0,W2>0];

solution = solvesdp(F,gamma,'solver','sedumi');

#=================================================================#

And we get the information of the solver function:

#=================================================================#

But how can you quickly assert such an "If you look at the solution, you'll see that the optimal solution tends to a solution where a row/column of Q is zero", can you please just teach me to avoid numerical error of sedumi solver in future.

Thank you for your consideration.

Johan Löfberg

unread,

Oct 23, 2018, 3:45:18 PM10/23/18

to YALMIP

First, make sure you have reasonable numbers (i.e. not 10^-15 and 10^12 etc) . Crap in crap out. If you use the latest version of YALMIP, coefficients are displayed for you

Here, the first constraint is a numerical disaster

>> Model = [10^30*x + 1e-12*y >= 1, 2*x + 3*y>=0]
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|                    Constraint|   Coefficient range|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Element-wise inequality 1x1|      1e-12 to 1e+30|
|   #2|   Element-wise inequality 1x1|              2 to 3|
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Second, accept that sedumi runs into numerical problems sometimes. Mosek is typically more numerically robust

and note that strict inequalities are not supported. YALMIP screams at you if you use this

Reply all

Reply to author

Forward