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

Superstructures and discrete choices¶

A superstructure represents the candidate units, streams and configurations considered in a design problem. Discrete variables select units or connections; continuous variables size and operate them. A structure omitted from this search space cannot be selected by optimization. Including an option does not guarantee its physical feasibility: balances, performance and specifications still have to be imposed.

BGW (1997), Example 15.2, pp. 500–502, considers three ammonia-process choices: tubular versus multibed-quench reactor, flash condensation versus absorption/distillation for product separation, and membrane versus no hydrogen recovery. Choosing exactly one option at each of three independent stages gives $2^3=8$ combinations. A tree lists those combinations; a network superstructure shares units and can represent recycle and simultaneous routes. BGW’s network may represent combinations beyond the eight tree leaves if multiple paths are allowed. Exactly-one assumptions must therefore be written as constraints, rather than inferred from the drawing.

A general mixed-integer nonlinear program is

\min_{\mathbf{x},\mathbf{y}} f(\mathbf{x},\mathbf{y})\quad\text{s.t. }\mathbf{h}(\mathbf{x},\mathbf{y})=\mathbf{0},\quad \mathbf{g}(\mathbf{x},\mathbf{y})\le\mathbf{0},\quad \mathbf{x}\in X\subseteq\mathbb R^n,\quad \mathbf{y}\in Y\subseteq\mathbb Z^t.

(1)

Here $\mathbf{x}$ is continuous and $\mathbf{y}$ is integer. Binary variables restrict $Y$ to a subset of $\{0,1\}^t$ . Both variable types can enter every function, including the equality constraints. In the linear special case the model is a MILP; if all its decision variables are integer it is an IP. An integer domain alone does not make nonlinear functions linear. Biegler (2010), Eq. (1.1), p. 3 and p. 4, gives the binary form and the class distinctions.

Three ammonia-process technology choices.

The reactor problem below is a compact superstructure with two parallel candidates. The two-reactor power-law NLP permits zero flows instead of explicit build switches; fixed charges will make selection a separate binary decision.

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 feed flowrate into reactor $r$ . Its cost contribution, including feed, is

C_r=a_r(x_r/x_{\mathrm{ref}})^{0.6}+c_f x_r,\quad x_{\mathrm{ref}}=1~\mathrm{kmol/hr},\quad c_f=5~\mathrm{\$/kmol}.

(2)

The reactor cost coefficients are rates in dollars per hour; the reference scale makes the power-law base dimensionless.

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

Count declared components before solving¶

With $n=|\mathcal R|$ reactors, the expanded model declares $1+2n$ continuous variables and $n+2$ linear equalities: one split, $n$ conversions and one mixer. The original NLP has $1+2n$ lower-bounded variables and no upper bounds. Its objective is nonlinear even though its constraints are linear.

These equality rows are independent when every $\theta_r>0$ . In a linear combination of rows, the coefficient of $x_0$ first forces the splitter multiplier to zero; the coefficient of each $x_r$ then forces its conversion multiplier to zero, and the $z_r$ coefficients force the mixer multiplier to zero. Thus the affine equality solution space has dimension $(1+2n)-(n+2)=n-1$ . For $P>0$ it contains strictly positive reactor flows, so nonnegativity leaves a feasible set of that relative dimension. At its endpoints the bounds restrict the feasible directions. Counting rows alone would not establish rank or feasibility in a different model.

For two reactors this gives 5 variables, 4 independent equalities and one independent split choice.

The MILP below retains those 5 continuous variables and 4 equalities, and adds 2 binaries and 2 big- $M$ inequalities. Its code also declares upper bounds of 20 on $x_0,x_1,x_2$ : 2 lower-only continuous variables and 5 variables bounded on both sides, including the binaries. The formulation without those explicit upper bounds has the same feasible set here: positive conversions and the production equality already imply $x_1\le12.5$ , $x_2\le15$ and $x_0\le15$ . A declared bound is still counted even when redundant.

The NLP demonstration later replaces the two reactor-feed lower bounds by 10^-6 to avoid singular derivatives. This slightly changes its feasible set; its output approximates the zero-flow endpoints.

Eliminate the expanded reactor balances using one independent flow.

Solve with the power-law cost model in Pyomo¶

The complete model uses the same sets, data and reference flow as the equations above.

