Primal solution computed cost did not match solver-returned cost

[pgkirsch]

I have gotten the above error a number of times for several different GPs and SPs.

Example (from the beam example discussed in #269):

from gpkit.shortcuts import *
from numpy import ones

N = 6
L = 5.
dx = L/(N-1)

EI = Var("EI",10)

P = 100
p = Vec(N, "p", descr="Distributed load")
p = p.sub(p, P*ones(N))

V  = Vec(N, "V", label="Internal shear")
M  = Vec(N, "M", label="Internal moment")
th = Vec(N, "th", label="Slope")
w  = Vec(N, "w", label="Displacement")

eps = 1E-6 #something small

substitutions = {var: eps for var in [V[-1], M[-1], th[0], w[0]]}

objective = w[-1]

constraints = [EI*V.left[1:N]     >= EI*V[1:N]    + 0.5*dx*p.left[1:N]     + 0.5*dx*p[1:N],
                       EI*M.left[1:N]     >= EI*M[1:N]    + 0.5*dx*V.left[1:N]     + 0.5*dx*V[1:N],
                                      EI*th.right[0:N-1] >= EI*th[0:N-1] + 0.5*dx*M.right[0:N-1]  + 0.5*dx*M[0:N-1],
                                                     EI*w.right[0:N-1]  >= EI*w[0:N-1]  + 0.5*dx*th.right[0:N-1] + 0.5*dx*th[0:N-1]]

m = Model(objective, constraints, substitutions)

sol = m.solve()

Output: (long traceback followed by..)

RuntimeWarning: Primal solution computed cost did not match solver-returned cost: 8.0000014 vs 0.78125166625

[bqpd]

What's m.program.solver_log and m.program.result? (so glad we have those now)

[pgkirsch]

['Number of Hessian non-zeros: 97',
 '* Solving exponential optimization problem on dual form. *',
 '* The following log information refers to the solution of the dual problem. *',
 'Computer',
 '  Platform               : MACOSX/64-X86   ',
 '  Cores                  : 4               ',
 '',
 'Problem',
 '  Name                   :                 ',
 '  Objective sense        : max             ',
 '  Type                   : GECO (general convex optimization problem)',
 '  Constraints            : 21              ',
 '  Cones                  : 0               ',
 '  Scalar variables       : 53              ',
 '  Matrix variables       : 0               ',
 '  Integer variables      : 0               ',
 '',
 'Optimizer started.',
 'Interior-point optimizer started.',
 'Presolve started.',
 'Linear dependency checker started.',
 'Linear dependency checker terminated.',
 'Eliminator - tries                  : 0                 time                   : 0.00            ',
 "Eliminator - elim's                 : 1               ",
 'Lin. dep.  - tries                  : 1                 time                   : 0.00            ',
 'Lin. dep.  - number                 : 0               ',
 'Presolve terminated. Time: 0.01    ',
 'Matrix reordering started.',
 'Local matrix reordering started.',
 'Local matrix reordering terminated.',
 'Matrix reordering terminated.',
 'Optimizer  - threads                : 4               ',
 'Optimizer  - solved problem         : the primal      ',
 'Optimizer  - Constraints            : 19',
 'Optimizer  - Cones                  : 0',
 'Optimizer  - Scalar variables       : 51                conic                  : 0               ',
 'Optimizer  - Semi-definite variables: 0                 scalarized             : 0               ',
 'Factor     - setup time             : 0.02              dense det. time        : 0.00            ',
 'Factor     - ML order time          : 0.00              GP order time          : 0.00            ',
 'Factor     - nonzeros before factor : 206               after factor           : 325             ',
 'Factor     - dense dim.             : 0                 flops                  : 0.00e+00        ',
 'ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  ',
 '0   5.5e+00  2.8e+01  7.3e+01  0.00e+00   -7.176322419e+01  0.000000000e+00   1.0e+00  0.03  ',
 '1   1.1e+00  5.7e+00  1.6e+01  -1.25e+00  -9.374628227e+01  3.723222950e+01   2.1e-01  0.04  ',
 '2   5.5e-01  2.6e+00  8.3e+00  -7.83e-01  -2.750022939e+01  3.290595624e+01   1.1e-01  0.04  ',
 '3   3.0e-01  1.2e+00  4.6e+00  3.42e+00   -6.304231453e+00  5.954818600e+00   5.9e-02  0.04  ',
 '4   3.5e-02  2.1e+00  7.0e-01  2.08e+00   -5.497130106e-01  5.101754563e-01   8.3e-03  0.04  ',
 '5   3.3e-03  1.9e+00  9.1e-02  1.21e+00   -2.692431110e-01  -1.362583019e-01  1.3e-03  0.04  ',
 '6   1.9e-04  1.7e+00  7.0e-03  1.07e+00   -2.483346936e-01  -2.383062906e-01  1.1e-04  0.04  ',
 '7   9.8e-06  1.5e+00  5.0e-04  1.05e+00   -2.469639799e-01  -2.462656833e-01  8.3e-06  0.04  ',
 '8   1.1e-07  1.5e+00  1.8e-05  1.06e+00   -2.468632866e-01  -2.468380828e-01  3.2e-07  0.04  ',
 '9   2.1e-09  1.2e+00  2.0e-06  1.19e+00   -2.468580903e-01  -2.468555987e-01  4.4e-08  0.04  ',
 '10  1.4e-11  5.2e-01  9.0e-08  1.12e+00   -2.468579488e-01  -2.468578470e-01  2.5e-09  0.04  ',
 '11  1.5e-13  2.1e-01  2.8e-09  1.04e+00   -2.468579452e-01  -2.468579422e-01  7.9e-11  0.04  ',
 '12  2.9e-15  2.3e-02  1.0e-10  1.06e+00   -2.468579451e-01  -2.468579450e-01  3.3e-12  0.04  ',
 '13  2.1e-15  1.2e-04  5.1e-13  1.00e+00   -2.468579451e-01  -2.468579451e-01  1.6e-14  0.04  ',
 '14  2.9e-15  5.9e-07  1.2e-15  1.00e+00   -2.468579451e-01  -2.468579451e-01  8.3e-17  0.04  ',
 'Interior-point optimizer terminated. Time: 0.05. ',
 '',
 'Optimizer terminated. Time: 0.08    ',
 'solsta = 1, prosta = 1',
 '',
 'Interior-point solution summary',
 '  Problem status  : PRIMAL_AND_DUAL_FEASIBLE',
 '  Solution status : OPTIMAL',
 '  Primal.  obj: -2.4685794513e-01   Viol.  con: 2e-15    var: 0e+00  ',
 '  Dual.    obj: -2.4685794513e-01   Viol.  con: 0e+00    var: 6e-07  ',
 '* End solution on dual form. *',
 'Transforming to primal solution.']

[bqpd]

Hmm, so the solver's returning that objective of 0.78. What's gp.program print?

[pgkirsch]

Do you mean m.program?

if so:

gpkit.GeometricProgram(
  # minimize
    w_(5,),
[ # subject to
    10*V_(0,)**-1 + V_(0,)**-1*V_(1,) <= 1,
    10*V_(1,)**-1 + V_(1,)**-1*V_(2,) <= 1,
    10*V_(2,)**-1 + V_(2,)**-1*V_(3,) <= 1,
    10*V_(3,)**-1 + V_(3,)**-1*V_(4,) <= 1,
    10*V_(4,)**-1 <= 1,
    0.05*M_(0,)**-1*V_(0,) + 0.05*M_(0,)**-1*V_(1,) + M_(0,)**-1*M_(1,) <= 1,
    0.05*M_(1,)**-1*V_(1,) + 0.05*M_(1,)**-1*V_(2,) + M_(1,)**-1*M_(2,) <= 1,
    0.05*M_(2,)**-1*V_(2,) + 0.05*M_(2,)**-1*V_(3,) + M_(2,)**-1*M_(3,) <= 1,
    0.05*M_(3,)**-1*V_(3,) + 0.05*M_(3,)**-1*V_(4,) + M_(3,)**-1*M_(4,) <= 1,
    0.05*M_(4,)**-1*V_(4,) + 1.05e-06*M_(4,)**-1 <= 1,
    0.05*M_(0,)*th_(1,)**-1 + 0.05*M_(1,)*th_(1,)**-1 + 1e-06*th_(1,)**-1 <= 1,
    0.05*M_(1,)*th_(2,)**-1 + 0.05*M_(2,)*th_(2,)**-1 + th_(1,)*th_(2,)**-1 <= 1,
    0.05*M_(2,)*th_(3,)**-1 + 0.05*M_(3,)*th_(3,)**-1 + th_(2,)*th_(3,)**-1 <= 1,
    0.05*M_(3,)*th_(4,)**-1 + 0.05*M_(4,)*th_(4,)**-1 + th_(3,)*th_(4,)**-1 <= 1,
    0.05*M_(4,)*th_(5,)**-1 + 5e-08*th_(5,)**-1 + th_(4,)*th_(5,)**-1 <= 1,
    0.05*th_(1,)*w_(1,)**-1 + 1.05e-06*w_(1,)**-1 <= 1,
    0.05*th_(1,)*w_(2,)**-1 + 0.05*th_(2,)*w_(2,)**-1 + w_(1,)*w_(2,)**-1 <= 1,
    0.05*th_(2,)*w_(3,)**-1 + 0.05*th_(3,)*w_(3,)**-1 + w_(2,)*w_(3,)**-1 <= 1,
    0.05*th_(3,)*w_(4,)**-1 + 0.05*th_(4,)*w_(4,)**-1 + w_(3,)*w_(4,)**-1 <= 1,
    0.05*th_(4,)*w_(5,)**-1 + 0.05*th_(5,)*w_(5,)**-1 + w_(4,)*w_(5,)**-1 <= 1,
])

[pgkirsch]

Also m.program.result as requested:

In [6]: m.program.result
Out[6]: 
{'cost': 0.7812516662500036,
 'sensitivities': {'monomials': array([  1.00000000e+00,   1.43999693e-02,   5.75998786e-02,
           2.06399567e-01,   6.87998533e-02,   3.35999300e-01,
           1.67999642e-01,   2.99199362e-01,   2.99199392e-01,
           4.49599086e-01,   7.19998479e-02,   5.75998786e-02,
           2.30399549e-01,   1.59999663e-01,   1.19999748e-01,
           3.59999336e-01,   1.77599627e-01,   1.18399753e-01,
           2.36799637e-01,   7.19998536e-02,   1.43999700e-01,
           7.20000048e-02,   1.64639649e-07,   7.83998406e-02,
           3.59999275e-01,   2.30399549e-01,   5.75998771e-07,
           1.00799814e-01,   4.59199534e-01,   1.79199649e-01,
           3.19999510e-02,   5.27999283e-01,   7.19998672e-02,
           1.91999706e-02,   3.79199461e-01,   4.80000032e-03,
           1.34399811e-01,   3.19999323e-09,   1.60000011e-03,
           6.55999334e-02,   1.34399713e-06,   1.05599857e-01,
           6.56012774e-02,   6.55999334e-02,   1.26399820e-01,
           2.36801067e-01,   1.05599857e-01,   1.34399811e-01,
           1.26399820e-01,   4.68800744e-01,   1.35999814e-01,
           7.29600375e-01,   1.34399811e-01]),
  'posynomials': array([ 1.        ,  0.07199985,  0.27519942,  0.50399894,  0.59839875,
          0.44959909,  0.35999928,  0.63999875,  0.53279902,  0.28799956,
          0.07840001,  0.5903994 ,  0.739199  ,  0.6319991 ,  0.40319943,
          0.13599981,  0.06560128,  0.23680107,  0.46880074,  0.72960038,  1.        ])},
 'variables': {V_(4,): 10.000000999999997,
  V_(3,): 20.000001000000033,
  V_(2,): 30.000001000000061,
  V_(1,): 40.000001000000125,
  V_(0,): 50.000001000000246,
  M_(0,): 12.500001500000044,
  M_(1,): 8.0000014000000235,
  M_(2,): 4.500001300000009,
  M_(3,): 2.0000012000000025,
  M_(4,): 0.50000109999999987,
  th_(1,): 1.0250011450000049,
  th_(2,): 1.650001280000007,
  th_(3,): 1.9750014050000084,
  th_(4,): 2.1000015200000095,
  th_(5,): 2.1250016250000114,
  w_(1,): 0.051251107250000302,
  w_(2,): 0.18500122850000103,
  w_(3,): 0.36625136275000192,
  w_(4,): 0.5700015090000029,
  w_(5,): 0.7812516662500042}}

[bqpd]

@whoburg, I think that primal_exp_vals[self.m_idxs[0]].sum() in GeometricProgram.check_solution() is getting the M_(0,) variable instead of w_(5,)...

[bqpd]

@pgkirsch, in the meantime maybe comment out self.check_solution(result)?

[pgkirsch]

Original example does not fail if N < 5

[pgkirsch]

Closed by 1375c5ab3ea9738a1c26e8f44d6e7c7477e3512e

[whoburg]

Nice work tracking this down, @pgkirsch!

[pgkirsch]

Another instance of this error with Mosek

RuntimeWarning: Primal solution computed cost did not match solver-returned cost: 239.485793868 vs 0.0

cvxopt solves it ok. mosek_cli gives non-zero exit status 152, but I suspect that's related to the 0 cost.

(code at https://github.mit.edu/convex/GPjet/commit/a84544241154938cb16eec2589710e0c584b4609)

[whoburg]

Potentially related to #417, where exit status 152 also showed up.

[pgkirsch]

Potentially yes, but it's worth mentioning that mosek_cli returns exit status 152 for basically everything it can't solve. I'd say at least 95% of the unsuccessful mosek_cli messages are exit status 152. It seems to be an umbrella status under which both primal infeasible and unknown fall.

[whoburg]

This was reopened for a 0-cost MOSEK solve, similar to #454 and #453, so I'll change the milestone. (Ned)

[whoburg]

Here's a MWE for this issue, encountered while investigating #519:

from gpkit import Variable, Model

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

m = Model(F,
         [F >= D + 300,
          D == Ap])
sol = m.solve(solver="mosek")

output:

Using solver 'mosek'
Solving for 3 variables.
Solving took 0.0133 seconds.

---------------------------------------------------------------------------
RuntimeWarning                            Traceback (most recent call last)
<ipython-input-26-51b27ced75a1> in <module>()
      9          [F >= D + 300,
     10           D == Ap])
---> 11 sol = m.solve(solver="mosek")

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

/Users/whoburg/MIT/dev/gpkit/gpkit/model.pyc in _solve(self, programType, solver, verbosity, skipsweepfailures, *args, **kwargs)
    500             self.program, solvefn = form_program(programType, signomials,
    501                                                  verbosity=verbosity)
--> 502             result = solvefn(*args, **kwargs)
    503             solution.append(parse_result(result, constants, beforesubs))
    504         solution.program = self.program

/Users/whoburg/MIT/dev/gpkit/gpkit/geometric_program.pyc in solve(self, solver, verbosity, *args, **kwargs)
    200                                  " model.feasibility().")
    201 
