Debugging UNKNOWN

[pgkirsch]

It would be a huge help if we could figure out a way of debugging programs that return UNKNOWN. Specifically, it would be useful to know what needs to be constrained for the problem to be solved.

A MWE for a problem that returns unknown:

In [1]: from gpkit import Variable, Model

In [2]: x = Variable('x')

In [3]: y = Variable('y')

In [4]: m = Model(x, [x*y >= 1, x+y <= 2])

In [5]: m.solve()
Using solver 'mosek'
Solving for 2 variables.
Solving took 0.0249 seconds.
---------------------------------------------------------------------------
RuntimeWarning                            Traceback (most recent call last)
<ipython-input-5-1223da3220cb> in <module>()
----> 1 m.solve()

/Users/philippekirschen/Documents/MIT/Research/GPkit/gpkit/gpkit/model.pyc in solve(self, solver, verbosity, skipsweepfailures, *args, **kwargs)
    340         try:
    341             return self._solve("gp", solver, verbosity, skipsweepfailures,
--> 342                                *args, **kwargs)
    343         except ValueError as err:
    344             if err.message == ("GeometricPrograms cannot contain Signomials"):

/Users/philippekirschen/Documents/MIT/Research/GPkit/gpkit/gpkit/model.pyc in _solve(self, programType, solver, verbosity, skipsweepfailures, *args, **kwargs)
    498             self.program, solvefn = form_program(programType, signomials,
    499                                                  verbosity=verbosity)
--> 500             result = solvefn(*args, **kwargs)
    501             solution.append(parse_result(result, constants, beforesubs))
    502         solution.program = self.program

/Users/philippekirschen/Documents/MIT/Research/GPkit/gpkit/gpkit/geometric_program.pyc in solve(self, solver, verbosity, *args, **kwargs)
    191                                  "not 'optimal'." %
    192                                  (solver, solver_out.get("status", None)) +
--> 193                                  "\n\nThe infeasible solve's result is stored"
    194                                  " in the 'result' attribute"
    195                                  " (model.program.result)"

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

The infeasible solve's result is stored in the 'result' attribute (model.program.result) and its raw output in 'solver_out'. If the problem was Primal Infeasible, you can generate a feasibility-finding relaxation of your Model with model.feasibility().

[pgkirsch]

My first proposed solution is very much brute force. It works for the MWE above though, and perhaps could be a starting point for something more sophisticated. It is based on the assumption that a program returns unknown because the problem is "underconstrained" from the perspective of the solver, and that all it needs is for a certain variable to be fixed. This code could help to identify which variable(s) are at the root of the problem.

def test_unknown(m):
    for key in m.varkeys:
        if key not in m.constants and key not in m.substitutions:
            m.substitutions.update({key: 1})
            print "Solving with fixed", key
            try:
                m.solve(verbosity=0)
                print "Problem solves"

            except ValueError:
                try:
                    m.localsolve(verbosity=0)

                except RuntimeWarning:
                    print "Runtime Warning"

            del m.substitutions[key]

For the MWE above:

In [6]: from test_unknown import test_unknown

In [7]: test_unknown(m)
Solving with fixed y
Problem solves
Solving with fixed x
Problem solves

[bqpd]

@whoburg, is this related to / generalizable with your result on one class of UNKNOWN / singular solves?

[whoburg]

Hmm. @pgkirsch, I'm really glad you started this issue, but I'm a bit skeptical on the MWE and proposed solution, for two reasons:

1. Isn't the feasible set of the MWE the single point (1, 1)? We've seen issues with problems that have a single feasible point before. Thus one hypothesis is that solvers sometimes return UNKNOWN on problems with a single feasible point. It would be good to know whether there are other potential causes of UNKNOWN. (sidenote: it looks so scary in all caps).
2. The proposed fix substitutes 1, which just so happens to be the correct single feasible value for this problem. Does the code run with m.substitutions.update({key: 1}) changed to m.substitutions.update({key: 2})? Assuming not, it seems to me that it's going to be difficult to guess what to substitute.

