Scaling for a simple set of LMI inequalities

Hi, I have a simple set of LMI’s for a feedback control problem that I have trouble solving due to numerics. I can’t seem to find a good way to scale the problem to make it easier for the solver. I am using Mosek. I have tried using slack variables which seemed help for some combinations of k1, k2, c1, c2, but I run into problems with the k1, k2, c1, and c2 below (mostly due to k2 and c2). Thank you for your help.

Code:

alpha = 0.01; %decay rate design variable

A = @(ks) [zeros(2,3);0,-ks,0];
B = @(cs) [eye(2);[0,-cs]];

k1 = 300;
k2 = 2.7e6;
c1 = 0.35;
c2 = 1000;

cvx_begin sdp
    variable Q(3,3) symmetric
    variable Y(2,3)
    
    A(k1)*Q + Q*A(k1)' + B(c1)*Y + Y'*B(c1)' + 2*alpha*Q <= -eye(3);
    A(k1)*Q + Q*A(k1)' + B(c2)*Y + Y'*B(c2)' + 2*alpha*Q <= -eye(3);
    A(k2)*Q + Q*A(k2)' + B(c1)*Y + Y'*B(c1)' + 2*alpha*Q <= -eye(3);
    A(k2)*Q + Q*A(k2)' + B(c2)*Y + Y'*B(c2)' + 2*alpha*Q <= -eye(3);

    Q >= eye(3);
cvx_end

K = Y/Q %feedback controller

Mosek output:

Calling Mosek_2 9.2.40: 30 variables, 12 equality constraints
   For improved efficiency, Mosek_2 is solving the dual problem.
------------------------------------------------------------

MOSEK Version 9.2.40 (Build date: 2021-3-10 14:33:40)
Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com
Platform: Windows/64-X86

Problem
  Name                   :                 
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 12              
  Cones                  : 0               
  Scalar variables       : 0               
  Matrix variables       : 5               
  Integer variables      : 0               

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 1                 time                   : 0.00            
Lin. dep.  - tries                  : 1                 time                   : 0.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.03    
Problem
  Name                   :                 
  Objective sense        : min             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 12              
  Cones                  : 0               
  Scalar variables       : 0               
  Matrix variables       : 5               
  Integer variables      : 0               

Optimizer  - threads                : 6               
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 12
Optimizer  - Cones                  : 0
Optimizer  - Scalar variables       : 0                 conic                  : 0               
Optimizer  - Semi-definite variables: 5                 scalarized             : 30              
Factor     - setup time             : 0.00              dense det. time        : 0.00            
Factor     - ML order time          : 0.00              GP order time          : 0.00            
Factor     - nonzeros before factor : 72                after factor           : 72              
Factor     - dense dim.             : 0                 flops                  : 3.20e+03        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   4.0e+00  2.7e+06  7.2e+00  0.00e+00   -8.160000000e+00  0.000000000e+00   1.0e+00  0.03  
1   1.2e+00  8.1e+05  1.4e+01  -2.46e+00  -8.695714696e+01  0.000000000e+00   3.0e-01  0.08  
2   1.1e-01  7.2e+04  2.2e+00  -2.18e+00  -1.477237605e+02  0.000000000e+00   2.7e-02  0.08  
3   9.8e-03  6.6e+03  1.5e+00  -1.17e+00  -6.947140589e+03  0.000000000e+00   2.4e-03  0.09  
4   1.5e-03  9.9e+02  5.4e-01  -1.00e+00  -4.291752444e+04  0.000000000e+00   3.7e-04  0.09  
5   3.8e-04  2.5e+02  2.8e-01  -1.10e+00  -1.749550417e+05  0.000000000e+00   9.4e-05  0.09  
6   2.9e-05  1.9e+01  4.6e-02  -7.03e-01  -8.031846519e+05  0.000000000e+00   7.1e-06  0.09  
7   4.5e-06  3.1e+00  1.7e-02  -9.29e-01  -4.638722612e+06  0.000000000e+00   1.1e-06  0.09  
8   9.2e-07  6.2e-01  8.7e-03  -9.65e-01  -2.771301660e+07  0.000000000e+00   2.3e-07  0.09  
9   2.1e-07  1.4e-01  4.6e-03  -1.01e+00  -1.514640916e+08  0.000000000e+00   5.2e-08  0.09  
10  4.4e-08  3.0e-02  1.9e-03  -8.66e-01  -5.714168838e+08  0.000000000e+00   1.1e-08  0.11  
11  1.4e-08  9.8e-03  1.3e-03  -1.16e+00  -2.591403085e+09  0.000000000e+00   3.6e-09  0.11  
Optimizer terminated. Time: 0.13    


Interior-point solution summary
  Problem status  : DUAL_INFEASIBLE
  Solution status : DUAL_INFEASIBLE_CER
  Primal.  obj: -3.0762116129e+00   nrm: 1e+00    Viol.  con: 4e-02    barvar: 0e+00  
