opposite dual values from cyipopt

[ZedongPeng]

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

[Robbybp]

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).

[michaelbynum]

Well said, @Robbybp

[mrmundt]

@ZedongPeng - Can this be closed via #2830 ?

[Robbybp]

Assuming #2830 is localized to MindtPy, it will not fix this issue.

[ZedongPeng]

2830 will not fix this issue.

[ZedongPeng]

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))

[Robbybp]

@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.