Manipulating Pyomo Models

This section gives an overview of commonly used scripting commands when working with Pyomo models. These commands must be applied to a concrete model instance or in other words an instantiated model.

Repeated Solves

>>> import pyomo.environ as pyo
>>> from pyomo.opt import SolverFactory
>>> model = pyo.ConcreteModel()
>>> model.nVars = pyo.Param(initialize=4)
>>> model.N = pyo.RangeSet(model.nVars)
>>> model.x = pyo.Var(model.N, within=pyo.Binary)
>>> model.obj = pyo.Objective(expr=pyo.summation(model.x))
>>> model.cuts = pyo.ConstraintList()
>>> opt = SolverFactory('glpk')
>>> opt.solve(model) 

>>> # Iterate, adding a cut to exclude the previously found solution
>>> for i in range(5):
...    expr = 0
...    for j in model.x:
...        if pyo.value(model.x[j]) < 0.5:
...            expr += model.x[j]
...        else:
...            expr += (1 - model.x[j])
...    model.cuts.add( expr >= 1 )
...    results = opt.solve(model)
...    print ("\n===== iteration",i)
...    model.display()

To illustrate Python scripts for Pyomo we consider an example that is in the file iterative1.py and is executed using the command

python iterative1.py

Note

This is a Python script that contains elements of Pyomo, so it is executed using the python command. The pyomo command can be used, but then there will be some strange messages at the end when Pyomo finishes the script and attempts to send the results to a solver, which is what the pyomo command does.

This script creates a model, solves it, and then adds a constraint to preclude the solution just found. This process is repeated, so the script finds and prints multiple solutions. The particular model it creates is just the sum of four binary variables. One does not need a computer to solve the problem or even to iterate over solutions. This example is provided just to illustrate some elementary aspects of scripting.

#  ___________________________________________________________________________
#
#  Pyomo: Python Optimization Modeling Objects
#  Copyright (c) 2008-2024
#  National Technology and Engineering Solutions of Sandia, LLC
#  Under the terms of Contract DE-NA0003525 with National Technology and
#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
#  rights in this software.
#  This software is distributed under the 3-clause BSD License.
#  ___________________________________________________________________________

# iterative1.py
import pyomo.environ as pyo
from pyomo.opt import SolverFactory


# Create a solver
opt = pyo.SolverFactory('glpk')

#
# A simple model with binary variables and
# an empty constraint list.
#
model = pyo.AbstractModel()
model.n = pyo.Param(default=4)
model.x = pyo.Var(pyo.RangeSet(model.n), within=pyo.Binary)


def o_rule(model):
    return pyo.summation(model.x)


model.o = pyo.Objective(rule=o_rule)
model.c = pyo.ConstraintList()

# Create a model instance and optimize
instance = model.create_instance()
results = opt.solve(instance)
instance.display()

# Iterate to eliminate the previously found solution
for i in range(5):
    expr = 0
    for j in instance.x:
        if pyo.value(instance.x[j]) == 0:
            expr += instance.x[j]
        else:
            expr += 1 - instance.x[j]
    instance.c.add(expr >= 1)
    results = opt.solve(instance)
    print("\n===== iteration", i)
    instance.display()

Let us now analyze this script. The first line is a comment that happens to give the name of the file. This is followed by two lines that import symbols for Pyomo. The pyomo namespace is imported as pyo. Therefore, pyo. must precede each use of a Pyomo name.

# iterative1.py
import pyomo.environ as pyo
from pyomo.opt import SolverFactory

An object to perform optimization is created by calling SolverFactory with an argument giving the name of the solver. The argument would be 'gurobi' if, e.g., Gurobi was desired instead of glpk:

# Create a solver
opt = pyo.SolverFactory('glpk')

The next lines after a comment create a model. For our discussion here, we will refer to this as the base model because it will be extended by adding constraints later. (The words “base model” are not reserved words, they are just being introduced for the discussion of this example). There are no constraints in the base model, but that is just to keep it simple. Constraints could be present in the base model. Even though it is an abstract model, the base model is fully specified by these commands because it requires no external data:

model = pyo.AbstractModel()
model.n = pyo.Param(default=4)
model.x = pyo.Var(pyo.RangeSet(model.n), within=pyo.Binary)


def o_rule(model):
    return pyo.summation(model.x)


model.o = pyo.Objective(rule=o_rule)

The next line is not part of the base model specification. It creates an empty constraint list that the script will use to add constraints.

model.c = pyo.ConstraintList()

