pyomo
pyomo copied to clipboard
Published
20 hours ago •
Pyomo
Pyomo
opposite dual values from cyipopt
Summary
I found an issue with cyipopt. If the objective sense of the model is minimization, the optimal duals provided by ipopt and cyipopt are opposite. Does this make sense?
Steps to reproduce the issue
from pyomo.environ import *
from pyomo.contrib.mindtpy.tests.eight_process_problem import EightProcessFlowsheet
from pyomo.contrib.mindtpy.tests.MINLP_simple import SimpleMINLP
m = SimpleMINLP()
# m = EightProcessFlowsheet(convex=True)
# m.objective.sense = maximize
TransformationFactory('core.relax_integer_vars').apply_to(m)
m.pprint()
m.dual = Suffix(direction=Suffix.IMPORT)
opt = SolverFactory('cyipopt')
result = opt.solve(m)
m.dual.pprint()
opt = SolverFactory('ipopt')
result = opt.solve(m)
m.dual.pprint()
Error Message
Dual from cyipopt
dual : Direction=Suffix.IMPORT, Datatype=Suffix.FLOAT
Key : Value
const1 : 1.2652108928398336
const2 : -6.237267596073834e-09
const3 : 7.012918179422071e-10
const4 : -4.733532397755043e-09
const5 : -3.4888414554416926e-09
const6 : 0.16666666746603248
const7 : -1.0000000050889761
Dual from ipopt
dual : Direction=Suffix.IMPORT, Datatype=Suffix.FLOAT
Key : Value
const1 : -1.2652108928434915
const2 : 6.234279237166242e-09
const3 : -7.013476064517299e-10
const4 : 4.734977354006692e-09
const5 : 3.48921680082169e-09
const6 : -0.16666666746193418
const7 : 1.0000000050778406
Information on your system
Pyomo version: 6.5.1.dev0 Python version: 3.7.12 Operating system: macOS Monterey Version 12.2 How Pyomo was installed (PyPI, conda, source): source Solver (if applicable): cyipopt, ipopt
Additional information
I believe this is because Ipopt's internal data structures and AMPL interface use different multiplier conventions. I believe the rationale was for the AMPL interface to be consistent with CPLEX's and Gurobi's dual convention.
In our "CyIpopt solver", we send duals straight from CyIpopt to the PyomoNLP:
nlp.set_duals(info["mult_g"]) # line 399 in cyipopt_solver.py
In PyomoNLP, we update the "dual" suffix straight from the NLP:
model_suffixes['dual'].update(zip(constraints, duals)) # line 444 in pyomo_nlp.py
A quick "fix" is to apply the "correct" factors of -1 in the CyIpopt solver. This will need to be done carefully (and with lots of clear tests) as IIRC the convention changes for maximization problems (and I forget if there is a difference for bound/inequality multipliers). However, if somebody knows about Ipopt's internal convention and is using these multipliers, this will break their code.
I think updating to match the Ipopt AMPL interface is probably the right solution, as this is what I think most people are familiar with. Obviously, this should be well-documented. And ideally, we somehow enforce that all solvers follow this convention (or somehow define what convention they are using).
The dual of cyipopt gave me a hard time. I am now using cyipopt to solve the grey box model.
If the grey box exits in the NLP model and the sense of the objective function is minimization, the dual provided by ipopt and cyipopt is the same. This contradicts with the above result.
Do you have any ideas about this? @Robbybp @michaelbynum
Many thanks.
FYI @bernalde .
Code to reproduce
from pyomo.environ import *
from pyomo.contrib.mindtpy.tests.eight_process_problem import EightProcessFlowsheet
from pyomo.contrib.mindtpy.tests.MINLP_simple import SimpleMINLP
m = SimpleMINLP(grey_box=True)
TransformationFactory('core.relax_integer_vars').apply_to(m)
m.dual = Suffix(direction=Suffix.IMPORT)
opt = SolverFactory('cyipopt')
opt.config.options['hessian_approximation'] = 'limited-memory'
opt.config.options['linear_solver'] = 'mumps'
result = opt.solve(m)
m.dual.pprint()
print(value(m.objective))
m = SimpleMINLP()
TransformationFactory('core.relax_integer_vars').apply_to(m)
# m.objective.sense = maximize
m.dual = Suffix(direction=Suffix.IMPORT)
opt = SolverFactory('ipopt')
result = opt.solve(m)
m.dual.pprint()
print(value(m.objective))
Output
dual : Direction=Suffix.IMPORT, Datatype=Suffix.FLOAT
Key : Value
con_X_1 : -2.3077366781562283
con_X_2 : -1.4318775699349497
con_Y_1 : -1.0000000035172412
con_Y_2 : -1.2500000032757335
con_Y_3 : -0.5000000022459474
const1 : -1.2652109007743766
const2 : 2.0358765165068324e-09
const3 : -1.6181781523874018e-09
const4 : -4.149376754573334e-10
const5 : -1.4658355348774077e-10
const6 : -0.16666666855761839
const7 : 1.0000000088083274
my_block.egb : {'my_block.egb.output_constraints[z]': -1.0000000027128453}
2.5323459411442015
WARNING: Implicitly replacing the Component attribute objective (type=<class
'pyomo.core.base.objective.ScalarObjective'>) on block SimpleMINLP with a new
Component (type=<class 'pyomo.core.base.objective.ScalarObjective'>). This is
usually indicative of a modelling error. To avoid this warning, use
block.del_component() and block.add_component().
dual : Direction=Suffix.IMPORT, Datatype=Suffix.FLOAT
Key : Value
const1 : -1.2652108928398416
const2 : 6.23725758862012e-09
const3 : -7.012918728680585e-10
const4 : 4.73353230919436e-09
const5 : 3.488841821543732e-09
const6 : -0.1666666674660251
const7 : 1.0000000050889548
const8 : -1.0
2.5323459511417323
Related import files
# MINLP_simple.py
# ___________________________________________________________________________
#
# Pyomo: Python Optimization Modeling Objects
# Copyright (c) 2008-2022
# National Technology and Engineering Solutions of Sandia, LLC
# Under the terms of Contract DE-NA0003525 with National Technology and
# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
# rights in this software.
# This software is distributed under the 3-clause BSD License.
# ___________________________________________________________________________
"""Implementation of MINLP problem in Assignment 6 of the Advanced PSE lecture at CMU.
Author: David Bernal <https://github.com/bernalde>
The expected optimal solution is 3.5.
Ref:
IGNACIO GROSSMANN.
CARNEGIE-MELLON UNIVERSITY , PITTSBURGH , PA.
Problem type: convex MINLP
size: 3 binary variables
3 continuous variables
7 constraints
"""
from pyomo.environ import (
Binary,
ConcreteModel,
Constraint,
NonNegativeReals,
Objective,
RangeSet,
Var,
minimize,
Reals,
Block,
maximize,
)
from pyomo.common.collections import ComponentMap
from pyomo.contrib.mindtpy.tests.MINLP_simple_grey_box import GreyBoxModel
from pyomo.common.dependencies import attempt_import
egb = attempt_import('pyomo.contrib.pynumero.interfaces.external_grey_box')[0]
def build_model_external(m):
ex_model = GreyBoxModel(initial={"X1": 0, "X2": 0, "Y1": 0, "Y2": 1, "Y3": 1})
m.egb = egb.ExternalGreyBoxBlock()
m.egb.set_external_model(ex_model)
class SimpleMINLP(ConcreteModel):
"""Convex MINLP problem Assignment 6 APSE."""
def __init__(self, grey_box=False, *args, **kwargs):
"""Create the problem."""
kwargs.setdefault('name', 'SimpleMINLP')
super(SimpleMINLP, self).__init__(*args, **kwargs)
m = self
"""Set declarations"""
I = m.I = RangeSet(1, 2, doc='continuous variables')
J = m.J = RangeSet(1, 3, doc='discrete variables')
# initial point information for discrete variables
initY = {
'sub1': {1: 1, 2: 1, 3: 1},
'sub2': {1: 0, 2: 1, 3: 1},
'sub3': {1: 1, 2: 0, 3: 1},
'sub4': {1: 1, 2: 1, 3: 0},
'sub5': {1: 0, 2: 0, 3: 0},
}
# initial point information for continuous variables
initX = {1: 0, 2: 0}
"""Variable declarations"""
# DISCRETE VARIABLES
Y = m.Y = Var(J, domain=Binary, initialize=initY['sub2'])
# CONTINUOUS VARIABLES
X = m.X = Var(I, domain=NonNegativeReals, initialize=initX)
"""Constraint definitions"""
# CONSTRAINTS
m.const1 = Constraint(expr=(m.X[1] - 2) ** 2 - m.X[2] <= 0)
m.const2 = Constraint(expr=m.X[1] - 2 * m.Y[1] >= 0)
m.const3 = Constraint(expr=m.X[1] - m.X[2] - 4 * (1 - m.Y[2]) <= 0)
m.const4 = Constraint(expr=m.X[1] - (1 - m.Y[1]) >= 0)
m.const5 = Constraint(expr=m.X[2] - m.Y[2] >= 0)
m.const6 = Constraint(expr=m.X[1] + m.X[2] >= 3 * m.Y[3])
m.const7 = Constraint(expr=m.Y[1] + m.Y[2] + m.Y[3] >= 1)
"""Cost (objective) function definition"""
m.objective = Objective(
expr=Y[1] + 1.5 * Y[2] + 0.5 * Y[3] + X[1] ** 2 + X[2] ** 2, sense=minimize
)
if not grey_box:
m.obj = Var(domain=Reals)
m.const8 = Constraint(expr=Y[1] + 1.5 * Y[2] + 0.5 * Y[3] + X[1] ** 2 + X[2] ** 2==m.obj)
m.objective = Objective(expr=m.obj,sense=minimize)
# m.objective = Objective(
# expr=Y[1] + 1.5 * Y[2] + 0.5 * Y[3] + X[1] ** 2 + X[2] ** 2,
# sense=minimize,
# )
else:
def _model_i(b):
build_model_external(b)
m.my_block = Block(rule=_model_i)
for i in m.I:
def eq_inputX(m):
return m.X[i] == m.my_block.egb.inputs["X" + str(i)]
con_name = "con_X_" + str(i)
m.add_component(con_name, Constraint(expr=eq_inputX))
for j in m.J:
def eq_inputY(m):
return m.Y[j] == m.my_block.egb.inputs["Y" + str(j)]
con_name = "con_Y_" + str(j)
m.add_component(con_name, Constraint(expr=eq_inputY))
# add objective
m.objective = Objective(expr=m.my_block.egb.outputs['z'], sense=minimize)
"""Bound definitions"""
# x (continuous) upper bounds
x_ubs = {1: 4, 2: 4}
for i, x_ub in x_ubs.items():
X[i].setub(x_ub)
m.optimal_value = 3.5
m.optimal_solution = ComponentMap()
m.optimal_solution[m.X[1]] = 1.0
m.optimal_solution[m.X[2]] = 1.0
m.optimal_solution[m.Y[1]] = 0.0
m.optimal_solution[m.Y[2]] = 1.0
m.optimal_solution[m.Y[3]] = 0.0
# MINLP_simple_grey_box.py
from pyomo.common.dependencies import numpy as np
import pyomo.common.dependencies.scipy.sparse as scipy_sparse
from pyomo.common.dependencies import attempt_import
egb = attempt_import('pyomo.contrib.pynumero.interfaces.external_grey_box')[0]
class GreyBoxModel(egb.ExternalGreyBoxModel):
"""Greybox model to compute the example OF."""
def __init__(self, initial, use_exact_derivatives=True, verbose=True):
"""
Parameters
use_exact_derivatives: bool
If True, the exact derivatives are used. If False, the finite difference
approximation is used.
verbose: bool
If True, print information about the model.
"""
self._use_exact_derivatives = use_exact_derivatives
self.verbose = verbose
self.initial = initial
# For use with exact Hessian
self._output_con_mult_values = np.zeros(1)
if not use_exact_derivatives:
raise NotImplementedError("use_exact_derivatives == False not supported")
def input_names(self):
"""Return the names of the inputs."""
self.input_name_list = ["X1", "X2", "Y1", "Y2", "Y3"]
return self.input_name_list
def equality_constraint_names(self):
"""Return the names of the equality constraints."""
# no equality constraints
return []
def output_names(self):
"""Return the names of the outputs."""
return ['z']
def set_output_constraint_multipliers(self, output_con_multiplier_values):
"""Set the values of the output constraint multipliers."""
# because we only have one output constraint
assert len(output_con_multiplier_values) == 1
np.copyto(self._output_con_mult_values, output_con_multiplier_values)
def finalize_block_construction(self, pyomo_block):
"""Finalize the construction of the ExternalGreyBoxBlock."""
if self.initial is not None:
print("initialized")
pyomo_block.inputs["X1"].value = self.initial["X1"]
pyomo_block.inputs["X2"].value = self.initial["X2"]
pyomo_block.inputs["Y1"].value = self.initial["Y1"]
pyomo_block.inputs["Y2"].value = self.initial["Y2"]
pyomo_block.inputs["Y3"].value = self.initial["Y3"]
else:
print("uninitialized")
for n in self.input_name_list:
pyomo_block.inputs[n].value = 1
pyomo_block.inputs["X1"].setub(4)
pyomo_block.inputs["X1"].setlb(0)
pyomo_block.inputs["X2"].setub(4)
pyomo_block.inputs["X2"].setlb(0)
pyomo_block.inputs["Y1"].setub(1)
pyomo_block.inputs["Y1"].setlb(0)
pyomo_block.inputs["Y2"].setub(1)
pyomo_block.inputs["Y2"].setlb(0)
pyomo_block.inputs["Y3"].setub(1)
pyomo_block.inputs["Y3"].setlb(0)
def set_input_values(self, input_values):
"""Set the values of the inputs."""
self._input_values = list(input_values)
def evaluate_equality_constraints(self):
"""Evaluate the equality constraints."""
# Not sure what this function should return with no equality constraints
return None
def evaluate_outputs(self):
"""Evaluate the output of the model."""
# form matrix as a list of lists
# M = self._extract_and_assemble_fim()
x1 = self._input_values[0]
x2 = self._input_values[1]
y1 = self._input_values[2]
y2 = self._input_values[3]
y3 = self._input_values[4]
# z
z = x1**2 + x2**2 + y1 + 1.5 * y2 + 0.5 * y3
if self.verbose:
pass
# print("\n Consider inputs [x1,x2,y1,y2,y3] =\n",x1, x2, y1, y2, y3)
# print(" z = ",z,"\n")
return np.asarray([z], dtype=np.float64)
def evaluate_jacobian_equality_constraints(self):
"""Evaluate the Jacobian of the equality constraints."""
return None
'''
def _extract_and_assemble_fim(self):
M = np.zeros((self.n_parameters, self.n_parameters))
for i in range(self.n_parameters):
for k in range(self.n_parameters):
M[i,k] = self._input_values[self.ele_to_order[(i,k)]]
return M
'''
def evaluate_jacobian_outputs(self):
"""Evaluate the Jacobian of the outputs."""
if self._use_exact_derivatives:
# compute gradient of log determinant
row = np.zeros(5) # to store row index
col = np.zeros(5) # to store column index
data = np.zeros(5) # to store data
row[0], col[0], data[0] = (0, 0, 2 * self._input_values[0]) # x1
row[0], col[1], data[1] = (0, 1, 2 * self._input_values[1]) # x2
row[0], col[2], data[2] = (0, 2, 1) # y1
row[0], col[3], data[3] = (0, 3, 1.5) # y2
row[0], col[4], data[4] = (0, 4, 0.5) # y3
# sparse matrix
return scipy_sparse.coo_matrix((data, (row, col)), shape=(1, 5))
@ZedongPeng Here is the load_state_into_pyomo function of PyomoNLPwithGreyBoxBlocks, which is called by our CyIpopt interface:
def load_state_into_pyomo(self, bound_multipliers=None):
# load the values of the primals into the pyomo
primals = self.get_primals()
for value, vardata in zip(primals, self._pyomo_model_var_datas):
vardata.set_value(value)
# get the active suffixes
m = self._pyomo_model
model_suffixes = dict(pyo.suffix.active_import_suffix_generator(m))
# we need to correct the sign of the multipliers based on whether or
# not we are minimizing or maximizing - this is done in the ASL interface
# for ipopt, but does not appear to be done in cyipopt.
obj_sign = 1.0
objs = list(m.component_data_objects(ctype=pyo.Objective, descend_into=True))
assert len(objs) == 1
if objs[0].sense == pyo.maximize:
obj_sign = -1.0
if 'dual' in model_suffixes:
model_suffixes['dual'].clear()
dual_values = self._dual_values_blockvector.flatten()
for value, t in zip(dual_values, self._constraint_datas):
if type(t) is tuple:
model_suffixes['dual'].setdefault(t[0], {})[t[1]] = (
-obj_sign * value
)
else:
# t is a constraint data
model_suffixes['dual'][t] = -obj_sign * value
if 'ipopt_zL_out' in model_suffixes:
model_suffixes['ipopt_zL_out'].clear()
if bound_multipliers is not None:
model_suffixes['ipopt_zL_out'].update(
zip(self._pyomo_model_var_datas, obj_sign * bound_multipliers[0])
)
if 'ipopt_zU_out' in model_suffixes:
model_suffixes['ipopt_zU_out'].clear()
if bound_multipliers is not None:
model_suffixes['ipopt_zU_out'].update(
zip(self._pyomo_model_var_datas, -obj_sign * bound_multipliers[1])
)
This is in contrast to the same method in PyomoNLP:
def load_state_into_pyomo(self, bound_multipliers=None):
primals = self.get_primals()
variables = self.get_pyomo_variables()
for var, val in zip(variables, primals):
var.set_value(val)
m = self.pyomo_model()
model_suffixes = dict(pyo.suffix.active_import_suffix_generator(m))
if 'dual' in model_suffixes:
duals = self.get_duals()
constraints = self.get_pyomo_constraints()
model_suffixes['dual'].clear()
model_suffixes['dual'].update(zip(constraints, duals))
if 'ipopt_zL_out' in model_suffixes:
model_suffixes['ipopt_zL_out'].clear()
if bound_multipliers is not None:
model_suffixes['ipopt_zL_out'].update(
zip(variables, bound_multipliers[0])
)
if 'ipopt_zU_out' in model_suffixes:
model_suffixes['ipopt_zU_out'].clear()
if bound_multipliers is not None:
model_suffixes['ipopt_zU_out'].update(
zip(variables, bound_multipliers[1])
)
The difference seems to be that PyomoNLPwithGreyBoxBlocks anticipates the difference in convention with the AMPL interface and attempts to adjust the multipliers to match.