Integer Programs - Optimization for Decision Science

# This code cell installs packages on Colab

import sys

if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    import helper

    helper.easy_install()
else:
    sys.path.insert(0, "../")
    import helper
helper.set_plotting_style()

Optimizing Across Process Alternatives¶

Reference: Example 15.3 from Biegler, Grossmann, and Westerberg (1997). See handout

Assume that we have the choice of selecting two reactors (shown below) for the reaction $A \rightarrow B$ . Reactor I has a higher conversion (80%) but it is more expensive; reactor II has a lower conversion (66.7%) but is cheaper. The cost of feed $A$ is $5/kmol. Which process alternative (reactor I, reactor II, or both) has the minimum costs to make 10 kmol/hr of product B?

Let $x_r$ be the size (flowrate into) of reactor $r$ .

Continuous cost model: $C = \underbrace{a_r (x_r)^{0.6}}_{\text{reactor}} + \underbrace{5.0 x_r}_{\text{feed}}$

Develop the Optimization Model¶

Draw a Picture

Reactor selection superstructure: feed A splits between two reactors, products recombine to give B.

Sets

Parameters

Variables

Objective

Constraints

Degree of Freedom Analysis

Solve with Continuous Cost Model in Pyomo¶

We start by defining the model in Pyomo.

Units expose a problem in the cost model¶

Declaring units on this model does not just document it --- it fails. The reactor capital cost $a_r x_r^{0.6}$ has $a_r$ in $/hr and $x_r$ in kmol/hr, so the product comes out in $\mathrm{\$ \cdot kmol^{0.6} \, hr^{-1.6}}$ , not $/hr. Pyomo raises UnitsError and it is right to: a power law of a dimensioned quantity is not dimensionally homogeneous.

Economy-of-scale correlations always carry an implicit reference scale. Writing it down,

C_r = a_r \left( \frac{x_r}{x_{\mathrm{ref}}} \right)^{0.6}, \qquad x_{\mathrm{ref}} = 1~\mathrm{kmol/hr},

(1)

makes the base of the power dimensionless and the cost $/hr. The numbers do not change, because $x_{\mathrm{ref}} = 1$ --- but the units now check, and the reader can see which flowrate unit the fitted coefficients 5.5 and 4.0 belong to. That information was previously nowhere in the model.

import pyomo.environ as pyo
from pyomo.environ import units as u
from pyomo.util.check_units import assert_units_consistent

# Pyomo's unit library has kmol and hr but no money, so declare a currency.
u.load_definitions_from_strings(["USD = [currency]"])

nlp = pyo.ConcreteModel()

## Define sets
nlp.REACTORS = pyo.Set(initialize=range(1, 3))

## Define parameters (data)

# Cost coefficient of reactor r, a_r [USD/hr]
cost_coefficient = {1: 5.5, 2: 4.0}
nlp.reactor_cost = pyo.Param(
    nlp.REACTORS, initialize=cost_coefficient, units=u.USD / u.hr
)

# Required product rate, P [kmol B/hr]
nlp.product_flowrate = pyo.Param(initialize=10.0, units=u.kmol / u.hr)

# Conversion of reactor r, theta_r [kmol B/kmol A]
reactor_conversion = {1: 0.8, 2: 2 / 3}
nlp.conversion = pyo.Param(nlp.REACTORS, initialize=reactor_conversion)

# Price of feed A, c_f [USD/kmol]
nlp.feed_cost = pyo.Param(initialize=5.0, units=u.USD / u.kmol)

# Reference flowrate for the economy-of-scale term, x_ref [kmol/hr]
nlp.flowrate_ref = pyo.Param(initialize=1.0, units=u.kmol / u.hr)


## Define variables

# Total feed of A, x_0 in the notes [kmol/hr]
nlp.feed_flowrate = pyo.Var(domain=pyo.NonNegativeReals, units=u.kmol / u.hr)

# Feed into reactor r, x_r in the notes [kmol/hr]
nlp.reactor_feed = pyo.Var(
    nlp.REACTORS, domain=pyo.NonNegativeReals, units=u.kmol / u.hr
)

# Effluent leaving reactor r, z_r in the notes [kmol B/hr]
nlp.reactor_effluent = pyo.Var(
    nlp.REACTORS, domain=pyo.NonNegativeReals, units=u.kmol / u.hr
)

## Define constraints

# Add your solution here

## Define objective


# Total cost [USD/hr]. The exponent 0.6 is the economy-of-scale factor beta,
# applied to the flowrate scaled by x_ref so that the power is dimensionless.
@nlp.Objective()
def cost(b):
    return (
        sum(
            b.reactor_cost[r] * (b.reactor_feed[r] / b.flowrate_ref) ** (0.6)
            for r in b.REACTORS
        )
        + b.feed_cost * b.feed_flowrate
    )

# Inspect the model
nlp.pprint()

# Raises UnitsError if any constraint or the objective is dimensionally inconsistent
assert_units_consistent(nlp)

print("Units are consistent.")

1 Set Declarations
    REACTORS : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    2 : {1, 2}

5 Param Declarations
    conversion : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=False
        Key : Value
          1 :                0.8
          2 : 0.6666666666666666
    feed_cost : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=USD/kmol
        Key  : Value
        None :   5.0
    flowrate_ref : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=kmol/h
        Key  : Value
        None :   1.0
    product_flowrate : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=kmol/h
        Key  : Value
        None :  10.0
    reactor_cost : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=True, Units=USD/h
        Key : Value
          1 :   5.5
          2 :   4.0