The next non-comment line creates the instantiated model and refers to the instance object with a Python variable instance. Models run using the pyomo script do not typically contain this line because model instantiation is done by the pyomo script. In this example, the create function is called without arguments because none are needed; however, the name of a file with data commands is given as an argument in many scripts.

instance = model.create_instance()

The next line invokes the solver and refers to the object contain results with the Python variable results.

results = opt.solve(instance)

The solve function loads the results into the instance, so the next line writes out the updated values.

instance.display()

The next non-comment line is a Python iteration command that will successively assign the integers from 0 to 4 to the Python variable i, although that variable is not used in script. This loop is what causes the script to generate five more solutions:

for i in range(5):

An expression is built up in the Python variable named expr. The Python variable j will be iteratively assigned all of the indexes of the variable x. For each index, the value of the variable (which was loaded by the load method just described) is tested to see if it is zero and the expression in expr is augmented accordingly. Although expr is initialized to 0 (an integer), its type will change to be a Pyomo expression when it is assigned expressions involving Pyomo variable objects:

    expr = 0
    for j in instance.x:
        if pyo.value(instance.x[j]) == 0:
            expr += instance.x[j]
        else:
            expr += 1 - instance.x[j]

During the first iteration (when i is 0), we know that all values of x will be 0, so we can anticipate what the expression will look like. We know that x is indexed by the integers from 1 to 4 so we know that j will take on the values from 1 to 4 and we also know that all value of x will be zero for all indexes so we know that the value of expr will be something like

0 + instance.x[1] + instance.x[2] + instance.x[3] + instance.x[4]

The value of j will be evaluated because it is a Python variable; however, because it is a Pyomo variable, the value of instance.x[j] not be used, instead the variable object will appear in the expression. That is exactly what we want in this case. When we wanted to use the current value in the if statement, we used the value function to get it.

The next line adds to the constraint list called c the requirement that the expression be greater than or equal to one:

    instance.c.add(expr >= 1)

The proof that this precludes the last solution is left as an exercise for the reader.

The final lines in the outer for loop find a solution and display it:

    results = opt.solve(instance)
    print("\n===== iteration", i)
    instance.display()

Note

The assignment of the solve output to a results object is somewhat anachronistic. Many scripts just use

>>> opt.solve(instance)

since the results are moved to the instance by default, leaving the results object with little of interest. If, for some reason, you want the results to stay in the results object and not be moved to the instance, you would use

>>> results = opt.solve(instance, load_solutions=False)

This approach can be useful if there is a concern that the solver did not terminate with an optimal solution. For example,

>>> results = opt.solve(instance, load_solutions=False) 
>>> if results.solver.termination_condition == TerminationCondition.optimal: 
...     instance.solutions.load_from(results)

Changing the Model or Data and Re-solving

The iterative1.py example above illustrates how a model can be changed and then re-solved. In that example, the model is changed by adding a constraint, but the model could also be changed by altering the values of parameters. Note, however, that in these examples, we make the changes to the concrete model instances. This is particularly important for AbstractModel users, as this implies working with the instance object rather than the model object, which allows us to avoid creating a new model object for each solve. Here is the basic idea for users of an AbstractModel:

  1. Create an AbstractModel (suppose it is called model)

  2. Call model.create_instance() to create an instance (suppose it is called instance)

  3. Solve instance

  4. Change something in instance

  5. Solve instance again

Note

Users of ConcreteModel typically name their models model, which can cause confusion to novice readers of documentation. Examples based on an AbstractModel will refer to instance where users of a ConcreteModel would typically use the name model.

If instance has a parameter whose name is Theta that was declared to be mutable (i.e., mutable=True) with an index that contains idx, then the value in NewVal can be assigned to it using

>>> instance.Theta[idx] = NewVal

For a singleton parameter named sigma (i.e., if it is not indexed), the assignment can be made using

>>> instance.sigma = NewVal

Note

If the Param is not declared to be mutable, an error will occur if an assignment to it is attempted.

For more information about access to Pyomo parameters, see the section in this document on Param access Accessing Parameter Values. Note that for concrete models, the model is the instance.

Fixing Variables and Re-solving

Instead of changing model data, scripts are often used to fix variable values. The following example illustrates this.

#  ___________________________________________________________________________
#
#  Pyomo: Python Optimization Modeling Objects
#  Copyright (c) 2008-2024
#  National Technology and Engineering Solutions of Sandia, LLC
#  Under the terms of Contract DE-NA0003525 with National Technology and
#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
#  rights in this software.
#  This software is distributed under the 3-clause BSD License.
#  ___________________________________________________________________________