[pgkirsch]

Yes, you're right, we have seen this before (specifically in issue #403). I reduced the example in issue #380 to another MWE:

from gpkit import Variable, Model

D  = Variable('D')
F  = Variable('F')
mi = Variable('m_i')

mb = Variable('m_b', 0.4)
V  = Variable('V', 300)

m = Model(F,
         [F >= D + V,
          V >= mi + mb,
          ])
sol = m.solve()

For this problem, the issue isn't that the feasible set is a single point. F can be anything greater than 300 + D, where D is completely unconstrained. mi is upper bounded, but there is no pressure on it, hence why I think it contributes to the solver returning UNKNOWN 👻 👹 😨 .

Glad it looks scary, it consistently ruins my day.

Edit by @whoburg: note that cvxopt solves the above problem -- only MOSEK returns UNKNOWN.

[bqpd]

We could lowercase the error messages? :p :p :p

The solver for example two does say "D was not lower-bounded" etc. More seriously we could remind the user what was unbounded and suggest that fixing those might fix the "UNKNOWN"

[bqpd]

I actually am surprised we don't remind the user of the unbounded vars on solve failures already. Should I just push that up?

[pgkirsch]

@whoburg you are also correct about point 2, that the proposed heuristic doesn't work for the case with one feasible point if the value is not exactly right.

[pgkirsch]

@bqpd For this example you do provide a warning message saying those two variables are unbounded, and again this might not be the greatest MWE, because the fix is quite obvious (and made clear by the message).

[bqpd]

Well, it's not entirely obvious what we would guess as bounds. And sometimes unbounded things are a model goof. But yeah, that's what we'd do for this example.

[pgkirsch]

A better MWE:

from gpkit import Variable, Model

Ap = Variable('A_p')
D  = Variable('D')
F  = Variable('F')
mi = Variable('m_i')
mf = Variable('m_f')
T  = Variable('T')

Fs = Variable('Fs', 0.9)
mb = Variable('m_b', 0.4)
rf = Variable('r_f', 0.01)
V  = Variable('V', 300)

m = Model(F,
         [F >= D + T,
          D == rf*V**2*Ap,
          T ==  mf*V,
          mf >= mi + mb,
          mf == rf*V,
          Fs <= mi,
          ])
sol = m.solve()

The problem here is that GPkit can't tell that Ap is not lower bounded because it is in an equality constraint. If you add a constraint, Ap >= 1, the problem solves.

[pgkirsch]

Simplified further:

from gpkit import Variable, Model

Ap = Variable('A_p')
D  = Variable('D')
F  = Variable('F')
mi = Variable('m_i')

Fs = Variable('Fs', 0.9)
mb = Variable('m_b', 0.4)
V  = Variable('V', 300)

m = Model(F,
         [F >= D + V**2,
          D == V**2*Ap,
          V >= mi + mb,
          Fs <= mi,
          ])
sol = m.solve()

edit by @whoburg: cvxopt solves this too -- mosek returns unknown.

[bqpd]

Yeah. I've been thinking for a while that we could figure this kind of thing out automatically. Any ideas for how to do that?

[pgkirsch]

Not exactly that, but my crappy proposed heuristic above at least works for this new MWE, regardless of choice of substitution value.

[pgkirsch]

Otherwise, maybe something like for every equality constraint, you could check with only one of the constituent inequality constraints at a time and see if you learn something?

[bqpd]

And check 2^(number of equality constraints) times? 👻 😨

[bqpd]

Everything else on the right side has a lower bound, everything on the left side has an upper bound...

[bqpd]

...so D/V**2 is upper-bounded, so D/V**2 >= Ap is a legit constraint

[bqpd]

let's divide each EQConstraint out to m == 1. I think that if one of the variables in m is unbounded (and all other variables are bounded), m_without_x == x is valid, because m_without_x >= x is clearly a good upper bound, and m_without_x <= x a good lower bound.

[bqpd]

If m_without_x is upper bounded, m_without_x >= x is legit, and similarly for lower-bounded and <=