Units in the power-law cost¶

The coefficients $a_r=5.5,4.0$ are specified in dollars per hour when the numerical flow is measured in kmol/hr. With those coefficient units, $a_r x_r^{0.6}$ does not have units of dollars per hour. A dimensional power is possible if its coefficient carries the corresponding compensating units; that is not how these coefficients are declared.

Write the implicit reference scale explicitly:

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

(3)

Now the exponent acts on a dimensionless ratio and $C_r$ has units of dollars per hour. The numerical objective is unchanged in the chosen units. If flow units change, convert $x_{\mathrm{ref}}$ too; changing only the numerical unit changes the fitted correlation. assert_units_consistent checks dimensional consistency, not physical adequacy, feasibility or optimality.

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 reactor I 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

Solve the current model and inspect its termination condition.

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.13.2: 

******************************************************************************
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 http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

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.

Expected solver failure: 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.13.2: 

******************************************************************************
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 http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

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 0.00e+00 5.56e+00  -1.0 2.22e-02    -  9.98e-01 5.00e-01f  2
   3  8.9541585e+01 2.78e-17 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 0.00e+00 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 3.39e-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
Complementarity.........:   9.0911618236617630e-10    5.4806318475633307e-09
Overall NLP error.......:   7.2654700034036068e-09    4.3800085248214534e-08


Number of objective function evaluations             = 30
Number of objective gradient evaluations             = 21
Number of equality constraint evaluations            = 30
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 CPU secs in IPOPT (w/o function evaluations)   =      0.002
Total CPU secs in NLP function evaluations           =      0.000

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] 		 1e-06
reactor_effluent[1] 		 9.999999339993963
reactor_effluent[2] 		 6.600060378065732e-07

Objective Name		Value
cost 		 87.53376237988246

# Full returned variable vector, reported compactly for the handout.
print(f"x0={pyo.value(nlp.feed_flowrate):.10f}, "
      f"x1={pyo.value(nlp.reactor_feed[1]):.10f}, "
      f"x2={pyo.value(nlp.reactor_feed[2]):.3e}")
print(f"z1={pyo.value(nlp.reactor_effluent[1]):.10f}, "
      f"z2={pyo.value(nlp.reactor_effluent[2]):.3e}; "
      f"C={pyo.value(nlp.cost):.8f} USD/hr")

x0=12.5000001650, x1=12.4999991750, x2=1.000e-06
z1=9.9999993400, z2=6.600e-07; C=87.53376238 USD/hr

Initialize to favor reactor II 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.13.2: 

******************************************************************************
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 http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

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 2.22e-16 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 2.39e-17 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 2.17e-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.70e-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.69e-20 9.28e+01  -3.8 4.85e-05   5.1 1.00e+00 4.60e-01f  2
  12  9.5315798e+01 1.86e-20 2.73e+02  -3.8 4.07e-04    -  1.00e+00 2.74e-02f  2
  13  9.5315678e+01 0.00e+00 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
Complementarity.........:   9.0911629856797057e-10    7.5358697536224235e-09
Overall NLP error.......:   7.2961370278790127e-09    6.0479323089120971e-08


Number of objective function evaluations             = 32
Number of objective gradient evaluations             = 24
Number of equality constraint evaluations            = 32
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 CPU secs in IPOPT (w/o function evaluations)   =      0.002
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.

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

Objective Name		Value
cost 		 95.31160515032713

Compare the two returned 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

Interpret the two returned solutions using the complete cost profile below.

The reactor NLP has two local minima¶

Cost along the production equality, with two endpoint minima and an interior maximum.

Eliminating $z_1,z_2,x_0,x_1$ gives

x_1(t)=\frac{10-(2/3)t}{0.8},\qquad 0\le t=x_2\le15,

(4)

C(t)=5.5\,x_1(t)^{0.6}+4t^{0.6}+5[x_1(t)+t],

(5)

where numerical flows use kmol/hr and the reference flow is 1 kmol/hr. On the open segment, each positive power-law term is strictly concave; the feed term is affine. Thus $C$ is strictly concave and continuous on the closed segment. For $0\le\lambda\le1$ and $t=15\lambda$ , concavity gives

C(t)\ge(1-\lambda)C(0)+\lambda C(15)\ge\min\{C(0),C(15)\}.

(6)