3 Var Declarations
    feed_flowrate : Size=1, Index=None, Units=kmol/h
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :  None :  None : False :  True : NonNegativeReals
    reactor_effluent : Size=2, Index=REACTORS, Units=kmol/h
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          1 :     0 :  None :  None : False :  True : NonNegativeReals
          2 :     0 :  None :  None : False :  True : NonNegativeReals
    reactor_feed : Size=2, Index=REACTORS, Units=kmol/h
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          1 :     0 :  None :  None : False :  True : NonNegativeReals
          2 :     0 :  None :  None : False :  True : NonNegativeReals

1 Objective Declarations
    cost : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : reactor_cost[1]*(1/flowrate_ref*reactor_feed[1])**0.6 + reactor_cost[2]*(1/flowrate_ref*reactor_feed[2])**0.6 + feed_cost*feed_flowrate

3 Constraint Declarations
    inlet_split : Size=1, Index=None, Active=True
        Key  : Lower : Body                                                : Upper : Active
        None :   0.0 : feed_flowrate - (reactor_feed[1] + reactor_feed[2]) :   0.0 :   True
    mixer : Size=1, Index=None, Active=True
        Key  : Lower            : Body                                      : Upper            : Active
        None : product_flowrate : reactor_effluent[1] + reactor_effluent[2] : product_flowrate :   True
    reactor_performance : Size=2, Index=REACTORS, Active=True
        Key : Lower : Body                                                     : Upper : Active
          1 :   0.0 :                reactor_effluent[1] - 0.8*reactor_feed[1] :   0.0 :   True
          2 :   0.0 : reactor_effluent[2] - 0.6666666666666666*reactor_feed[2] :   0.0 :   True

13 Declarations: REACTORS reactor_cost product_flowrate conversion feed_cost flowrate_ref feed_flowrate reactor_feed reactor_effluent inlet_split reactor_performance mixer cost
Units are consistent.

# Mass balance over the splitter [kmol/hr]
@nlp.Constraint()
def inlet_split(b):
    return b.feed_flowrate == sum(b.reactor_feed[r] for r in b.REACTORS)


# Reactor conversion, z_r = theta_r * x_r [kmol B/hr]
@nlp.Constraint(nlp.REACTORS)
def reactor_performance(b, r):
    return b.reactor_effluent[r] == b.conversion[r] * b.reactor_feed[r]


# Mass balance over the mixer meets product requirements [kmol B/hr]
@nlp.Constraint()
def mixer(b):
    return b.product_flowrate == sum(b.reactor_effluent[r] for r in b.REACTORS)

Initialize to Favor Reaction 1 and Solve¶

def initialize(model, reactor_choice=1):
    """Initialize all of the variables in the model to demonstrate local solutions

    Arguments:
        model: Pyomo model
        reactor_choice: 1 or 2

    Returns:
        nothing

    Action:
        initializes model

    """

    # Guess 20 kmol/hr feed of A
    model.feed_flowrate = 20.0

    # Either assign all of the feed to reactor 1 or 2
    if reactor_choice == 1:
        model.reactor_feed[1] = 20.0
        model.reactor_feed[2] = 0
    elif reactor_choice == 2:
        model.reactor_feed[1] = 0
        model.reactor_feed[2] = 20.0
    else:
        raise ValueError("Argument reactor_choice needs value 1 or 2.")

    # Based on the feed assignments, calculate effluent flowrate
    for r in model.REACTORS:
        model.reactor_effluent[r] = model.reactor_feed[r]() * model.conversion[r]


initialize(nlp, reactor_choice=1)
nlp.pprint()

1 Set Declarations
    REACTORS : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    2 : {1, 2}

5 Param Declarations
    conversion : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=False
        Key : Value
          1 :                0.8
          2 : 0.6666666666666666
    feed_cost : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=USD/kmol
        Key  : Value
        None :   5.0
    flowrate_ref : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=kmol/h
        Key  : Value
        None :   1.0
    product_flowrate : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=kmol/h
        Key  : Value
        None :  10.0
    reactor_cost : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=True, Units=USD/h
        Key : Value
          1 :   5.5
          2 :   4.0

3 Var Declarations
    feed_flowrate : Size=1, Index=None, Units=kmol/h
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :  20.0 :  None : False : False : NonNegativeReals
    reactor_effluent : Size=2, Index=REACTORS, Units=kmol/h
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          1 :     0 :  16.0 :  None : False : False : NonNegativeReals
          2 :     0 :   0.0 :  None : False : False : NonNegativeReals
    reactor_feed : Size=2, Index=REACTORS, Units=kmol/h
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          1 :     0 :  20.0 :  None : False : False : NonNegativeReals
          2 :     0 :     0 :  None : False : False : NonNegativeReals

1 Objective Declarations
    cost : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : reactor_cost[1]*(1/flowrate_ref*reactor_feed[1])**0.6 + reactor_cost[2]*(1/flowrate_ref*reactor_feed[2])**0.6 + feed_cost*feed_flowrate

3 Constraint Declarations
    inlet_split : Size=1, Index=None, Active=True
        Key  : Lower : Body                                                : Upper : Active
        None :   0.0 : feed_flowrate - (reactor_feed[1] + reactor_feed[2]) :   0.0 :   True
    mixer : Size=1, Index=None, Active=True
        Key  : Lower            : Body                                      : Upper            : Active
        None : product_flowrate : reactor_effluent[1] + reactor_effluent[2] : product_flowrate :   True
    reactor_performance : Size=2, Index=REACTORS, Active=True
        Key : Lower : Body                                                     : Upper : Active
          1 :   0.0 :                reactor_effluent[1] - 0.8*reactor_feed[1] :   0.0 :   True
          2 :   0.0 : reactor_effluent[2] - 0.6666666666666666*reactor_feed[2] :   0.0 :   True