[bqpd]

Now, if m_without_x is upper bounded and m_without_y is lower bounded, m_without_x >= x is legit, and it seems like m_without_y <= y is legit, so I'd expect m/x >= 1, m/y <= 1, m == 1 to be legit for the solvers. Note that it algebraically reduces to x <= 1, y >= 1.

[bqpd]

Given a monomial m, let's call a variable which is either in the numerator and upper-unbounded or in the denominator and lower-unbounded a "potential increaser" because if that variable's value has the potential to increase the value of m unboundedly. The opposite of that we'll call a "potential decreaser". I'm reasonably sure that if a monomial m has at most one potential increaser and one potential decreaser (which may be the same variable), then by the above m == 1 will be a dual-feasible constraint because m_without_increasers >= increaser and m_without_decreasers <= decreaser are dual-feasible constraints. (where by dual-feasible I really just mean 'solveable for some values of the constants in the problem')

[whoburg]

@bqpd, I've been trying to follow your comments above but getting caught up on a few details -- let's discuss in person.

Alternatively, and even better, I suggest writing up a short pdf to explain your point.

[whoburg]

This is an important ticket -- the UNKNOWN issue is a big problem.

I'd like us to start a working latex document to categorize the types of failures that result in UNKNOWN.

Here's a start.

Unbounded Variable(s)

