ENV590.05 - Econ of Modern Power Systems - Module 6 - Intro to LPs in Python with PYOMO

Learning outcomes:

  • Learn how to implement a simple LP in Python using PYOMO;
  • Learn hor to implement an energy storage system management for residential applications.

Simple Example

You have two machines A and B and you can produce two types of water heaters. Profit from producing type I is \$800 and type II \$600. The units produced are constrained by machine availability. You only have 60 hours in machine A and 48 in machine B.

  • Type I requires 4h in MA and 2h in MB
  • Type II requires 2h in MA and 4h in MB

Screen Shot 2021-09-23 at 10.23.03 PM.png

This is a small instance with 2 decision variables and 2 constraints. We will explore two ways to implement this model to get you familiar with the Pyomo environment.

In [ ]:
from google.colab import drive
Mounted at /content/drive
In [ ]:
import os
os.chdir('/content/drive/MyDrive/Colab Notebooks/')

Installing and Running Pyomo on Google Colab

To import/install a library that's not in Colaboratory by default, you can use !pip install. This needs to be done at the begining of you notebook. And you only need to run it once at the start of each Colab session.

In [ ]:
!pip install pyomo
Collecting pyomo
  Downloading Pyomo-6.1.2-cp37-cp37m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl (9.1 MB)
     |████████████████████████████████| 9.1 MB 5.2 MB/s 
Collecting ply
  Downloading ply-3.11-py2.py3-none-any.whl (49 kB)
     |████████████████████████████████| 49 kB 6.0 MB/s 
Installing collected packages: ply, pyomo
Successfully installed ply-3.11 pyomo-6.1.2

Installing optimization solver

Pyomo does not include any optimization solvers. Therefore, you will need to install third-party solvers to solve optimization models built with Pyomo. There are several solver options, for this class we will use glpk. We'll install glpk using apt-get.

In [ ]:
!apt-get install -y -qq glpk-utils
Selecting previously unselected package libsuitesparseconfig5:amd64.
(Reading database ... 155013 files and directories currently installed.)
Preparing to unpack .../libsuitesparseconfig5_1%3a5.1.2-2_amd64.deb ...
Unpacking libsuitesparseconfig5:amd64 (1:5.1.2-2) ...
Selecting previously unselected package libamd2:amd64.
Preparing to unpack .../libamd2_1%3a5.1.2-2_amd64.deb ...
Unpacking libamd2:amd64 (1:5.1.2-2) ...
Selecting previously unselected package libcolamd2:amd64.
Preparing to unpack .../libcolamd2_1%3a5.1.2-2_amd64.deb ...
Unpacking libcolamd2:amd64 (1:5.1.2-2) ...
Selecting previously unselected package libglpk40:amd64.
Preparing to unpack .../libglpk40_4.65-1_amd64.deb ...
Unpacking libglpk40:amd64 (4.65-1) ...
Selecting previously unselected package glpk-utils.
Preparing to unpack .../glpk-utils_4.65-1_amd64.deb ...
Unpacking glpk-utils (4.65-1) ...
Setting up libsuitesparseconfig5:amd64 (1:5.1.2-2) ...
Setting up libcolamd2:amd64 (1:5.1.2-2) ...
Setting up libamd2:amd64 (1:5.1.2-2) ...
Setting up libglpk40:amd64 (4.65-1) ...
Setting up glpk-utils (4.65-1) ...
Processing triggers for libc-bin (2.27-3ubuntu1.3) ...
/sbin/ldconfig.real: /usr/local/lib/python3.7/dist-packages/ideep4py/lib/libmkldnn.so.0 is not a symbolic link

Processing triggers for man-db (2.8.3-2ubuntu0.1) ...

Import Pyomo and solver

The first step in any Pyomo project is to import relevant components of the Pyomo library. This can be done with the following python statement 'from pyomo.environ import *'. \

We use the * symbol to elimate the need of using the expression pyomo.environ every time we need to use a pyomo function. \

You also need to load the solver. The Pyomo libary includes a SolverFactory() class used to specify a solver. Here we will use glpk which works for linear problems, put cbc from coin OR could be used for nonlinear applications.

In [ ]:
from pyomo.environ import *
#Import solver

Approach 1: Create Model using scalars only

The next thing you need to do is initialize the model. The are two options with Pyomo for model building: AbstractModel() and ConcreteModel(). The file “PyomoFundamentals.pdf” will explain the differences in details. In a nutshell concrete models are immediately constructed and data is presented when components like variables, constraints and objective function are created. The abstract model is useful when you will run the same optimization model for different data sets, or when you will import data that is stored in a excel spreadsheet. You can create an abstract model even before entering/importing the data set.

For our simple example, since we already have all the data we need, let’s build a concrete model. Let’s create the model and name it “model” using the following command.

In [ ]:
#Creating model
model = ConcreteModel()

The first components we will add to the model are the decision variables. We will use class Var() to specify the variable type using the domain argument. The number of units type I and II will be a continuous number greater than equal to zero.

In [ ]:
#Adding decision variables
model.x1 = Var(domain=NonNegativeReals)
model.x2 = Var(domain=NonNegativeReals)
#model.x2 = Var(domain=Reals,bounds=(0,None))  #alternative way to define variables

The domain argument could also be used to specify other type of variables like Real, Integers and Booleans. The class Var() has also an argument to set bounds on decision variables, for example Var(domain=Reals, bounds=(0,None)) is the same of defining domain equal NonNegativeReals. \