Comparing endpoints therefore proves the global optimum of this particular model: reactor I only, $x_1=12.5,x_2=0$ , costing about 87.53 dollars/hr. Reactor II alone costs about 95.31 dollars/hr. The interior maximum is at $x_2\approx11.26$ ; BGW’s sketch labels it 11.4. These values use the exact conversion $2/3$ , rather than the rounded 0.67 in the printed equations.

Concavity alone does not make both endpoints local minima. Here the incoming unused-reactor cost has a positive $h^{0.6}$ term that dominates the $O(h)$ change in the active reactor and feed. Near $t=0$ , $C(t)-C(0)=4t^{0.6}+O(t)>0$ ; near $t=15$ , $C(15-h)-C(15)=5.5(5h/6)^{0.6}+O(h)>0$ . Hence both endpoints are strict local minima relative to the feasible segment.

The derivative of $x^{0.6}$ is singular at zero. The usual differentiable KKT theorem cannot be applied at these exact endpoints. The notebook instead solves the slightly perturbed model with reactor-feed lower bounds 10^-6. Its two initializations demonstrate two local returns; there is no general rule that a local solver selects the nearest endpoint, and an optimal termination message alone does not prove the global result. The endpoint comparison above supplies that proof.

The analytic curve is owned by figures/plots/reactor-concave-local-minima.py; it is not a second implementation of solver output. Source: BGW (1997), Example 15.3, pp. 509–511.

Linearized Discrete Cost Model¶

Replace the power-law equipment cost by a fixed charge plus a constant marginal cost. This gives a different cost approximation whose resulting model is a MILP. A globally solved MILP certifies that approximation; it does not certify an arbitrary original nonlinear cost model. The on/off link below imposes an upper capacity when a reactor is selected. It does not impose a positive minimum operating flow.

Fixed charges and what the binary means¶

Power-law and fixed-charge equipment costs, with and without feed cost; the inset shows the jump at zero.

For each reactor introduce $y_r\in\{0,1\}$ and impose $0\le x_r\le M_r y_r$ . Its equipment cost is $\bar b_r y_r+\bar a_r x_r$ , with $(\bar b_1,\bar b_2)=(7.5,5.5)$ dollars/hr and $(\bar a_1,\bar a_2)=(1.4,1.0)$ dollars/kmol. Add the feed cost $5x_0$ . A one-time capital expense must be annualized or otherwise put on the same time basis before combining it with a rate in dollars/hr.

The link proves $y_r=0\Rightarrow x_r=0$ , or equivalently $x_r>0\Rightarrow y_r=1$ . It permits $y_r=1,x_r=0$ : a selected unit can be idle. With a positive fixed charge and no other reason to retain an idle unit, minimization switches it off. A positive minimum throughput would require the additional row $L_r y_r\le x_r$ for a justified $L_r>0$ .

Allowed nonnegative flows when a reactor is off or selected.

After minimizing over a free switch for one unit, its cost as a function of flow is piecewise:

\widehat C_r(x_r)=\begin{cases}0,&x_r=0,\\\bar b_r+\bar a_r x_r,&0<x_r\le M_r.\end{cases}

(7)

This is a discontinuous function of $x_r$ alone. The binary gives it a linear mixed-integer representation; it is not correct to say no such scalar function exists. The plotted zero-cost point represents an unbuilt unit, and the positive-flow branch represents a selected unit.

The shared analytic figure contains equipment-only and equipment-plus-feed panels and an origin inset. The separate exploratory cell below lets you change the two total-cost lines. It plots the selected-unit branch at zero as well; the shared figure distinguishes that branch from the unbuilt zero-cost point. Both implementations serve these stated, different roles.

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 valid upper bound (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*}

