1.4. Continuous Optimization: Linear Programming#
# 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()
import pandas as pd
import pyomo.environ as pyo
import idaes
1.4.1. Linear Programs: Student Diet Example#
Reference: https://docs.mosek.com/modeling-cookbook/linear.html
You want to save money eating while remaining healthy. A healthy diet requires at least P=6 units of protein, C=15 units of carbohydrates, F=5 units of fats and V=7 units of vitamins. Due to compounding factors (blizzard during Lent), our campus only has these options:
# Load data from file, use the first column (0, recall Python starts counting at 0) as the index
food_options = pd.read_csv('https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/data/student_diet.csv',index_col=0)
# Print up the the first 10 rows of data
food_options.head(10)
| P | C | F | V | price | |
|---|---|---|---|---|---|
| takeaway | 3.0 | 3 | 2 | 1 | 5 |
| vegetables | 1.0 | 2 | 0 | 4 | 1 |
| bread | 0.5 | 4 | 1 | 0 | 2 |
Let’s build a Python dictionary to store the nutrient requirements. (I strongly recommend not touching Python until we write the model on paper. I am including this here to avoid scrolling between the problem description and this cell.)
# Uncomment and fill in with all of the data
# nutrient_requirements = {'P':6, 'C':15 }
# Add your solution here
1.4.1.1. Propose an Optimization Model#
Sets
Click to expand
\[
\mathcal{F} = \{\text{take away, veggies, bread}\}
\]
\[
\mathcal{N} = \{P, C, F, V\}
\]
Parameters
Click to expand
\(\pi_i\): price for food option \(i \in \mathcal{F}\)
\(r_j\): nutrition requirement for nutrient \(j \in \mathcal{N}\)
\(D_{ij}\): nutrition info for food \(i \in \mathcal{F}\) and nutrient \(j \in \mathcal{N}\)
Variables
Click to expand
\(x_i\): amount of food \(i \in \mathcal{F}\) eaten/purchased
Objective and Constraints
Click to expand
Minimize: \( \sum_{i \in \mathcal{F}} x_i \cdot \pi_i \quad \text{(cost)} \)
Subject to: \( \sum_{i \in \mathcal{F}} D_{ij} \cdot x_i \geq r_j \quad \forall j \in \mathcal{N} \quad \text{(be healthy)} \)
Bounds: \( x_i \geq 0 \quad \forall i \in \mathcal{F} \)
Full Model
Let’s put all of the components together into the full model written in compact notation.
\[\begin{split}
\begin{align*}
\min_{x} \quad & \sum_{i \in \mathcal{F}} x_i \cdot \pi_i & \text{(cost)} \\
\text{s.t.} \quad & \sum_{i \in \mathcal{F}} D_{ij} \cdot x_i \geq r_j, \quad \forall j \in \mathcal{N} & \text{(be healthy)} \\
& x_i \geq 0, \quad \forall i \in \mathcal{F} & \text{(non-negative bounds)}
\end{align*}
\end{split}\]
Degree of Freedom Analysis
We will later learn more about how to factor inequality constraints into degree of freedom analysis. For now, please count the number of equality and inequality constraints separately.
Click to expand
Continuous variables: 3
Inequality constraints: 4
Bounds: 3
Question: Why is it okay to have more inequality constraints than variables?
Answer: Not all inequality constraints may be active at a time.
1.4.1.2. Solve in Pyomo#
With our optimization model written on paper, we can proceed to solve in Pyomo. Before we start, let’s review a few pandas tricks.
# Extract the column names, convert to a list
food_options.columns.to_list()
['P', 'C', 'F', 'V', 'price']
# Same as above, but drop the last entry
nutrients = food_options.columns.to_list()[0:4]
nutrients
['P', 'C', 'F', 'V']
# Extract the index names, convert to a list
foods = food_options.index.to_list()
foods
['takeaway', 'vegetables', 'bread']
# Create a dictionary with keys such as ('takeaway', 'P')
# Do not include 'price'
food_info = food_options[nutrients].stack().to_dict()
food_info
{('takeaway', 'P'): 3.0,
('takeaway', 'C'): 3.0,
('takeaway', 'F'): 2.0,
('takeaway', 'V'): 1.0,
('vegetables', 'P'): 1.0,
('vegetables', 'C'): 2.0,
('vegetables', 'F'): 0.0,
('vegetables', 'V'): 4.0,
('bread', 'P'): 0.5,
('bread', 'C'): 4.0,
('bread', 'F'): 1.0,
('bread', 'V'): 0.0}
# Create dictionary of only prices
price = food_options['price'].to_dict()
price
{'takeaway': 5, 'vegetables': 1, 'bread': 2}
Now let’s build our Pyomo model!
# Add your solution here
2 Set Declarations
FOOD : Size=1, Index=None, Ordered=Insertion
Key : Dimen : Domain : Size : Members
None : 1 : Any : 3 : {'takeaway', 'vegetables', 'bread'}
NUTRIENTS : Size=1, Index=None, Ordered=Insertion
Key : Dimen : Domain : Size : Members
None : 1 : Any : 4 : {'P', 'C', 'F', 'V'}
3 Param Declarations
food_info : Size=12, Index=FOOD*NUTRIENTS, Domain=Any, Default=None, Mutable=False
Key : Value
('bread', 'C') : 4.0
('bread', 'F') : 1.0
('bread', 'P') : 0.5
('bread', 'V') : 0.0
('takeaway', 'C') : 3.0
('takeaway', 'F') : 2.0
('takeaway', 'P') : 3.0
('takeaway', 'V') : 1.0
('vegetables', 'C') : 2.0
('vegetables', 'F') : 0.0
('vegetables', 'P') : 1.0
('vegetables', 'V') : 4.0
needs : Size=4, Index=NUTRIENTS, Domain=Any, Default=None, Mutable=False
Key : Value
C : 15
F : 5
P : 6
V : 7
price : Size=3, Index=FOOD, Domain=Any, Default=None, Mutable=False
Key : Value
bread : 2
takeaway : 5
vegetables : 1
1 Var Declarations
food_eaten : Size=3, Index=FOOD
Key : Lower : Value : Upper : Fixed : Stale : Domain
bread : 0 : 1.0 : None : False : False : NonNegativeReals
takeaway : 0 : 1.0 : None : False : False : NonNegativeReals
vegetables : 0 : 1.0 : None : False : False : NonNegativeReals
1 Objective Declarations
cost : Size=1, Index=None, Active=True
Key : Active : Sense : Expression
None : True : minimize : 5*food_eaten[takeaway] + food_eaten[vegetables] + 2*food_eaten[bread]
1 Constraint Declarations
diet_min : Size=4, Index=NUTRIENTS, Active=True
Key : Lower : Body : Upper : Active
C : 15.0 : 3.0*food_eaten[takeaway] + 2.0*food_eaten[vegetables] + 4.0*food_eaten[bread] : +Inf : True
F : 5.0 : 2.0*food_eaten[takeaway] + 0.0*food_eaten[vegetables] + food_eaten[bread] : +Inf : True
P : 6.0 : 3.0*food_eaten[takeaway] + food_eaten[vegetables] + 0.5*food_eaten[bread] : +Inf : True
V : 7.0 : food_eaten[takeaway] + 4.0*food_eaten[vegetables] + 0.0*food_eaten[bread] : +Inf : True
8 Declarations: FOOD NUTRIENTS needs food_info price food_eaten diet_min cost
Click to see the solution to the activity
m = pyo.ConcreteModel()
## Define sets
m.FOOD = pyo.Set(initialize = foods)
m.NUTRIENTS = pyo.Set(initialize = nutrients)
## Define parameters
m.needs = pyo.Param(m.NUTRIENTS, initialize=nutrient_requirements)
m.food_info = pyo.Param(m.FOOD, m.NUTRIENTS, initialize=food_info)
m.price = pyo.Param(m.FOOD, initialize=price)
## Define variables
m.food_eaten = pyo.Var(m.FOOD,initialize=1.0,domain=pyo.NonNegativeReals)
## Define constraints
def diet_min_rule(m, n):
return sum(m.food_info[f,n] * m.food_eaten[f] for f in m.FOOD) >= m.needs[n]
m.diet_min = pyo.Constraint(m.NUTRIENTS, rule=diet_min_rule)
## Define objective
m.cost = pyo.Objective(expr=sum(m.food_eaten[f]*m.price[f] for f in m.FOOD), sense=pyo.minimize)
# Print model
m.pprint()
Activity
Check the Pyomo model. Specifically, are the input (parameter) data correct? Do the equations match our model written on paper?
# Specify the solver
solver = pyo.SolverFactory('ipopt')
# Solve
results = solver.solve(m, tee=True)
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...: 0
Number of nonzeros in inequality constraint Jacobian.: 10
Number of nonzeros in Lagrangian Hessian.............: 0
Total number of variables............................: 3
variables with only lower bounds: 3
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 4
inequality constraints with only lower bounds: 4
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 8.0000000e+00 6.00e+00 1.10e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 8.6444860e+00 5.02e+00 9.73e-01 -1.0 1.16e+00 - 2.18e-01 1.54e-01h 1
2 1.3016656e+01 0.00e+00 5.33e-01 -1.0 8.85e-01 - 6.53e-01 1.00e+00h 1
3 1.2884986e+01 0.00e+00 6.91e-02 -1.7 1.53e-01 - 7.53e-01 8.71e-01f 1
4 1.2512801e+01 0.00e+00 1.19e-01 -2.5 4.08e+00 - 1.17e-01 6.52e-01f 1
5 1.2513485e+01 0.00e+00 2.83e-08 -2.5 5.33e-02 - 1.00e+00 1.00e+00f 1
6 1.2500398e+01 0.00e+00 1.50e-09 -3.8 4.46e-02 - 1.00e+00 1.00e+00f 1
7 1.2500005e+01 0.00e+00 1.84e-11 -5.7 5.47e-04 - 1.00e+00 1.00e+00f 1
8 1.2500000e+01 0.00e+00 2.52e-14 -8.6 1.25e-05 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 8
(scaled) (unscaled)
Objective...............: 1.2499999882508366e+01 1.2499999882508366e+01
Dual infeasibility......: 2.5153648267429247e-14 2.5153648267429247e-14
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 2.5136445446423307e-09 2.5136445446423307e-09
Overall NLP error.......: 2.5136445446423307e-09 2.5136445446423307e-09
Number of objective function evaluations = 9
Number of objective gradient evaluations = 9
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 9
Number of equality constraint Jacobian evaluations = 0
Number of inequality constraint Jacobian evaluations = 9
Number of Lagrangian Hessian evaluations = 8
Total CPU secs in IPOPT (w/o function evaluations) = 0.004
Total CPU secs in NLP function evaluations = 0.000
EXIT: Optimal Solution Found.
Activity
Does your degree of freedom analysis match Ipopt?
Solution
Add Solution HERE
1.4.1.3. Analyze Results#
Let’s extract the solution.
Activity
Write Python code to extract and print the solution.
# Add your solution here
Units of takeaway eaten = 1.0
Units of vegetables eaten = 1.5
Units of bread eaten = 3.0
Click here to see the activity solution.
for f in m.FOOD:
print("Units of",f,"eaten =",round(m.food_eaten[f](),2))
Next, let us extract the KKT multipliers.
### Declare all suffixes
# https://pyomo.readthedocs.io/en/stable/pyomo_modeling_components/Suffixes.html#exporting-suffix-data
# Ipopt bound multipliers
m.ipopt_zL_out = pyo.Suffix(direction=pyo.Suffix.IMPORT)
m.ipopt_zU_out = pyo.Suffix(direction=pyo.Suffix.IMPORT)
# Ipopt constraint multipliers
m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT_EXPORT)
# Resolve the model
results = solver.solve(m, tee=True)
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...: 0
Number of nonzeros in inequality constraint Jacobian.: 10
Number of nonzeros in Lagrangian Hessian.............: 0
Total number of variables............................: 3
variables with only lower bounds: 3
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 0
Total number of inequality constraints...............: 4
inequality constraints with only lower bounds: 4
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.2500000e+01 0.00e+00 1.10e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 1.2689467e+01 0.00e+00 1.72e-02 -1.0 1.22e-01 - 9.83e-01 1.00e+00h 1
2 1.2509157e+01 0.00e+00 2.83e-08 -2.5 1.30e-01 - 1.00e+00 1.00e+00f 1
3 1.2500405e+01 0.00e+00 1.50e-09 -3.8 2.95e-02 - 1.00e+00 1.00e+00f 1
4 1.2500005e+01 0.00e+00 1.84e-11 -5.7 7.04e-04 - 1.00e+00 1.00e+00f 1
5 1.2500000e+01 0.00e+00 2.57e-14 -8.6 1.24e-05 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 5
(scaled) (unscaled)
Objective...............: 1.2499999882508405e+01 1.2499999882508405e+01
Dual infeasibility......: 2.5730637052582628e-14 2.5730637052582628e-14
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 2.5137374227288735e-09 2.5137374227288735e-09
Overall NLP error.......: 2.5137374227288735e-09 2.5137374227288735e-09
Number of objective function evaluations = 6
Number of objective gradient evaluations = 6
Number of equality constraint evaluations = 0
Number of inequality constraint evaluations = 6
Number of equality constraint Jacobian evaluations = 0
Number of inequality constraint Jacobian evaluations = 6
Number of Lagrangian Hessian evaluations = 5
Total CPU secs in IPOPT (w/o function evaluations) = 0.002
Total CPU secs in NLP function evaluations = 0.000
EXIT: Optimal Solution Found.
#Inspect dual variables for lower bound
m.ipopt_zL_out.display()
ipopt_zL_out : Direction=IMPORT, Datatype=FLOAT
Key : Value
food_eaten[bread] : 8.35310647704211e-10
food_eaten[takeaway] : 2.5068018506518723e-09
food_eaten[vegetables] : 1.6730627698360893e-09
#Inspect dual variables for upper bound
m.ipopt_zU_out.display()
ipopt_zU_out : Direction=IMPORT, Datatype=FLOAT
Key : Value
# Inspect duals
m.dual.display()
dual : Direction=IMPORT_EXPORT, Datatype=FLOAT
Key : Value
diet_min[C] : 8.368398117061357e-10
diet_min[F] : 1.78571428033602
diet_min[P] : 0.42857143096267053
diet_min[V] : 0.14285714142265304
1.4.2. Take Away Messages#
Linear programs are convex. We will learn this means all local optima are global optima.
There are specialize solves for linear programs, quadratic programs, and convex programs. In this class, we will focus on more general algorithms for (non)convex nonlinear programs including the algorithms used by the
ipoptsolver.