--> 202         self.check_solution(result["cost"], primal, nu, la)
    203 
    204         if verbosity > 1:

/Users/whoburg/MIT/dev/gpkit/gpkit/geometric_program.pyc in check_solution(self, cost, primal, nu, la, tol)
    234             raise RuntimeWarning("Primal solution computed cost did not match"
    235                                  " solver-returned cost: %s vs %s" %
--> 236                                  (primal_exp_vals[self.m_idxs[0]].sum(), cost))
    237         for mi in self.m_idxs[1:]:
    238             if primal_exp_vals[mi].sum() > 1 + tol:

RuntimeWarning: Primal solution computed cost did not match solver-returned cost: 300.0 vs 0.0

[whoburg]

Here's an even more minimal MWE:

from gpkit import Variable, Model

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

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

After a lot of traceback, the output error is:

RuntimeWarning: Primal solution computed cost did not match solver-returned cost: 1.0 vs 0.0

Solved with mosek_cli instead, the output is

CalledProcessError: Command '['mskexpopt', '/var/folders/6r/98btd1dj215001_h77p4mx3m0000gn/T/tmpIL2Mwf/gpkit_mosek']' returned non-zero exit status 152

The corresponding input .eo file for mosek is

1
2
3

*c
1.00000000000000000000e+00
1.00000000000000000000e+00
1.00000000000000000000e+00

