Practical Applications of Pyomo: Optimization Problems and Solutions

This assignment provides a comprehensive exploration of optimization problems using the Pyomo library in Python. It begins with an introduction to Pyomo and its basic functionalities, including importing necessary libraries and defining decision variables, objective functions, and constraints. The assignment then delves into several real-world applications, starting with a product mix problem where the goal is to maximize profit by determining the optimal number of plants of different types, subject to fertilizer constraints. The solutions are presented with Python code, demonstrating how to define and solve these problems using Pyomo. The assignment extends the product mix problem by adding additional constraints, such as storage capacity. It further explores the concept of the dual problem, providing the primal and dual formulations and their solutions. The assignment also includes product mix challenges, golf bag manufacturing, baseball glove production, bike frame design, investment portfolio optimization, production scheduling, assembly line equipment, rock crushing, hospital scheduling, and butter production problems. For each problem, the assignment provides the linear program (LP) formulation, the corresponding Pyomo code, and the business decisions derived from the solutions.

Statistics: Optimization in Pyomo
Introduction
The first step in solving the optimization problems in Pyomo is to import libraries using the
code:
from pylab import *
import shutil
import sys
import os.path
from pyomo.environ import *
Product Mix Problem
Let the number of plant type A be A and type be B. Then the linear program (LP) is
max profit
A , B
: 2.25 A+ 2.60 B
Subject to:
2 A +B≤ 4000 (fertilizer 1 constraint)
A+2 B≤ 5000 (fertilizer 2 constraint)
A , B ≥ 0 (Non-negativity constraint)
The python code for solving the problem is:
model = ConcreteModel(name="(Product_Mix)")
# declare decision variables
model.A = Var(domain=NonNegativeReals)
model.B = Var(domain=NonNegativeReals)
# declare objective
model.profit = Objective(
expr = 2.25*model.A + 2.6*model.B,
sense = maximize)
# declare constraints
model.fertilzer1 = Constraint(expr = 2*model.A + model.B <= 4000)
model.fertilizer2 = Constraint(expr = model.A + 2*model.B <= 5000)
# solve
SolverFactory('glpk').solve(model)
# Print the business decision
print("The business decision is to plant ", int (model.A()), " of Roses
and ", int (model.B()),
" of Begonias for a profit of $", model.profit(), " per season")

The business decision is to plant 1000 of Roses and 2000 of Begonias for a profit of $ 7450.0 per
season.
Product Mix Continued
The additional constraint in the product mix problem is on the available fertilizer storage
capacity. Then, the LP becomes:
max profit
A , B
: 2.25 A+ 2.60 B
Subject to:
2 A +B≤ 4000 (fertilizer 1 constraint)
A+2 B≤ 5000 (fertilizer 2 constraint)
3 A +3 B≤ 8000 (capacity constraint)
A , B ≥ 0 (Non-negativity constraint)
model = ConcreteModel(name="(Product Mix Continued)")
# declare decision variables
model.A = Var(domain=NonNegativeReals)
model.B = Var(domain=NonNegativeReals)
# declare objective
model.profit = Objective(
expr = 2.25*model.A + 2.6*model.B,
sense = maximize)
# declare constraints
model.fertilzer1 = Constraint(expr = 2*model.A + model.B <= 4000)
model.fertilizer2 = Constraint(expr = model.A + 2*model.B <= 5000)
model.storage = Constraint(expr = 3*model.A + 3*model.B <= 8000)
# solve
SolverFactory('glpk').solve(model)
#Business decision
print("To maximize profit, the decision is to plant ", int (model.A()), "
of roses and ",
int (model.B()), " of begonias. Therefore, making a profit of $",
int(model.profit())"." )
To maximize profit, the decision is to plant 333 of roses and 2333 of begonias. Therefore,
making a profit of $ 6816.
The Dual Problem
The primal problem above is written as dual.