13 Declarations: REACTORS reactor_cost product_flowrate conversion feed_cost flowrate_ref feed_flowrate reactor_feed reactor_effluent inlet_split reactor_performance mixer cost

Now let’s solve the model.

from pyomo.common.errors import ApplicationError

solver = pyo.SolverFactory("ipopt")
try:
    results = solver.solve(nlp, tee=True)
except ApplicationError as error:
    # This failure is intentional: the next cell explains the singular derivative.
    print(f"Expected solver failure: {type(error).__name__}: {error}")

Error evaluating

"var =" definition -1: can't evaluate pow'(0,0.6).

Ipopt 3.14.19: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:        9
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        2

ERROR: Solver (ipopt) returned non-zero return code (1)

ERROR: See the solver log above for diagnostic information.

---------------------------------------------------------------------------
ApplicationError                          Traceback (most recent call last)
Cell In[5], line 2
      1 solver = pyo.SolverFactory("ipopt")
----> 2 results = solver.solve(nlp, tee=True)

File ~/opt/anaconda3/envs/optimization_fall2026/lib/python3.13/site-packages/pyomo/opt/base/solvers.py:657, in OptSolver.solve(self, *args, **kwds)
    655     elif hasattr(_status, 'log') and _status.log:
    656         logger.error("Solver log:\n" + str(_status.log))
--> 657     raise ApplicationError("Solver (%s) did not exit normally" % self.name)
    658 solve_completion_time = time.time()
    659 if self._report_timing:

ApplicationError: Solver (ipopt) did not exit normally

What happened? 0^0.6 is fine --- it is 0. The derivative is not: $\frac{d}{dx}x^{0.6} = 0.6\,x^{-0.4} \to \infty$ as $x \to 0$ , and Ipopt needs it. Read the message again: can't evaluate pow'(0,0.6) --- note the prime. Work around? Let’s set the lower bound to something really small:

small_number = 1e-6
for r in nlp.REACTORS:

    # Set lower bound
    nlp.reactor_feed[r].setlb(small_number)

    # Adjust initial point if needed
    nlp.reactor_feed[r] = max(nlp.reactor_feed[r](), small_number)

nlp.pprint()

1 Set Declarations
    REACTORS : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    2 : {1, 2}

5 Param Declarations
    conversion : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=False
        Key : Value
          1 :                0.8
          2 : 0.6666666666666666
    feed_cost : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=USD/kmol
        Key  : Value
        None :   5.0
    flowrate_ref : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=kmol/h
        Key  : Value
        None :   1.0
    product_flowrate : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=kmol/h
        Key  : Value
        None :  10.0
    reactor_cost : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=True, Units=USD/h
        Key : Value
          1 :   5.5
          2 :   4.0

3 Var Declarations
    feed_flowrate : Size=1, Index=None, Units=kmol/h
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :  20.0 :  None : False : False : NonNegativeReals
    reactor_effluent : Size=2, Index=REACTORS, Units=kmol/h
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          1 :     0 :  16.0 :  None : False : False : NonNegativeReals
          2 :     0 :   0.0 :  None : False : False : NonNegativeReals
    reactor_feed : Size=2, Index=REACTORS, Units=kmol/h
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          1 : 1e-06 :  20.0 :  None : False : False : NonNegativeReals
          2 : 1e-06 : 1e-06 :  None : False : False : NonNegativeReals

1 Objective Declarations
    cost : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : reactor_cost[1]*(1/flowrate_ref*reactor_feed[1])**0.6 + reactor_cost[2]*(1/flowrate_ref*reactor_feed[2])**0.6 + feed_cost*feed_flowrate

3 Constraint Declarations
    inlet_split : Size=1, Index=None, Active=True
        Key  : Lower : Body                                                : Upper : Active
        None :   0.0 : feed_flowrate - (reactor_feed[1] + reactor_feed[2]) :   0.0 :   True
    mixer : Size=1, Index=None, Active=True
        Key  : Lower            : Body                                      : Upper            : Active
        None : product_flowrate : reactor_effluent[1] + reactor_effluent[2] : product_flowrate :   True
    reactor_performance : Size=2, Index=REACTORS, Active=True
        Key : Lower : Body                                                     : Upper : Active
          1 :   0.0 :                reactor_effluent[1] - 0.8*reactor_feed[1] :   0.0 :   True
          2 :   0.0 : reactor_effluent[2] - 0.6666666666666666*reactor_feed[2] :   0.0 :   True

13 Declarations: REACTORS reactor_cost product_flowrate conversion feed_cost flowrate_ref feed_flowrate reactor_feed reactor_effluent inlet_split reactor_performance mixer cost

Now let’s resolve:

solver = pyo.SolverFactory("ipopt")
results = solver.solve(nlp, tee=True)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

Ipopt 3.14.19: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:        9
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        2