(8)

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_{\mathbf{x},\mathbf{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*}

(9)

Enumerate the four binary assignments¶

For each fixed $(y_1,y_2)$ , solve the remaining LP. The same production target, conversions and costs apply in every row.

$y_1$	$y_2$	$x_1$	$x_2$	Cost ($/hr)	Interpretation
0	0	0	0	infeasible	$0.8x_1+(2/3)x_2=0\ne10$
1	0	12.5	0	87.5	reactor I alone
0	1	0	15	95.5	reactor II alone
1	1	12.5	0	93.0	reactor II is selected but idle and still costs 5.5

The marginal costs per unit of product, including feed, are $6.4/0.8=8$ for reactor I and $6/(2/3)=9$ for reactor II. Thus, when both are selected, all production still uses I. The cheapest feasible row proves the global optimum of this two-reactor MILP. Its agreement in reactor choice with the nonlinear model is specific to these data; their objective values differ.

With $t$ binary choices, exhaustive enumeration has $2^t$ assignments: 1,024 at $t=10$ and 1,073,741,824 at $t=30$ . General integer variables need their own finite bounds before exhaustive enumeration is finite. Practical methods discard groups using valid bounds rather than solving every assignment. The complete branch-and-bound algorithm is taught in the course’s MINLP Algorithms guide; this section supplies the modeling example. Source: BGW (1997), pp. 512–513.

Solve with Pyomo¶

Create and inspect the model. The explicit upper bounds on feed variables are redundant for these data, but remain part of its declared component count.

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)

# On/off link x_r <= M y_r [kmol/hr].
# Add your solution here

# Equipment and feed 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)

# A compact report of the actual solved MILP for the handout.
print(f"MILP flows: x0={pyo.value(milp.feed_flowrate):.4f}, "
      f"x1={pyo.value(milp.reactor_feed[1]):.4f}, "
      f"x2={pyo.value(milp.reactor_feed[2]):.4f} kmol/hr")
print(f"Switches: y1={pyo.value(milp.reactor_boolean[1]):.0f}, "
      f"y2={pyo.value(milp.reactor_boolean[2]):.0f}; "
      f"C={pyo.value(milp.cost):.4f} USD/hr")

MILP flows: x0=12.5000, x1=12.5000, x2=0.0000 kmol/hr
Switches: y1=1, y2=0; C=87.5000 USD/hr

Valid bounds and the LP relaxation¶

For a nonnegative flow, $x_r\le M_r y_r$ is valid only when $M_r$ is at least the largest feasible value of $x_r$ in the intended physical model. With $\theta_r>0$ , $x\ge0$ and $\sum_r\theta_r x_r=P$ , every term is at most $P$ , so

x_r\le P/\theta_r;\qquad M_1=12.5,\quad M_2=15

(10)

are valid individual bounds here. The common value 20 is valid but weaker. A bound of 10 excludes the feasible reactor-I-only design. Bounds must be re-derived if conversion, production requirements, additional reactions or flow signs change.

Relaxing a binary switch gives $0\le y_r\le1$ . At $x_1=12.5$ , the link requires $y_1\ge12.5/M_1$ . With common $M=20$ , the smallest switch is 0.625; with 20,000 it is 0.000625. The physical production equality already bounds both flows, so increasing this valid logical bound does not add new integer designs in this particular example. It does allow the LP to pay a smaller fraction of the fixed charge. The scalar switch picture below omits the reactor production equality and explicit physical caps so its mechanism is visible.

The on/off link and its LP relaxation for two choices of a valid upper bound.

For fixed $M_r$ , the single-unit set $\{(0,0)\}\cup\{(x,1):0\le x\le M_r\}$ has convex hull $0\le y\le1$ , $0\le x\le M_r y$ . If $y>0$ , any relaxed point equals $(1-y)(0,0)+y(x/y,1)$ ; the $y=0$ case is $(0,0)$ . This proof concerns one switch. Intersecting many such hulls with balances does not generally give the convex hull of the whole process.

For minimization the relaxed feasible set contains every integer-feasible point, so its globally optimal value is a lower bound on the MILP optimum, provided the same objective is used. The table computed below solves the complete reactor model, not just the scalar picture. The common bound of 20 gives 84.6875 dollars/hr; 20,000 gives 80.0046875; the individual production bounds give 87.5 and close the objective gap for these data. This is not a general integrality guarantee.

Source: BGW (1997), Eqs. (15.10)–(15.12), p. 512, and Appendix A, §A.3.2, pp. 758, 760–761. The production-bound and convex-combination arguments are supplied here in full.