min profit
y1 , y1
: 4000 y1+5000 y2 +8000 y3
Subject to:
2 y1+ y2+3 y1 ≤2.25 (plant A)
y1 + y2+3 y1 ≤ 2.60 (plant B)
y1 , y2 , y1 ≥ 0 (Non-negativity constraint)
The langrage optimal are obtained as follows:
model = ConcreteModel(name="(Product Mix Continued Dual)")
# declare decision variables
model.y1 = Var(domain=NonNegativeReals)
model.y2 = Var(domain=NonNegativeReals)
model.y3 = Var(domain=NonNegativeReals)
# declare objective
model.profit = Objective(
expr = 4000*model.y1 + 5000*model.y2 + 8000*model.y3,
sense = minimize)
# declare constraints
model.Aconst = Constraint(expr = 2*model.y1 + model.y2 + 3*model.y3 >=
2.25)
model.Bconst = Constraint(expr = model.y1 + model.y2 + 3*model.y3 >= 2.6)
# solve
SolverFactory('glpk').solve(model)
# optiamal values of dual problem
print("y1 = ", model.y1(), "y2 = ", model.y2(), "y3 = ", model.y3())
y1=0 , y2=0 , and y3=0.87
Using these values to get the number of plants A and B, and the profits:
model = ConcreteModel(name="(Product Mix Continued dual continued)")
# declare decision variables
model.A = Var(domain=NonNegativeReals)
model.B = Var(domain=NonNegativeReals)
# declare objective
model.profit = Objective(
expr = -2.25*model.A - 2.6*model.B,
sense = minimize)
# declare constraints
model.dual = Constraint(expr = 0.87*(3*model.A + 3*model.B - 8000)<=0)

# solve
SolverFactory('glpk').solve(model)
#Business decision
print("To maximize profit, the decision is to plant ", int (model.A()), "
of roses and ",
int (model.B()), " of begonias. Therefore, making a profit of $", -
int(model.profit()) )
To maximize profit, the decision is to plant 0 of roses and 2666 of begonias. Therefore, making a
profit of $ 6933.
Product Mix Challenge
The LP problem is defined as follows:
max profit
A , B
: A +B
Subject to:
5 A +3 B≤ 120 (blobs constraint)
3 A +5 B≤ 120 (globs constraint)
A , B ≥ 0 (Non-negativity constraint)
model = ConcreteModel(name="(Product Mix Challenge)")
# declare decision variables
model.A = Var(domain=NonNegativeReals)
model.B = Var(domain=NonNegativeReals)
# declare objective
model.profit = Objective(
expr = model.A + model.B,
sense = maximize)
# declare constraints
model.blobs = Constraint(expr = 5*model.A + 3*model.B <= 120)
model.globs = Constraint(expr = 3*model.A + 5*model.B <= 120)
# solve
SolverFactory('glpk').solve(model)
# print the results
print("To maximize profit, produce ", int (model.A()), " of product A and
",
int (model.B()), " of product B. Therefore, making a profit of $",
model.profit() )