Total number of variables............................:        5
                     variables with only lower bounds:        5
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        4
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.3344037e+02 6.01e+00 8.35e-01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  8.9515896e+01 1.78e-15 1.12e+01  -1.0 7.77e+00    -  4.39e-02 1.00e+00f  1
   2  8.9555393e+01 2.78e-17 5.56e+00  -1.0 2.22e-02    -  9.98e-01 5.00e-01f  2
   3  8.9541585e+01 1.78e-15 2.24e-05  -1.0 3.90e-03    -  1.00e+00 1.00e+00f  1
   4  8.7593746e+01 1.78e-17 5.53e+00  -2.5 5.40e-01    -  9.89e-01 6.15e-01f  1

   5  8.7593068e+01 1.08e-19 1.75e-01  -2.5 1.74e-05   4.0 1.00e+00 1.00e+00f  1
   6  8.7589190e+01 1.78e-15 2.20e-02  -2.5 9.70e-05    -  1.00e+00 1.00e+00f  1
   7  8.7533830e+01 8.16e-20 8.54e+01  -3.8 1.31e-03    -  1.00e+00 6.29e-01f  1
   8  8.7542212e+01 1.78e-15 2.05e+02  -3.8 8.17e-05    -  2.27e-03 5.00e-01f  2
   9  8.7536537e+01 1.69e-21 1.89e+01  -3.8 3.29e-05   5.3 1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  8.7535629e+01 1.78e-15 7.28e+00  -3.8 6.69e-06   5.8 1.00e+00 5.00e-01f  2
  11  8.7536094e+01 1.78e-15 9.18e-01  -3.8 1.64e-06    -  1.00e+00 1.00e+00f  1
  12  8.7536019e+01 0.00e+00 1.80e-02  -3.8 2.76e-07    -  1.00e+00 1.00e+00f  1
  13  8.7536022e+01 1.78e-15 2.74e-05  -3.8 1.04e-08    -  1.00e+00 1.00e+00f  1
  14  8.7533756e+01 1.78e-15 3.24e+01  -5.7 1.03e-05    -  1.00e+00 5.98e-01f  1
  15  8.7533783e+01 1.78e-15 5.52e-02  -5.7 4.47e-08    -  1.00e+00 1.00e+00f  1
  16  8.7533768e+01 1.78e-15 1.74e-02  -5.7 2.57e-08    -  1.00e+00 1.00e+00f  1
  17  8.7533768e+01 1.78e-15 1.13e-07  -5.7 6.42e-11    -  1.00e+00 1.00e+00h  1
  18  8.7533756e+01 1.78e-15 6.33e-02  -8.6 1.93e-08    -  1.00e+00 9.76e-01f  1
  19  8.7533756e+01 1.06e-22 1.79e-08  -8.6 2.50e-11    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  8.7533756e+01 1.78e-15 7.27e-09  -9.0 1.59e-11    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 20

                                   (scaled)                 (unscaled)
Objective...............:   1.4519923360965409e+01    8.7533756344723201e+01
Dual infeasibility......:   7.2654700034036068e-09    4.3800085248214534e-08
Constraint violation....:   1.7763568394002505e-15    1.7763568394002505e-15
Variable bound violation:   9.9909432901401031e-09    9.9909432901401031e-09
Complementarity.........:   9.0911618236617692e-10    5.4806318475633340e-09
Overall NLP error.......:   7.2654700034036068e-09    4.3800085248214534e-08


Number of objective function evaluations             = 28
Number of objective gradient evaluations             = 21
Number of equality constraint evaluations            = 28
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 21
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 20
Total seconds in IPOPT                               = 0.116

EXIT: Optimal Solution Found.

Now we can print the variable names and values:

def print_solution(model):
    """Print variable names and values

    Arguments:
        model: Pyomo model

    """

    print("Variable Names\t\tValue")
    for c in model.component_data_objects(pyo.Var):
        print(c.name, "\t\t", pyo.value(c))

    print("\nObjective Name\t\tValue")
    for c in model.component_data_objects(pyo.Objective):
        print(c.name, "\t\t", pyo.value(c))


print_solution(nlp)

Variable Names		Value
feed_flowrate 		 12.50000016500151
reactor_feed[1] 		 12.499999174992452
reactor_feed[2] 		 9.900090567098599e-07
reactor_effluent[1] 		 9.999999339993963
reactor_effluent[2] 		 6.600060378065732e-07

Objective Name		Value
cost 		 87.53375634472322

Initialize to Favor Reaction 2 and Solve¶

# Initialize
initialize(nlp, reactor_choice=2)

# Correct for bound
# Note: I would have put this in the initialize function but I wanted to show
# the error in class
for r in nlp.REACTORS:
    # Adjust initial point if needed
    nlp.reactor_feed[r] = max(nlp.reactor_feed[r](), small_number)

results = solver.solve(nlp, tee=True)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)
print_solution(nlp)

WARNING (W1002): Setting Var 'reactor_feed[1]' to a numeric value `0` outside
the bounds (1e-06, None).
    See also https://pyomo.readthedocs.io/en/stable/errors.html#w1002

Ipopt 3.14.19: 

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:        9
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        2

