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.

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.

In [1]:

```
from google.colab import drive
drive.mount('/content/drive')
```

Mounted at /content/drive

In [2]:

```
import os
os.chdir('/content/drive/MyDrive/Colab Notebooks/')
```

In [3]:

```
!pip install pyomo
```

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 [4]:

```
!apt-get install -y -qq glpk-utils
```

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 [5]:

```
from pyomo.environ import *
#Import solver
opt=SolverFactory('glpk')
```

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 [6]:

```
#Creating model
model = ConcreteModel()
```

In [7]:

```
#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 [8]:

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

In [9]:

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

In [10]:

```
#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 [11]:

```
model.x1()
```

Out[11]:

12.0

In [12]:

```
print('\nConstraints')
print('Mach A = ', model.machA())
print('Mach B = ', model.machB())
```

Constraints Mach A = 60.0 Mach B = 48.0

In [13]:

```
#results.write()
if results.solver.status == 'ok':
model.pprint()
```

In [14]:

```
#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.

In [15]:

```
#Solve Model
results = opt.solve(model)
print("Shadow Prices")
model.dual.pprint()
```

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$.

In [16]:

```
#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='')
```

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))
```

In [24]:

```
#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})
```

In [25]:

```
#add dec variables
model.X=Var(model.P,domain=NonNegativeReals)
#add obj func
model.profit=Objective(sense=maximize,expr=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)
```

In [26]:

```
#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)
```

Profit = 13200.0 Decision Variables X[TypeI] 12.0 X[TypeII] 6.0