Optimizer summary
  Optimizer                 -                        time: 0.13    
    Interior-point          - iterations : 11        time: 0.11    
      Basis identification  -                        time: 0.00    
        Primal              - iterations : 0         time: 0.00    
        Dual                - iterations : 0         time: 0.00    
        Clean primal        - iterations : 0         time: 0.00    
        Clean dual          - iterations : 0         time: 0.00    
    Simplex                 -                        time: 0.00    
      Primal simplex        - iterations : 0         time: 0.00    
      Dual simplex          - iterations : 0         time: 0.00    
    Mixed integer           - relaxations: 0         time: 0.00    

------------------------------------------------------------
Status: Infeasible
Optimal value (cvx_optval): +Inf

You haven’'t told us what alpha is. But it clashes with MATLAB’s alpha.

You haven’t shown us the CVX output. So we don’t know whether Mosek was provided the dual, and what status CVX reported. If CVX sent the dual, Mosek’s status means the problem is primal infeasible. If so, read https://yalmip.github.io/debugginginfeasible . if not, read https://yalmip.github.io/debuggingunbounded .

And try installing the most recent Mosek 10., and select that version with cvx_solver.

As for the scaling, can you divide all the input numbers by 10000, without getting into trouble with the 1 elements in` eye? Not obvious to me you can, but I leave any reformulation to your determination.

You haven’'t told us what alpha is. But it clashes with MATLAB’s alpha.

alpha is given in the first line of the code. It is set as 0.01. It is a scalar>0. It is a design parameter that sets the slowest decay rate of the closed loop system A+B*K for all combinations of A and B within the values k1-k2 and c1-c2

You haven’t shown us the CVX output.

Not sure what you mean, the output in the second code block is what cvx put in the Matlab command window. Is there a command I should run to display more information?

And try installing the most recent Mosek 10., and select that version with cvx_solver.

Tried that and I get effectively the same thing

Calling Mosek_3 10.0.26: 30 variables, 12 equality constraints
   For improved efficiency, Mosek_3 is solving the dual problem.
------------------------------------------------------------

MOSEK Version 10.0.26 (Build date: 2022-10-24 12:03:48)
Copyright (c) MOSEK ApS, Denmark WWW: mosek.com
Platform: Windows/64-X86

Problem
  Name                   :                 
  Objective sense        : minimize        
  Type                   : CONIC (conic optimization problem)
  Constraints            : 12              
  Affine conic cons.     : 0               
  Disjunctive cons.      : 0               
  Cones                  : 0               
  Scalar variables       : 0               
  Matrix variables       : 5               
  Integer variables      : 0               

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 1                 time                   : 0.00            
Lin. dep.  - tries                  : 1                 time                   : 0.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.01    
Problem
  Name                   :                 
  Objective sense        : minimize        
  Type                   : CONIC (conic optimization problem)
  Constraints            : 12              
  Affine conic cons.     : 0               
  Disjunctive cons.      : 0               
  Cones                  : 0               
  Scalar variables       : 0               
  Matrix variables       : 5               
  Integer variables      : 0               

Optimizer  - threads                : 6               
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 12
Optimizer  - Cones                  : 0
Optimizer  - Scalar variables       : 0                 conic                  : 0               
Optimizer  - Semi-definite variables: 5                 scalarized             : 30              
Factor     - setup time             : 0.00              dense det. time        : 0.00            
Factor     - ML order time          : 0.00              GP order time          : 0.00            
Factor     - nonzeros before factor : 72                after factor           : 72              
Factor     - dense dim.             : 0                 flops                  : 2.64e+03        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   4.0e+00  2.7e+06  7.2e+00  0.00e+00   -8.160000000e+00  0.000000000e+00   1.0e+00  0.03  
1   1.2e+00  8.1e+05  1.4e+01  -2.46e+00  -8.690561756e+01  0.000000000e+00   3.0e-01  0.09  
2   1.1e-01  7.2e+04  2.2e+00  -2.18e+00  -1.469136633e+02  0.000000000e+00   2.7e-02  0.11  
3   9.8e-03  6.6e+03  1.5e+00  -1.17e+00  -6.891501352e+03  0.000000000e+00   2.4e-03  0.11  
4   1.5e-03  9.9e+02  5.4e-01  -1.00e+00  -4.252466407e+04  0.000000000e+00   3.7e-04  0.11  
5   3.8e-04  2.5e+02  2.8e-01  -1.10e+00  -1.738778952e+05  0.000000000e+00   9.4e-05  0.11  
6   2.8e-05  1.9e+01  4.5e-02  -7.03e-01  -8.004366305e+05  0.000000000e+00   7.1e-06  0.11  
7   4.6e-06  3.1e+00  1.7e-02  -9.27e-01  -4.551400384e+06  0.000000000e+00   1.1e-06  0.11  
8   9.5e-07  6.4e-01  8.9e-03  -9.73e-01  -2.749332029e+07  0.000000000e+00   2.4e-07  0.11  
9   1.5e-07  9.9e-02  3.8e-03  -9.71e-01  -2.134236352e+08  0.000000000e+00   3.7e-08  0.12  
10  6.2e-08  4.2e-02  2.7e-03  -1.15e+00  -5.874594096e+08  0.000000000e+00   1.6e-08  0.12  
11  1.1e-08  7.6e-03  8.2e-04  -7.55e-01  -1.645255744e+09  0.000000000e+00   2.8e-09  0.12  
Optimizer terminated. Time: 0.16    


Interior-point solution summary
  Problem status  : DUAL_INFEASIBLE
  Solution status : DUAL_INFEASIBLE_CER
  Primal.  obj: -2.4901412395e+00   nrm: 1e+00    Viol.  con: 3e-02    barvar: 0e+00  
Optimizer summary
  Optimizer                 -                        time: 0.16    
    Interior-point          - iterations : 11        time: 0.12    
      Basis identification  -                        time: 0.00    
        Primal              - iterations : 0         time: 0.00    
        Dual                - iterations : 0         time: 0.00    
        Clean primal        - iterations : 0         time: 0.00    
        Clean dual          - iterations : 0         time: 0.00    
    Simplex                 -                        time: 0.00    
      Primal simplex        - iterations : 0         time: 0.00    
      Dual simplex          - iterations : 0         time: 0.00    
    Mixed integer           - relaxations: 0         time: 0.00    

------------------------------------------------------------
Status: Infeasible
Optimal value (cvx_optval): +Inf

As for the scaling, can you divide all the input numbers by 10000, without getting into trouble with the 1 elements in` eye?