Total number of variables............................:        5
                     variables with only lower bounds:        5
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        4
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.2448375e+02 3.34e+00 8.57e-01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  9.7212017e+01 1.78e-15 9.39e+00  -1.0 5.30e+00    -  6.57e-02 1.00e+00f  1
   2  9.9691054e+01 1.78e-15 2.07e+00  -1.0 2.28e+00    -  1.39e-01 1.00e+00f  1
   3  9.9758192e+01 1.78e-15 6.79e-04  -1.0 2.43e-01    -  1.00e+00 1.00e+00f  1
   4  9.9364492e+01 1.78e-15 2.03e-02  -1.7 9.59e-01  -2.0 9.78e-01 1.00e+00f  1
   5  9.5513279e+01 1.78e-15 3.07e+00  -2.5 2.56e+00  -1.6 1.00e+00 7.26e-01f  1
   6  9.5434966e+01 1.78e-15 1.31e+00  -2.5 2.97e-03   2.5 1.00e+00 1.00e+00f  1

   7  9.5393290e+01 1.08e-19 1.79e+01  -2.5 4.95e-02   2.0 1.00e+00 2.31e-02f  2
   8  9.5392828e+01 0.00e+00 1.60e-04  -2.5 1.08e-05    -  1.00e+00 1.00e+00f  1
   9  9.5311712e+01 3.68e-20 8.53e+01  -3.8 1.85e-03    -  1.00e+00 6.16e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  9.5322431e+01 1.78e-15 1.70e+02  -3.8 1.78e-04    -  1.20e-03 2.50e-01f  3

  11  9.5318416e+01 1.78e-15 9.28e+01  -3.8 4.85e-05   5.1 1.00e+00 4.60e-01f  2
  12  9.5315798e+01 1.78e-15 2.73e+02  -3.8 4.07e-04    -  1.00e+00 2.74e-02f  2
  13  9.5315678e+01 1.69e-21 1.44e-02  -3.8 4.39e-07    -  1.00e+00 1.00e+00f  1

  14  9.5314637e+01 0.00e+00 1.31e+00  -3.8 3.55e-06    -  1.00e+00 1.00e+00f  1
  15  9.5314724e+01 1.78e-15 1.37e-02  -3.8 2.76e-07    -  1.00e+00 1.00e+00f  1
  16  9.5314720e+01 1.78e-15 2.30e-05  -3.8 1.16e-08    -  1.00e+00 1.00e+00h  1

  17  9.5311597e+01 1.78e-15 3.25e+01  -5.7 1.25e-05    -  1.00e+00 5.96e-01f  1
  18  9.5311634e+01 1.78e-15 5.58e-02  -5.7 5.39e-08    -  1.00e+00 1.00e+00f  1
  19  9.5311613e+01 1.06e-22 1.77e-02  -5.7 3.11e-08    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  9.5311612e+01 1.06e-22 1.14e-07  -5.7 7.73e-11    -  1.00e+00 1.00e+00h  1

  21  9.5311597e+01 1.78e-15 5.31e-02  -8.6 2.32e-08    -  1.00e+00 9.76e-01f  1
  22  9.5311597e+01 1.78e-15 1.80e-08  -8.6 3.00e-11    -  1.00e+00 1.00e+00f  1
  23  9.5311597e+01 1.78e-15 7.30e-09  -9.0 1.91e-11    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 23

                                   (scaled)                 (unscaled)
Objective...............:   1.1498251558692399e+01    9.5311596851998971e+01
Dual infeasibility......:   7.2961370278790127e-09    6.0479323089120971e-08
Constraint violation....:   1.7763568394002505e-15    1.7763568394002505e-15
Variable bound violation:   9.9909242686312207e-09    9.9909242686312207e-09
Complementarity.........:   9.0911629856797057e-10    7.5358697536224235e-09
Overall NLP error.......:   7.2961370278790127e-09    6.0479323089120971e-08


Number of objective function evaluations             = 29
Number of objective gradient evaluations             = 24
Number of equality constraint evaluations            = 29
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 24
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 23
Total seconds in IPOPT                               = 0.261

EXIT: Optimal Solution Found.

Variable Names		Value
feed_flowrate 		 14.999999801998188
reactor_feed[1] 		 9.900090757313687e-07
reactor_feed[2] 		 14.99999881198911
reactor_effluent[1] 		 7.92007260585095e-07
reactor_effluent[2] 		 9.99999920799274

Objective Name		Value
cost 		 95.31159685199898

Compare Solutions¶

	Solution 1	Solution 2
Feed Flowrate (kmol/hr)	12.50	15.00
Reactor 1 Feed (kmol/hr)	12.50	0
Reactor 2 Feed (kmol/hr)	0	15.00
Reactor 1 Effluent (kmol/hr)	10.00	0
Reactor 2 Effluent (kmol/hr)	0	10.00
Cost ($/hr)	87.53	95.31

Which solution is better? Why are there multiple solutions?

Why there are two solutions, as one picture¶

Cost along the mass balance as a function of the feed to reactor II: a concave curve on a segment, with local minima at both endpoints and a maximum in the interior.

Eliminating $x_1$ through $0.8 x_1 + (2/3) x_2 = 10$ leaves the cost a function of $x_2$ alone on the segment $0 \le x_2 \le 15$ . Because each $x^{0.6}$ term is concave, the cost is concave, so its minima sit at the two endpoints — use reactor I only, or use reactor II only — and its maximum is in the interior. Both endpoints satisfy the KKT conditions, so a local solver returns whichever one the initial guess was nearer. That is precisely what the two runs above did.

Rendered from figures/plots/reactor-concave-local-minima.py, and the same image printed in the course pack. It re-derives the curve with numpy rather than Pyomo, so regenerating it needs no solver binary. It recomputes Biegler, Grossmann & Westerberg (1997) Figure 15.12, p. 510. With $\theta_2 = 2/3$ the right-hand endpoint is $x_2 = 15$ exactly and $C = \$95.31$ /hr, matching their printed 95.3; the interior maximum recomputes to $x_2 = 11.26$ against their hand-drawn 11.4.

Linearized Discrete Cost Model¶

We want to modify the model to:

Easily find the best global solution (reduce impacts of initialization)
Account for the fact there is a minimum reactor size we can purchase

