Pyomo basics

7. Pyomo basics#

Note

This chapter is a Jupyter notebook. You can download and run it yourself (click on the “” icon in the toolbar above and choose “.ipynb”).

Pyomo extends Python to provide a mathematical modelling language with which you can formulate optimisation problems. Like Python itself, it is free and open-source.

At any point, you can refer to the official Pyomo documentation to learn more about the many advanced features of Pyomo that we do not cover here.

Before using Pyomo, it needs to be imported. By saying import as pyo we are defining pyo as a shorthand under which to access everything that Pyomo provides. This is best practice: it means that everything coming from Python itself or from other packages is cleanly separated from Pyomo - anytime we need Pyomo, we access it through pyo, as you’ll see below.

import pyomo.environ as pyo

Next, we define a variable to hold the model. By convention, we call this model. The model itself is created by calling ConcreteModel() (if you know about object-oriented programming - this creates an object which is a member of the ConcreteModel class).

Later on, we will fill this ConcreteModel with data such as variables and constraints.

model = pyo.ConcreteModel(name="Economic dispatch")

Giving the model a name is optional but leads to nicer output later on.

Remember, an optimisation problem consists of:

variables
the objective function, and
the constraints.

We now set these up in turn, using the economic dispatch example from the first lecture.

7.1. Variables#

Variables, are defined with Var objects. We “attach” them to the model as follows: model.variable_name = .... Here, variable_name is the name of the variable inside the model, and must therefore be unique to the model. For example:

model.coal_power = pyo.Var(within=pyo.NonNegativeReals)
model.gas_power = pyo.Var(within=pyo.NonNegativeReals)

within= is optional and specifies the domain within which the variable can vary.

We could also call the variables p1 and p2 as in the lecture. By using more expressive names we are making our model easier to understand.

Commonly used domains for within include:

pyo.NonNegativeReals: any real number greater than or equal to zero
pyo.Reals: any real number, including negative ones
pyo.Binary: a variable that can be either 0 or 1
pyo.Integer: a variable that can be any integer number (…, -2, -1, 0, 1, 2, …)

By using pyo.NonNegativeReals, we are making our life a little easier as we do not need to add extra “greater than or equal to zero” constraints on those variables.

Note

You can define a large number of variables by first defining sets, i.e. dimensions, over which the variables are valid. For example, you might want to have a variable called hourly_power_generation and you might want to model power generation over a longer period of time - say a year, which would be 8760 hours. You can first define a set of timesteps as a list, say [0, 1, ..., 8759], and then tell Pyomo that your variable can take on a value for each timestep by defining the variable over the set.

In our example, instead of defining two separate variables for coal_power and gas_power we could do it as follows:

power_plant_types = ["coal", "gas"]
model.power_generation = pyo.Var(power_plant_types, within=pyo.NonNegativeReals)

We can then access the two “component variables” as model.power_generation["coal"] and model.power_generation["gas"].

In the future, we could expand the set of power_plant_types by simply expanding the list, say, to ["coal", "gas", "wind", "pv", ...].

The Pyomo documentation provides more detail on this functionality, including the use of pyo.Set objects which provide advanced functionality when dealing with sets.

7.2. Objective function#

Just like we did for the variables, we create a pyo.Objective and attach it to our model. We’re calling it generation_cost but any name is possible (for example, you might want to always call the objective the same - obj or objective).

We give the objective two arguments:

expr is the mathematical expression of the objective function
sense is used to specify whether we are looking for a maximum (pyo.maximize) or a minimum (pyo.minimize) of this function. If sense is omitted, minimisation is assumed.

model.generation_cost = pyo.Objective(
    expr = 3 * model.coal_power + 4 * model.gas_power,
    sense = pyo.minimize
)

Note

An alternative way to formulate the objective if we use sets (see the box above) would look something like this:

prices = {"coal": 3, "gas": 4}

model.generation_cost = pyo.Objective(
    expr = sum(prices[i] * model.power_generation[i] for i in power_plant_types),
    sense = pyo.minimize
)

It makes use of additional features of Python, like the sum() function, and a so-called “comprehension” that iterates over all elements in the set (for i in power_plant_types). With this kind of approach, you can easily build a model with thousands of variables, which can easily be dealt with in mostly automated ways, like here summing up in the objective function.