The second component we will add is objective function. Recall an objective function is an expression involving decision variables. We will store it in model.profit and use the class Objective() to define the expression and sense of the optimization model.

In [ ]:
#Adding objective function
model.profit = Objective(expr = 800*model.x1 + 600*model.x2, sense=maximize)

The third component we will add is the constraints using the class Constraint(). The expression is specified using the argument expr.

In [ ]:
#Adding constraints
model.machA = Constraint(expr = 4*model.x1 + 2*model.x2 <= 60)
model.machB = Constraint(expr = 2*model.x1 + 4*model.x2 <= 48)

We already specified the solver now all we need to do is call solve(). The solve() method attempts to solve the model using the specified solver. And you can add a few lines to print the results as below.

In [ ]:
#Solve Model
results = opt.solve(model)

#Print results
print('Profit = ', model.profit())
print('\nDecision Variables')
print('x1 = ', model.x1())
print('x2 = ', model.x2())
Profit =  13200.0

Decision Variables
x1 =  12.0
x2 =  6.0
In [ ]:
Out[ ]:
In [ ]:
print('Mach A = ', model.machA())
print('Mach B = ', model.machB())
Mach A =  60.0
Mach B =  48.0
In [ ]:
if results.solver.status == 'ok':
2 Var Declarations
    x1 : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :  12.0 :  None : False : False : NonNegativeReals
    x2 : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :   6.0 :  None : False : False : NonNegativeReals

1 Objective Declarations
    profit : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : maximize : 800*x1 + 600*x2

2 Constraint Declarations
    machA : Size=1, Index=None, Active=True
        Key  : Lower : Body        : Upper : Active
        None :  -Inf : 4*x1 + 2*x2 :  60.0 :   True
    machB : Size=1, Index=None, Active=True
        Key  : Lower : Body        : Upper : Active
        None :  -Inf : 2*x1 + 4*x2 :  48.0 :   True

5 Declarations: x1 x2 profit machA machB

Optional: How to get shadow prices?

Access to shadow prices values is similar to accessing decision variables values, except that shadow prices are not captured by default so additional specifications are needed before solving the model to signal that shadow prices are desired. The dual suffix is where you can access shadow prices.

In [ ]:
#calculating shadow prices

Then all you need to do is solve the model and print results plus shadow prices.

In [ ]:
#Solve Model
results = opt.solve(model)

print("Shadow Prices")
Shadow Prices
dual : Direction=Suffix.IMPORT, Datatype=Suffix.FLOAT
    Key   : Value
    machA : 166.666666666667
    machB : 66.6666666666667

Approach 2: Create model using data in matrix/vector format and numpy

So far we used scalar modeling components, model.x1 = Var() and model.x2 = Var(), and each constraint was added as a separate line in the model. Now I want to show you how you could create the same model using numpy arrays to enter the data and use indexed variables and constraints. This will come in handy when you create larger instances of LP models. And that will be the approach needed for Assignment 2.

Recall the original representation with scalar, let’s build a matrix representation. Screen Shot 2021-09-23 at 10.05.11 PM.png

Note that the matrix representation can be used for instances with more machines and more water heater types. The size of matrix A and vector c, b and x will depend on the sets. \

$I=\{Type I,Type II\}$ \


The indices are represented by the range() statement. We have two sets, one for the chemical solution types I and II and one for machines. Entering the numpy arrays is straight forward. The only part that is not familiar to you is the sum() used to define the expression for constraints and objective function. You can see this $sum(A[j,i]*x[i] \ \ for \ \ i \ \ in \ \ I)$ as equivalent to the expression \


In [ ]:
#Indexed variables and constraints
import numpy as np
#LP matrix and vector - standard form
A = np.array([[4, 2], [2, 4]])
b = np.array([60, 48])
c = np.array([800,600])

I=range(2)  #set of water heater type I and II
J=range(2)  #set of machines A and B

model = ConcreteModel()

#Index decision variables by water heater type set

#Adding constraints
model.constraints = ConstraintList()
for j in J:
    model.constraints.add(sum(A[j,i]*model.x[i] for i in I) <= b[j])

#Adding objective
model.profit = Objective(expr = sum(c[i]*model.x[i] for i in I), sense=maximize)

#Solvign model

#Print results
print('Profit = ', model.profit())
print('\nDecision Variables')
for j in J:
    print('x',(j+1),' = ', model.x[j].value,sep='')

Profit =  13200.0

Decision Variables
x1 = 12.0
x2 = 6.0
In [ ]:
#Create index sets
#model.I = Set(initialize=range(2))  
#model.J = Set(initialize=range(2))  

Approach 3 - Optional - Create model using Set() and Param() objects

In [ ]:
#using sets and parameter

model.M=Set(initialize=['MA','MB'])  #set of machines
model.P=Set(initialize=['TypeI','TypeII']) #set of solution types

In [ ]:
#add dec variables

#add obj func
model.profit=Objective(sense=maximize,exp=sum(model.c[p]*model.X[p] for p in model.P)

#add const
def mach_hours(model,m):
    return sum(model.a[m,p]*model.X[p] for p in model.P) <= model.H[m]
In [ ]:
#Solve model

#Print results
print("Profit =",model.profit())
print("Decision Variables")
for p in model.P: