None

This notebook contains material from cbe30338-2021; content is available on Github.

< 5.6 Design of a Cold Weather Fuel | Contents | Tag Index | 5.99 Pyomo Examples >

5.98 Linear Production Model in Pyomo

Pyomo is a state-of-the-art language for solving optimization problems embedded within Python. Using Pyomo, a user can describe optimization model by specifying decision variables, constraints, and an optimization objective. Pyomo includes a rich set of features to enable modeling of complex systems, specifying a solver, and displaying the solution.

5.98.1 Installation

5.98.1.1 Anaconda

These instructions assume you have previously installed a suitable Python development environment, in particular the Anaconda package for Python 3.x.

The following commands will install Pyomo and extra files plus the glpk (MILP) and ipopt (nonlinear) solvers. These commands should be executed one at a time from the terminal window on MacOS, or the Anaconda command prompt on Windows 10.

conda install -c conda-forge pyomo
conda install -c conda-forge pyomo.extras
conda install -c conda-forge glpk
conda install -c conda-forge ipopt

Optionally, you may also wish to install the COIN-OR cbc solvers.

conda install -c conda-forge coincbc

This installation provides a capable open-source optimization suite with multiple solvers well suited a wide range of process applications. The solvers are

• glpk. "GNU Linear Programming Kit" is a software package written in highly portable C for the solution of mixed integer linear programming and related problems.
• ipopt. "Interior Point Optimizer" for large scale nonlinear optimization of continuous systems.
• cbc. "COIN-OR Branch and Cut" is a mixed integer linear programming solver written in C++. It generally solves the same problems as glpk, but can run multiple threads, and is often faster and more robust.

This suite can be further augmented by installing additional solvers from open-source and commercial sources.

NOTE: If command window procedure fails on your laptop, you might try running these commands in a Jupyter notebook. In that case, each command needs to be prefixed with the shell escape !, and followed by the -y option to handle the y/n interaction. For example, the basic pyomo install command would read

!conda install -c conda-forge pyomo -y

To provide timely feedback, each install command should be executed in a separate cell.

5.98.1.2 Google Colab

%%capture
!pip install -q pyomo
!apt-get install -y -qq glpk-utils

5.98.2 Algebraic Modeling Languages

One of the difficult aspects of applying optimization methods to large scale applications is creating and maintaining the underlying model. An application may include tens, hundreds, even thousdands of variables and constraints, and may incorporate data extracted in real time from sensor networks and corporate data bases. The resulting models can be very complex, and essentially impossible to create and maintain without use of software tools.

