PAO Documentation¶
PAO is a Python-based package for Adversarial Optimization. The goal of this package is to provide a general modeling and analysis capability for bilevel, trilevel and other multilevel optimization forms that express adversarial dynamics. PAO integrates two different modeling abstractions:
- Algebraic models extend the modeling concepts in the Pyomo algebraic modeling language to express problems with an intuitive algebraic syntax. Thus, we expect that this modeling abstraction will commonly be used by PAO end-users.
- Compact models express objective and constraints in a manner that is typically used to express the mathematical form of these problems (e.g. using vector and matrix data types). PAO defines custom Multilevel Problem Representations (MPRs) that simplify the implementation of solvers for bilevel, trilevel and other multilevel optimization problems.
Indices and tables¶
PAO Resources¶
PAO development is hosted at GitHub:
The OR-Fusion GitHub organization is used to coordinate installation of PAO with other OR-related capabilities:
Ask a question on StackOverflow:
Installation¶
PAO currently supports the following versions of Python:
- CPython: 3.7, 3.8, 3.9
Using PIP¶
The standard utility for installing Python packages is pip. You can install the latest release of PAO by executing the following:
pip install pao
You can also use pip to install from the PAO software repository. For example, the master branch can be installed as follows:
python -m pip install https://github.com/or-fusion/pao.git
Using CONDA¶
The conda utility can be used to install the latest release of PAO
using the conda-forge
channel:
conda install -c conda-forge pao
Using GIT¶
PAO can be installed by cloning the PAO software repostory and then directly installing the software. For example, the master branch can be installed as follows:
git clone https://github.com/or-fusion/pao.git
cd pao
python setup.py develop
Conditional Dependencies¶
Both conda and pip can be used to install the third-party packages that are needed to model problems with PAO. We recommend conda because it has better support for optimization solver packages.
PAO intrinsically depends on Pyomo, both for the representation of algebraic problems but also for interfaces to numerical optimizers used by PAO solvers. Pyomo is installed with PAO, but the Pyomo website [PyomoWeb] and GitHub site [PyomoGithub] provide additional resources for installing Pyomo and related software.
PAO and Pyomo have conditional dependencies on a variety of third-party packages, including Python packages like scipy, numpy and optimization solvers. Optimization solvers are particularly important, and a commercial optimizer may be needed to analyze complex, real-world applications.
The following optimizers are used to test the PAO solvers:
- glpk - An open-source mixed-integer linear programming solver
- cbc - An open-source mixed-integer linear programming solver
- ipopt - An open-source interior point optimizer for continuous problems
Additionally, PAO can interface with optimization solvers at NEOS.
Overview¶
Many planning situations involve the analysis of a hierarchy of decision-makers with competing objectives. For example, policy decisions are made at different levels of a government, each of which has a different objective and decision space. Similarly, robust planning against adversaries is often modeled with a 2-level hierarchy, where the defensive planner makes decisions that account for adversarial response. Multilevel optimization techniques partition control over decision variables amongst the levels. Decisions at each level of the hierarchy may be constrained by decisions at other levels, and the objectives for each level may account for decisions made at other levels. The PAO library is designed to express this structure in a manner that is intuitive to users, and which facilitates the application of appropriate optimization solvers. In particular, PAO extends the modeling concepts in the Pyomo algebraic modeling language to express problems with an intuitive algebraic syntax.
However, data structures used to represent algebraic models are often poorly suited for complex numerical algorithms. Consequently, PAO supports distinct representations for algebraic models (using Pyomo) and compact problem representations that express objective and constraints using vector and matrix data types. Currently, PAO includes several Multilevel Problem Representations (MPRs) that represent multilevel optimization problems with an explicit, compact representation that simplifies the implementation of solvers for bilevel, trilevel and other multilevel optimization problems.
For example, PAO includes a compact representation for linear bilevel
problems, LinearMultilevelProblem
. Several solvers have been
developed for problems expressed as a LinearMultilevelProblem
,
including the big-M method proposed by Fortuny-Amat and McCarl
[FortunyMcCarl]. PAO can similarly express linear bilevel problems in
Pyomo using a Submodel
component, which was previously introduced in
the pyomo.bilevel package [PyomoBookII]. Further, these Pyomo models
can be automatically converted to the compact LinearMultilevelProblem
representation and optimized with solvers tailored for that representation.
The use of independent problem representations in this manner has
several implications for PAO. First, this design facilitates the
development of solvers for algebraic modeling languages like Pyomo
that are intrinsically more robust. Compact representations like
LinearMultilevelProblem
enable the development of solvers using
natural operations (e.g. matrix-vector multiplication). Thus, we expect
these solvers to be more robust and easier to maintain when compared to
solvers developed using Pyomo data structures (e.g. expression trees).
Additionally, the conversion of a Pyomo representation to a compact
representation provides a context for verifying that the Pyomo model
represents the intended problem form.
Second, this design facilitates the development of different problem representations that may or may not be inter-operable. Although PAO is derived from initial efforts with pyomo.bilevel, it has evolved from an extension of Pyomo’s modeling capability to be a library that is synergistic with Pyomo but not strictly dependent on it.
The following sections provide detailed descriptions of these representations that are illustrated with increasingly complex examples, including
- bilevel
- bilevel with multiple lower-levels
- trilevel
Additionally, the following sections describe the setup of linear and quadratic problems, and the transformations that can be applied to them in PAO.
Note
We do not restrict the description of PAO representations to models that PAO can solve. Rather, a goal of this documentation to illustrate the breadth of the adversarial optimization problems that can be expressed with PAO, thereby motivating the development of new solvers for new classes of problems.
Note
Although the modeling abstractions in PAO are suitable for both pessimistic and optimistic multilevel optimization formulations, PAO currently assumes that all problems are expressed and solved in an an optimistic form. Thus, when there multiple solutions in lower-level problems, the choice of the solution is determined by the follower not the leader. This convention reflects the fact that solution techniques for pessimistic forms are not well-developed.
Simple Examples¶
We illustrate the PAO algebraic and compact representations with a series of simple examples. Consider the following bilevel problem introduced by Bard [Bard98] (example 5.1.1):
This problem has linear upper- and lower-level problems with different
objectives in each level. Thus, this problem can be represented in
PAO using a Pyomo model representation or the LinearMultilevelProblem
representation.
Using Pyomo¶
The following python script defines a bilevel problem in Pyomo:
>>> import pyomo.environ as pe
>>> from pao.pyomo import *
# Create a model object
>>> M = pe.ConcreteModel()
# Define decision variables
>>> M.x = pe.Var(bounds=(0,None))
>>> M.y = pe.Var(bounds=(0,None))
# Define the upper-level objective
>>> M.o = pe.Objective(expr=M.x - 4*M.y)
# Create a SubModel component to declare a lower-level problem
# The variable M.x is fixed in this lower-level problem
>>> M.L = SubModel(fixed=M.x)
# Define the lower-level objective
>>> M.L.o = pe.Objective(expr=M.y)
# Define lower-level constraints
>>> M.L.c1 = pe.Constraint(expr= -M.x - M.y <= -3)
>>> M.L.c2 = pe.Constraint(expr= -2*M.x + M.y <= 0)
>>> M.L.c3 = pe.Constraint(expr= 2*M.x + M.y <= 12)
>>> M.L.c4 = pe.Constraint(expr= 3*M.x - 2*M.y <= 4)
# Create a solver and apply it
>>> with Solver('pao.pyomo.FA') as solver:
... results = solver.solve(M)
# The final solution is loaded into the model
>>> print(M.x.value)
4.0
>>> print(M.y.value)
4.0
The SubModel
component defines a Pyomo block object within which the
lower-level problem is declared. The fixed
option is used to specify
the upper-level variables whose value is fixed in the lower-level problem.
The pao.pyomo.FA
uses the method for solving linear bilevel
programs with big-M relaxations described by Fortuny-Amat and McCarl
[FortunyMcCarl]. In this example, the default big-M value is used, but
in general a user may need to explore suitable values for their problem.
By default, the final values of the upper- and lower-level variables
are loaded back into the Pyomo model. The results
object contains
information about the problem, the solver and the solver status.
Using LinearMultilevelProblem¶
Bilevel linear problems can also be represented using the
LinearMultilevelProblem
class. This class provides a simple mechanism
for organizing data for variables, objectives and linear constraints. The following
examples illustrate the use of LinearMultilevelProblem
for Bard’s example 5.1.1 described
above, but this representation can naturally be used to express multi-level problems as well
as problems with multiple lower-levels.
Using Numpy and Scipy Data¶
The following python script defines a bilevel problem using LinearMultilevelProblem
with
numpy and scipy data:
>>> import numpy as np
>>> from scipy.sparse import coo_matrix
>>> from pao.mpr import *
# Create a model object
>>> M = LinearMultilevelProblem()
# Declare the upper- and lower-levels, including the number of decision-variables
# nxR=1 means there will be 1 real-valued decision variable
>>> U = M.add_upper(nxR=1)
>>> L = U.add_lower(nxR=1)
# Declare the bounds on the decision variables
>>> U.x.lower_bounds = np.array([0])
>>> L.x.lower_bounds = np.array([0])
# Declare the upper-level objective
# U.c[X] is the array of coefficients in the objective for variables in level X
>>> U.c[U] = np.array([1])
>>> U.c[L] = np.array([-4])
# Declare the lower-level objective, which has no upper-level decision-variables
>>> L.c[L] = np.array([1])
# Declare the lower-level constraints
# L.A[X] is the matrix coefficients in the constraints for variables in level X
>>> L.A[U] = coo_matrix((np.array([-1, -2, 2, 3]), # Coefficients
... (np.array([0, 1, 2, 3]), # Row indices of matrix entries
... np.array([0, 0, 0, 0])))) # Column indices of matrix entries
>>> L.A[L] = coo_matrix((np.array([-1, 1, 1, -2]),
... (np.array([0, 1, 2, 3]),
... np.array([0, 0, 0, 0]))))
# Declare the constraint right-hand-side
# By default, constraints are inequalities, so these are upper-bounds
>>> L.b = np.array([-3, 0, 12, 4])
# Create a solver and apply it
>>> with Solver('pao.mpr.FA') as solver:
... results = solver.solve(M)
# The final solution is loaded into the model
>>> print(U.x.values[0])
4.0
>>> print(L.x.values[0])
4.0
The U
and L
objects represent the upper- and lower-level
respectively. When declaring these objects, the user specifies the number
of real, integer and binary variables. The remaining declarations assume
that these variables are used in that order. Thus, there is a single
declaration for the objective coefficients, c
, which is an array
with values for each of the declared variables. However, the upper-
and lower-level objective coefficients are separately declared for
the upper- and lower-level variables by indexing c
with U
and
L
respectively. This example includes declarations for the upper-
and lower-level variable bounds and objective coefficients. There are no
upper-level constraints, so only the lower-level constriants are declared.
Note that the syntax for specifying solvers is analogous to that used with Pyomo models. The same solver options are available. The principle difference is the specification of the solver name that indicates the expected type of the model that will be solved.
Using Python Lists and Dictionaries¶
Although the constraint matrices are dense in this example, the
coo_matrix
is used to illustrate the general support for sparse data.
The LinearMultilevelProblem
class also supports a simpler syntax
where dense arrays can be specified and Python lists and sparse matrices
can be specified with Python tuple and dictionary objects:
>>> from pao.mpr import *
>>> M = LinearMultilevelProblem()
>>> U = M.add_upper(nxR=1)
>>> L = U.add_lower(nxR=1)
>>> U.x.lower_bounds = [0]
>>> L.x.lower_bounds = [0]
>>> U.c[U] = [1]
>>> U.c[L] = [-4]
>>> L.c[L] = [1]
>>> L.A[U] = (4,1), {(0,0):-1, (1,0):-2, (2,0):2, (3,0): 3}
>>> L.A[L] = (4,1), {(0,0):-1, (1,0): 1, (2,0):1, (3,0):-2}
>>> L.b = [-3, 0, 12, 4]
>>> with Solver('pao.mpr.FA') as solver:
... results = solver.solve(M)
>>> print(U.x.values[0])
4.0
>>> print(L.x.values[0])
4.0
When specifying a sparse matrix, a tuple is provided (e.g. for
L.A[U]
). The first element is a 2-tuple that defines the shape
of the matrix, and the second element is a dictionary that defines the
non-zero values in the sparse matrix.
Similarly, a list-of-lists syntax can be used to specify dense matrices:
>>> from pao.mpr import *
>>> M = LinearMultilevelProblem()
>>> U = M.add_upper(nxR=1)
>>> L = U.add_lower(nxR=1)
>>> U.x.lower_bounds = [0]
>>> L.x.lower_bounds = [0]
>>> U.c[U] = [1]
>>> U.c[L] = [-4]
>>> L.c[L] = [1]
>>> L.A[U] = [[-1], [-2], [2], [3]]
>>> L.A[L] = [[-1], [1], [1], [-2]]
>>> L.b = [-3, 0, 12, 4]
>>> with Solver('pao.mpr.FA') as solver:
... results = solver.solve(M)
>>> print(U.x.values[0])
4.0
>>> print(L.x.values[0])
4.0
When native Python data values are used to initialize a
LinearMultilevelProblem
, they are converted into numpy and scipy
data types. This facilitates the use of LinearMultilevelProblem
objects for defining numerical solvers using a consistent, convenient
API for numerical operations (e.g. matrix-vector multiplication).
Pyomo Models¶
PAO can be used to express linear and quadratic problems in Pyomo using a SubModel
component, which was previously introduced in the pyomo.bilevel
package [PyomoBookII]. Pyomo
represents optimization models using an model objects that are
annotated with modeling component objects. Thus, SubModel
component is a simple extension of the modeling concepts in Pyomo.
Hint
Advanced Pyomo users will realize that the PAO SubModel
component
is a special case of the Pyomo Block
component, which is used to
structure the expression of Pyomo models.
A SubModel
component creates a context for expressing the
objective and constraints in a lower-level model. Pyomo models can
include nested and parallel SubModel
components to express
complex multilevel problems.
Consider the following bilevel problem:
This problem has has linear upper- and lower-level problems with different objectives in each level.
>>> M = pe.ConcreteModel()
>>> M.x = pe.Var(bounds=(2,6))
>>> M.y = pe.Var()
>>> M.L = pao.pyomo.SubModel(fixed=[M.x,M.y])
>>> M.L.z = pe.Var(bounds=(0,None))
>>> M.o = pe.Objective(expr=M.x + 3*M.L.z, sense=pe.minimize)
>>> M.c = pe.Constraint(expr= M.x + M.y == 10)
>>> M.L.o = pe.Objective(expr=M.L.z, sense=pe.maximize)
>>> M.L.c1 = pe.Constraint(expr= M.x + M.L.z <= 8)
>>> M.L.c2 = pe.Constraint(expr= M.x + 4*M.L.z >= 8)
>>> M.L.c3 = pe.Constraint(expr= M.x + 2*M.L.z <= 13)
>>> opt = pao.Solver("pao.pyomo.FA")
>>> results = opt.solve(M)
>>> print(M.x.value, M.y.value, M.L.z.value)
6.0 4.0 2.0
This example illustrates the flexibility of Pyomo representations in PAO:
- Each level can express different objectives with different senses
- Variables can be bounded or unbounded
- Equality and inequality constraints can be expressed
The SubModel
component is used to define a logically separate
optimization model that includes variables that are dynamically fixed
by upper-level problems. All of the Pyomo objective and constraint
declarations contained in the SubModel
declaration are included
in the sub-problem that it defines, even if they are nested in Pyomo
Block
components. The SubModel
component also declares
which variables are fixed in a lower-level problem. The value of the
fixed argument is a Pyomo variable or a list of variables. For example,
the following model expresses the upper-level variables with a single
variable, M.x, which is fixed in the SubModel
declaration:
>>> M = pe.ConcreteModel()
>>> M.x = pe.Var([0,1])
>>> M.x[0].setlb(2)
>>> M.x[0].setub(6)
>>> M.L = pao.pyomo.SubModel(fixed=M.x)
>>> M.L.z = pe.Var(bounds=(0,None))
>>> M.o = pe.Objective(expr=M.x[0] + 3*M.L.z, sense=pe.minimize)
>>> M.c = pe.Constraint(expr= M.x[0] + M.x[1] == 10)
>>> M.L.o = pe.Objective(expr=M.L.z, sense=pe.maximize)
>>> M.L.c1 = pe.Constraint(expr= M.x[0] + M.L.z <= 8)
>>> M.L.c2 = pe.Constraint(expr= M.x[0] + 4*M.L.z >= 8)
>>> M.L.c3 = pe.Constraint(expr= M.x[0] + 2*M.L.z <= 13)
>>> opt = pao.Solver("pao.pyomo.FA")
>>> results = opt.solve(M)
>>> print(M.x[0].value, M.x[1].value, M.L.z.value)
6.0 4.0 2.0
Although a lower-level problem is logically a separate optimization model,
you cannot use a SubModel
that is defined with a separate Pyomo
model object. Pyomo implicitly requires that all variables used in
objective and constraint expressions are attributes of the same Pyomo model.
However, the location of variable declarations in a Pyomo model does not denote their
use in upper- or lower-level problems. For example, consider the following
model that re-expresses the previous problem:
>>> M = pe.ConcreteModel()
>>> M.x = pe.Var(bounds=(2,6))
>>> M.y = pe.Var()
>>> M.z = pe.Var(bounds=(0,None))
>>> M.o = pe.Objective(expr=M.x + 3*M.z, sense=pe.minimize)
>>> M.c = pe.Constraint(expr= M.x + M.y == 10)
>>> M.L = pao.pyomo.SubModel(fixed=[M.x,M.y])
>>> M.L.o = pe.Objective(expr=M.z, sense=pe.maximize)
>>> M.L.c1 = pe.Constraint(expr= M.x + M.z <= 8)
>>> M.L.c2 = pe.Constraint(expr= M.x + 4*M.z >= 8)
>>> M.L.c3 = pe.Constraint(expr= M.x + 2*M.z <= 13)
>>> opt = pao.Solver("pao.pyomo.FA")
>>> results = opt.solve(M)
>>> print(M.x.value, M.y.value, M.z.value)
6.0 4.0 2.0
Note that all of the decision variables are declared outside of the
SubModel
component, even though the variable M.z
is a
lower-level variable. The declarations of SubModel
components
defines the mathematical role of all decision variables in a Pyomo model.
As this example illustrates, the specification of a bilevel problem can
be simplified if all variables are expressed at once.
Finally, we observe that PAO’s Pyomo representation only works with a subset of the many different modeling components that are supported in Pyomo:
Set
- Set declarationsParam
- Parameter declarationsVar
- Variable declarationsBlock
- Defines a subset of a modelObjective
- Define a model objectiveConstraint
- Define model constraints
Additional Pyomo modeling components will be added to PAO as motivating applications arise and as suitable solvers become available.
Multilevel problems can be easily expressed with Pyomo using multiple declarations
of SubModel
.
Consider the following bilevel problem that extends the PAO1 model to include two equivalent lower-levels:
The PAO2 model can be expressed in Pyomo as follows:
>>> M = pe.ConcreteModel()
>>> M.x = pe.Var(bounds=(2,6))
>>> M.y = pe.Var()
>>> M.z = pe.Var([1,2], bounds=(0,None))
>>> M.o = pe.Objective(expr=M.x + 3*M.z[1]+3*M.z[2], sense=pe.minimize)
>>> M.c = pe.Constraint(expr= M.x + M.y == 10)
>>> M.L1 = pao.pyomo.SubModel(fixed=[M.x])
>>> M.L1.o = pe.Objective(expr=M.z[1], sense=pe.maximize)
>>> M.L1.c1 = pe.Constraint(expr= M.x + M.z[1] <= 8)
>>> M.L1.c2 = pe.Constraint(expr= M.x + 4*M.z[1] >= 8)
>>> M.L1.c3 = pe.Constraint(expr= M.x + 2*M.z[1] <= 13)
>>> M.L2 = pao.pyomo.SubModel(fixed=[M.y])
>>> M.L2.o = pe.Objective(expr=M.z[2], sense=pe.maximize)
>>> M.L2.c1 = pe.Constraint(expr= M.y + M.z[2] <= 8)
>>> M.L2.c2 = pe.Constraint(expr= M.y + 4*M.z[2] >= 8)
>>> M.L2.c3 = pe.Constraint(expr= M.y + 2*M.z[2] <= 13)
>>> opt = pao.Solver("pao.pyomo.FA")
>>> results = opt.solve(M)
>>> print(M.x.value, M.y.value, M.z[1].value, M.z[2].value)
2.0 8.0 5.5 0.0
Trilevel problems can be described with nested declarations of
SubModel
components. Consider the following trilevel continuous
linear problem described by Anadalingam [Anadalingam]:
This model can be expressed in Pyomo as follows:
>>> M = pe.ConcreteModel()
>>> M.x1 = pe.Var(bounds=(0,None))
>>> M.x2 = pe.Var(bounds=(0,None))
>>> M.x3 = pe.Var(bounds=(0,0.5))
>>> M.L = pao.pyomo.SubModel(fixed=M.x1)
>>> M.L.B = pao.pyomo.SubModel(fixed=M.x2)
>>> M.o = pe.Objective(expr=-7*M.x1 - 3*M.x2 + 4*M.x3)
>>> M.L.o = pe.Objective(expr=-M.x2)
>>> M.L.B.o = pe.Objective(expr=-M.x3)
>>> M.L.B.c1 = pe.Constraint(expr= M.x1 + M.x2 + M.x3 <= 3)
>>> M.L.B.c2 = pe.Constraint(expr= M.x1 + M.x2 - M.x3 <= 1)
>>> M.L.B.c3 = pe.Constraint(expr= M.x1 + M.x2 + M.x3 >= 1)
>>> M.L.B.c4 = pe.Constraint(expr= - M.x1 + M.x2 + M.x3 <= 1)
Note
PAO solvers cannot currently solve trilevel problems like this, but an issue has been submitted to add this functionality.
PAO models using Pyomo represent general quadratic problems with quadratic terms in the objective and constraints at each level. The special case where bilinear terms arise with an upper-level binary variable multiplied with a lower-level variable is common in many applications. For this case, the PAO solvers for Pyomo models include an option to linearize these bilinear terms.
The following models considers a variation of the PAO1 model where binary variables control the expression of lower-level constraints:
The PAO3 model can be expressed in Pyomo as follows:
>>> M = pe.ConcreteModel()
>>> M.w = pe.Var([1,2], within=pe.Binary)
>>> M.x = pe.Var(bounds=(2,6))
>>> M.y = pe.Var()
>>> M.z = pe.Var(bounds=(0,None))
>>> M.o = pe.Objective(expr=M.x + 3*M.z+5*M.w[1], sense=pe.minimize)
>>> M.c1 = pe.Constraint(expr= M.x + M.y == 10)
>>> M.c2 = pe.Constraint(expr= M.w[1] + M.w[2] >= 1)
>>> M.L = pao.pyomo.SubModel(fixed=[M.x,M.y,M.w])
>>> M.L.o = pe.Objective(expr=M.z, sense=pe.maximize)
>>> M.L.c1 = pe.Constraint(expr= M.x + M.w[1]*M.z <= 8)
>>> M.L.c2 = pe.Constraint(expr= M.x + 4*M.z >= 8)
>>> M.L.c3 = pe.Constraint(expr= M.x + 2*M.w[2]*M.z <= 13)
>>> opt = pao.Solver("pao.pyomo.FA", linearize_bigm=100)
>>> results = opt.solve(M)
>>> print(M.x.value, M.y.value, M.z.value, M.w[1].value, M.w[2].value)
6.0 4.0 3.5 0 1
In general, it may be difficult to determine a valid value of $M$. Thus, this transformation may result in a restriction or relaxation of the original problem (depending on where the big-M values are introduced).
Multilevel Representations¶
PAO includes several Multilevel Problem Representations (MPRs) that represent multilevel optimization problems with an explicit, compact representation that simplifies the implementation of solvers for bilevel, trilevel and other multilevel optimization problems. These MPRs express objective and constraints using vector and matrix data types. However, they organize these data types to provide a semantically clear organization of multilevel problems. Additionally, the MPRs provide checks to ensure the consistency of the data within and across levels.
The classes LinearMultilevelProblem
and QuadraticMultilevelProblem
respectively
represent linear and quadratic multilevel problems. Although
QuadraticMultilevelProblem
is a generalization of the
representation in LinearMultilevelProblem
, the use of tailored
representations for different classes of problems clarifies the semantic
context when using them. For example, this allows for simple error checking
to confirm that a problem is linear.
Currently, all PAO solvers for MPRs support only linear problems, so the
following sections focus on LinearMultilevelProblem
. However,
we conclude with an example of model transformations that enable the
solution of quadratic problems using QuadraticMultilevelProblem
.
Note
We do not expect many users to directly employ a MPR data representation for their applications. Perhaps this would be desirable if their problem was already represented with matrix and vector data. In general, the algebraic representation supported by Pyomo will be more convenient for large, complex applications.
We expect this representation to be more useful for researchers developing multilevel solvers, since the MPR representations provide structure that simplifies the expression of necessary mathematical operations for these problems.
We consider again the bilevel problem PAO1 (1). This problem has linear upper- and lower-level problems with different objectives in each level.
>>> M = pao.mpr.LinearMultilevelProblem()
>>> U = M.add_upper(nxR=2)
>>> L = U.add_lower(nxR=1)
>>> U.x.lower_bounds = [2, np.NINF]
>>> U.x.upper_bounds = [6, np.PINF]
>>> L.x.lower_bounds = [0]
>>> L.x.upper_bounds = [np.PINF]
>>> U.c[U] = [1, 0]
>>> U.c[L] = [3]
>>> L.c[L] = [1]
>>> L.maximize = True
>>> U.equalities = True
>>> U.A[U] = [[1, 1]]
>>> U.b = [10]
>>> L.A[U] = [[ 1, 0],
... [-1, 0],
... [ 1, 0]]
>>> L.A[L] = [[ 1],
... [-4],
... [ 2]]
>>> L.b = [8, -8, 13]
>>> M.check()
>>> opt = pao.Solver("pao.mpr.FA")
>>> results = opt.solve(M)
>>> print(M.U.x.values)
[6.0, 4.0]
>>> print(M.U.LL.x.values)
[2.0]
The example illustrates both the flexibility of the MPR representions
in PAO but also the structure they enforce on the multilevel problem
representation. The upper-level problem is created by calling the
add_upper()
method, which takes arguments that specify the
variables at that level:
nxR
- The number of real variables (Default: 0)nxZ
- The number of general integer variables (Default: 0)nxB
- The number of binary variables (Default: 0)
In each level, the variables are represented as a vector of values, ordered in this manner.
Similarly, the add_lower()
method is used to generate a lower-level
problem from a given level. Note that this allows for the specification
of arbitrary nesting of levels, since a lower-level can be defined
relative to any other level in the model. Additionally, multiple
lower-levels can be specified for relative to a single level (see below).
The add_upper()
and add_lower()
methods return the
corresponding level object, which is used to specify data in the model
later.
For a given level object, Z
, the data Z.x
contains
information about the decision variables. In particular, the values
Z.x.lower_bounds
and Z.x.upper_bounds
can be set with arrays
of numeric values to specify lower- and upper-bounds on the decision
variables. Note that missing lower- and upper-bounds are specified with
numpy.NINF
and numpy.PINF
respectively.
The Z.c
data specifies coefficients of the objective function for
this level. This data is indexed by a level object B
to indicate
the data associated with the variables in B
. In the example above:
U.c[U]
is the array of coefficients of the upper-level objective for the variables in the upper-level,U.c[L]
is the array of coefficients of the upper-level objective for the variables in the lower-level, andL.c[L]
is the array of coefficients of the lower-level objective for the variables in the lower-level.
Since L.c[U]
is not specified, it has a value None
that
indicates that no upper-level variables have non-zero coefficients in the
lower-level objective. The Z.A
data specifies the matrix coefficients
for the constraints using a similar indexing notation and semantics.
The values Z.minimize
and Z.maximize
can be set to True
to indicate whether the objective in Z
minimizes or maximizes.
(The default is minimize.) Similarly the value Z.inequalities
and Z.equalities
can be set to True
to indicate whether the
constraints in Z
are inequalities or equalities. (The default
is inequalities.) Finally, the value Z.b
defines the array of
constraint right-hand-side values.
The :meth:check
method provides a convenient sanity check that the
data is defined consistently within each level and between levels.
Note that PAO supports a consistent interface for creating a solver
interface and for applying solvers. In fact, the user should be
aware that Pyomo and MPR solvers are named in a consistent fashion.
For example, the Pyomo solver pao.pyomo.FA calls the MPR solver
pao.mpr.FA after automatically converting the Pyomo representation
to a :class:LinearMultilevelProblem
representation. This example
illustrates that values Z.x.values
contains the values of each level
Z
after optimization.
Multilevel problems can be easily expressed using the same MPR data representation.
We consider again the bilevel problem PAO2 (2), which has multiple lower-level problems. The PAO2 model can be expressed as a linear multilevel problem as follows:
>>> M = pao.mpr.LinearMultilevelProblem()
>>> U = M.add_upper(nxR=2)
>>> L1 = U.add_lower(nxR=1)
>>> L2 = U.add_lower(nxR=1)
>>> U.x.lower_bounds = [2, np.NINF]
>>> U.x.upper_bounds = [6, np.PINF]
>>> U.c[U] = [1, 0]
>>> U.c[L1] = [3]
>>> U.c[L2] = [3]
>>> U.equalities = True
>>> U.A[U] = [[1, 1]]
>>> U.b = [10]
>>> L1.x.lower_bounds = [0]
>>> L1.x.upper_bounds = [np.PINF]
>>> L1.c[L1] = [1]
>>> L1.maximize = True
>>> L1.A[U] = [[ 1, 0],
... [-1, 0],
... [ 1, 0]]
>>> L1.A[L1] = [[ 1],
... [-4],
... [ 2]]
>>> L1.b = [8, -8, 13]
>>> L2.x.lower_bounds = [0]
>>> L2.x.upper_bounds = [np.PINF]
>>> L2.c[L2] = [1]
>>> L2.maximize = True
>>> L2.A[U] = [[0, 1],
... [0, -1],
... [0, 1]]
>>> L2.A[L2] = [[ 1],
... [-4],
... [ 2]]
>>> L2.b = [8, -8, 13]
>>> opt = pao.Solver("pao.mpr.FA")
>>> results = opt.solve(M)
>>> print(U.x.values)
[2.0, 8.0]
>>> print(L1.x.values)
[5.5]
>>> print(L2.x.values)
[0.0]
The declarataion of the two lower level problems is naturally contained
within the data of the L1
and L2
objects. Further, the
cross-level interactions are intuitively represented using the index
notation for the objective and constraint data objects.
Note that this more explicit representation clarifies some ambiguity in the expression of lower-levels in the Pyomo representation. The Pyomo representation of PAO2 only specifies the fixed variables that are used in each of the two lower-level problems. PAO analyzes the use of decision variables in Pyomo models, and treats unused variables as fixed. Thus, the Pyomo and MPR representations generate a consistent interpretation of the variable specifications. However, the MPR representation is more explicit in this regard.
We consider again the trilevel problem described by Anadalingam (3), which can be expressed as a trilevel linear problem as follows:
>>> M = pao.mpr.LinearMultilevelProblem()
>>> U = M.add_upper(nxR=1)
>>> U.x.lower_bounds = [0]
>>> L = U.add_lower(nxR=1)
>>> L.x.lower_bounds = [0]
>>> B = L.add_lower(nxR=1)
>>> B.x.lower_bounds = [0]
>>> B.x.upper_bounds = [0.5]
>>> U.minimize = True
>>> U.c[U] = [-7]
>>> U.c[L] = [-3]
>>> U.c[B] = [4]
>>> L.minimize = True
>>> L.c[L] = [-1]
>>> B.minimize = True
>>> B.c[B] = [-1]
>>> B.inequalities = True
>>> B.A[U] = [[1], [ 1], [-1], [-1]]
>>> B.A[L] = [[1], [ 1], [-1], [ 1]]
>>> B.A[B] = [[1], [-1], [-1], [ 1]]
>>> B.b = [3,1,-1,1]
The QuadraticMultilevelProblem
class can represent general
quadratic problems with quadratic terms in the objective and constraints
at each level. The special case where bilinear terms arise with an
upper-level binary variable multiplied with a lower-level variable is
common in many applications. For this case, PAO provides a function to
linearize these bilinear terms.
We consider again the bilevel problem PAO3 (4), which is represented and solved as follows:
>>> M = pao.mpr.QuadraticMultilevelProblem()
>>> U = M.add_upper(nxR=2, nxB=2)
>>> L = U.add_lower(nxR=1)
>>> U.x.lower_bounds = [2, np.NINF, 0, 0]
>>> U.x.upper_bounds = [6, np.PINF, 1, 1]
>>> L.x.lower_bounds = [0]
>>> L.x.upper_bounds = [np.PINF]
>>> U.c[U] = [1, 0, 5, 0]
>>> U.c[L] = [3]
>>> L.c[L] = [1]
>>> L.maximize = True
>>> U.A[U] = [[ 1, 1, 0, 0],
... [-1, -1, 0, 0],
... [ 0, 0, -1, -1]
... ]
>>> U.b = [10, -10, -1]
>>> L.A[U] = [[ 1, 0, 0, 0],
... [-1, 0, 0, 0],
... [ 1, 0, 0, 0]]
>>> L.A[L] = [[ 0],
... [-4],
... [ 0]]
>>> L.Q[U,L] = (3,4,1), {(0,2,0):1, (2,3,0):2}
>>> L.b = [8, -8, 13]
>>> lmr, soln = pao.mpr.linearize_bilinear_terms(M, 100)
>>> opt = pao.Solver("pao.mpr.FA")
>>> results = opt.solve(lmr)
>>> soln.copy(From=lmr, To=M)
>>> print(U.x.values)
[6.0, 4.0, 0, 1]
>>> print(L.x.values)
[3.5]
The data L.Q[U,L]
specifies the bilinear terms multiplying
variables from level U
with variables from level L
, which
are included in the constraints in level L
. Note that Q
is a
tensor, which is indexed over the constraints, upper-level variables
and lower-level variables. A similar syntax is used to define bilinear
terms in objectives, P
, though that is represented as a sparse matrix.
Quadratic terms can be specified simply by using the same levels to index
Q
or P
.
Model transformations like :func:.linearize_bilinear_terms
are
described in further detail in the next section. Note that this function
returns both the transformed model as well as a helper class that maps
solutions back to the original model. This logic facilitates the
automation of model transformations within PAO.
Solvers¶
After formulating a multilevel problem, PAO users will generally need to (1) transform the model to a standard form, and (2) apply an optimizer to solve the problem. The examples in the previous sections illustrate that step (1) is often optional; PAO automates the applications of several model transformations, particularly for problems formulated with Pyomo. The following section summarizes the solvers available in PAO, and describes how PAO manages solvers. Section Model Transformations describes model transformations in PAO.
Summary of PAO Solvers¶
The following summarizes the current solvers available in PAO:
pao.mpr.FA, pao.pyomo.FA
PAO solver for Multilevel Problem Representations that define linear bilevel problems. Solver uses big-M relaxations discussed by Fortuny- Amat and McCarl (1981).
pao.mpr.MIBS, pao.pyomo.MIBS
PAO solver for Multilevel Problem Representations using the COIN-OR MibS solver by Tahernejad, Ralphs, and DeNegre (2020).
pao.mpr.PCCG, pao.pyomo.PCCG
PAO solver for Multilevel Problem Representations that define linear bilevel problems. Solver uses projected column constraint generation algorithm described by Yue et al. (2017).
pao.mpr.REG, pao.pyomo.REG
PAO solver for Multilevel Problem Representations that define linear bilevel problems. Solver uses regularization discussed by Scheel and Scholtes (2000) and Ralph and Wright (2004).
The following table summarize key features of the problems these solvers can be applied to:
Solver | |||||
Problem Feature | FA | REG | PCCG | MibS | |
Equation Structure | Linear | Y | Y | Y | Y |
Bilinear | Y | Y | Y | Y | |
Nonlinear | |||||
Upper-Level Variables | Integer | Y | Y | Y | |
Real | Y | Y | Y | Y | |
Lower-Level Variables | Integer | Y | Y | ||
Real | Y | Y | Y | Y | |
Multilevel Representation | Bilevel | Y | Y | Y | Y |
Trilevel | |||||
k-Bilevel | Y | Y |
Note
The iterface to MibS is a prototype that has not been well-tested. This interface will be documented and finalized in an upcoming release of PAO.
The Solver Interface¶
The Solver
object provides a single interface for setting up
an interface to optimizers in PAO. This includes both PAO solvers for
multilevel optimization problems, but also interfaces to conventional
numerical solvers that are used by PAO solvers. We illustrate this
distinction with the following example, which optimizes the PAO1
(1) example:
>>> # Create an interface to the PAO FA solver
>>> opt = pao.Solver("pao.pyomo.FA")
>>> # Optimize the model
>>> # By default, FA uses the glpk MIP solver
>>> results = opt.solve(M)
>>> print(M.x.value, M.y.value, M.L.z.value)
6.0 4.0 2.0
>>> # Create an interface to the PAO FA solver, using cbc
>>> opt = pao.Solver("pao.pyomo.FA", mip_solver="cbc")
>>> # Optimize the model using cbc
>>> results = opt.solve(M)
>>> print(M.x.value, M.y.value, M.L.z.value)
6.0 4.0 2.0
The Solver
object is initialized using the solver name followed
by solver-specific options. In this case, the FA algorithm accepts
the mip_solver
option that specifies the mixed-integer programming
(MIP) solver that is used to solve the MIP that is generated by FA after
reformulating the bilevel problem. The value of mip_solver
is itself
an optimizer. As illustrated here, this option can simply be the string
name of the MIP solver that will be used. However, the Solver
object can be used to define a MIP solver interface as well:
>>> # Create an interface to the cbc MIP solver
>>> mip = pao.Solver("cbc")
>>> # Create an interface to the PAO FA solver, using cbc
>>> opt = pao.Solver("pao.pyomo.FA", mip_solver=mip)
>>> # Optimize the model using cbc
>>> results = opt.solve(M)
>>> print(M.x.value, M.y.value, M.L.z.value)
6.0 4.0 2.0
This enables the customization of the MIP solver used by FA. Note that
the solve()
method accepts the same options as Solve
.
This allows for more dynamic specification of solver options:
>>> # Create an interface to the cbc MIP solver
>>> cbc = pao.Solver("cbc")
>>> # Create an interface to the glpk MIP solver
>>> glpk = pao.Solver("glpk")
>>> # Create an interface to the PAO FA solver
>>> opt = pao.Solver("pao.pyomo.FA")
>>> # Optimize the model using cbc
>>> results = opt.solve(M, mip_solver=cbc)
>>> print(M.x.value, M.y.value, M.L.z.value)
6.0 4.0 2.0
>>> # Optimize the model using glpk
>>> results = opt.solve(M, mip_solver=glpk)
>>> print(M.x.value, M.y.value, M.L.z.value)
6.0 4.0 2.0
Warning
The solve()
current passes unknown keyword arguments to the
optimizer used by PAO solvers, but this feature will be disabled.
PAO Solvers¶
Solvers developed in PAO have names that begin with pao.
.
The current set of available PAO solvers can be queried using the
Solver
object:
>>> for name in pao.Solver:
... print(name)
pao.mpr.FA
pao.mpr.MIBS
pao.mpr.PCCG
pao.mpr.REG
pao.pyomo.FA
pao.pyomo.MIBS
pao.pyomo.PCCG
pao.pyomo.REG
>>> pao.Solver.summary()
pao.mpr.FA
PAO solver for Multilevel Problem Representations that define linear
bilevel problems. Solver uses big-M relaxations discussed by Fortuny-
Amat and McCarl (1981).
pao.mpr.MIBS
PAO solver for Multilevel Problem Representations using the COIN-OR
MibS solver by Tahernejad, Ralphs, and DeNegre (2020).
pao.mpr.PCCG
PAO solver for Multilevel Problem Representations that define linear
bilevel problems. Solver uses projected column constraint generation
algorithm described by Yue et al. (2017).
pao.mpr.REG
PAO solver for Multilevel Problem Representations that define linear
bilevel problems. Solver uses regularization discussed by Scheel and
Scholtes (2000) and Ralph and Wright (2004).
pao.pyomo.FA
PAO solver for Pyomo models that define linear and bilinear bilevel
problems. Solver uses big-M relaxations discussed by Fortuny-Amat and
McCarl (1981).
pao.pyomo.MIBS
PAO solver for Multilevel Problem Representations using the COIN-OR
MibS solver by Tahernejad, Ralphs, and DeNegre (2020).
pao.pyomo.PCCG
PAO solver for Pyomo models that define linear and bilinear bilevel
problems. Solver uses projected column constraint generation
algorithm described by Yue et al. (2017)
pao.pyomo.REG
PAO solver for Pyomo models that define linear and bilinear bilevel
problems. Solver uses regularization discussed by Scheel and Scholtes
(2000) and Ralph and Wright (2004).
The solve()
method includes documentation describing the keyword
arguments for a specific solver. For example:
>>> opt = pao.Solver("pao.pyomo.FA")
>>> help(opt.solve)
Help on method solve in module pao.pyomo.solvers.mpr_solvers:
solve(model, **options) method of pao.pyomo.solvers.mpr_solvers.PyomoSubmodelSolver_FA instance
Executes the solver and loads the solution into the model.
Parameters
----------
model
The model that is being optimized.
options
Keyword options that are used to configure the solver.
Keyword Arguments
-----------------
tee
If True, then solver output is streamed to stdout. (default is False)
load_solutions
If True, then the finale solution is loaded into the model. (default is True)
linearize_bigm
The name of the big-M value used to linearize bilinear terms. If this is not specified, then the solver will throw an error if bilinear terms exist in the model.
mip_solver
The MIP solver used by FA. (default is glpk)
Returns
-------
Results
A summary of the optimization results.
The solve()
method returns a results object that contains
data about the optimization process. In particular, this object
contains information about the termination conditions for the solver.
The check_optimal_termination()
method can be used confirm that the
termination condition indicates that an optimal solution was found. For example:
>>> nlp = pao.Solver('ipopt', print_level=3)
>>> opt = pao.Solver('pao.pyomo.REG', nlp_solver=nlp)
>>> results = opt.solve(M)
>>> print(results.solver.termination_condition)
TerminationCondition.optimal
>>> results.check_optimal_termination()
True
Pyomo Solvers¶
The Solver
object also provides a convenient interface to
conventional numerical solvers. Currently, solver objects constructed
by Solver
are simple wrappers around Pyomo optimization
solver objects. This interface supports two types of solver
interfaces: (1) solvers that execute locally, and (2) solvers that execute
on remote servers.
When optimizating a Pyomo model, solver parameters can be setup both when the solver interface is created and when a model is optimized. For example:
>>> # This is a nonlinear toy problem modeled with Pyomo
>>> NLP = pe.ConcreteModel()
>>> A = list(range(10))
>>> NLP.x = pe.Var(A, bounds=(0,None), initialize=1)
>>> NLP.o = pe.Objective(expr=sum(pe.sin((i+1)*NLP.x[i]) for i in A))
>>> NLP.c = pe.Constraint(expr=sum(NLP.x[i] for i in A) >= 1)
>>> nlp = pao.Solver('ipopt', print_level=3)
>>> # Apply ipopt with print level 3
>>> results = nlp.solve(NLP)
>>> # Override the default print level to using 5
>>> results = nlp.solve(NLP, print_level=5)
However, PAO users will typically setup solver parameters when the Pyomo solver is initially created:
>>> nlp = pao.Solver('ipopt', print_level=3)
>>> opt = pao.Solver('pao.pyomo.REG', nlp_solver=nlp)
>>> results = opt.solve(M)
When executing locally, the executable
option can be used
to explicitly specify the path to the executable that is used by this solver.
This is helpful in contexts where Pyomo is not automatically finding the correct
optimizer executable in a user’s shell environment.
When executing on a remote server, the server
is used to
specify the server that is used. Currently, only the neos
server is
supported, which allows the user to perform optimization at NEOS [NEOS].
The NEOS server requires a user to specify a valid email address:
Warning
There is no common reference for solver-specific parameters for the solvers available in Pyomo. These are generally documented with solver documentation, and users should expect to contact solver developers to learn about these.
Model Transformations¶
PAO includes a variety of functions that transform models, which generally are applied as follows:
Here, the function transform_func
generates the model new_model
from the model old_model
. The
object soln
is used to map a solution back to the old model:
The following transformation functions are documented and suitable for use by end-users:
pao.pyomo.convert.convert_pyomo2MultilevelProblem()
This function generates a
LinearMultilevelProblem
orQuadraticMultilevelProblem
from a Pyomo model. By default, all constraints in the MPR representation are inequalities.pao.mpr.convert_repn.linearize_bilinear_terms()
This function generates a
LinearMultilevelProblem
from aQuadraticMultilevelProblem
that only contains bilinear terms. This transformation currently is limited to MPRs that only contain inequality constraints.pao.mpr.convert_repn.convert_to_standard_form()
This function generates an equivalent linear multilevel representation for which all variables are non-negative and all constraints have the same form (inequalities or equalities). This simplifies the implementation of solvers, which typically assume a standard form for subproblems.
Library Reference¶
The following classes and functions represent the core functionality in PAO:
Warning
The logic in pao.duality
is currently disabled. There are known errors in this code
that will be resolved by re-implementing it using the logic in pao.mpr
.
Bibliography¶
[Anadalingam] | G. Anandalingam, A mathematical programming model of decentralized multi-level systems. J. Oper.Res. Soc.39(11), 1021–1033 (1988) |
[Bard98] | J. F. Bard, Practical bilevel optimization: Algorithms and applications. Kluwer Academic Publishers, 1998. |
[FortunyMcCarl] | (1, 2) J. Fortuny-Amat and B. McCarl, A representation and economic interpretation of a two-level programming problem, The Journal of the Operations Research Society. 32(9), 783-792, 1981. |
[NEOS] | https://neos-server.org/neos/ |
[PyomoBookII] | (1, 2) W. E. Hart, C. D. Laird, J.-P. Watson, D. L. Woodruff, G. A. Hackebeil, B. L. Nicholson, J. D. Siirola. Pyomo - Optimization Modeling in Python, 2nd Edition. Springer Optimization and Its Applications, Vol 67. Springer, 2017. |
[PyomoWeb] | http://www.pyomo.org/ |
[PyomoGithub] | https://github.com/Pyomo/pyomo |