MOSEK DUAL INFEASIBILITY

[p4wp4w]

Hello everyone,

I am using Yalmip (solver: MOSEK) to solve a conic optimization problem. However, the solver throws the error of DUAL INFEASIBILITY. Interestingly, under the same conditions (Hurwitz minimal controllable & observable) algorithm works fine for another system with the same dimensions. Also, there exists a proof of solution existence!

Below, I have provided the code:

clear 

close all

clc

%%

load plant Ap Bp Hp Cp Dp

np = sqrt(numel(Ap));

nu = numel(Bp)/np;

nw = numel(Hp)/np;

ny = numel(Cp)/np;

%% Controller

load controller

nk = sqrt(numel(Ak));

%%

Delta = (eye(ny)-Dp*Dk)^-1;

Deltat = (eye(nu)-Dk*Dp)^-1;

Ab = [Ap+Bp*Deltat*Dk*Cp Bp*Deltat*Ck;Bk*Delta*Cp Ak+Bk*Delta*Dp*Ck];

B0 = [Bp*Deltat;Bk*Delta*Dp];

Bb = [Bp*Deltat -Bp*Deltat*Dk;Bk*Delta*Dp -Bk*Delta];

Cb1 = [Deltat*Dk*Cp Deltat*Ck];

D01 = Deltat*Dk*Dp;

Db1 = [Deltat -Deltat*Dk];

Cb2 = [Delta*Cp Delta*Dp*Ck];

D02 = Delta*Dp;

Db2 = [Delta*Dp -Delta*Dp*Dk];

%% LMI

% Define variables

Q = sdpvar(np+nk,np+nk);

L = sdpvar(nu+ny,nu);

U_dum = sdpvar(nu,1);

U = diag(U_dum);

gama = sdpvar(1);

% Define constraints and objective

Pi1 = [Q*Ab'+Ab*Q B0*U+Bb*L-Q*Cb1' zeros(np+nk,nu) Q*Cb2';(B0*U+Bb*L-Q*Cb1')' -2*U-U*D01'-D01*U-L'*Db1'-Db1*L,eye(nu,nu) U*D02'+L'*Db2';zeros(nu,np+nk) eye(nu,nu) -gama*eye(nu,nu) zeros(nu,ny);(Q*Cb2')' (U*D02'+L'*Db2')' zeros(ny,nu) -gama*eye(ny,ny)];

Constraints = [Pi1 <= -1e-3*eye(size(Pi1)),gama >= 0,Q >= 0,U_dum(:) >= 1e4];

Objective = [gama];

% Set some options for YALMIP and solver

options = sdpsettings('verbose',1,'solver','mosek');

sol = optimize(Constraints,Objective,options);

if sol.problem == 0

    % Extract and display value

    Q = value(Q);

    U_dum = value(U_dum);

    U = diag(U_dum);

    gama = value(gama);

    L = value(L);

    Pi1 = [Q*Ab'+Ab*Q B0*U+Bb*L-Q*Cb1' zeros(np+nk,nu) Q*Cb2';(B0*U+Bb*L-Q*Cb1')' -2*U-U*D01'-D01*U-L'*Db1'-Db1*L,eye(nu,nu) U*D02'+L'*Db2';zeros(nu,np+nk) eye(nu,nu) -gama*eye(nu,nu) zeros(nu,ny);(Q*Cb2')' (U*D02'+L'*Db2')' zeros(ny,nu) -gama*eye(ny,ny)];

    eig(Pi1)

    

    theta = L*inv(U);

else

    

     display('Hmm, something went wrong!');

     sol.info

     yalmiperror(sol.problem)

end

I use MATLAB 2016B and I get the following as output:

MOSEK Version 8.0.0.46 (Build date: 2016-12-2 14:03:28)

Copyright (c) MOSEK ApS, Denmark. WWW: mosek.com

Platform: Windows/64-X86

Problem

  Name                   :                 

  Objective sense        : min             

  Type                   : CONIC (conic optimization problem)

  Constraints            : 89              

  Cones                  : 0               

  Scalar variables       : 3               

  Matrix variables       : 2               

  Integer variables      : 0               

Optimizer started.

Conic interior-point optimizer started.

Presolve started.

Linear dependency checker started.

Linear dependency checker terminated.

Eliminator - tries                  : 0                 time                   : 0.00            

Lin. dep.  - tries                  : 1                 time                   : 0.01            

Lin. dep.  - number                 : 0               

Presolve terminated. Time: 0.01    

Optimizer  - threads                : 2               

Optimizer  - solved problem         : the primal      

Optimizer  - Constraints            : 89

Optimizer  - Cones                  : 0

Optimizer  - Scalar variables       : 3                 conic                  : 0               

Optimizer  - Semi-definite variables: 2                 scalarized             : 249             

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 : 4005              after factor           : 4005            