7.3. Constraints#

Finally, we add the constraints, again giving each constraint an expr which is its mathematical formulation:

model.coal_min = pyo.Constraint(expr=model.coal_power >= 50)
model.coal_max = pyo.Constraint(expr=model.coal_power <= 300)
model.gas_min = pyo.Constraint(expr=model.gas_power >= 100)
model.gas_max = pyo.Constraint(expr=model.gas_power <= 400)
model.demand = pyo.Constraint(expr=model.gas_power + model.coal_power == 500)

7.4. Solving the model#

By default, Pyomo captures only a limited amount of information from the solver. Because we later want to look at dual variables, we need to specify this before solving with the following statement (what exactly the Suffix, direction, etc are, is not relevant here - simply always use this line of code to ensure that you can later access dual variables):

model.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)

Next, in order to solve the model, we first need to let Pyomo make the connection to a solver - a separate program that knows how to solve a specific kind of optimisation problem. In our case, we are solving a linear (LP) problem. Gurobi is a commercial solver that uses the Simplex algorithm to solve such problems, and is free to use for problems with up to 2000 variables and constraints, so we can use it here.

solver = pyo.SolverFactory('gurobi')

Note

To use Gurobi, you need to ensure that it is installed. If you follow the reader’s installation instructions, your setup will include Gurobi. For more details on solvers and alternative solver choices, see the “Solver choice” section below.

Finally, we can pass our model to the solver, which does its thing and then passes back the solution to Pyomo (or tells us if something went wrong and no optimal solution was found):

solver.solve(model)

{'Problem': [{'Name': 'x1', 'Lower bound': 1700.0, 'Upper bound': 1700.0, 'Number of objectives': 1, 'Number of constraints': 5, 'Number of variables': 2, 'Number of binary variables': 0, 'Number of integer variables': 0, 'Number of continuous variables': 2, 'Number of nonzeros': 6, 'Sense': 'minimize'}], 'Solver': [{'Status': 'ok', 'Return code': 0, 'Message': 'Model was solved to optimality (subject to tolerances), and an optimal solution is available.', 'Termination condition': 'optimal', 'Termination message': 'Model was solved to optimality (subject to tolerances), and an optimal solution is available.', 'Wall time': 0.001322031021118164, 'Error rc': 0}], 'Solution': [OrderedDict({'number of solutions': 0, 'number of solutions displayed': 0})]}

You can see that solver.solve(model) returns some summary information about the model solution. In this case, you can see 'Termination condition': 'optimal' – that’s good, we found an optimal solution.

If we do not want to display this information but capture it for later use, we could assign it to a variable, i.e.:

result_object = solver.solve(model)

If you want to see the output from the solver while it is solving the model, you can set the tee argument to True - in the above case, you would see the output from Gurobi, giving more detail on what it’s doing while solving the problem (e.g. what exact algorithm it is using, etc):

solver.solve(model, tee=True)

Now that the model is solved, we can display detailed information about it, including the state of the variables, the objective, and the constraints, with model.display():

model.display()

Model Economic dispatch

  Variables:
    coal_power : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 : 300.0 :  None : False : False : NonNegativeReals
    gas_power : Size=1, Index=None
        Key  : Lower : Value : Upper : Fixed : Stale : Domain
        None :     0 : 200.0 :  None : False : False : NonNegativeReals

  Objectives:
    generation_cost : Size=1, Index=None, Active=True
        Key  : Active : Value
        None :   True : 1700.0

  Constraints:
    coal_min : Size=1
        Key  : Lower : Body  : Upper
        None :  50.0 : 300.0 :  None
    coal_max : Size=1
        Key  : Lower : Body  : Upper
        None :  None : 300.0 : 300.0
    gas_min : Size=1
        Key  : Lower : Body  : Upper
        None : 100.0 : 200.0 :  None
    gas_max : Size=1
        Key  : Lower : Body  : Upper
        None :  None : 200.0 : 400.0
    demand : Size=1
        Key  : Lower : Body  : Upper
        None : 500.0 : 500.0 : 500.0

We can use the Python print function to set up a nicer overview of the results we care about most, for example:

print(f"Generation costs (objective) = {model.generation_cost()} EUR")
print(f"Gas plant generation = {model.gas_power()} MW")
print(f"Coal plant generation = {model.coal_power()} MW")