# iterative2.py

import pyomo.environ as pyo
from pyomo.opt import SolverFactory

# Create a solver
opt = pyo.SolverFactory('cplex')

#
# A simple model with binary variables and
# an empty constraint list.
#
model = pyo.AbstractModel()
model.n = pyo.Param(default=4)
model.x = pyo.Var(pyo.RangeSet(model.n), within=pyo.Binary)


def o_rule(model):
    return pyo.summation(model.x)


model.o = pyo.Objective(rule=o_rule)
model.c = pyo.ConstraintList()

# Create a model instance and optimize
instance = model.create_instance()
results = opt.solve(instance)
instance.display()

# "flip" the value of x[2] (it is binary)
# then solve again

if pyo.value(instance.x[2]) == 0:
    instance.x[2].fix(1)
else:
    instance.x[2].fix(0)

results = opt.solve(instance)
instance.display()

In this example, the variables are binary. The model is solved and then the value of model.x[2] is flipped to the opposite value before solving the model again. The main lines of interest are:


if pyo.value(instance.x[2]) == 0:
    instance.x[2].fix(1)
else:
    instance.x[2].fix(0)

results = opt.solve(instance)

This could also have been accomplished by setting the upper and lower bounds:

>>> if instance.x[2].value == 0:
...     instance.x[2].setlb(1)
...     instance.x[2].setub(1)
... else:
...     instance.x[2].setlb(0)
...     instance.x[2].setub(0)

Notice that when using the bounds, we do not set fixed to True because that would fix the variable at whatever value it presently has and then the bounds would be ignored by the solver.

For more information about access to Pyomo variables, see the section in this document on Var access Accessing Variable Values.

Note that

>>> instance.x.fix(1)

is equivalent to

>>> instance.x.value = 1
>>> instance.x.fixed = True
and
>>> instance.x.fix()

is equivalent to

>>> instance.x.fixed = True

Extending the Objective Function

One can add terms to an objective function of a ConcreteModel (or and instantiated AbstractModel) using the expr attribute of the objective function object. Here is a simple example:

>>> import pyomo.environ as pyo
>>> from pyomo.opt import SolverFactory

>>> model = pyo.ConcreteModel()

>>> model.x = pyo.Var(within=pyo.PositiveReals)
>>> model.y = pyo.Var(within=pyo.PositiveReals)

>>> model.sillybound = pyo.Constraint(expr = model.x + model.y <= 2)

>>> model.obj = pyo.Objective(expr = 20 * model.x)

>>> opt = SolverFactory('glpk') 
>>> opt.solve(model) 

>>> model.pprint() 

>>> print ("------------- extend obj --------------") 
>>> model.obj.expr += 10 * model.y

>>> opt.solve(model) 
>>> model.pprint()

Activating and Deactivating Objectives

Multiple objectives can be declared, but only one can be active at a time (at present, Pyomo does not support any solvers that can be given more than one objective). If both model.obj1 and model.obj2 have been declared using Objective, then one can ensure that model.obj2 is passed to the solver as shown in this simple example:

>>> model = pyo.ConcreteModel()
>>> model.obj1 = pyo.Objective(expr = 0)
>>> model.obj2 = pyo.Objective(expr = 0)

>>> model.obj1.deactivate()
>>> model.obj2.activate()

For abstract models this would be done prior to instantiation or else the activate and deactivate calls would be on the instance rather than the model.

Activating and Deactivating Constraints

Constraints can be temporarily disabled using the deactivate() method. When the model is sent to a solver inactive constraints are not included. Disabled constraints can be re-enabled using the activate() method.

>>> model = pyo.ConcreteModel()
>>> model.v = pyo.Var()
>>> model.con = pyo.Constraint(expr=model.v**2 + model.v >= 3)
>>> model.con.deactivate()
>>> model.con.activate()

Indexed constraints can be deactivated/activated as a whole or by individual index:

>>> model = pyo.ConcreteModel()
>>> model.s = pyo.Set(initialize=[1,2,3])
>>> model.v = pyo.Var(model.s)
>>> def _con(m, s):
...    return m.v[s]**2 + m.v[s] >= 3
>>> model.con = pyo.Constraint(model.s, rule=_con)
>>> model.con.deactivate()   # Deactivate all indices
>>> model.con[1].activate()  # Activate single index