Cvxpy: Infeasible QP reported for an MPC problem with input constraints only

[gbanjac]

I tried to run a variant of the MPC example from this link. I modified the code so that the terminal constraints are not considered.

# Generate data for the control problem
import numpy as np
np.random.seed(2)
n = 8
m = 2
T = 50
alpha = 0.2
beta = 5
A = np.eye(n) + alpha*np.random.randn(n,n)
B = np.random.randn(n,m)
x_0 = beta*np.random.randn(n)

# Form and solve the control problem
from cvxpy import *
x = Variable((n, T+1))
u = Variable((m, T))
cost = 0
constraints = [x[:,0] == x_0]
for t in range(T):
    cost += sum_squares(x[:,t+1]) + sum_squares(u[:,t])
    constraints += [x[:,t+1] == A*x[:,t] + B*u[:,t], norm(u[:,t], 'inf') <= 1]
#constraints += [x[:,T] == 0]
prob = Problem(Minimize(cost), constraints)
prob.solve(verbose=True)

OSQP reports that the problem is primal infeasible. However, this does not seem to be possible since we only have constraints on the control inputs u. One feasible solution would be u[:,t]=0 for all t.

There seems to be a problem with the canonicalization.

Related issues: https://github.com/cvxgrp/cvxpy/issues/476 and https://github.com/cvxgrp/cvxpy/issues/540.

[bstellato]

ECOS has also numerical issues and MOSEK takes more than 1000 iterations to declare it infeasible.

[rileyjmurray]

To determine where the issue is, let's go ahead and turn it into a feasibility problem. Instead of minimizing a long sum of squares, just minimize "zero". When I did that with OSQP I got

status:               solved
solution polish:      unsuccessful
number of iterations: 471375
optimal objective:    0.0000
run time:             2.30e+01s
optimal rho estimate: 1.00e-06

This suggests digging deeper into sum_squares canonicalization (how many new variables are introduced?).

Edit: if I set cost = sum_squares(u), then the result is

status:               solved
solution polish:      unsuccessful
number of iterations: 175
optimal objective:    16.3602
run time:             9.04e-03s
optimal rho estimate: 7.69e-01

This result seems wrong, since u = 0 should be feasible (but the "optimal" objective is > 16).
On a related note, if cost = sum_squares(x), then OSQP again returns "infeasible."

[aszekMosek]

If I change to T=25 then the solution is

Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
  Primal.  obj: 1.9253654826e+06    nrm: 5e+05    Viol.  con: 1e-06    var: 0e+00    cones: 0e+00  
  Dual.    obj: 1.9253654901e+06    nrm: 1e+06    Viol.  con: 0e+00    var: 5e-08    cones: 2e-11

so you might say the solution norm is rather large (both primal and dual). For T=26 Mosek says unknown, and for some large T it starts claiming infeasibility.

If you instead make a feasibility problem then everything is perfectly fine.

Here is what Mosek gets for n=1, m=1, T=1. It looks reasonably good to me. (Except maybe for the unnecessary empty constraints).

[objective minimize]
   x0000 + x0001
[/objective]

[constraints]
  [con c0000]  <= 1e+00 [/con]
  [con c0001]  <= 1e+00 [/con]
  [con c0002]  x0004 - x0005 <= 0e+00 [/con]
  [con c0003]  - x0004 - x0005 <= 0e+00 [/con]
  [con c0004]  x0005 <= 1e+00 [/con]
  [con c0005]  x0002 = -1.068098047834227e+01 [/con]
  [con c0006]  - 9.166484305189059e-01 x0002 + x0003 + 5.626682722632947e-02 x0004 = 0e+00 [/con]
  [con c0007]  - x0000 + x0006 = 1e+00 [/con]
  [con c0008]  x0000 + x0007 = 1e+00 [/con]
  [con c0009]  - 2e+00 x0003 + x0008 = 0e+00 [/con]
  [con c0010]  - x0001 + x0009 = 1e+00 [/con]
  [con c0011]  x0001 + x0010 = 1e+00 [/con]
  [con c0012]  - 2e+00 x0004 + x0011 = 0e+00 [/con]
[/constraints]

[bounds]
  [b]               x0000,x0001,x0002,x0003,x0004,x0005 free [/b]
  [b]               x0006,x0007,x0008,x0009,x0010,x0011 free [/b]
  [cone quad k0000] x0006, x0007, x0008 [/cone]
  [cone quad k0001] x0009, x0010, x0011 [/cone]
[/bounds]

So possibly the problem is numerically nasty rather than anything wrong with cvxpy.

[joachimdahl]

ECOS has also numerical issues and MOSEK takes more than 1000 iterations to declare it infeasible.

I think you are misinterpreting the MOSEK log output. You may have looked at the infeasibility report (which for some reason is being printed by default) and misinterpreted the variable numbers as iterations counts - MOSEK has a default iteration limit of 400. I tried to run the example posted above, and MOSEK converged in less than 30 iterations with an infeasibility certificate. Not a very strong certificate, but as mentioned by Michal, the problem is numerically challenging, perhaps close to ill-posed.

[bstellato]

@joachimdahl I apologize for having misinterpreted the numbers but the output was so long it did not fit in my terminal emulator. The log file is here. The main point of this issue is not to undermine how fast MOSEK can detect infeasible problems but rather to understand why the problem is so badly scaled. In addition, the canonicalization for MOSEK and ECOS is in conic form and is different than the new QP one introduced in 1.0.