To maximize profit, produce 15 of product A and 15 of product B. Therefore, making a profit of
$ 30.0
Golf Bags
Let the number of standard bags manufactured be S and the deluxe bags be D, then the LP is:
max profit
S , D
:10 S+ 9 D
Subject to:
42
60 S+D ≤ 630 (cutting time constraint)
0.5 S+ 50
60 D ≤600 (sewing time constraint)
S+ 40
60 D≤ 708 (finishing constraint)
6
60 S+ 15
60 D ≤135 (finishing constraint)
S , D≥ 0 (Non-negativity constraint)
model = ConcreteModel(name="(Golf Bags)")
# declare decision variables
model.S = Var(domain=NonNegativeReals)
model.D = Var(domain=NonNegativeReals)
# declare objective
model.profit = Objective(
expr = 10*model.S + 9*model.D,
sense = maximize)
# declare constraints
model.cutting = Constraint(expr = (42/60)*model.S + model.D <= 630)
model.sewing = Constraint(expr = 0.5*model.S + (50/60)*model.D <= 600)
model.finishing = Constraint(expr = model.S + (40/60)*model.D <= 708)
model.Inspection_packaging = Constraint(expr = (6/60)*model.S +
(15/60)*model.D <= 135)
# solve
SolverFactory('glpk').solve(model)
# Print the results
print("To maximize profit, manufacture ", int (model.S()), " standard bags
and ",
int (model.D()), " delux bags over the next three months. Therefore,
making a profit of $", model.profit() )

To maximize profit, manufacture 540 standard bags and 252 deluxe bags over the next three
months. Therefore, making a profit of $ 7668.
Baseball Gloves
Let the number of fielder gloves manufactured be F and catcher be C, then the LP is:
max profit
A , B
:5 F +8 C
Subject to:
F+ 1.5C ≤ 900 (cutting and sewing constraint)
0.5 F+ 1
3 C ≤300 (finishing constraint)
12.5
60 F+ 15
60 C ≤135 (packaging and shipping constraint)
F , C ≥ 0 (Non-negativity constraint)
model = ConcreteModel(name="(Baseball Gloves)")
# declare decision variables
model.F = Var(domain=NonNegativeReals)
model.C = Var(domain=NonNegativeReals)
# declare objective
model.profit = Objective(
expr = 5*model.F + 8*model.C,
sense = maximize)
# declare constraints
model.cutting_sewing = Constraint(expr = model.F + 1.5*model.C <= 900)
model.Finishing = Constraint(expr = 0.5*model.F + (1/3)*model.C <= 300)
model.packaging_shipping = Constraint(expr = (12.5/60)*model.F +
(15/60)*model.C <= 135)
# solve
SolverFactory('glpk').solve(model)
print("To maximize profit, manufacture ", int (model.F()), " Fielder mitts
and ",
int (model.C()), " Catcher over the next month. Therefore, making a
profit of $", model.profit() )
To maximize profit, manufacture 0 Fielder mitts and 540 Catcher over the next month.
Therefore, making a profit of $ 4320.
Bike Frames

Let the yards of standard material used be S and of professional be P. The LP is
min cost
S , D
:7.5 S+9 P
Subject to:
0.84 S+ 0.58 P ≤ 0.7(30) (fiber glass constraint)
0. 10 S+0.30 P ≤ 0.2(30) (carbon fiber constraint)
0.06 S+0.12 P ≤0.1(30) (Kevlar constraint)
S+ P ≤3 0 (length constraint)
S , P ≥ 0 (Non-negativity constraint)
model = ConcreteModel(name="(Bike Frames)")
# declare decision variables
model.S = Var(domain=NonNegativeReals)
model.P = Var(domain=NonNegativeReals)
# declare objective
model.cost = Objective(
expr = 7.5*model.S + 9*model.P,
sense = minimize)
# declare constraints
model.fiberglass = Constraint(expr = 0.84*model.S + 0.58*model.P <=
0.70*30)
model.carbonfiber = Constraint(expr = 0.10*model.S + 0.30*model.P >=
0.20*30)
model.kevlar = Constraint(expr = 0.06*model.S + 0.12*model.P <= 0.10*30)
model.length = Constraint(expr = model.S + model.P == 30)
# solve
SolverFactory('glpk').solve(model)
# print the results
print("To minimize cost, blend ", int (model.S()), " yards of standard
material and ",
int (model.P()),
" yards of professional material in making each frame. Therefore,
using a minimum cost of $", int (model.cost()) )
To minimize cost, blend 13 yards of standard material and 16 yards of professional material in
making each frame. Therefore, using a minimum cost of $ 249.
Investment Portfolio

This is an allocation problem which can be solved using optimization. Let the amount invested in
growth be X1, income be X2 and money market fund be X3. The liner program becomes:
min risk
X1 , X 2 , X 3
=0.1 ( X1
1000000 ) + 0.05 ( X2
1000000 )+ 0.0 1 ( X3
1000000 )
Subject to:
X1
1000000 ≥ 0.1 (growth constraint)
X2
1000000 ≥ 0.1 (income constraint)
X3
1000000 ≥ 0 (market fund constraint)
X1 + X 2+ X3=1000000 (investment constraint)
0.1 ( X1
1000000 )+0. 05 ( X2
1000000 )+ 0.0 1 ( X3
1000000 )≤ 0.05 (maximum risk constraint)
X1 , X2 , X3 ≥ 0 (Non-negativity constraint)
model = ConcreteModel(name="(Investment Portfolio)")
# declare decision variables
model.growth = Var(domain=NonNegativeReals)
model.income = Var(domain=NonNegativeReals)
model.market_fund = Var(domain=NonNegativeReals)
# declare objective
model.risks = Objective(
expr = 0.1*(model.growth/1000000) + 0.05*(model.income/1000000) +
0.01*(model.market_fund/1000000),
sense = minimize)
# declare constraints
model.growthconst = Constraint(expr = model.growth/1000000 >= 0.1)
model.incomeconst = Constraint(expr = model.income/1000000 >= 0.1)
model.market_fundconst = Constraint(expr = model.market_fund/1000000 >=
0.2)
model.investmentconst = Constraint(expr = model.growth + model.income +
model.market_fund == 1000000)
model.riskconst = Constraint(expr = 0.1*(model.growth/1000000) +
0.05*(model.income/1000000) + 0.01*(model.market_fund/1000000)
<= 0.05)

# solve
SolverFactory('glpk').solve(model)
print("To minimize risks, the client's $1,000,000 should be invested as
follows: $", int(model.growth()),
" on growth, $", int (model.income()), " income and $",
int(model.market_fund()),
" on money market fund at a minimum overall risk of ",
100*model.risks(), "%.")
To minimize risks, the client's $1,000,000 should be invested as follows: $ 100000 on growth, $
100000 income and $ 800000 on money market fund at a minimum overall risk of 2.3 %.
Production Scheduling
model = ConcreteModel(name="(Production scheduling)")
# declare decision variables
model.hours_M100 = Var(domain=NonNegativeReals)
model.hours_M200 = Var(domain=NonNegativeReals)
# declare objective
model.profit = Objective(
expr = 160*model.hours_M100 + 345*model.hours_M200,
sense = maximize)
# declare constraints
model.chemical = Constraint(expr = 40*model.hours_M100 +
50*model.hours_M200 <= 1000)
model.Maxhours_M100 = Constraint(expr = model.hours_M100 <= 15)
model.Maxhours_M200 = Constraint(expr = model.hours_M200 <= 10)
model.Minhours_M100 = Constraint(expr = model.hours_M100 >= 5)
model.Minhours_M200 = Constraint(expr = model.hours_M200 >= 5)
# solve
SolverFactory('glpk').solve(model)
# print the results
print("The business decision for next week is to operate M-100 moulding
machine for ", int (model.hours_M100()),
" hours and M-200 machine for ", int (model.hours_M200())," hours
per week, and generate a profit of $",
int (model.profit()), ".")
The business decision for next week is to operate M-100 moulding machine for 12 hours and M-
200 machine for 10 hours per week, and generate a profit of $ 5450.
Assembly Line Equipment

model = ConcreteModel(name="(Assembly Line Equipment)")
# declare decision variables
model.M1 = Var(domain=NonNegativeReals)
model.M2 = Var(domain=NonNegativeReals)
model.M3 = Var(domain=NonNegativeReals)
# declare objective
model.cost = Objective(
expr = 18500*model.M1 + 25000*model.M2 + 35000*model.M3,
sense = minimize)
# declare constraints
model.type1 = Constraint(expr = 100*model.M1 + 265*model.M2 + 200*model.M3
>= 3200)
model.type2 = Constraint(expr = 130*model.M1 + 235*model.M2 + 160*model.M3
>= 2500)
model.type3 = Constraint(expr = 140*model.M1 + 170*model.M2 + 260*model.M3
>= 3500)
model.type4 = Constraint(expr = 210*model.M1 + 220*model.M2 + 180*model.M3
>= 3000)
model.type5 = Constraint(expr = 80*model.M1 + 120*model.M2 + 220*model.M3
>= 2500)
# solve
SolverFactory('glpk').solve(model)
# print the results
print("In order to most efficiently meet the production needs, the company
should buy ", int (model.M1()),
" type I machine, ", int (model.M2())," type II machine, and ", int
(model.M3()),
" type III machine at a total cost of $", int (model.cost()), ".")
In order to most efficiently meet the production needs, the company should buy 9 type I
machine, 4 type II machine, and 5 type III machine at a total cost of $ 477057 .
Crushing More Rocks
model = ConcreteModel(name="(Crushing More Rocks)")
# declare decision variables
model.fine = Var(domain=NonNegativeReals)
model.medium = Var(domain=NonNegativeReals)
model.coarse = Var(domain=NonNegativeReals)
# declare objective
model.cost = Objective(
expr = 8*model.fine + 5*model.medium + 3*model.coarse,
sense = minimize)
# declare constraints

model.limestone = Constraint(expr = 0.5*model.fine + 0.2*model.medium +
0.05*model.coarse >= 50)
model.chat = Constraint(expr = 0.3*model.fine + 0.4*model.medium +
0.2*model.coarse >= 60)
model.redimix = Constraint(expr = 0.2*model.fine + 0.3*model.medium +
0.35*model.coarse >= 70)
model.rough = Constraint(expr = 0*model.fine + 0.1*model.medium +
0.4*model.coarse >= 30)
model.capacity = Constraint(expr = model.fine + model.medium +
model.coarse >= 50)
# solve
SolverFactory('glpk').solve(model)
# print the results
print("The company should process ", int (model.fine()), "tones,", int
(model.medium()), "tones, and ", int (model.coarse()),
" tones of fine, medium, and coarse respectively on each setting.",
"The total processing cost for the coming week will be $", int
(model.cost()),
".")
The company should process 76 tones, 25 tones, and 134 tons of fine, medium, and coarse
respectively on each setting. The total processing cost for the coming week will be $ 1141.
Hospital Scheduling
model = ConcreteModel(name="(Hospital Scheduling)")
# declare decision variables
model.facelift = Var(domain=NonNegativeReals)
model.lipo = Var(domain=NonNegativeReals)
model.implants = Var(domain=NonNegativeReals)
# declare objective
model.profit = Objective(
expr = 240*model.facelift + 225*model.lipo + 425*model.implants,
sense = maximize)
# declare constraints
model.days = Constraint(expr = 3*model.facelift + 5*model.lipo +
6*model.implants <= 70)
model.hours = Constraint(expr = 2*model.facelift + 1.5*model.lipo +
3*model.implants <= 165)
# solve
SolverFactory('glpk').solve(model)
# print the results

print("The hospital should admit ", int (model.facelift()), " face lift
patients,", int (model.lipo()),
" lipo patients, and ", int (model.implants()), " implant patients,
each week.",
"This will give the hospital a profit of $", int (model.profit()), "
per week.")
The hospital should admit 23 face lift patients, 0 lipo patients, and 0 implant patients, each
week. This will give the hospital a profit of $ 5599 per week.
Butter
model = ConcreteModel(name="(Butter)")
# declare decision variables
model.peanut = Var(domain=NonNegativeReals)
model.apple = Var(domain=NonNegativeReals)
# declare objective
model.profit = Objective(
expr = 1100*model.peanut + 1300*model.apple,
sense = maximize)
# declare constraints
model.sterilization = Constraint(expr = 4*model.peanut + 5*model.apple <=
40)
model.packaging = Constraint(expr = 6*model.peanut + 4*model.apple <= 40)
model.appleconst = Constraint(expr = model.apple >= 5)
# solve
SolverFactory('glpk').solve(model)
# print the results
print("Quart Industries should produce ", int(model.peanut()), " batches
of 1000 jars of peanut butter, and ", int(model.apple())
, " batches of 1000 jars of apple butter this week. This will give
the Quart's a profit of $",
int (model.profit()), " this week.")
Quart Industries should produce 2 batches of 1000 jars of peanut butter, and 5 batches of 1000
jars of apple butter this week. This will give the Quart's a profit of $ 10571 this week.