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
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.
from google.colab import drive
drive.mount('/content/drive')
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.
!pip install pyomo
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.
!apt-get install -y -qq glpk-utils
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.
from pyomo.environ import *
#Import solver
opt=SolverFactory('glpk')
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.
#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.
#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.
#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.
#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.
#Solve Model
results = opt.solve(model)
#Print results
print('Profit = ', model.profit())
print('\nDecision Variables')
print('x1 = ', model.x1())
print('x2 = ', model.x2())
model.x1()
print('\nConstraints')
print('Mach A = ', model.machA())
print('Mach B = ', model.machB())
#results.write()
if results.solver.status == 'ok':
model.pprint()
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.
#calculating shadow prices
model.dual=Suffix(direction=Suffix.IMPORT)
Then all you need to do is solve the model and print results plus shadow prices.
#Solve Model
results = opt.solve(model)
print("Shadow Prices")
model.dual.pprint()
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.
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\}$ \
$J=\{MA,MB\}$
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 \
$\sum_{i∈I}a_{j,i}*x_i$.
#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
model.x=Var(I)
#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
results=opt.solve(model)
#Print results
print('Profit = ', model.profit())
print('\nDecision Variables')
for j in J:
print('x',(j+1),' = ', model.x[j].value,sep='')
#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¶
#using sets and parameter
model=ConcreteModel()
#Sets
model.M=Set(initialize=['MA','MB']) #set of machines
model.P=Set(initialize=['TypeI','TypeII']) #set of solution types
#Parameters
model.c=Param(model.P,initialize={'TypeI':800,'TypeII':600})
model.H=Param(model.M,initialize={'MA':60,'MB':48})
model.a=Param(model.M,model.P,initialize={
('MA','TypeI'):4,
('MA','TypeII'):2,
('MB','TypeI'):2,
('MB','TypeII'):4})
#add dec variables
model.X=Var(model.P,domain=NonNegativeReals)
#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]
model.mach=Constraint(model.M,rule=mach_hours)
#Solve model
results=opt.solve(model)
#Print results
print("Profit =",model.profit())
print("Decision Variables")
for p in model.P:
print(model.X[p],model.X[p].value)