Typical MWE: a lower-unbounded variable inside a posynomial constraint
Typical cvxopt behavior: returns unknown
Typical mosek behavior: computed cost mismatch (issue #361)

MWE:

from gpkit import Variable, Model

x = Variable('x')
t = Variable('t')

m = Model(t, [t >= 1 + x])
sol = m.solve()

On this MWE, cvxopt returns unknown and mosek has an issue #361 error.

Note that there are a number of strange tweaks that can cause cvxopt to actually return solutions (with very small values for the unbounded variable). Examples are listed in a comment on issue #460. Mosek, on the other hand, continues raising an issue #361 error.

Note this appears to be the failure mode associated with @pgkirsch's MWE above adapted from #380.

Nearly dual-infeasible models

Typical MWE: specific, precise exponents are required on a constraint to achieve dual feasibility
Typical cvxopt behavior: returns unknown
Typical mosek behavior: solves for correctly-chosen exponents; returns dual-infeasibile for slightly-tweaked exponents.

MWE:

from gpkit import Variable, Model
import numpy as np

CL = Variable("CL")
CD = Variable("CD")
V = Variable("V", "m/s", "cruise speed")
W = Variable("W", 200, "lbf", "weight")
rho = Variable(r"\rho", 1.2, "kg/m^3", "air density")
S = Variable("S", 190, "ft^2", "wing area")
m = Model(V*(W/(CL/CD)),
          [0.5*rho*CL*S*V**2 == W,
           CL**1.5/CD <= 20])
sol = m.solve("cvxopt")

With a little math (to eliminate V), it becomes clear that the objective value here is proportional to CD/CL**1.5. Thus, the bound required on CL and CD is a bound on CL**1.5/CD. Change the 1.5 to 1.51, or any other value, and the problem becomes dual infeasible according to mosek.

Note that when the == is changed to >=, cvxopt returns a rank error instead of unknown.

In my opinion there is a big opportunity with this class of problems. In particular, I think we might be able to use some linear algebra tricks to identify problems that are nearly dual infeasible and provide other useful information to modelers. The is still an idea being formed; I'll elaborate in a separate ticket.

[bqpd]

Instead of being clever with math we've recently been just using resolves.

[pgkirsch]

@bqpd what do you mean by "resolves"?

[bqpd]

re-solving the model with boundedconstraintsets / realaxed models

[pgkirsch]

Do boundedconstraintsets and relaxed models help with the nearly-dual-infeasible case? p.s. what do you mean by relaxed models?

[bqpd]

Relaxed models are my new name for primal-feasibility models. It may or may not stick.

[pgkirsch]

I'm wondering if there is a lazy approach to handling the nearly dual infeasible case, for example relaxing the convergence criterion of the solver.

[whoburg]

😨

[pgkirsch]

Ha that went over about as well as I was expecting it to......

[bqpd]

This issue is partly covered by the new .debug features in feasmodtest. I think the comments on this issue might usefully enhance the new debugging documentation, but otherwise this issue can be closed when that gets merged?

[whoburg]

before we close this based on the existence of .debug, let's at least go through all the MWEs in this issue and see what .debug says about them.

[bqpd]

1st example: does not solve with cvxopt debug
2nd example: does not solve with cvxopt debug, prints bounds warnings.
3rd example: does not solve with cvxopt debug
4th example: solves with cvxopt debug just as it solves with cvxopt

:-/

Should be tried again with MOSEK

[mayork]

@1ozturkbe this is useful

[whoburg]

related: JuliaOpt/MathProgBase.jl#164

[pgkirsch]

I have a model that solves perfectly well for a given set of inputs, however, when slightly higher performance is asked of the model (higher speed, heavier payload) it returns unknown, and the solver iterations look completely different (reaches max iterations after plateauing pretty early and making no convergence progress. When I say a slightly higher performance, we're talking about less than 1% increases. Plotting a sweep of this variable, this point typically occurs as the objective starts growing exponentially, so it's understandable why the solver may have issues with it. It's possible that the model is genuinely becoming infeasible, but I am skeptical that this is the case.

When I "diff" the solution of the original model with the solver output at the end of max iterations with the higher performance model there are substantial differences, most noticeably in variables I wouldn't expect to care. From what I recall, these solvers typically move through infeasible regions as they optimize the dual function (is that right?) but is there anything meaningful to be gained from seeing how the solution evolves iteration-to-iteration? And does GPkit have an interface that allows this kind of solution stepping?

Coming at it from a different angle, are there exposed solver parameters that control how a solver responds to perceived near-dual-infeasibility?

Finally, does anyone from the modeling side of things (@mayork, @1ozturkbe, @mjburton11) have any good debugging tips they have discovered for this specific type of unknown problem? .debug() didn't seem to tell me much (see #1143), other than the fact that no variables are unbounded.

[pgkirsch]

Eesh. cvxopt doesn't seem to be giving very credible results around this feasibility threshold. As I increment one of the "performance" variables (e.g. speed) in very small steps (0.04% increase) the result changes dramatically (50% increase in objective value). Based on conversation with @bqpd on #1143 it seems that the problem is becoming infeasible. This raises two questions for me: (1) is the dramatically different threshold solution truly globally optimal? (2) why isn't the solver returning primal infeasible instead of unknown once that threshold is crossed?

[pgkirsch]

This is most likely specific to my particular model, but I have actually had some success by increasing the resolution of discretization used. i.e. More elements in key vector variables take this model from returning unknown to returning a feasible solution.

[bqpd]

.debug now debugs all of the above, if you also check the output and note that certain variables are getting quite small...

@pgkirsch do you mind if I close this thread? I like that it's been a place to bring new weird UNKNOWNS, but it's gotten a bit unwieldy...

[pgkirsch]

Sure go for it. On Nov 9, 2017, at 22:07, bqpd <notifications@github.com<mailto:notifications@github.com>> wrote: .debug now debugs all of the above, if you also check the output and note that certain variables are getting quite small... @pgkirsch<https://github.com/pgkirsch> do you mind if I close this thread? I like that it's been a place to bring new weird UNKNOWNS, but it's gotten a bit unwieldy... — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<#519 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AITKc6IolrWShSjIc7cPdNSei9NSA9e-ks5s0-f6gaJpZM4HaHod>.

[bqpd]

Excitingly, we can now catch all of @pgkirsch's examples before even solving, by implementing an algorithm similar to the one I describe above...

[bqpd]

Catching @whoburg's marginally feasible examples will probably require the SVD analysis he discussed!

[whoburg]

Cool! If anyone wants to work on that let me know, I'd be excited to discuss the math.