The model predictive control problem is strange because of various reasons:

• It should not be infeasible by construction (the terminal constraint is commented out).
• Gurobi solver (with the new QP canonicalization) finds the optimal solution in around 40 iterations. I have attached the output in case there are any confusions.
• OSQP solver (with the new QP canonicalization) report infeasibility.

I check the conditioning of the problem in the OSQP form and it seems quite strange for a toy MPC instance:

data = prob.get_problem_data(cvx.OSQP)[0]
A_osqp = spa.vstack([data['A'], data['F']]).tocsc()
print("Condition number of P: ", np.linalg.cond(data['P'].todense()))
print("Condition number of A: ", np.linalg.cond(A_osqp.todense()))

giving

Condition number of P:  inf
Condition number of A:  31027974.93709921

I wrote the same problem in YALMIP and it is easily solved without any infeasibility issues. Also the condition number of matrix A is much lower:

Condition number of P: Inf
Condition number of A: 4.0194843824

Note that the problems in YALMIP and CVXPY can be slightly different because of the random number generators. However, this should not cause any infeasibility issues.

Edit: minor typo in generating matrix A for OSQP in CVXPY. The condition number got even worse.

[SteveDiamond]

Gurobi succeeds due to its presolve, and the final objective value is enormous. I'm not sure what YALMIP is doing. I tried reformulating the problem so the number of variables was minimal, and it still failed. I don't know why the condition number is so high.

[gbanjac]

I think there are some numerical issues because the system is unstable (check the eigenvalues of A) and constraints on the inputs are too tight to stabilize the system. This is why the objective value is huge.

• If we set the constraints to norm(u[:,t], 'inf') <= 10, then OSQP converges in 200 iterations. Also, Gurobi returns a solution with the objective value by almost 4 orders of magnitude smaller than in the original case.

• If the prediction horizon is reduced to T=20, then OSQP finds a solution again.

• I tried to reproduce the same example with YALMIP and Gurobi still returns a solution with a huge objective value, while OSQP still reports infeasibility.

% Problem data
nx = 8;
nu = 2;
T = 50;
alpha = 0.2;
beta = 5;
A = [ 9.16648431e-01, -1.12533654e-02, -4.27239219e-01, ...
      3.28054162e-01, -3.58687117e-01, -1.68349473e-01, ...
      1.00576283e-01, -2.49057617e-01;
     -2.11590444e-01,  8.18198477e-01,  1.10290809e-01, ...
      4.58441603e-01,  8.30787860e-03, -2.23585089e-01, ...
      1.07811664e-01, -1.19231940e-01;
     -3.82609930e-03,  2.35000244e-01,  8.50425810e-01, ...
      1.80505019e-03, -1.75621579e-01, -3.12868341e-02, ...
      5.13140904e-02, -1.97755810e-01;
     -6.77643932e-02, -4.72368062e-02, -1.27531002e-01, ...
      7.62477543e-01, -2.84243445e-01, -3.06990391e-02, ...
     -5.38113920e-02,  4.46273358e-01;
     -4.86953515e-01,  2.25453010e-02,  7.40889073e-02, ...
      2.71926773e-01,  1.10037144e+00, -1.68842741e-01, ...
      1.95229432e-06,  1.08470514e-01;
     -6.27016394e-02,  1.54202348e-01, -3.73618131e-01, ...
      3.46236933e-01,  2.93535602e-01,  9.32864532e-01, ...
      1.22268156e-01,  9.59411837e-03;
     -1.65827058e-01,  1.75420437e-02,  2.00073177e-01, ...
     -7.62185035e-02, -7.51338846e-02, -1.48941526e-02, ...
      1.08669927e+00,  2.55675846e-01;
     -1.26935861e-01,  1.01679249e-01,  4.32232013e-02, ...
     -3.71722477e-01, -8.38632964e-02, -2.64657797e-02, ...
     -7.91404794e-03,  1.06520069e+00];
B = [-2.04032305,  0.04625552;
     -0.67767558, -1.43943903;
      0.52429643,  0.73527958;
     -0.65325027,  0.84245628;
     -0.38151648,  0.06648901;
     -1.09873895,  1.58448706;
     -2.65944946, -0.09145262;
      0.69511961, -2.03346655];
x_0 = [-0.94734632, -0.38609333,  4.12351503,  6.2410646 , ...
       -2.01946135, -6.92259333,  6.83617712,  6.08942817]';


% MPC
x = sdpvar(repmat(nx,1,T+1),ones(1,T+1));
u = sdpvar(repmat(nu,1,T),ones(1,T));

objective = 0;
constraints = [x{1} == x_0];
for k = 1:T
    objective = objective + x{k}'*x{k} + u{k}'*u{k};
    constraints = [constraints, x{k+1} == A*x{k} + B*u{k}];
    constraints = [constraints, -1 <= u{k} <= 1];
end
optimize(constraints,objective,sdpsettings('solver','gurobi'));
value(objective)

I think we can close the issue, but one should definitely make a more reasonable MPC example here.

[bstellato]

I agree in closing the issue. I suspect the YALMIP generated problem had stable dynamics thanks to the seed used. The control example needs to be fixed but there is already this issue https://github.com/cvxgrp/cvxpy/issues/540.