The linearized cost model, as one picture¶

Fixed-charge linearization against the concave power law, both reactors: left, vessel cost only; right, with the $5/kmol feed added. The linearized model is discontinuous at zero feed.

Rendered from figures/plots/reactor-cost-linearization.py, and the same image printed in the course pack, so the handout and the website cannot drift apart. The left panel is Biegler, Grossmann & Westerberg, Systematic Methods of Chemical Process Design (1997), Figure 15.13, p. 511; the right panel is the cell below.

Note the jump at $x_r = 0$ : the fixed charge is paid only if the unit is built. That discontinuity is exactly what a binary variable is for.

Run the cell below to reproduce the right panel and change the constants yourself. (It plots only the total cost, so the two models nearly coincide — the feed term is linear in both.)

import matplotlib.pyplot as plt
import numpy as np

feed = np.linspace(0, 20, 101)

# Price of feed A, c_f [USD/kmol] -- the same for both reactors
c_feed = 5

for reactor in [1, 2]:

    if reactor == 1:
        c = 5.5
        c_linear = 6.4
        c_fixed = 7.5
        color = "r"
    elif reactor == 2:
        c = 4.0
        c_linear = 6.0
        c_fixed = 5.5
        color = "b"
    else:
        break

    cost_continuous = c * (feed) ** (0.6) + c_feed * feed
    cost_linearized = c_linear * feed + c_fixed

    plt.plot(
        feed,
        cost_continuous,
        label="Reactor " + str(reactor) + " (continuous)",
        color=color,
        linestyle="-",
    )
    plt.plot(
        feed,
        cost_linearized,
        label="Reactor " + str(reactor) + " (linearized)",
        color=color,
        linestyle="--",
    )

plt.legend()
plt.xlabel("Feed Flowrate (kmol/hr)")
plt.ylabel(r"Cost (\$/hr)")
plt.show()

Define binary variable:

$y_r = \begin{cases} 1 & \text{if reactor } r \text{ is selected} \\ 0 & \text{otherwise} \end{cases}$

Add logical constraints to include the maximum flowrate (20 kmol/hr):

$x_r \leq 20 y_r, \quad \forall r \in \mathcal{R}$

Thus if $y_r = 0$ , then $x_r = 0$

Finally, we can define the new linearized reactor cost model:

\begin{align*} \text{Reactor I:} \quad & 7.5 y_1 + 1.4 x_1 \\ \text{Reactor II:} \quad & 5.5 y_2 + 1.0 x_2 \end{align*}

(2)

Recall that $x_0 = x_1 + x_2$ . Thus, the feed cost is $5x_0 = 5 x_1 + 5 x_2$ .

Putting this all together gives the following MILP optimization problem:

\begin{align*} \min_{x,y} \quad & C = 7.5 y_1 + 6.4 x_1 + 5.5 y_2 + 6.0 x_2 \\ \text{s.t.} \quad & 0.8 x_1 + (2/3) x_2 = 10 \\ & x_1 - 20 y_1 \leq 0 \\ & x_2 - 20 y_2 \leq 0 \\ & x_1, x_2 \geq 0, \quad y_1, y_2 \in \{0, 1\} \end{align*}

(3)

Enumerate the solutions¶

As an illustration, enumerate through the following four options:

No reactor
Reactor I only
Reactor II only
Reactor I and II only

For each option, ask:

Are the constraints feasible?
What is the objective?

Reactor 1	Reactor 2	Cost
0	0	infeasible
1	0	87.5
0	1	95.5
1	1	93.0

Solve with Pyomo¶

Create and inspect the model.

milp = pyo.ConcreteModel()

## Define sets
milp.REACTORS = pyo.Set(initialize=range(1, 3))

## Define parameters (data)

# Big-M throughput limit, M [kmol/hr]
milp.max_flowrate = pyo.Param(initialize=20.0, units=u.kmol / u.hr)

# Marginal (linear) cost of reactor r, abar_r [USD/kmol]
cost_coefficient1 = {1: 1.4, 2: 1.0}
milp.reactor_cost_linear = pyo.Param(
    milp.REACTORS, initialize=cost_coefficient1, units=u.USD / u.kmol
)

# Fixed charge for selecting reactor r, bbar_r [USD/hr]
cost_coefficient2 = {1: 7.5, 2: 5.5}
milp.reactor_cost_fixed = pyo.Param(
    milp.REACTORS, initialize=cost_coefficient2, units=u.USD / u.hr
)

# Required product rate, P [kmol B/hr]
milp.product_flowrate = pyo.Param(initialize=10.0, units=u.kmol / u.hr)

# Conversion of reactor r, theta_r [kmol B/kmol A]
reactor_conversion = {1: 0.8, 2: 2 / 3}
milp.conversion = pyo.Param(milp.REACTORS, initialize=reactor_conversion)

# Price of feed A, c_f [USD/kmol]
milp.feed_cost = pyo.Param(initialize=5.0, units=u.USD / u.kmol)


## Define variables

# Total feed of A, x_0 in the notes [kmol/hr]
milp.feed_flowrate = pyo.Var(
    domain=pyo.NonNegativeReals, bounds=(0, milp.max_flowrate), units=u.kmol / u.hr
)

# Feed into reactor r, x_r in the notes [kmol/hr]
milp.reactor_feed = pyo.Var(
    milp.REACTORS,
    domain=pyo.NonNegativeReals,
    bounds=(0, milp.max_flowrate),
    units=u.kmol / u.hr,
)