Factor     - dense dim.             : 0                 flops                  : 8.12e+005       

ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  

0   4.0e+000 2.0e+004 1.0e+000 0.00e+000  -1.152000000e+000 0.000000000e+000  1.0e+000 0.01  

1   5.2e-001 2.6e+003 3.8e-001 1.03e+000  -4.201467791e+000 -5.351723879e+001 1.3e-001 0.05  

2   1.1e-001 5.6e+002 2.5e-001 1.23e+000  -7.028350534e+000 -4.618250341e+001 2.8e-002 0.05  

3   3.3e-002 1.6e+002 1.8e-001 1.91e+000  -1.424900864e+001 -2.237265998e+001 8.1e-003 0.05  

4   8.5e-003 4.2e+001 4.2e-002 4.78e-001  -4.010920949e+001 -3.891769814e+001 2.1e-003 0.05  

5   2.3e-003 1.2e+001 1.1e-002 8.87e-002  -6.837484226e+001 -6.279696916e+001 5.8e-004 0.06  

6   6.1e-004 3.1e+000 2.6e-003 -1.92e-001 -1.193555651e+002 -1.082321495e+002 1.5e-004 0.06  

7   1.2e-004 5.8e-001 4.3e-004 -1.88e-001 -2.143005469e+002 -1.972093275e+002 2.9e-005 0.06  

8   2.0e-005 1.0e-001 5.7e-005 -2.64e-001 -3.951409717e+002 -3.640944928e+002 5.0e-006 0.06  

9   3.4e-006 1.7e-002 6.1e-006 -4.00e-001 -7.988941643e+002 -7.209643049e+002 8.4e-007 0.06  

10  4.4e-007 2.2e-003 4.5e-007 -4.96e-001 -2.029401325e+003 -1.776605875e+003 1.1e-007 0.06  

11  9.1e-008 4.5e-004 5.7e-008 -5.00e-001 -4.748096623e+003 -4.093553719e+003 2.2e-008 0.08  

12  5.8e-008 2.9e-004 2.9e-008 5.47e-001  -7.221039724e+003 -6.214942003e+003 1.7e-008 0.08  

13  1.7e-008 8.3e-005 2.8e-009 -9.07e-001 -5.437096766e+004 -4.511366774e+004 6.2e-009 0.08  

14  4.3e-010 2.1e-006 1.1e-011 -1.02e+000 -2.153917593e+006 -1.784016654e+006 1.6e-010 0.08  

15  9.8e-012 4.9e-008 4.0e-014 -9.92e-001 -8.811637075e+007 -7.306154430e+007 3.8e-012 0.08  

16  1.9e-013 9.7e-010 1.2e-016 -9.89e-001 -4.103041443e+009 -3.408505242e+009 7.5e-014 0.08  

17  1.4e-013 1.6e-011 2.6e-019 -9.86e-001 -2.263629472e+011 -1.884431244e+011 1.2e-015 0.09  

18  1.2e-013 1.4e-011 2.2e-019 -9.91e-001 -2.551020110e+011 -2.123968985e+011 1.1e-015 0.09  

19  1.2e-013 1.4e-011 2.2e-019 -1.06e+000 -2.551020110e+011 -2.123968985e+011 1.1e-015 0.09  

Interior-point optimizer terminated. Time: 0.11. 

MOSEK DUAL INFEASIBILITY REPORT.

Problem status: The problem is dual infeasible

Optimizer terminated. Time: 0.13    

Interior-point solution summary

  Problem status  : DUAL_INFEASIBLE

  Solution status : NEAR_DUAL_INFEASIBLE_CER

  Primal.  obj: -2.1758478422e+000  nrm: 1e+003   Viol.  con: 9e-009   var: 7e-012   barvar: 0e+000 

Optimizer summary

  Optimizer                 -                        time: 0.13    

    Interior-point          - iterations : 20        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    

Mosek error: MSK_RES_TRM_STALL ()

Hmm, something went wrong!

Can anyone help me in this regards?

Best,

Atta

[Johan Löfberg]

The fact that already this fails

optimize(Pi1 <= 0)

should tell you that there is something bad with your model.

Heck, already this fails

optimize(Q*Ab'+Ab*Q <= 0)

Although Ab is stable, it appears to be a challenging matrix for some reason. Combining that with your requirement on huge elements U_dum(:) >= 1e4 leads to even worse numerics

The infeasibility though is most likely due to your randomly picked value 1e-3 when trying to force Pi1 to be negative definite. If you remove that it does not lead to infeasibility, but just numerical issues of varying degree

My guess is that beyond the ad numerics of the dynamics gama has to be insanely large for the problem to be feasible, so the solver fails.

[Erling D. Andersen]

Just to make sure. The problem  is dual infeasible and hence unbounded if the problem has a feasible solution.