Generation costs (objective) = 1700.0 EUR
Gas plant generation = 200.0 MW
Coal plant generation = 300.0 MW

7.5. Dual variables#

Recall that constraints in the primal problem are associated with variables in the dual problem.

The solver solves the dual problem for us at the same time as solving the primal problem. Because we set up the model to extract this information earlier (model.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)), we can now get the shadow prices of the primal problem by looking up the optimal values of the variables of the dual problem like so:

model.dual[model.my_constraint_name]

Let’s print out the shadow price of some of our constraints:

print(f"Shadow price of `demand` = {model.dual[model.demand]}")
print(f"Shadow price of `coal_max` = {model.dual[model.coal_max]}")
print(f"Shadow price of `coal_min` = {model.dual[model.coal_min]}")

Shadow price of `demand` = 4.0
Shadow price of `coal_max` = -1.0
Shadow price of `coal_min` = 0.0

7.6. An easier way to access results#

Note

To make things easier for you when getting started with Pyomo, I have written a function for you to print a table with summary information about a solved model inside the notebook environment.

This function is part of a custom package that you need to import first by running from optimutils import summarise_results. Usually we import everything we will need at the very start of the notebook, but here, we only import it below, before using it the first time. If you are knowledgeable about Python and interested in the details, you can find the code here.

from optimutils import summarise_results

Then, you can run the function summarise_results, giving it the name of your model, in this case, simply model, to get an overview of the objective function value, the variables and their value in the solution, as well as the constraints, including the shadow price extracted from the solution of the dual problem:

summarise_results(model)

Objective
	Name	Value
0	generation_cost	1700.000000

Variables
	Name	Value
0	coal_power	300.000000
1	gas_power	200.000000

Constraints
	Name	Expression	Value	Shadow price	Binding
0	coal_min	50 <= coal_power	300.000000	0.000000	False
1	coal_max	coal_power <= 300	300.000000	-1.000000	True
2	gas_min	100 <= gas_power	200.000000	0.000000	False
3	gas_max	gas_power <= 400	200.000000	0.000000	False
4	demand	gas_power + coal_power == 500	500.000000	4.000000	True

7.7. Advanced topics#

7.7.1. Modifying and re-running a model, then saving results#

Note

The following uses Python to a larger extent than the basic overview of Pyomo above. If you are comfortable with using some Python, comfortable with basic programming concepts from other programming languages, or willing to experiment a little bit, then you will be able to exploit the full power of automating tasks to generate a larger set of results quickly by following the recipes below.

Let’s set up the exact same model again (“overwriting” the old model variable with a new version of the model). The one change is that we replace the value 4 for gas cost in the objective function with a new kind of object - a mutable (changeable) parameter represented by pyo.Param(mutable=True). We call this parameters gas_cost and attach it to the model. The rest of the model setup is exactly as above, but all combined in one notebook cell.

Parameters are fixed while the model is optimised - unlike variables that the solver algorithm can change to find an optimal solution.

However, by creating a mutable parameter we are telling Pyomo that we would like to change this parameter’s value in between several optimisation runs. This does not change anything about the mathematical formulation of the problem: mathematically, a mutable parameter is still fixed in the optimisation problem. However, as you will see below, it makes our lives easier, because we can change this one parameter and re-run the optimisation without creating a whole new model. This is helpful, for example, to rapidly scan a whole range of values for this parameter and see how it affects the solution.

model = pyo.ConcreteModel(name="Economic dispatch, version 2")

##
# This is the key bit: we define a `pyo.Param` which is mutable,
# allowing us to set (and change) its value later
##
model.gas_cost = pyo.Param(mutable=True)

model.coal_power = pyo.Var(within=pyo.NonNegativeReals)
model.gas_power = pyo.Var(within=pyo.NonNegativeReals)

model.generation_cost = pyo.Objective(
    expr = 3 * model.coal_power + model.gas_cost * model.gas_power,
    sense = pyo.minimize
)

model.coal_min = pyo.Constraint(expr=model.coal_power >= 50)
model.coal_max = pyo.Constraint(expr=model.coal_power <= 300)
model.gas_min = pyo.Constraint(expr=model.gas_power >= 100)
model.gas_max = pyo.Constraint(expr=model.gas_power <= 400)
model.demand = pyo.Constraint(expr=model.gas_power + model.coal_power == 500)