*p_idxs
0
1
1

*t j A_tj
2 0 1.00000000000000000000e+00
0 1 1.00000000000000000000e+00
1 1 -1.00000000000000000000e+00
2 1 -1.00000000000000000000e+00

Running the compiled mskexpopt program on this file at the command line, the output lists a more precise Return code 1432 at the bottom:

Reading 'gpkit_mosek'
Data setup for expopt begin.
Warning: The variable with index '0' has only positive coefficients akj.
 The problem is possibly ill-posed.
.
Number of Hessian non-zeros: 4
Data setup for expopt end.
* Solving exponential optimization problem on dual form. *
* The following log information refers to the solution of the dual problem. *
Computer
  Platform               : MACOSX/64-X86   
  Cores                  : 4               

Problem
  Name                   :                 
  Objective sense        : max             
  Type                   : GECO (general convex optimization problem)
  Constraints            : 3               
  Cones                  : 0               
  Scalar variables       : 3               
  Matrix variables       : 0               
  Integer variables      : 0               

Optimizer started.
Interior-point optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator - tries                  : 0                 time                   : 0.00            
Eliminator - elim's                 : 0               
Lin. dep.  - tries                  : 1                 time                   : 0.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.00    
Interior-point optimizer terminated. Time: 0.00. 

Optimizer terminated. Time: 0.01    
solsta = 1, prosta = 1

Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
* End solution on dual form. *
Transforming to primal solution.

Interior-point solution summary
  Problem status  : PRIMAL_AND_DUAL_FEASIBLE
  Solution status : OPTIMAL
Return code: 1432
Description: MSK_RES_ERR_USER_NLO_FUNC [The user-defined nonlinear function reported an error.]

Not sure how best to proceed -- seems like we could

1. try to catch the 1432 return code or MSK_RES_ERR_USER_NLO_FUNC somehow
2. complain to mosek that solutions that are PRIMAL_AND_DUAL_FEASIBLE and OPTIMAL should not also have error return codes ? I also inspected the result in m.program.result, and it's not a valid solution -- the cost check is just the first check that it fails.