2.7. Non-Continuous Objectives#
Keywords: linear programming, cbc usage, production models
This notebook demonstrates the use of linear programming to maximize profit for a simple model of a multiproduct production facility. The notebook uses Pyomo to represent the model with the COINOR-CBC solver to calculate solutions.
2.7.1. Imports#
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import shutil
import sys
import os.path
if not shutil.which("pyomo"):
!pip install -q pyomo
assert(shutil.which("pyomo"))
if not (shutil.which("cbc") or os.path.isfile("cbc")):
if "google.colab" in sys.modules:
!apt-get install -y -qq coinor-cbc
else:
try:
!conda install -c conda-forge coincbc
except:
pass
assert(shutil.which("cbc") or os.path.isfile("cbc"))
from pyomo.environ import *
2.7.2. Example: Production Bonuses#
A plant produces three products in the amounts \(x\), \(y\), and \(z\) with unit profit of \\(40, \\\)30, and \$50, respectively. There are several constraints imposed by product demand and the availability of specialized labor
In addition, the plant receives bonuses for meeting certain production targets
If the plant produces more than 20 \(y\) items, then the unit profit for \(y\) will be \\(50 plus a fixed bonus of \\\)200.
If the plant produces more than 30 \(z\) items, then the unit profit for \(z\) will be \\(60 plus a fixed bonus of \\\)300.
Find the optimal production targets.
2.7.3. Model 1. Without bonuses#
from pyomo.environ import *
model = ConcreteModel()
model.x = Var(domain=NonNegativeReals)
model.y = Var(domain=NonNegativeReals)
model.z = Var(domain=NonNegativeReals)
model.profit = Objective(expr = 40*model.x + 30*model.y + 50*model.z, sense=maximize)
model.demand = Constraint(expr = model.x <= 40)
model.laborA = Constraint(expr = model.x + model.y <= 80)
model.laborB = Constraint(expr = 2*model.x + model.z <= 100)
model.laborC = Constraint(expr = model.z <= 50)
# solve
SolverFactory('cbc').solve(model).write()
# ==========================================================
# = Solver Results =
# ==========================================================
# ----------------------------------------------------------
# Problem Information
# ----------------------------------------------------------
Problem:
- Name: unknown
Lower bound: -5150.0
Upper bound: -5150.0
Number of objectives: 1
Number of constraints: 5
Number of variables: 4
Number of nonzeros: 2
Sense: maximize
# ----------------------------------------------------------
# Solver Information
# ----------------------------------------------------------
Solver:
- Status: ok
User time: -1.0
System time: 0.0
Wallclock time: 0.0
Termination condition: optimal
Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
Statistics:
Branch and bound:
Number of bounded subproblems: None
Number of created subproblems: None
Black box:
Number of iterations: 1
Error rc: 0
Time: 0.02187371253967285
# ----------------------------------------------------------
# Solution Information
# ----------------------------------------------------------
Solution:
- number of solutions: 0
number of solutions displayed: 0
The results of the solution step show the solver has converged to an optimal solution. Next we display the particular components of the model of interest to us.
print(f"Profit = ${model.profit()}")
print(f"X = {model.x()} units")
print(f"Y = {model.y()} units")
print(f"Z = {model.z()} units")
Profit = $5150.0
X = 25.0 units
Y = 55.0 units
Z = 50.0 units
2.7.4. Model 2. Include bonuses#
To incorporate the bonus structure into the objective function, we introduce new decision variables
where \(y\) and \(z\) denote total production, and \(y_1\) and \(z_1\) denote production above bonus levels. In addition, we introduce two binary variables, \(b_y\) and \(b_z\), to indicate when \(y\) or \(z\), respectively, are in a bonus condition.
The profit objective becomes (after some arithmetic)
For a sufficiently large value of \(M\), the constraints
enforce the bonus conditions.
Let’s verify this formulation and new objective before proceeding further.
y = np.linspace(0, 100, 1000)
profit_y = [30*_ if _ <= 20 else (50*_ + 200) for _ in y]
y0 = np.array([min(20, _) for _ in y])
y1 = y - y0
by = np.array([1 if _ > 0 else 0 for _ in y1])
bonus_y = 30*y0 + 50*y1 + 600*by
fig, ax = plt.subplots(2, 2)
ax[0,0].plot(y, y0, label='y0')
ax[0,0].plot(y, y1, label='y1')
ax[0,0].plot(y, y0 + y1, '--', label='y')
ax[0,0].legend()
ax[1,0].plot(y, profit_y)
ax[1,0].plot(y, bonus_y)
ax[1,0].set_title("profit_y")
z = np.linspace(0, 100, 1000)
profit_z = [50*_ if _ <= 30 else (60*_ + 300) for _ in z]
z0 = np.array([min(30, _) for _ in z])
z1 = z - z0
bz = np.array([1 if _ > 0 else 0 for _ in z1])
bonus_z = 50*z0 + 60*z1 + 600*bz
ax[0,1].plot(z, z0, label='z0')
ax[0,1].plot(z, z1, label='z1')
ax[0,1].plot(z, z0 + z1, '--', label='z')
ax[0,1].legend()
ax[1,1].plot(y, profit_z)
ax[1,1].plot(y, bonus_z)
ax[1,1].set_title("profit_z")
fig.tight_layout()
The corresponding Pyomo model follows. Note that this could be streamlined in various ways.
from pyomo.environ import *
M = 100
model = ConcreteModel()
model.x = Var(domain=NonNegativeReals)
model.y = Var(domain=NonNegativeReals)
model.y0 = Var(bounds=(0, 20))
model.y1 = Var(domain=NonNegativeReals)
model.by = Var(domain=Binary)
model.z = Var(domain=NonNegativeReals)
model.z0 = Var(bounds=(0, 30))
model.z1 = Var(domain=NonNegativeReals)
model.bz = Var(domain=Binary)
model.profit = Objective(sense=maximize, expr =
+ 40*model.x \
+ 30*model.y0 + 50*model.y1 + 600*model.by \
+ 50*model.z0 + 60*model.z1 + 600*model.bz)
model.dy = Constraint(expr = model.y == model.y0 + model.y1)
model.dz = Constraint(expr = model.z == model.z0 + model.z1)
model.demand = Constraint(expr = model.x <= 40)
model.laborA = Constraint(expr = model.x + model.y <= 80)
model.laborB = Constraint(expr = 2*model.x + model.z <= 100)
model.laborC = Constraint(expr = model.z <= 50)
model.bonus_y0 = Constraint(expr = model.y0 >= 20 - M*(1 - model.by))
model.bonus_z0 = Constraint(expr = model.z0 >= 30 - M*(1 - model.bz))
model.bonus_y1 = Constraint(expr = model.y1 <= M*model.by)
model.bonus_z1 = Constraint(expr = model.z1 <= M*model.bz)
# solve
SolverFactory('cbc').solve(model).write()
# ==========================================================
# = Solver Results =
# ==========================================================
# ----------------------------------------------------------
# Problem Information
# ----------------------------------------------------------
Problem:
- Name: unknown
Lower bound: -7500.0
Upper bound: -7500.0
Number of objectives: 1
Number of constraints: 7
Number of variables: 8
Number of binary variables: 2
Number of integer variables: 2
Number of nonzeros: 7
Sense: maximize
# ----------------------------------------------------------
# Solver Information
# ----------------------------------------------------------
Solver:
- Status: ok
User time: -1.0
System time: 0.05
Wallclock time: 0.06
Termination condition: optimal
Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
Statistics:
Branch and bound:
Number of bounded subproblems: 0
Number of created subproblems: 0
Black box:
Number of iterations: 0
Error rc: 0
Time: 0.08483672142028809
# ----------------------------------------------------------
# Solution Information
# ----------------------------------------------------------
Solution:
- number of solutions: 0
number of solutions displayed: 0
print(f"Profit = ${model.profit()}")
print(f"X = {model.x()} units")
print(f"Y = {model.y0()} + {model.y1()} = {model.y()} units {model.by()}")
print(f"Z = {model.z0()} + {model.z1()} = {model.z()} units {model.bz()}")
Profit = $7500.0
X = 0.0 units
Y = 20.0 + 60.0 = 80.0 units 1.0
Z = 30.0 + 20.0 = 50.0 units 1.0
2.7.6. Exercise#
Reformulate this model using the Pyomo.GDP package. Parameterize values for the profit objective, bonus amounts, and bonus levels.
2.7.5. Comments#
This particular formulation leaves much to be desired.
The main thing is a lack of numerical robustness on the constraints for \(z_0\) and \(y_0\). * * It would also be useful to parameterize and general the bonus structure.
The big M method often yields redundant constraints and numerical performance that depends on the choice of M. Generalized disjunctive constraints provide more solution options.