Algebraic Modeling Languages (AMLs) were developed to assist with the creation and maintainence of optimization models. GAMS (General Algebraic Modeling System, https://www.gams.com/), first proposed in 1976, was among the first and still widely used. Other notable examples include AIMMS, AMPL, and FICO XPRESS.

In recent years, modeling for optimization has become more tightly integrated with scripting languages commonly used in science and engineering. Open-source examples include YALMIP and CVX which are tightly integrated with Matlab, JuMP which works with Julia, and a variety of systems integrated with Python.

Of the Python options, the open-source Pyomo is perhaps the most ambitious and certainly the most aligned with the needs of process systems engineering.

5.98.3 Key Concepts in Pyomo

A typical Pyomo application involves the creation of at least two Pyomo objects from the following classes:

• ConcreteModel() A python object representing the optimization problem to be solved.
• SolverFactory() A python object representing the mathematical progamming software to be used for calculating a solution.

There are a number of of open-source and commercial solvers that can be used with Pyomo. This simple structure allows one to test and find solvers suited to a particular application.

A model, in turn, is composed of additional objects used to specify a problem. A minimal set of classes is needed to create useful models. These classes are:

• Var() Objects representing variables determined in the course of solving a particular problem.
• Objective() An object denoting the problem objective function that is to be either minimized or maximized.
• Contraint() Objects representing problem constraints.

The following cell presents a typical example.

5.98.4 Example: Linear Production Model with Constraints

The mathematical formulation of a simple linear production model is given by

\begin{align} \max_{x,y \geq 0} &\ 40\ x + 30\ y & \mbox{objective}\\ \mbox{subject to:}\qquad \\ x & \leq 40 & \mbox{demand constraint} \\ x + y & \leq 80 & \mbox{labor A constraint} \\ 2x + y & \leq 100 & \mbox{labor B constraint} \end{align}

where $x$ and $y$ are the rates of production (in units per week) for two products.

5.98.4.1 Step 1. Import Pyomo.

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 *

5.98.4.2 Step 2. Create the model object.

Pyomo provides two distinct methods to create models. Here we use ConcreteModel() to create a model instance which is appropriate when all of the data needed to complete the model is avaiable at the current time.

model = ConcreteModel()

The Python variable model stores the model instance. This could be any valid Python variable name.

5.98.4.3 Step 3. Add the Decision Variables, Objective, and Constraints to the model object.

The first major component of a Pyomo model are decision variables which are added as fields to model. In the case we name the fields model.x and model.y corresponding to $x$ and $y$ in the process model. The Python class Var() is used to specify these as real numbers that must be greater than or equal to zero.

model.x = Var(domain=NonNegativeReals)
model.y = Var(domain=NonNegativeReals)

As we will see in other examples, the domain can specify other types of decision variables including reals, integers, and booleans.

The objective is specified as an algebraic expression involving the decision variables. Here we store the objective in model.profit, and use the optional keyword sense to specify a maximization problem.

model.profit = Objective(expr = 40*model.x + 30*model.y, sense=maximize)

Constraints are added as fields to the model, each constraint created using the Constraint() class. The constraints are specified using the expr keywork in the form of linear algebraic expressions of the decision variables.

model.demand = Constraint(expr = model.x <= 40)
model.laborA = Constraint(expr = model.x + model.y <= 80)
model.laborB = Constraint(expr = 2*model.x + model.y <= 100)

5.98.4.4 Step 4. Create a solver object and solve.

The Pyomo libary includes a SolverFactory() class used to specify a solver. In this case, because the problem is a linear programming problem, we use the glpk solver.

results = SolverFactory('glpk').solve(model)
results.write()
if results.solver.status == 'ok':
    model.pprint()

The solve() method attempts to solve the model using the specified solver, and returns results which can be inspected to see if, in fact, a solution was found. If a solution is found, then model will have a pretty-print method pprint() that creates a summary of the problem solution.

5.98.4.5 Step 5. Display the solution.

Most applications will require access to specific components of the model. If a solution is found, the model will include methods with the same name as the fields created above, and which can be used to access solution values.

print('Profit = ', model.profit())

print('\nDecision Variables')
print('x = ', model.x())
print('y = ', model.y())

print('\nConstraints')
print('Demand  = ', model.demand())
print('Labor A = ', model.laborA())
print('Labor B = ', model.laborB())

from pyomo.environ import *

# create a model
model = ConcreteModel()

# declare decision variables
model.x = Var(domain=NonNegativeReals)
model.y = Var(domain=NonNegativeReals)

# declare objective
model.profit = Objective(expr = 40*model.x + 30*model.y, sense=maximize)

# declare constraints
model.demand = Constraint(expr = model.x <= 40)
model.laborA = Constraint(expr = model.x + model.y <= 80)
model.laborB = Constraint(expr = 2*model.x + model.y <= 100)

# solve
results = SolverFactory('glpk').solve(model)
results.write()
if results.solver.status:
    model.pprint()

# display solution
print('\nProfit = ', model.profit())

print('\nDecision Variables')
print('x = ', model.x())
print('y = ', model.y())

print('\nConstraints')
print('Demand  = ', model.demand())
print('Labor A = ', model.laborA())
print('Labor B = ', model.laborB())

# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 2600.0
  Upper bound: 2600.0
  Number of objectives: 1
  Number of constraints: 4
  Number of variables: 3
  Number of nonzeros: 6
  Sense: maximize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 0
      Number of created subproblems: 0
  Error rc: 0
  Time: 0.01840686798095703
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 0
  number of solutions displayed: 0
2 Var Declarations
    x : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :  20.0 :  None : False : False : NonNegativeReals
    y : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 :  60.0 :  None : False : False : NonNegativeReals

1 Objective Declarations
    profit : Size=1, Index=None, Active=True
        Key  : Active : Sense    : Expression
        None :   True : maximize : 40*x + 30*y

3 Constraint Declarations
    demand : Size=1, Index=None, Active=True
        Key  : Lower : Body : Upper : Active
        None :  -Inf :    x :  40.0 :   True
    laborA : Size=1, Index=None, Active=True
        Key  : Lower : Body  : Upper : Active
        None :  -Inf : x + y :  80.0 :   True
    laborB : Size=1, Index=None, Active=True
        Key  : Lower : Body    : Upper : Active
        None :  -Inf : 2*x + y : 100.0 :   True

6 Declarations: x y profit demand laborA laborB

Profit =  2600.0

Decision Variables
x =  20.0
y =  60.0

Constraints
Demand  =  20.0
Labor A =  80.0
Labor B =  100.0

5.98.5 Python Lists, Sets, Dictionaries, and Iterators

Pyomo is integrated with the Python language, and inherits significant functionality by leveraging the basic data structures of Python. To make the best use of Pyomo, it is important to understand the basics of Python lists, sets, and dictionaries.

5.98.5.1 Data in Matrix/Vector Format

The example above used scalar modeling components, model.x = Var() and model.y = Var(), and each constraint was added as a separate line in the model. This is fine for small problems with a just a few unknowns, but becomes impractical for larger problems.

Here we repeat the example above, this time using numpy arrays to enter the data. An indexed variable model.x represents the unknowns, and the constraints are indexed as well. The indices are represented by the Python range() statement.

from pyomo.environ import *
import numpy as np

# enter data as numpy arrays
A = np.array([[1, 0], [1, 1],[2,1]])
b = np.array([40, 80,100])
c = np.array([40,30])

# set of row indices
I = range(len(A))

# set of column indices
J = range(len(A.T))

# create a model instance
model = ConcreteModel()

# create x and y variables in the model
model.x = Var(J)

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

# add a model objective
model.objective = Objective(expr = sum(c[j]*model.x[j] for j in J), sense=maximize)

# create a solver
solver = SolverFactory('glpk')

# solve
solver.solve(model)

# print solutions
for j in J:
    print("x[", j, "] =", model.x[j].value)

model.constraints.display()

x[ 0 ] = 20.0
x[ 1 ] = 60.0
constraints : Size=3
    Key : Lower : Body  : Upper
      1 :  None :  20.0 :  40.0
      2 :  None :  80.0 :  80.0
      3 :  None : 100.0 : 100.0

< 5.6 Design of a Cold Weather Fuel | Contents | Tag Index | 5.99 Pyomo Examples >