Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

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
from pyomo.environ import units as u
from pyomo.util.check_units import assert_units_consistent

# Pyomo ships the SI and US customary units, but not money, servings of food, or
# "units of protein". Declare the three this problem needs.
u.load_definitions_from_strings(
    ["USD = [currency]", "serving = [serving]", "nutrient = [nutrient]"]
)

Linear Programs: Student Diet Example

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 the first 10 rows of data
food_options.head(10)
Loading...

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

Propose an Optimization Model

Sets

Click to expand
F={take away, veggies, bread}\mathcal{F} = \{\text{take away, veggies, bread}\}
(1)
N={P,C,F,V}\mathcal{N} = \{P, C, F, V\}
(2)

Parameters

Click to expand

πi\pi_i: price for food option i∈Fi \in \mathcal{F}

rjr_j: nutrition requirement for nutrient j∈Nj \in \mathcal{N}

DijD_{ij}: nutrition info for food i∈Fi \in \mathcal{F} and nutrient j∈Nj \in \mathcal{N}

Variables

Click to expand

xix_i: amount of food i∈Fi \in \mathcal{F} eaten/purchased

Objective and Constraints

Click to expand

Minimize: ∑i∈Fxi⋅πi(cost) \sum_{i \in \mathcal{F}} x_i \cdot \pi_i \quad \text{(cost)}

Subject to: ∑i∈FDij⋅xi≥rj∀j∈N(be healthy) \sum_{i \in \mathcal{F}} D_{ij} \cdot x_i \geq r_j \quad \forall j \in \mathcal{N} \quad \text{(be healthy)}

Bounds: xi≥0∀i∈F 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.

min⁡x∑i∈Fxi⋅πi(cost)s.t.∑i∈FDij⋅xi≥rj,∀j∈N(be healthy)xi≥0,∀i∈F(non-negative bounds)\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*}
(3)

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.

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!

Aside: two ways to write a constraint rule

The “be healthy” constraint is our first indexed constraint: one inequality for every nutrient j∈Nj \in \mathcal{N}. Pyomo builds an indexed constraint from a rule --- a function of the model and the index --- which it calls once per index. There are two ways to hand Pyomo that rule, and they build exactly the same model.

The older form defines the function, then attaches it with the rule= keyword:

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)

The decorator form attaches the function as it is defined:

@m.Constraint(m.NUTRIENTS)
def diet_min(b, n):
    return sum(b.food_info[f, n] * b.food_eaten[f] for f in b.FOOD) >= b.needs[n]

This course uses the decorator from here on, and so does our Pyomo style guide. Two reasons: the name is written once instead of three times (function name, _rule suffix, component name), so the component name cannot drift away from the function that defines it; and the lines are shorter, which matters when the model goes into a handout. Note the first argument is the block Pyomo is building, so we call it b rather than m.

You still need to be able to read the rule= form. It is what the Pyomo textbook uses and what you will meet in most existing research code. When you see rule=some_function, translate it into the decorator in your head. This is the only place in the course notes where we write a constraint both ways.

# Add your solution here
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, units=u.nutrient)
m.food_info = pyo.Param(
    m.FOOD, m.NUTRIENTS, initialize=food_info, units=u.nutrient / u.serving
)
m.price = pyo.Param(m.FOOD, initialize=price, units=u.USD / u.serving)

## Define variables
m.food_eaten = pyo.Var(
    m.FOOD, initialize=1.0, domain=pyo.NonNegativeReals, units=u.serving
)


## Define constraints
# Be healthy: sum_i D_ij * x_i >= r_j [nutrient]
@m.Constraint(m.NUTRIENTS)
def diet_min(b, n):
    return sum(b.food_info[f, n] * b.food_eaten[f] for f in b.FOOD) >= b.needs[n]


## Define objective
# Total cost of the diet, sum_i pi_i * x_i [USD]
@m.Objective(sense=pyo.minimize)
def cost(b):
    return sum(b.food_eaten[f] * b.price[f] for f in b.FOOD)


# Print model
m.pprint()

Check the units

Every Var and Param above carries a units= keyword. That is a declaration, and a declaration by itself proves nothing --- Pyomo will happily build a model that adds servings to dollars. assert_units_consistent walks every constraint and the objective and raises UnitsError if the two sides disagree.

Declare and then check. Units nobody checks only look verified.

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

print("Units are consistent.")
Units are consistent.

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)
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...:        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.59e-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.5871698107970097e-14    2.5871698107970097e-14
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.5136445446423332e-09    2.5136445446423332e-09
Overall NLP error.......:   2.5136445446423332e-09    2.5136445446423332e-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.001
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

Analyze Results

Let’s extract the solution.

Activity

Write Python code to extract and print the solution.
# Add your solution here
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)
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...:        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.5137374227288727e-09    2.5137374227288727e-09
Overall NLP error.......:   2.5137374227288727e-09    2.5137374227288727e-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.000
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.368398117055004e-10
    diet_min[F] :      1.78571428033602
    diet_min[P] :   0.42857143096267036
    diet_min[V] :   0.14285714142265307

Take Away Messages

  • Linear programs are convex. We will learn this means all local optima are global optima.

  • There are specialized solvers 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 ipopt solver.

Optimization Modeling in Pyomo
Your First Optimization Problem
Optimization Modeling in Pyomo
Continuous Optimization: Nonlinear Programming