Tried that without much success

Sorry. I guess I had some scrolling issues. Your problem was determined to be (primal) infeasible, so read the corresponding ink from my previous post.

Id certain combinations of deleting two of the first 4 constraints are made, the problem becomes feasible. Presuming you want to keep tall those constraints:

You can add slack*eye(3) to the LHS of the last constraint. For values of slack up to about 1e7, that problem is reported as infeasible. Above that, it gets ill-posed as Mosek has numerical difficulties. Mosek 10.0.26 claims to have solved For slack = 1e7, but the solution is huge, and perhaps is just meaningless garbage. Perhaps Mosek 9.2.40 would report differently in this region. And it fails for slack = 1e-7. Mosek failed with primal ill-posed when II made slack a variable and minimize it, so that means the slack minimization problem is dual ill-posed. Mosek runs into never never land in the 1e7 region.

So basically, `Q doesn’t “want” to be PSD, even if given a huge eigenvalue increase. I’'m not a control theory person, so I’ll defer further diagnosis to someone else.

When slack = 5e-7

value(Q)

ans =
   1.0e+07 *
   8.255901620180961                   0                   0
                   0   0.000000000370846   0.000004219064422
                   0   0.000004219064422  -4.910555297931870

Edit: The problem is reported feasible for values of k2 up to just under 2e4, and infeasible for k2 >= 2e4.

Perhaps bad scaling is contributing to the difficulties. Is the problem inherently badly scaled, corresponding to vastly different time scales in the same problem? Maybe inherently ill-posed, at least for a double precision solver? Maybe this set of LMI’s isn’t so simple?

That is exactly what I was noticing as well. It does seem to be inherently badly scaled due to different time scales. The design goal is to create a controller that can handle a range of spring constants k1 to k2. The larger spring constant k2 corresponds to a much smaller time scale compared to the smaller constant k1. Also, a larger alpha corresponds to a more aggressive controller response (lower alpha should help with convergence but is not ideal for the actual problem). The problem might be solved if I could use a higher precision solver.

Using the ‘MSK_DPAR_INTPNT_CO_TOL_PFEAS’ option does seem to help at least for k2 = 2.7e5 (not 2.7e6) when using a slack of 1e-9. I can go up to alpha = 0.03 using this which is still lower than what I would like (ideally alpha > 0.1)

This is what I get for:
alpha = 0.03
k1 = 300;
k2 = 2.7e5;
c1 = 0.35;
c2 = 1000;
‘MSK_DPAR_INTPNT_CO_TOL_PFEAS’: 1e-12

Q = 
      544.344091573019      0	                       0
      0	                    8.93344183234658e-07	   0.00549554973744356
      0	                    0.00549554973744356	       33.8580384942786

Y = 
      -288.502368533700	     0	                        0
      0	                     -9.44005060501605e-06	    0.00481277501969532

K = 
      -0.530000000000001     0                          0
      0                      -7553.67764674437          1.22619105623424|

Maybe there’s a reason your design goal has not previously been achieved?

There are high precision LMI solvers, but not available under CVX. And I can’t vouch for any of them.

For instance:

http://www.math.umbc.edu/~potra/sdpham6/sdpham6Doc.pdf

https://www.tuhh.de/ti3/jansson/vsdp_cj.html (not really higher precision, but uses interval (and radial) arithmetic in INTLAB to do verified calculations - perhaps it will conclude the problem is in no man’s land where it can’t make a reliable feasibility determination one way or another. If you are going to get involved in these kinds of hairy situations, where you have to know for sure worst case error bounds, INTLAB (a MATLAB add-on) might be useful anyway.

1 Like