# Count declared rows and bounds, rather than subtracting rows from variables.
def declared_counts(model):
    variables = list(model.component_data_objects(pyo.Var))
    rows = list(model.component_data_objects(pyo.Constraint, active=True))
    return dict(
        continuous=sum(v.is_continuous() for v in variables),
        binary=sum(v.is_binary() for v in variables),
        linear_equalities=sum(c.equality and c.body.polynomial_degree() in (0, 1) for c in rows),
        linear_inequalities=sum(not c.equality and c.body.polynomial_degree() in (0, 1) for c in rows),
        lower_only=sum(v.lb is not None and v.ub is None for v in variables),
        both_bounds=sum(v.lb is not None and v.ub is not None for v in variables),
    )

print("NLP declarations:", declared_counts(nlp))
print("MILP declarations:", declared_counts(milp))
highs = pyo.SolverFactory("appsi_highs")
print("Fixed binaries | x1 | x2 | cost [USD/hr]")
for first, second in [(0, 0), (1, 0), (0, 1), (1, 1)]:
    fixed = milp.clone()
    fixed.reactor_boolean[1].fix(first)
    fixed.reactor_boolean[2].fix(second)
    result = highs.solve(fixed, load_solutions=False)
    if result.solver.termination_condition == pyo.TerminationCondition.infeasible:
        print(f"{first}{second}: infeasible")
    else:
        assert pyo.check_optimal_termination(result)
        fixed.solutions.load_from(result)
        print(f"{first}{second}: {pyo.value(fixed.reactor_feed[1]):.4f} | "
              f"{pyo.value(fixed.reactor_feed[2]):.4f} | {pyo.value(fixed.cost):.4f}")

print("LP bound | y1 | x1 | cost [USD/hr]")
for label, bound in [("M=20", 20.0), ("M=20000", 20000.0), ("production bounds", None)]:
    relaxed = milp.clone()
    pyo.TransformationFactory("core.relax_integer_vars").apply_to(relaxed)
    if bound is not None:
        relaxed.max_flowrate.set_value(bound)
    else:
        for r in relaxed.REACTORS:
            capacity = relaxed.product_flowrate / relaxed.conversion[r]
            relaxed.toggle_reactor[r].set_value(
                relaxed.reactor_feed[r] <= capacity * relaxed.reactor_boolean[r])
    result = highs.solve(relaxed)
    assert pyo.check_optimal_termination(result)
    print(f"{label}: {pyo.value(relaxed.reactor_boolean[1]):.8f} | "
          f"{pyo.value(relaxed.reactor_feed[1]):.4f} | {pyo.value(relaxed.cost):.7f}")

NLP declarations: {'continuous': 5, 'binary': 0, 'linear_equalities': 4, 'linear_inequalities': 0, 'lower_only': 5, 'both_bounds': 0}
MILP declarations: {'continuous': 5, 'binary': 2, 'linear_equalities': 4, 'linear_inequalities': 2, 'lower_only': 2, 'both_bounds': 5}
Fixed binaries | x1 | x2 | cost [USD/hr]

00: infeasible

10: 12.5000 | 0.0000 | 87.5000

01: 0.0000 | 15.0000 | 95.5000

11: 12.5000 | 0.0000 | 93.0000
LP bound | y1 | x1 | cost [USD/hr]

M=20: 0.62500000 | 12.5000 | 84.6875000

M=20000: 0.00062500 | 12.5000 | 80.0046875

production bounds: 1.00000000 | 12.5000 | 87.5000000

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}

(11)

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¶

Check the following, then explain your reasoning:

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}

(12)

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

Special forms and modeling applications¶

A common special MINLP form is

\min_{\mathbf{x},\mathbf{y}} \mathbf{c}^T\mathbf{y}+f(\mathbf{x})\quad\text{s.t. }\mathbf{h}(\mathbf{x})=\mathbf{0},\quad \mathbf{g}(\mathbf{x})+\mathbf{A}\mathbf{y}\le\mathbf{0},\quad \mathbf{x}\in X,\quad \mathbf{y}\in\{0,1\}^t.

(13)

Here fixed charges are linear in $\mathbf{y}$ and binaries enter the coupling inequalities linearly. The matrix $\mathbf{A}$ contains signed coefficients; it is not the positive scalar upper bound $M$ . This is BGW problem (P1), p. 513, a specialization of the general form above. Binaries can enter nonlinear equalities in the general form, so $\mathbf{h}(\mathbf{x})$ is not a restriction to impose on all MINLPs.

