Unbounded variables masked by constraints pushing in both directions

[pgkirsch]

I was adamant that the GP below is feasible so I solved it by hand. It is indeed feasible, but when I run the code, it gives me the following message:

RuntimeWarning: final status of solver 'mosek' was 'UNKNOWN', not 'optimal'.

To generate a feasibility-finding relaxed version of your
Geometric Program, run model.program.feasibility_search().

The infeasible solve's result is stored in the 'result' attribute
(model.program.result) and its raw output in 'solver_out'.

This suggests that my problem is infeasible, which is not true because I can tell you a value for each variable that will satisfy all of the constraints.

Original code:

from gpkit.shortcuts import *
import numpy as np

A_bypass = Var('A_bypass') #  bypass cross sectional area 
A_pod = Var('A_pod') #pod cross sectional area 
D = Var('D') # total drag 
d = Var('d') # deceleration 
F = Var('F') # net force on pod (positive defined backwards) 
F_ski = Var('F_ski') #  force generated by air bearings 
M = Var('M') #  cruise Mach number 
M_bypass = Var('M_bypass') #  external bypass Mach number 
dotm_bypass = Var('dotm_bypass') # bypass mass flow 
dotm_inlet = Var('dotm_inlet') # compressor mass flow 
dotm_ski = Var('dotm_ski') # air bearing mass flow 
dotm_infty = Var('dotm_infty') # freestream mass flow 
T = Var('T') # total thrust 
V = Var('V') # cruise velocity 

A_tube = Var('A_tube', value=2) # cross sectional tube area 
m = Var('m', 1000) # total pod mass 
C_D = Var('C_D', value=0.04) #  Drag coefficient 
g = 9.81 # gravitational acceleration
rho_infty = Var('rho_infty', value=0.00116) #  freestream air density 
K = Var('K', value=50000)

m = Model(d,
         [d == F/m, # Newton's 2nd law
          F >= D + T, # Force balance
          D == 0.5*rho_infty*V**2*A_pod*C_D, # Simplistic drag model 
          T == dotm_infty*V, # anti-thrust
          dotm_infty >= dotm_inlet + dotm_bypass, # Conservation of mass
          dotm_infty == rho_infty*V*A_tube, # Freestream mass flow
          dotm_inlet == dotm_ski, # Conservation of mass
          F_ski <= dotm_ski*K, # Force generated by skis
          F_ski == m*g, # Force equilibrium
          ],
         {V: 300, dotm_bypass: 0.4})
sol = m.solve()

It turns out the issue is that A_pod is not bounded. When I add a lower bound on A_pod, e.g. A_pod >= 1, I get a solution.

The tricky thing is that it looks like A_pod is bounded because it has an equality constraint, however, the solver is quite happy to let D go to zero along with A_pod.

tl;dr: I think the bound checking heuristic might need to be enhanced.

[pgkirsch]

Also the To generate a feasibility-finding... message should only be displayed when we can confidently say that a GP is infeasible. Otherwise it misleads a user into thinking their GP is infeasible.

[whoburg]

Agreed -- I think the short-term fix (for 0.3.2) should be to only print the feasibility-finding suggestion for primal-infeasible GPs (and not for UNKNOWN).

Once that's fixed, we should leave this ticket open and move to FUTURE to work on better diagnostics/error messages for tricky GPs like this where a variable is unbounded but masked by constraints in both directions.

edit: the issue of better ways to handle UNKNOWN is captured in #519.

[bqpd]

We don't have standard status messages besides "optimal" right now, so I just reworded the error message in the PR. Should we standardize status messages?

[whoburg]

Ok, looks good for now. Standardizing / parsing status messages from various solvers sounds too complicated for 0.3.2.