# Effluent leaving reactor r, z_r in the notes [kmol B/hr]
milp.reactor_effluent = pyo.Var(
    milp.REACTORS, domain=pyo.NonNegativeReals, units=u.kmol / u.hr
)

# Binary variables: y_r = 1 if reactor r is selected
# Add your solution here

## Define constraints


# Mass balance over the splitter [kmol/hr]
@milp.Constraint()
def inlet_split(b):
    return b.feed_flowrate == sum(b.reactor_feed[r] for r in b.REACTORS)


# Reactor conversion, z_r = theta_r * x_r [kmol B/hr]
@milp.Constraint(milp.REACTORS)
def reactor_performance(b, r):
    return b.reactor_effluent[r] == b.conversion[r] * b.reactor_feed[r]


# Mass balance over the mixer meets product requirements [kmol B/hr]
@milp.Constraint()
def mixer(b):
    return b.product_flowrate == sum(b.reactor_effluent[r] for r in b.REACTORS)


# Big-M logical constraint, x_r <= M * y_r [kmol/hr]
# Add your solution here

## Define objective


# BGW fold the feed price into each flow coefficient, abar_r + c_f = 6.4
# and 6.0; here the two terms are kept separate.
# Total cost [USD/hr]
# Add your solution here

# Inspect the model
milp.pprint()

# Raises UnitsError if any constraint or the objective is dimensionally inconsistent
assert_units_consistent(milp)

print("Units are consistent.")

1 Set Declarations
    REACTORS : Size=1, Index=None, Ordered=Insertion
        Key  : Dimen : Domain : Size : Members
        None :     1 :    Any :    2 : {1, 2}

6 Param Declarations
    conversion : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=False
        Key : Value
          1 :                0.8
          2 : 0.6666666666666666
    feed_cost : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=USD/kmol
        Key  : Value
        None :   5.0
    max_flowrate : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=kmol/h
        Key  : Value
        None :  20.0
    product_flowrate : Size=1, Index=None, Domain=Any, Default=None, Mutable=True, Units=kmol/h
        Key  : Value
        None :  10.0
    reactor_cost_fixed : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=True, Units=USD/h
        Key : Value
          1 :   7.5
          2 :   5.5
    reactor_cost_linear : Size=2, Index=REACTORS, Domain=Any, Default=None, Mutable=True, Units=USD/kmol
        Key : Value
          1 :   1.4
          2 :   1.0

4 Var Declarations
    feed_flowrate : Size=1, Index=None, Units=kmol/h
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :  None :  20.0 : False :  True : NonNegativeReals
    reactor_boolean : Size=2, Index=REACTORS
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          1 :     0 :  None :     1 : False :  True : Binary
          2 :     0 :  None :     1 : False :  True : Binary
    reactor_effluent : Size=2, Index=REACTORS, Units=kmol/h
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          1 :     0 :  None :  None : False :  True : NonNegativeReals
          2 :     0 :  None :  None : False :  True : NonNegativeReals
    reactor_feed : Size=2, Index=REACTORS, Units=kmol/h
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          1 :     0 :  None :  20.0 : False :  True : NonNegativeReals
          2 :     0 :  None :  20.0 : False :  True : NonNegativeReals

1 Objective Declarations
    cost : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize : reactor_cost_linear[1]*reactor_feed[1] + reactor_cost_fixed[1]*reactor_boolean[1] + reactor_cost_linear[2]*reactor_feed[2] + reactor_cost_fixed[2]*reactor_boolean[2] + feed_cost*feed_flowrate

4 Constraint Declarations
    inlet_split : Size=1, Index=None, Active=True
        Key  : Lower : Body                                                : Upper : Active
        None :   0.0 : feed_flowrate - (reactor_feed[1] + reactor_feed[2]) :   0.0 :   True
    mixer : Size=1, Index=None, Active=True
        Key  : Lower            : Body                                      : Upper            : Active
        None : product_flowrate : reactor_effluent[1] + reactor_effluent[2] : product_flowrate :   True
    reactor_performance : Size=2, Index=REACTORS, Active=True
        Key : Lower : Body                                                     : Upper : Active
          1 :   0.0 :                reactor_effluent[1] - 0.8*reactor_feed[1] :   0.0 :   True
          2 :   0.0 : reactor_effluent[2] - 0.6666666666666666*reactor_feed[2] :   0.0 :   True
    toggle_reactor : Size=2, Index=REACTORS, Active=True
        Key : Lower : Body                                              : Upper : Active
          1 :  -Inf : reactor_feed[1] - max_flowrate*reactor_boolean[1] :   0.0 :   True
          2 :  -Inf : reactor_feed[2] - max_flowrate*reactor_boolean[2] :   0.0 :   True

16 Declarations: REACTORS max_flowrate reactor_cost_linear reactor_cost_fixed product_flowrate conversion feed_cost feed_flowrate reactor_feed reactor_effluent reactor_boolean inlet_split reactor_performance mixer toggle_reactor cost
Units are consistent.

# Binary variables: y_r = 1 if reactor r is selected
milp.reactor_boolean = pyo.Var(milp.REACTORS, domain=pyo.Binary)


# Big-M logical constraint, x_r <= M * y_r [kmol/hr]
@milp.Constraint(milp.REACTORS)
def toggle_reactor(b, r):
    return b.reactor_feed[r] <= b.reactor_boolean[r] * b.max_flowrate