Operation	Result
Fix all binaries to a specified assignment	continuous NLP subproblem; it can be infeasible
Relax binaries to their bounds $0\le\mathbf{y}\le1$	continuous NLP relaxation
Use affine objective and constraint expressions, retaining binaries	MILP
Start with a MILP and fix or relax every integer variable	LP subproblem or LP relaxation

Replacing a nonlinear cost by an affine approximation changes the model; it is not an exact identity. For a general nonlinear relaxation, a local returned objective need not be a lower bound: the lower-bound statement requires the globally minimized relaxation.

Scheduling and planning use binaries for assignments, operating campaigns and facility locations. Unit commitment uses start/on decisions with continuous output; nonlinear heat-rate functions lead to MINLP, while suitable linear data/formulations lead to MILP. Discrete equipment design can select tray numbers, feed-tray location, exchanger connections or catalog sizes. Biegler (2010), Table 1.1, p. 6 and Exercise 6, p. 16, gives engineering examples. Further worked MO-book models: workforce shift scheduling, machine scheduling, facility location, dinner seating, and exam-room scheduling.

For a signed regression coefficient $\gamma_j\in[L_j,U_j]$ , with valid $L_j\le0\le U_j$ , a selection binary can impose

L_j y_j\le\gamma_j\le U_j y_j.

(14)

Required signed coefficient values for an omitted or selected term.

Both sides matter: $y_j=0$ then forces $\gamma_j=0$ . The nonnegative reactor’s one-sided switch does not apply unchanged to a coefficient that can be negative. A cardinality budget $\sum_jy_j\le k$ limits selected terms. The bound validity and choice of $k$ are additional modeling decisions.

Some nonsmooth functions have continuous reformulations: positive absolute-value penalties and max epigraphs can be LPs, as developed in the continuous formulation reference. Exact nonconvex graph relations or logical alternatives may require binaries. It is incorrect to claim that every absolute value or maximum introduces integer variables.

Broader modeling and theory references: H. P. Williams, Model Building in Mathematical Programming, 5th ed., Wiley, 2013; G. L. Nemhauser and L. A. Wolsey, Integer and Combinatorial Optimization, Wiley, 1988.

A bounded integer $q\in\{L,\ldots,U\}$ can be represented by $q=L+\sum_{k=0}^{b-1}2^k z_k$ , $z_k\in\{0,1\}$ , with $b=\lceil\log_2(U-L+1)\rceil$ and the additional restriction $\sum_k2^kz_k\le U-L$ . The restriction excludes unused binary codes when the range size is not a power of two. The case $U=L$ needs no binary. This represents a bounded integer domain; it does not linearize a nonlinear function of that integer.

Further literature on superstructure optimization¶

The lecture introduces the representation and a two-unit model. The extended reading list is retained here for independent study: reviews and frameworks first, then carbon capture/utilization, membrane separation, water and exchanger networks, energy systems, intensification and waste systems. These are the existing course references; the short lecture does not teach each application.

Naonori Nishida et al. (1981), A Review of Process Synthesis.
Qi Chen et al. (2017), Recent Developments and Challenges in Optimization-Based Process Synthesis.
Luca Mencarelli et al. (2020), A Review on Superstructure Optimization Approaches in Process System Engineering.
Efstratios N. Pistikopoulos et al. (2024), Advanced Modeling and Optimization Strategies for Process Synthesis.
J. Rafael Alcántara-Avila (2023), 50 Years of Optimization in Japan Using Superstructures.
Hector Yeomans et al. (1999), A Systematic Modeling Framework of Superstructure Optimization in Process Synthesis.
Zev Friedman et al. (2013), Block-Oriented Modeling of Superstructure Optimization Problems.
Sunwoo Kim et al. (2026), Decarbonization as a Systems and Control Challenge: From Superstructure Optimization to Multi-Period Decision-Making under Uncertainty.
Xuechong Ding et al. (2024), Superstructure-Based Carbon Capture and Utilization Process Design.
Natsayi Chiwaye et al. (2025), Superstructure-Based Optimization of Membrane Gas Separation Processes: A Review.
Rouzbeh Ramezani et al. (2023), Using a Superstructure Approach for Techno-Economic Analysis of Membrane Processes.
Ignacio E. Grossmann et al. (2014), Review of Optimization Models for Integrated Process Water Networks and Their Application to Biofuel Processes.
Elvis Ahmetović et al. (2011), Global Superstructure Optimization for the Design of Integrated Process Water Networks.
Ramkumar Karuppiah et al. (2006), Global Optimization for the Synthesis of Integrated Water Systems in Chemical Processes.
Siqi Liu et al. (2026), Advancements in Superstructures Design and Solver Approaches for Heat and Mass Exchange Network Synthesis.
Jian Gong et al. (2015), Sustainable Design and Synthesis of Energy Systems.
Philip Voll et al. (2013), Automated Superstructure-Based Synthesis and Optimization of Distributed Energy Supply Systems.
Mirko Skiborowski (2018), Process Synthesis and Design Methods for Process Intensification.
Channarong Puchongkawarin et al. (2020), Development of a Superstructure Optimization Framework for the Design of Municipal Solid Waste Facilities.