model.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)

Now that we have a parameter in the model, we can solve the model several times with different values for this parameter. First, we define a list of values we want to explore:

gas_costs = [1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0]

Next, we are going to iterate over each value in this list, set the gas_cost parameter to the current value from the list, solve the model, and store three results we care about - the generation cost (objective function value), gas power and coal power generation, into a list that we call result_accumulator - since it accumulates results (this happens with the Python functionality to append a new entry to an existing list: result_accumulator.append(...)).

result_accumulator = []

for i in gas_costs:
    model.gas_cost = i
    solver.solve(model)
    result_accumulator.append([model.generation_cost(), model.gas_power(), model.coal_power()])

Finally - we are going to turn the accumulated results into a DataFrame - a data structure provided by the pandas package which represents a table of data that is easy to further process, analyse, or save to a file. For that, we import the pandas package and make it available as pd, then create a DataFrame from the results, giving names for the columns, and specifying that we want the index (first) column to be the gas_costs which we iterated over to build these results).

import pandas as pd

results = pd.DataFrame(result_accumulator, columns=["generation cost", "gas power", "coal power"], index=gas_costs)

We can look at this table:

results

	generation cost	gas power	coal power
1.0	700.0	400.0	100.0
1.5	900.0	400.0	100.0
2.0	1100.0	400.0	100.0
2.5	1300.0	400.0	100.0
3.0	1500.0	400.0	100.0
3.5	1600.0	200.0	300.0
4.0	1700.0	200.0	300.0
4.5	1800.0	200.0	300.0
5.0	1900.0	200.0	300.0

And we can save it to a tabular CSV file that can be opened in Excel for further analysis:

results.to_csv("results.csv")

Finally, we can also easily plot these results:

results.plot(xlabel="gas cost", ylabel="cost (EUR) / power (MW)", marker='o', linestyle=':')

<Axes: xlabel='gas cost', ylabel='cost (EUR) / power (MW)'>

../_images/e51773baacdbcdbfb010ca928cfcc770ca5e2e3983191acfa431debb062356dd.png

7.7.2. Solver choice#

Many different solvers exist which implement a range of algorithms to solve optimisation problems.

Above, we used Gurobi, which can handle continuous linear problems, mixed-integer linear problems, as well as quadratic problems and some other types of nonlinear problems. For larger problems (more than a few 10,000 variables or constraints, and/or for mixed-integer problems) Gurobi is significantly faster than open-source alternatives.

To use it for problems with more than 2000 variables or constraints, you need to obtain a licens, which you can obtain for free as a student or academic researcher (see the Gurobi website).

To solve general nonlinear problems, you can use the IPOPT solver:

solver = pyo.SolverFactory('ipopt')

You could also solve a linear problem with a nonlinear solver like IPOPT, but that is not usually a good idea: the nonlinear solver is not guaranteed to find the global optimum. So remember to pick the right solver for the job!

Finally, if you prefer to use a free and open-source solver also for linear problems, you can instead use GLPK:

solver = pyo.SolverFactory('glpk')

If you have followed the reader’s installation instructions, then the following solvers are already installed and available:

Solver	Problem type(s)	Free?	Website	How to install	Use in Pyomo
GLPK	LP, MILP	Yes	https://www.gnu.org/software/glpk/	Easiest with Anaconda	`SolverFactory('glpk')`
Gurobi	LP, MILP, QP, MIQP	No	https://www.gurobi.com/	Easiest with Anaconda	`SolverFactory('gurobi')`
IPOPT	NLP	Yes	https://coin-or.github.io/Ipopt/	Easiest with Anaconda	`SolverFactory('ipopt')`
MindtPy	MINLP	Yes	https://pyomo.readthedocs.io/en/stable/explanation/solvers/mindtpy.html	Included in Pyomo	`SolverFactory('mindtpy')`

Pyomo basics

Contents

7. Pyomo basics#

7.1. Variables#

7.2. Objective function#

7.3. Constraints#

7.4. Solving the model#

7.5. Dual variables#

7.6. An easier way to access results#

7.7. Advanced topics#

7.7.1. Modifying and re-running a model, then saving results#

7.7.2. Solver choice#