# Total cost [USD/hr]
@milp.Objective()
def cost(b):
    return (
        sum(
            b.reactor_cost_linear[r] * b.reactor_feed[r]
            + b.reactor_cost_fixed[r] * b.reactor_boolean[r]
            for r in b.REACTORS
        )
        + b.feed_cost * b.feed_flowrate
    )

Solve the model using appsi_highs (the course default for LPs and MILPs), bonmin, gurobi (need license), or cplex (need license).

# Set solver
# solver = pyo.SolverFactory('gurobi')
solver = pyo.SolverFactory("appsi_highs")
# solver = pyo.SolverFactory('bonmin')

# Add your solution here

print_solution(milp)

Is rounding good enough?¶

Linear Program (Relaxation)¶

\begin{align}\min_{x_1,x_2} \quad & x_2 \\ \mathrm{s.t.} \quad & 2 x_1 + x_2 \geq 13 \\ & 5 x_1 + 2 x_2 \leq 30 \\ & -x_1 + x_2 \geq 5 \\ & x_1, x_2 \geq 0 \end{align}

(4)

import pyomo.environ as pyo

m = pyo.ConcreteModel()

# Declare variables with bounds
m.x1 = pyo.Var(domain=pyo.NonNegativeReals)
m.x2 = pyo.Var(domain=pyo.NonNegativeReals)

# Constraint 1
m.con1 = pyo.Constraint(expr=2 * m.x1 + m.x2 >= 13)

# Constraint 2
m.con2 = pyo.Constraint(expr=5 * m.x1 + 2 * m.x2 <= 30)

# Constraint 3
m.con3 = pyo.Constraint(expr=-m.x1 + m.x2 >= 5)

# Objective
m.obj = pyo.Objective(expr=m.x2)

# Print model
m.pprint()

2 Var Declarations
    x1 : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :  None :  None : False :  True : NonNegativeReals
    x2 : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :  None :  None : False :  True : NonNegativeReals

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize :         x2

3 Constraint Declarations
    con1 : Size=1, Index=None, Active=True
        Key  : Lower : Body      : Upper : Active
        None :  13.0 : 2*x1 + x2 :  +Inf :   True
    con2 : Size=1, Index=None, Active=True
        Key  : Lower : Body        : Upper : Active
        None :  -Inf : 5*x1 + 2*x2 :  30.0 :   True
    con3 : Size=1, Index=None, Active=True
        Key  : Lower : Body      : Upper : Active
        None :   5.0 : - x1 + x2 :  +Inf :   True

6 Declarations: x1 x2 con1 con2 con3 obj

# Set solver
# solver = pyo.SolverFactory('gurobi')
# solver = pyo.SolverFactory('ipopt')
solver = pyo.SolverFactory("appsi_highs")


# Solve
results = solver.solve(m, tee=True)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

# Print solution
print(" ")
print("x1 = ", pyo.value(m.x1))
print("x2 = ", pyo.value(m.x2))

 
x1 =  2.6666666666666665
x2 =  7.666666666666667

Rounding¶

With your neighbor, discuss:

If you round $x_1$ and $x_2$ to an integer, are all of the constraints feasible?
How would you go about checking the optimality of a feasible integer solution?

Integer Program¶

Consider the following integer program:

\begin{align}\min_{x_1,x_2} \quad & x_2 \\ \mathrm{s.t.} \quad & 2 x_1 + x_2 \geq 13 \\ & 5 x_1 + 2 x_2 \leq 30 \\ & -x_1 + x_2 \geq 5 \\ & x_1, x_2 \in \mathcal{Z} := \{0,1,2,...\} \end{align}

(5)

m2 = pyo.ConcreteModel()

# Declare variables as non-negative integers
m2.x1 = pyo.Var(domain=pyo.NonNegativeIntegers)
m2.x2 = pyo.Var(domain=pyo.NonNegativeIntegers)

# Constraint 1
m2.con1 = pyo.Constraint(expr=2 * m2.x1 + m2.x2 >= 13)

# Constraint 2
m2.con2 = pyo.Constraint(expr=5 * m2.x1 + 2 * m2.x2 <= 30)

# Constraint 3
m2.con3 = pyo.Constraint(expr=-m2.x1 + m2.x2 >= 5)

# Objective
m2.obj = pyo.Objective(expr=m2.x2)

m2.pprint()

2 Var Declarations
    x1 : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :  None :  None : False :  True : NonNegativeIntegers
    x2 : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :  None :  None : False :  True : NonNegativeIntegers

1 Objective Declarations
    obj : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : minimize :         x2

3 Constraint Declarations
    con1 : Size=1, Index=None, Active=True
        Key  : Lower : Body      : Upper : Active
        None :  13.0 : 2*x1 + x2 :  +Inf :   True
    con2 : Size=1, Index=None, Active=True
        Key  : Lower : Body        : Upper : Active
        None :  -Inf : 5*x1 + 2*x2 :  30.0 :   True
    con3 : Size=1, Index=None, Active=True
        Key  : Lower : Body      : Upper : Active
        None :   5.0 : - x1 + x2 :  +Inf :   True

6 Declarations: x1 x2 con1 con2 con3 obj

# Set solver
# solver = pyo.SolverFactory('gurobi')
solver = pyo.SolverFactory("appsi_highs")
# solver = pyo.SolverFactory('bonmin')

# Solve
results = solver.solve(m2, tee=True)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

# Print solution
print(" ")
print("x1 = ", m2.x1())
print("x2 = ", m2.x2())

 
x1 =  2.0
x2 =  9.0

Why rounding does not always work¶

Reference: http://personal.lse.ac.uk/williahp/publications/wp 10-118.pdf