References¶

Nishida, N., Stephanopoulos, G., & Westerberg, A. W. (1981). A review of process synthesis. AIChE Journal, 27(3), 321–351. 10.1002/aic.690270302
Chen, Q., & Grossmann, I. E. (2017). Recent Developments and Challenges in Optimization-Based Process Synthesis. Annual Review of Chemical and Biomolecular Engineering, 8(1), 249–283. 10.1146/annurev-chembioeng-080615-033546
Mencarelli, L., Chen, Q., Pagot, A., & Grossmann, I. E. (2020). A review on superstructure optimization approaches in process system engineering. Computers & Chemical Engineering, 136, 106808. 10.1016/j.compchemeng.2020.106808
Pistikopoulos, E. N., & Tian, Y. (2024). Advanced Modeling and Optimization Strategies for Process Synthesis. Annual Review of Chemical and Biomolecular Engineering, 15(1), 81–103. 10.1146/annurev-chembioeng-100522-112139
Alcántara-Avila, J. R. (2023). 50 Years of Optimization in Japan Using Superstructures. Journal of Chemical Engineering of Japan, 56(1). 10.1080/00219592.2023.2188112
Yeomans, H., & Grossmann, I. E. (1999). A systematic modeling framework of superstructure optimization in process synthesis. Computers & Chemical Engineering, 23(6), 709–731. 10.1016/s0098-1354(99)00003-4
Friedman, Z., Ingalls, J., Siirola, J. D., & Watson, J.-P. (2013). Block-oriented modeling of superstructure optimization problems. Computers & Chemical Engineering, 57, 10–23. 10.1016/j.compchemeng.2013.04.008
Kim, S., & Lee, J. H. (2026). Decarbonization as a systems and control challenge: From superstructure optimization to multi-period decision-making under uncertainty. Journal of Process Control, 164, 103775. 10.1016/j.jprocont.2026.103775
Ding, X., Li, J., Chen, H., & Zhou, T. (2024). Superstructure-based carbon capture and utilization process design. Current Opinion in Chemical Engineering, 43, 100995. 10.1016/j.coche.2023.100995
Chiwaye, N., Majozi, T., & Daramola, M. O. (2025). Superstructure-Based Optimization of Membrane Gas Separation Processes: A Review. Industrial & Engineering Chemistry Research, 64(28), 13905–13919. 10.1021/acs.iecr.4c03693
Ramezani, R., Randon, A., Felice, L. D., & Gallucci, F. (2023). Using a superstructure approach for techno-economic analysis of membrane processes. Chemical Engineering Research and Design, 199, 296–311. 10.1016/j.cherd.2023.10.007
Grossmann, I. E., Martín, M., & Yang, L. (2014). Review of optimization models for integrated process water networks and their application to biofuel processes. Current Opinion in Chemical Engineering, 5, 101–109. 10.1016/j.coche.2014.07.003
Karuppiah, R., & Grossmann, I. E. (2006). Global optimization for the synthesis of integrated water systems in chemical processes. Computers & Chemical Engineering, 30(4), 650–673. 10.1016/j.compchemeng.2005.11.005
Skiborowski, M. (2018). Process synthesis and design methods for process intensification. Current Opinion in Chemical Engineering, 22, 216–225. 10.1016/j.coche.2018.11.004
Puchongkawarin, C., & Mattaraj, S. (2020). Development of a superstructure optimization framework for the design of municipal solid waste facilities. Sustainable Environment Research, 30(1). 10.1186/s42834-020-00071-7