Description
Summary
Ipopt appears to use the following convention for the sign of its bound multipliers:
| objective type | zL | zU |
|---|---|---|
| maximize | <= 0 | >= 0 |
| minimize | >= 0 | <= 0 |
I believe that, for an NLP with equality constraint bodies c, inequality constraint bodies d (here these are just, e.g. x - xL), and a Lagrangian written
f + c^T * y + d^T * z
that the objective type, multiplier sign, and form of inequality constraint (i.e. d <= 0 vs d >= 0) must belong to the same row of the following table:
| objective | multiplier sign | inequality form |
|---|---|---|
| max | - | >= 0 |
| max | + | <= 0 |
| min | + | >= 0 |
| min | - | <= 0 |
Ipopt's sign convention for the multipliers suggests that it considers variable bounds to be inequalities of the following forms, regardless of objective type
x - xL >= 0
x - xU <= 0
However, the interior point interface appears to write the barrier Lagrangian (in evaluate_primal_dual_kkt_rhs) under the assumption that the Lagrangian of the original problem is:
f + c^T * y - (x - xL)^T * zL - (x - xU)^T * zU
When I compute the primal gradient of this Lagrangian after solving with Ipopt, with duals returned by Ipopt, I do not get a zero vector. However, when I compute the primal gradient of the Lagrangian that I would expect Ipopt to use,
f + c^T * y + (x - xL)^T * zL + (x - xU)^T * zU
I also do not get a zero vector. To get a primal gradient that is zero at the solution of a simple NLP, I need to use the following Lagrangian:
-f + c^T * y + (x - xL)^T * zL + (x - xU)^T * zU
This is thoroughly confusing to me. So I have two issues, only one of which relates to Pyomo. The first, unrelated, is that I do not know how Ipopt constructs its Lagrangian. The Pyomo issue is that, however Ipopt constructs its Lagrangian, the Pyomo interior point interface does not seem to be consistent. My questions are: Are others aware of this? Have I missed anything? Why is my primal Lagrangian gradient (grad_lag_vec below) not zero?
Steps to reproduce the issue
I have been using the following script to test this. grad_lag_map, which contains the gradient of the last Lagrangian above, seems to consistently be zero. evaluate_grad_lag_primals, which is my attempt to mirror the interior point interface's barrier Lagrangian gradient, for the original problem, does not return zero.
import pyomo.environ as pyo
from pyomo.common.collections import ComponentMap
from pyomo.contrib.interior_point.interface import InteriorPointInterface
m = pyo.ConcreteModel()
m.ipopt_zL_out = pyo.Suffix(direction=pyo.Suffix.IMPORT)
m.ipopt_zU_out = pyo.Suffix(direction=pyo.Suffix.IMPORT)
m.ipopt_zL_in = pyo.Suffix(direction=pyo.Suffix.EXPORT)
m.ipopt_zU_in = pyo.Suffix(direction=pyo.Suffix.EXPORT)
m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT_EXPORT)
m.v1 = pyo.Var(initialize=2.0)
m.v2 = pyo.Var(initialize=2.0)
m.v3 = pyo.Var(initialize=2.0)
m.v1.setlb(1.5)
m.v2.setlb(1.5)
m.v1.setub(10.0)
m.v2.setub(10.0)
m.eq_con = pyo.Constraint(expr=m.v1*m.v2*m.v3 - 2.0 == 0)
obj_factor = 1
m.obj = pyo.Objective(
expr=obj_factor*(m.v1**2 + m.v2**2 + m.v3**2),
sense=pyo.minimize,
)
solver = pyo.SolverFactory("ipopt")
solver.solve(m, tee=True)
def evaluate_grad_lag_primals(ip):
nlp = ip._nlp
grad_obj = nlp.evaluate_grad_objective()
jac_eq = nlp.evaluate_jacobian_eq()
jac_ineq = nlp.evaluate_jacobian_ineq()
duals_eq = nlp.get_duals_eq()
duals_ineq = nlp.get_duals_ineq()
duals_primals_lb = ip.get_duals_primals_lb()
duals_primals_ub = ip.get_duals_primals_ub()
grad_lag = (
grad_obj +
jac_eq.transpose() * duals_eq +
jac_ineq.transpose() * duals_ineq -
duals_primals_lb -
duals_primals_ub
)
return grad_lag
ip = InteriorPointInterface(m)
ip.set_barrier_parameter(0)
nlp = ip._nlp
rhs = ip.evaluate_primal_dual_kkt_rhs()
var_vec = nlp.get_pyomo_variables()
#grad_lag_vec = -rhs.get_block(0)
# Evaluating grad_lag via kkt rhs is unreliable at the solution
# if variables are exactly at their bounds as it will divide by zero.
# This is why I define a function to compute grad_lag of the original
# problem, which is more what I'm interested in anyway.
grad_lag_vec = evaluate_grad_lag_primals(ip)
grad_lag_map = ComponentMap()
grad_lag_map[m.v1] = (
-(obj_factor*2*m.v1) + m.dual[m.eq_con]*m.v2*m.v3 +
m.ipopt_zL_out[m.v1] + m.ipopt_zU_out[m.v1]
)
grad_lag_map[m.v2] = (
-(obj_factor*2*m.v2) + m.dual[m.eq_con]*m.v1*m.v3 +
m.ipopt_zL_out[m.v2] + m.ipopt_zU_out[m.v2]
)
grad_lag_map[m.v3] = (
-(obj_factor*2*m.v3) + m.dual[m.eq_con]*m.v1*m.v2
)
print("grad lag constructed manually")
for var, expr in grad_lag_map.items():
#print(var.name, expr)
print(var.name, pyo.value(expr))
print("grad lag constructed via ip interface")
for var, val in zip(var_vec, grad_lag_vec):
print(var.name, val)Information on your system
Pyomo version: 6.0.dev0
Python version: 3.7.8
Ipopt verson: 3.12