Getting Started with Coopr 3.1

1. Preface

This book provides a quick introduction to Coopr. Coopr is a collection of Python software packages that supports a diverse set of optimization capabilities for formulating and analyzing optimization models. A central component of Coopr is Pyomo, which supports the formulation and analysis of mathematical models for complex optimization applications. This capability is commonly associated with algebraic modeling languages (AMLs), which support the description and analysis of mathematical models with a high-level language. Although most AMLs are implemented in custom modeling languages, Pyomo’s modeling objects are embedded within Python, a full-featured high-level programming language that contains a rich set of supporting libraries.

Coopr has also proven an effective framework for developing high-level optimization and analysis tools. For example, the PySP package provides generic solvers for stochastic programming. PySP leverages the fact that Pyomo’s modeling objects are embedded within a full-featured high-level programming language, which allows for transparent parallelization of subproblems using Python parallel communication libraries.

1.1. Goals of the Book

In this book, we provide a broad overview of different components of the Coopr software. There are roughly two main goals for this book:

Help users get started with different Coopr capabilities: Our goal is not to provide a comprehensive reference, but rather to provide a tutorial with simple and illustrative examples. Also, we aim to provide explanations behind the design and philosophy of Coopr.
Provide preliminary documentation of new features and capabilities in Coopr. We know that a new feature or capability probably will not be used unless it is documented. As Coopr evolves, we plan to use this book to document these features. This provides users some context concerning the focus of Coopr development, and it also provides an opportunity to get early feedback on new features before they are documented in other contexts.

1.2. Who Should Read This Book

This book is intended to be a reference for students, academic researchers and practitioners. Coopr has been effectively used in the classroom with undergraduate and graduate students. However, we assume that the reader is generally familiar with optimization and mathematical modeling. Although this book does not contain a glossary, we recommend the Mathematical Programming Glossary [MPG] as a reference for the reader. We also assume that the reader is generally familiar with the Python programming language. There are a variety of books describing Python, as well as excellent documentation of the Python language and the software packages bundled with Python distributions.

1.3. Comments and Questions

Further information about Pyomo and Coopr is available on the Coopr wiki:

https://software.sandia.gov/trac/coopr

Coopr is also hosted at COIN-OR:

https://projects.coin-or.org/Coopr

We strongly encourage feedback from readers about the software on the Coopr Forum:

coopr-forum@googlegroups.com

If you have comments about this pre-alpha rough draft of the book, please send them to DLWoodruff@UCDavis.edu. We hope this will include feedback on typos and errors in our examples. Additionally, we welcome comments on the presentation of this material, and suggestions for material that we should develop in the other book chapters.

Good Luck!

2. Overview

This overview chapter outlines the book.

3. Pyomo Overview

3.1. Mathematical Modeling

This chapter provides an introduction to Pyomo: Python Optimization Modeling Objects. A more complete description is contained in <[PyomoBook]>. Pyomo supports the formulation and analysis of mathematical models for complex optimization applications. This capability is commonly associated with algebraic modeling languages (AMLs) such as AMPL [AMPL] AIMMS [AIMMS] and GAMS [GAMS]. Pyomo’s modeling objects are embedded within Python, a full-featured high-level programming language that contains a rich set of supporting libraries.

Modeling is a fundamental process in many aspects of scientific research, engineering and business. Modeling involves the formulation of a simplified representation of a system or real-world object. Thus, modeling tools like Pyomo can be used in a variety of ways[Sch04]:

Explain phenomena that arise in a system,
Make predictions about future states of a system,
Assess key factors that influence phenomena in a system,
Identify extreme states in a system, that might represent worst-case scenarios or minimal cost plans, and
Analyze trade-offs to support human decision makers.

Mathematical models represent system knowledge with a formalized mathematical language. The following mathematical concepts are central to modern modeling activities:

variables: Variables represent unknown or changing parts of a model (e.g. whether or not to make a decision, or the characteristic of a system outcome). The values taken by the variables are often referred to as a solution and are usually an output of the optimization process.
parameters: Parameters represents the data that must be supplied to perform the optimization. In fact, in some settings the word data is used in place of the word parameters.
relations: These are equations, inequalities or other mathematical relationships that define how different parts of a model are connected to each other.
goals: These are functions that reflect goals and objectives for the system being modeled.

The widespread availability of computing resources has made the numerical analysis of mathematical models a commonplace activity. Without a modeling language, the process of setting up input files, executing a solver and extracting the final results from the solver output is tedious and error prone. This difficulty is compounded in complex, large-scale real-world applications which are difficult to debug when errors occur. Additionally, there are many different formats used by optimization software packages, and few formats are recognized by many optimizers. Thus the application of multiple optimization solvers to analyze a model introduces additional complexities.

Pyomo is an AML that extends Python to include objects for mathematical modeling. Hart et al. [] compare Pyomo with other AMLs. Although many good AMLs have been developed for optimization models, the following are motivating factors for the development of Pyomo:

Open Source: Pyomo is developed within Coopr’s open source project to promote transparency of the modeling framework and encourage community development of Pyomo capabilities.
Customizable Capability: Pyomo supports a customizable capability through the extensive use of plug-ins to modularize software components.
Solver Integration: Pyomo models can be optimized with solvers that are written either in Python or in compiled, low-level languages.
Programming Language: Pyomo leverages a high-level programming language, which has several advantages over custom AMLs: a very robust language, extensive documentation, a rich set of standard libraries, support for modern programming features like classes and functions, and portability to many platforms.

3.2. Overview of Modeling Components and Processes

Pyomo supports an object-oriented design for the definition of optimization models. The basic steps of a simple modeling process are:

Create model and declare components
Instantiate the model
Apply solver
Interrogate solver results

In practice, these steps may be applied repeatedly with different data or with different constraints applied to the model. However, we focus on this simple modeling process to illustrate different strategies for modeling with Pyomo.

A Pyomo model consists of a collection of modeling components that define different aspects of the model. Pyomo includes the modeling components that are commonly supported by modern AMLs: index sets, symbolic parameters, decision variables, objectives, and constraints. These modeling components are defined in Pyomo through the following Python classes:

Set: set data that is used to define a model instance
Param: parameter data that is used to define a model instance
Var: decision variables in a model
Objective: expressions that are minimized or maximized in a model
Constraint: constraint expressions that impose restrictions on variable values in a model

3.3. Abstract Versus Concrete Models

A mathematical model can be defined using symbols that represent data values. For example, the following equations represent a linear program (LP) to find optimal values for the vector $x$ with parameters $n$ and $b$, and parameter vectors $a$ and $c$:

$\begin{array}{lll} \min & \sum_{j=1}^n c_j x_j &\\ \mathrm{s.t.} & \sum_{j=1}^n a_{ij} x_j \geq b_i & \forall i = 1 \ldots m\\ & x_j \geq 0 & \forall j = 1 \ldots n \end{array}$

Note	As a convenience, we use the symbol $\forall$ to mean “for all” or “for each.”

We call this an abstract or symbolic mathematical model since it relies on unspecified parameter values. Data values can be used to specify a model instance. The AbstractModel class provides a context for defining and initializing abstract optimization models in Pyomo when the data values will be supplied at the time a solution is to be obtained.

In some contexts a mathematical model can be directly defined with the data values supplied at the time of the model definition and built into the model. We call these concrete mathematical models. For example, the following LP model is a concrete instance of the previous abstract model:

$\begin{array}{ll} \min & 2 x_1 + 3 x_2\\ \mathrm{s.t.} & 3 x_1 + 4 x_2 \geq 1\\ & x_1, x_2 \geq 0 \end{array}$

The ConcreteModel class is used to define concrete optimization models in Pyomo.

3.4. A Simple Abstract Pyomo Model

We repeat the abstract model already given:

$\begin{array}{lll} \min & \sum_{j=1}^n c_j x_j &\\ \mathrm{s.t.} & \sum_{j=1}^n a_{ij} x_j \geq b_i & \forall i = 1 \ldots m\\ & x_j \geq 0 & \forall j = 1 \ldots n \end{array}$

One way to implement this in Pyomo is as follows:

from coopr.pyomo import *

model = AbstractModel()

model.m = Param(within=NonNegativeIntegers)
model.n = Param(within=NonNegativeIntegers)

model.I = RangeSet(1, model.m)
model.J = RangeSet(1, model.n)

model.a = Param(model.I, model.J)
model.b = Param(model.I)
model.c = Param(model.J)

# the next line declares a variable indexed by the set J
model.x = Var(model.J, domain=NonNegativeReals)

def obj_expression(model):
    return summation(model.c, model.x)

model.OBJ = Objective(rule=obj_expression)

def ax_constraint_rule(model, i):
    # return the expression for the constraint for i
    return sum(model.a[i,j] * model.x[j] for j in model.J) >= model.b[i]

# the next line creates one constraint for each member of the set model.I
model.AxbConstraint = Constraint(model.I, rule=ax_constraint_rule)

Note

Python is interpreted one line at a time and if statements need to span a line for some reason then a line continuation character, which is backslash, must be used. In Python, indentation has meaning and must be consistent. Lines inside a function definition, for example, must be indented and the end of the indentation is used by Python to signal the end of the definition.

We will now examine the lines in this example beginning with the import line that is required in every Pyomo model. Its purpose is to make the symbols used by Pyomo known to Python.

from coopr.pyomo import *

The declaration of a model is also required. The use of the name model is not required. Almost any name could be used, but we will use the name model most of the time in this book. In this example, we are declaring that it will be an abstract model.

model = AbstractModel()

We declare the parameters $m$ and $n$ using the Pyomo Param function. This function can take a variety of arguments; this example illustrates use of the within option that is used by Pyomo to validate the data value that is assigned to the parameter. If this option were not given, then Pyomo would not object to any type of data being assigned to these parameters. As it is, assignment of a value that is not a non-negative integer will result in an error.

model.m = Param(within=NonNegativeIntegers)
model.n = Param(within=NonNegativeIntegers)

Although not required, it is convenient to define index sets. In this example we use the RangeSet function to declare that the sets will be a sequence of integers starting at 1 and ending at a value specified by the the parameters model.m and model.n.

model.I = RangeSet(1, model.m)
model.J = RangeSet(1, model.n)

The coefficient and right-hand-side data are defined as indexed parameters. When sets are given as arguments to the Param function, they indicate that the set will index the parameter.

model.a = Param(model.I, model.J)
model.b = Param(model.I)
model.c = Param(model.J)

Note	In Python, and therefore in Pyomo, any text after pound sign is considered to be a comment.

The next line interpreted by Python as part of the model declares the variable $x$. The first argument to the Var function is a set, so it is defined as an index set for the variable. In this case the variable has only one index set, but multiple sets could be used as was the case for the declaration of the parameter model.a. The second argument specifies a domain for the variable. This information is part of the model and will passed to the solver when data is provided and the model is solved. Specification of the NonNegativeReals domain implements the requirement that the variables be greater than or equal to zero.

# the next line declares a variable indexed by the set J
model.x = Var(model.J, domain=NonNegativeReals)

In abstract models, Pyomo expressions are usually provided to objective function and constraint declarations via a function defined with a Python def statement. The def statement establishes a name for a function along with its arguments. When Pyomo uses a function to get objective function or constraint expressions, it always passes in the model (i.e., itself) as the the first argument so the model is always the first formal argument when declaring such functions in Pyomo. Additional arguments, if needed, follow. Since summation is an extremely common part of optimization models, Pyomo provides a flexible function to accommodate it. When given two arguments, the summation function returns an expression for the sum of the product of the two arguments over their indexes. This only works, of course, if the two arguments have the same indexes. If it is given only one argument it returns an expression for the sum over all indexes of that argument. So in this example, when summation is passed the arguments model.c, model.x it returns an internal representation of the expression $\sum_{j=1}^{n}c_{j} x_{j}$.

def obj_expression(model):
    return summation(model.c, model.x)

To declare an objective function, the Pyomo function called Objective is used. The rule argument gives the name of a function that returns the expression to be used. The default sense is minimization. For maximization, the sense=maximize argument must be used. The name that is declared, which is OBJ in this case, appears in some reports and can be almost any name.

model.OBJ = Objective(rule=obj_expression)

Declaration of constraints is similar. A function is declared to deliver the constraint expression. In this case, there can be multiple constraints of the same form because we index the constraints by $i$ in the expression $\sum_{j=1}^n a_{ij} x_j \geq b_i \;\;\forall i = 1 \ldots m$, which states that we need a constraint for each value of $i$ from one to $m$. In order to parametrize the expression by $i$ we include it as a formal parameter to the function that declares the constraint expression. Technically, we could have used anything for this argument, but that might be confusing. Using an i for an $i$ seems sensible in this situation.

def ax_constraint_rule(model, i):
    # return the expression for the constraint for i
    return sum(model.a[i,j] * model.x[j] for j in model.J) >= model.b[i]

Note	In Python, indexes are in square brackets and function arguments are in parentheses.

In order to declare constraints that use this expression, we use the Pyomo Constraint function that takes a variety of arguments. In this case, our model specifies that we can have more than one constraint of the same form and we have created a set, model.I, over which these constraints can be indexed so that is the first argument to the constraint declaration function. The next argument gives the rule that will be used to generate expressions for the constraints. Taken as a whole, this constraint declaration says that a list of constraints indexed by the set model.I will be created and for each member of model.I, the function ax_constraint_rule will be called and it will be passed the model object as well as the member of model.I.

# the next line creates one constraint for each member of the set model.I
model.AxbConstraint = Constraint(model.I, rule=ax_constraint_rule)

In the object oriented view of all of this, we would say that model object is a class instance of the AbstractModel class, and model.J is a Set object that is contained by this model. Many modeling components in Pyomo can be optionally specified as indexed components: collections of components that are referenced using one or more values. In this example, the parameter model.c is indexed with set model.J.

In order to use this model, data must be given for the values of the parameters. Here is one file that provides data.

# one way to input the data in AMPL format
# for indexed parameters, the indexes are given before the value

param m := 1 ;
param n := 2 ;

param a :=
1 1 3
1 2 4
;

param c:=
1 2
2 3
;

param b := 1 1 ;

There are multiple formats that can be used to provide data to a Pyomo model, but the AMPL format works well for our purposes because it contains the names of the data elements together with the data. In AMPL data files, text after a pound sign is treated as a comment. Lines generally do not matter, but statements must be terminated with a semi-colon.

For this particular data file, there is one constraint, so the value of model.m will be one and there are two variables (i.e., the vector model.x is two elements long) so the value of model.n will be two. These two assignments are accomplished with standard assignments. Notice that in AMPL format input, the name of the model is omitted.

param m := 1 ;
param n := 2 ;

There is only one constraint, so only two values are needed for model.a. When assigning values to arrays and vectors in AMPL format, one way to do it is to give the index(es) and the the value. The line 1 2 4 causes model.a[1,2] to get the value 4. Since model.c has only one index, only one index value is needed so, for example, the line 1 2 causes model.c[1] to get the value 2. Line breaks generally do not matter in AMPL format data files, so the assignment of the value for the single index of model.b is given on one line since that is easy to read.

param a :=
1 1 3
1 2 4
;

param c:=
1 2
2 3
;

param b := 1 1 ;

When working with Pyomo (or any other AML), it is convenient to write abstract models in a somewhat more abstract way by using index sets that contain strings rather than index sets that are implied by $1,\ldots,m$ or the summation from 1 to $n$. When this is done, the size of the set is implied by the input, rather than specified directly. Furthermore, the index entries may have no real order. Often, a mixture of integers and indexes and strings as indexes is needed in the same model. To start with an illustration of general indexes, consider a slightly different Pyomo implementation of the model we just presented.

# abstract2.py

from coopr.pyomo import *

model = AbstractModel()

model.I = Set()
model.J = Set()

model.a = Param(model.I, model.J)
model.b = Param(model.I)
model.c = Param(model.J)

# the next line declares a variable indexed by the set J
model.x = Var(model.J, domain=NonNegativeReals)

def obj_expression(model):
    return summation(model.c, model.x)

model.OBJ = Objective(rule=obj_expression)

def ax_constraint_rule(model, i):
    # return the expression for the constraint for i
    return sum(model.a[i,j] * model.x[j] for j in model.J) >= model.b[i]

# the next line creates one constraint for each member of the set model.I
model.AxbConstraint = Constraint(model.I, rule=ax_constraint_rule)

To get the same instantiated model, the following data file can be used.

# abstract2a.dat AMPL format

set I := 1 ;
set J := 1 2 ;

param a :=
1 1 3
1 2 4
;

param c:=
1 2
2 3
;

param b := 1 1 ;

However, this model can also be fed different data for problems of the same general form using meaningful indexes.

# abstract2.dat AMPL data format

set I := TV Film ;
set J := Graham John Carol ;

param a :=
TV  Graham 3
TV John 4.4
TV Carol 4.9
Film Graham 1
Film John 2.4
Film Carol 1.1
;

param c := [*]
 Graham 2.2
 John 3.1416
 Carol 3
;

param b := TV 1 Film 1 ;

3.5. A Simple Concrete Pyomo Model

It is possible to get nearly the same flexible behavior from models declared to be abstract and models declared to be concrete in Pyomo; however, we will focus on a straightforward concrete example here where the data is hard-wired into the model file. Python programmers will quickly realize that the data could have come from other sources.

We repeat the concrete model already given:

$\begin{array}{ll} \min & 2 x_1 + 3 x_2\\ \mathrm{s.t.} & 3 x_1 + 4 x_2 \geq 1\\ & x_1, x_2 \geq 0 \end{array}$

This is implemented as a concrete model as follows:

from coopr.pyomo import *

model = ConcreteModel()

model.x = Var([1,2], domain=NonNegativeReals)

model.OBJ = Objective(expr = 2*model.x[1] + 3*model.x[2])

model.Constraint1 = Constraint(expr = 3*model.x[1] + 4*model.x[2] >= 1)

Although rule functions can also be used to specify constraints and objectives, in this example we use the expr option that is available only in concrete models. This option gives a direct specification of the expression.

3.6. Solving the Simple Examples

Pyomo supports modeling and scripting but does not install a solver automatically. In order to solve a model, there must be a solver installed on the computer to be used. If there is a solver, then the pyomo command can be used to solve a problem instance.

Suppose that the solver named glpk (also known as glpsol) is installed on the computer. Suppose further that an abstract model is in the file named abstract1.py and a data file for it is in the file named abstract1.dat. From the command prompt, with both files in the current directory, a solution can be obtained with the command:

pyomo abstract1.py abstract1.dat --solver=glpk

Since glpk is the default solver, there really is no need specify it so the --solver option can be dropped.

Note	There are two dashes before the command line option names such as `solver`.

To continue the example, if CPLEX is installed then it can be listed as the solver. The command to solve with CPLEX is

pyomo abstract1.py abstract1.dat --solver=cplex

This yields the following output on the screen:

[    0.00] Setting up Pyomo environment
[    0.00] Applying Pyomo preprocessing actions
[    0.07] Creating model
[    0.15] Applying solver
[    0.37] Processing results
    Number of solutions: 1
    Solution Information
      Gap: 0.0
      Status: optimal
      Function Value: 0.666666666667
    Solver results file: results.json
[    0.39] Applying Pyomo postprocessing actions
[    0.39] Pyomo Finished

The numbers is square brackets indicate how much time was required for each step. Results are written to the file named results.json, which has a special structure that makes it useful for post-processing. To see a summary of results written to the screen, use the --summary option:

pyomo abstract1.py abstract1.dat --solver=cplex --summary

To see a list of Pyomo command line options, use:

pyomo --help

Note	There are two dashes before `help`.

For a concrete model, no data file is specified on the Pyomo command line.

4. Sets

Sets can be declared using the Set and RangeSet functions or by assigning set expressions. The simplest set declaration creates a set and postpones creation of its members:

model.A = Set()

The Set function takes optional arguments such as:

doc = String describing the set
dimen = Dimension of the members of the set
filter = A boolean function used during construction to indicate if a potential new member should be assigned to the set
initialize = A function that returns the members to initialize the set. ordered = A boolean indicator that the set is ordered; the default is False
validate = A boolean function that validates new member data
virtual = A boolean indicator that the set will never have elements; it is unusual for a modeler to create a virtual set; they are typically used as domains for sets, parameters and variables
within = Set used for validation; it is a super-set of the set being declared.

To create a set whose members will be two dimensional, use

model.B = Set(dimen=2)

To create a set of all the numbers in set model.A doubled, one could use

def doubleA_init(model):
    return (i*2 for i in model.A)
model.C = Set(initialize=DoubleA_init)

As an aside we note that as always in Python, there are lot of ways to accomplish the same thing. Also, note that this will generate an error if model.A contains elements for which multiplication times two is not defined.

The initialize option can refer to a Python set, which can be returned by a function or given directly as in

model.D = Set(initialize=['red', 'green', 'blue'])

The initialize option can also specify a function that is applied sequentially to generate set members. Consider the case of a simple set. In this case, the initialization function accepts a set element number and model and returns the set element associated with that number:

def Z_init(model, i):
    if i > 10:
        return Set.End
    return 2*i+1
model.Z = Set(initialize=Z_init)

The Set.End return value terminates input to the set. Additional information about iterators for set initialization is in [PyomoBook].

Note	Data specified in an input file will override the data specified by the initialize options.

If sets are given as arguments to Set without keywords, they are interpreted as indexes for an array of sets. For example, to create an array of sets that is indexed by the members of the set model.A, use

model.E = Set(model.A)

Arguments can be combined. For example, to create an array of sets with three dimensional members indexed by set model.A, use

model.F = Set(model.A, dimen=3)

The initialize option can be used to create a set that contains a sequence of numbers, but the RangeSet function provides a concise mechanism for simple sequences. This function takes as its arguments a start value, a final value, and a step size. If the RangeSet has only a single argument, then that value defines the final value in the sequence; the first value and step size default to one. If two values given, they are the first and last value in the sequence and the step size defaults to one. For example, the following declaration creates a set with the numbers 1.5, 5 and 8.5:

model.G = RangeSet(1.5, 10, 3.5)

Sets may also be created by assigning other Pyomo sets as in these examples that also illustrate the set operators union, intersection, difference, and exclusive-or:

model.H = model.A
model.I = model.A | model.D # union
model.J = model.A & model.D # intersection
model.K = model.A - model.D # difference
model.L = model.A ^ model.D # exclusive-or

4.1. Predefined Virtual Sets

For use in specifying domains for sets, parameters and variables, Pyomo provides the following pre-defined virtual sets:

Any: all possible values
Reals : floating point values
PositiveReals: strictly positive floating point values
NonPositiveReals: non-positive floating point values
NegativeReals: strictly negative floating point values
NonNegativeReals: non-negative floating point values
PercentFraction: floating point values in the interval [0,1]
Integers: integer values
PositiveIntegers: positive integer values
NonPositiveIntegers: non-positive integer values
NegativeIntegers: negative integer values
NonNegativeIntegers: non-negative integer values
Boolean: boolean values, which can be represented as False/True, 0/1, ’False’/’True’ and ’F’/’T’
Binary: same as boolean

For example, if the set model.M is declared to be within the virtual set NegativeIntegers then an attempt to add anything other than a negative integer will result in an error. Here is the declaration:

model.M = Set(within=NegativeIntegers)

4.2. Sparse Index Sets

It is common for an index set to be used to index variables and parameters in a constraint as well as to index other index sets. For example, one may want to have a constraint that holds

for i in model.I, k in model.K, v in model.V[k]

There are many ways to accomplish this, but one good way is to create a set of tuples composed of all of model.k, model.V[k] pairs. This can be done as follows:

def kv_init(model):
    return ((k,v) for k in model.K for v in model.V[k])
model.KV=Set(dimen=2, initialize=kv_init)

So then if there was a constraint defining rule such as

def MyC_rule(model, i, k, v):
   return ...

Then a constraint could be declared using

model.MyConstraint = Constraint(model.I,model.KV,rule=c1Rule)

Here is the first few lines of a model that illustrates this:

from coopr.pyomo import *

model = AbstractModel()

model.I=Set()
model.K=Set()
model.V=Set(model.K)

def kv_init(model):
    return ((k,v) for k in model.K for v in model.V[k])
model.KV=Set(dimen=2, initialize=kv_init)

model.a = Param(model.I, model.K)

model.y = Var(model.I)
model.x = Var(model.I, model.KV)

#include a constraint
#x[i,k,v] <= a[i,k]*y[i], for i in model.I, k in model.K, v in model.V[k]

def c1Rule(model,i,k,v):
   return model.x[i,k,v] <= model.a[i,k]*model.y[i]
model.c1 = Constraint(model.I,model.KV,rule=c1Rule)

5. Parameters

The word "parameters" is used in many settings. When discussing a Pyomo model, we use the word to refer to data that must be provided in order to find an optimal (or good) assignment of values to the decision variables. Parameters are declared with the Param function, which takes arguments that are very similar to the Set function. For example, the following code snippet declares sets model.A, model.B and then a parameter array model.P that is indexed by model.A:

model.A = Set()
model.B = Set()
model.P = Param(model.A, model.B)

In addition to sets that serve as indexes, the Param function takes the following command options:

default = The value absent any other specification.
doc = String describing the parameter
initialize = A function (or Python object) that returns the members to initialize the parameter values.
rule = (this is a synonym for initilize)
validate = A boolean function with arguments that are the prospective parameter value, the parameter indices and the model.
within = Set used for validation; it specifies the domain of the parameter values.

These options perform in the same way as they do for Set. For example, suppose that Model.A = RangeSet(1,3), then there are many ways to create a parameter that is a square matrix with 9, 16, 25 on the main diagonal zeros elsewhere, here are two ways to do it. First using a Python object to initialize:

v={}
v[1,1] = 9
v[2,2] = 16
v[3,3] = 25
model.S = Param(model.A, model.A, initialize=v, default=0)

And now using an initialization rule that is automatically called once for each index tuple (remember that we are assuming that model.A contains 1,2,3)

def s_init(model, i, j):
    if i == j:
        return i*i
    else:
        return 0.0
model.S = Param(model.A, model.A, rule=s_init)

In this example, the index set contained integers, but index sets need not be numeric. It is very common to use strings.

Note	Data specified in an input file will override the data specified by the initialize options.

6. Variables

Variables are intended to ultimately be given values by an optimization package. The are declared and optionally bounded, given initial values, and documented using the Pyomo Var function. If index sets are given as arguments to this function they are used to index the variable, other optional directives include:

bounds = A function (or Python object) that gives a (lower,upper) bound pair for the variable
domain = A set that is a super-set of the values the variable can take on.
initialize = A function (or Python object) that gives a starting value for the variable; this is particularly important for non-linear models
within = (synonym for domain)

The following code snippet illustrates some aspects of these options by declaring a singleton (i.e. unindexed) variable named model.LumberJack that will take on real values between zero and 6 and it initialized to be 1.5:

model.LumberJack = Var(within=NonNegativeReals, bounds=(0,6), initialize=1.5)

Instead of the initialize option, initialization is sometimes done with a Python assignment statement as in

model.LumberJack = 1.5

For indexed variables, bounds and initial values are often specified by a rule (a Python function) that itself may make reference to parameters or other data. The formal arguments to these rules begins with the model followed by the indexes. This is illustrated in the following code snippet that makes use of Python dictionaries declared as lb and ub that are used by a function to provide bounds:

model.A = Set(initialize=['Scones', 'Tea']
lb = {'Scones':2, 'Tea':4}
ub = {'Scones':5, 'Tea':7}
def fb(model, i):
   return (lower[i], upper[i])
model.PriceToCharge = Var(model.A, domain=PositiveInteger, bounds=fb)

Note	Many of the pre-defined virtual sets that are used as domains imply bounds. A strong example is the set `Boolean` that implies bounds of zero and one.

7. Objectives

An objective is a function of variables that returns a value that an optimization package attempts to maximize or minimize. The Objective function in Pyomo declares an objective. Although other mechanisms are possible, this function is typically passed the name of a another function that gives the expression. Here is a very simple version of such a function that assumes model.x has previously been declared as a Var:

def ObjRule(model):
    return 2*model.x[1] + 3*model.x[2]
model.g = Objective(rule=ObjRule)

It is more common for an objective function to refer to parameters as in this example that assumes that model.p has been declared as a parameters and that model.x has been declared with the same index set, while model.y has been declared as a singleton:

def profrul(model):
   return summation(model.p, model.x) + model.y
model.Obj = Objective(rule=ObjRule, sense=maximize)

This example uses the sense option to specify maximization. The default sense is minimize.

8. Constraints

Most constraints are specified using equality or inequality expressions that are created using a rule, which is a Python function. For example, if the variable model.x has the indexes butter and scones, then this constraint limits the sum for them to be exactly three:

def teaOKrule(model):
    return(model.x['butter'] + model.x['scones'] == 3)
model.TeaConst = Constraint, rule=teaOKrule)

Instead of expressions involving equality (==) or inequalities (<= or >=), constraints can also be expressed using a 3-tuple if the form (lb, expr, ub) where lb and ub can be None, which is interpreted as lb <= expr <= ub. Variables can appear only in the middle expr. For example, the following two constraint declarations have the same meaning:

model.x = Var()

def aRule(model):
   return model.x >= 2
Boundx = Constraint(rule=aRule)

def bRule(model):
   return (2, model.x, None)
Boundx = Constraint(rule=bRule)

For this simple example, it would also be possible to declare model.x with a bound option to accomplish the same thing.

Constraints (and objectives) can be indexed by lists or sets. When the declaration contains lists or sets as arguments, the elements are iteratively passed to the rule function. If there is more than one, then the cross product is sent. For example the following constraint could be interpreted as placing a budget of $i$ on the $i^{\mbox{th}}$ item to buy where the cost per item is given by the parameter model.a:

model.A = RangeSet(1,10)
model.a = Param(model.A, within=PostiveReals)
model.ToBuy = Var(model.A)
def bud_rule(model, i):
    return model.a[i]*model.ToBuy[i] <= i
aBudget = Constraint(model.A)

Note	Python and Pyomo are case sensitive so `model.a` is not the same as `model.A`.

9. Rules to Generate Expressions

Both objectives and constraints make use of rules to generate expressions. These are Python functions that return the appropriate expression. These are first-class functions that can access global data as well as data passed in, including the model object.

Operations on model elements results in expressions, which seems natural in expression like the constraints we have seen so far. It is also possible to build up expressions. The following example illustrates this along with a reference to global Pyton data in the form of a Python variable called switch:

switch = 3

model.A = RangeSet(1, 10)
model.c = Param(model.A)
model.d = Param()
model.x = Var(model.A, domain=Boolean)

def pi_rule(model)
    accexpr = summation(model.c, model.x)
    if switch >= 2:
        accexpr = accexpr - model.d
    return accexpr >= 0.5
PieSlice = Constraint(rule=pi_rule)

In this example, the constraint that is generated depends on the value of the Python variable called switch. If the value is 2 or greater, then the constraint is summation(model.c, model.x) - model.d >= 0.5; otherwise, the model.d term is not present.

Caution

Because model elements result in expressions, not values, the following does not work as expected in an abstract model!

model.A = RangeSet(1, 10)
model.c = Param(model.A)
model.d = Param()
model.x = Var(model.A, domain=Boolean)

def pi_rule(model)
    accexpr = summation(model.c, model.x)
    if model.d >= 2:  # NOT in an abstract model!!
        accexpr = accexpr - model.d
    return accexpr >= 0.5
PieSlice = Constraint(rule=pi_rule)

The trouble is that model.d >= 2 results in an expression, not its evaluated value. Instead use if value(model.d) >= 2

10. Non-linear Expressions

Pyomo supports non-linear expressions and can call non-linear solvers such as Ipopt.

11. Data Input

For abstract models, Pyomo supports data input from a data command file or using ModelData object.

11.1. Data Command Files

The syntax of a data command files is based on AMPL’s data commands. Often, however, the data command file refers to files in other formats.

Here are the commands of interest to us that can be used in data command files:

• set declares set data. • param declares a table of parameter data, which can also include the declaration of the set data used to index parameter data. • import defines import from external data sources such as ASCII table files, CSV files, ranges in spreadsheets, and database tables. • include specifies a data command file that is to be pro- cessed immediately. • The data and end commands do not perform any actions, but they provide compatibility with AMPL scripts that define data commands.

The namespace declaration allows data commands to be organized into named groups that can be enabled or disabled during model construction.

Note

Apart from the set and param commands, Pyomo’s data commands do not exactly correspond to AMPL data commands. In particular, the syntax of the AMPL table command allows the user to specify complex mappings from table data values to corresponding model parameters and sets. The corresponding Pyomo import command supports much simpler mappings. Complex mappings are accomplished in Pyomo via scripting.

11.2. ModelData Objects

The ModelData object reads in data sequentially. Any data previously read is overwritten.

Here is an excerpt from an example in [PyomoBook] that illustrates the use of ModelData objects:

from coopr.pyomo import ModelData
from coopr.opt import SolverFactory
from DiseaseEstimation import model

modeldata = ModelData()
modeldata.add(’DiseaseEstimation.dat’)
modeldata.add(’DiseasePop.dat’)
modeldata.read(model)
instance = model.create(modeldata)

12. PySP Overview

This chapter describes PySP: (Pyomo Stochastic Programming), where parameters are allowed to be uncertain.

12.1. Overview of Modeling Components and Processes

The sequence of activities is typically the following:

Create a deterministic model and declare components
Develop base-case data for the deterministic model
Test, verify and validate the deterministic model
Model the stochastic processes
Develop a way to generate scenarios (in the form of a tree if there are more than two stages)
Create the data files need to describe the stochastics
Use PySP to solve stochastic problem

When viewed from the standpoint of file creation, the process is

Create an abstract model for the deterministic problem in a file called ReferenceModel.py
Specify data for this model in a file called ReferenceModel.dat
Specify the stochastics in a file called ScenarioStructure.dat
Specify scenario data

12.2. Birge and Louveaux’s Farmer Problem

Birge and Louveaux [BirgeLouveauxBook] make use of the example of a farmer who has 500 acres that can be planted in wheat, corn or sugar beets, at a per acre cost of $150, $230 and $260, respectively. The farmer needs to have at least 200 tons of wheat and 240 tons of corn to use as feed, but if enough is not grown, those crops can be purchased for $238 and $210, respectively. Corn and wheat grown in excess of the feed requirements can be sold for $170 and $150, respectively. A price of $36 per ton is guaranteed for the first 6000 tons grown by any farmer, but beets in excess of that are sold for $10 per ton. The yield is 2.5, 3, and 20 tons per acre for wheat, corn and sugar beets, respectively.

12.2.1. ReferenceModel.py

So far, this is a deterministic problem because we are assuming that we know all the data. The Pyomo model for this problem shown here is in the file ReferenceModel.py in the sub-directory examples/pysp/farmer/models that is distributed with Coopr.

# Farmer: rent out version has a singleton root node var
# note: this will minimize
#
# Imports
#

from coopr.pyomo import *

#
# Model
#

model = AbstractModel()

#
# Parameters
#

model.CROPS = Set()

model.TOTAL_ACREAGE = Param(within=PositiveReals)

model.PriceQuota = Param(model.CROPS, within=PositiveReals)

model.SubQuotaSellingPrice = Param(model.CROPS, within=PositiveReals)

def super_quota_selling_price_validate (model, value, i):
    return model.SubQuotaSellingPrice[i] >= model.SuperQuotaSellingPrice[i]

model.SuperQuotaSellingPrice = Param(model.CROPS, validate=super_quota_selling_price_validate)

model.CattleFeedRequirement = Param(model.CROPS, within=NonNegativeReals)

model.PurchasePrice = Param(model.CROPS, within=PositiveReals)

model.PlantingCostPerAcre = Param(model.CROPS, within=PositiveReals)

model.Yield = Param(model.CROPS, within=NonNegativeReals)

#
# Variables
#

model.DevotedAcreage = Var(model.CROPS, bounds=(0.0, model.TOTAL_ACREAGE))

model.QuantitySubQuotaSold = Var(model.CROPS, bounds=(0.0, None))
model.QuantitySuperQuotaSold = Var(model.CROPS, bounds=(0.0, None))

model.QuantityPurchased = Var(model.CROPS, bounds=(0.0, None))

model.FirstStageCost = Var()
model.SecondStageCost = Var()

#
# Constraints
#

def ConstrainTotalAcreage_rule(model):
    return summation(model.DevotedAcreage) <= model.TOTAL_ACREAGE

model.ConstrainTotalAcreage = Constraint(rule=ConstrainTotalAcreage_rule)

def EnforceCattleFeedRequirement_rule(model, i):
    return model.CattleFeedRequirement[i] <= (model.Yield[i] * model.DevotedAcreage[i]) + model.QuantityPurchased[i] - model.QuantitySubQuotaSold[i] - model.QuantitySuperQuotaSold[i]

model.EnforceCattleFeedRequirement = Constraint(model.CROPS)

def LimitAmountSold_rule(model, i):
    return model.QuantitySubQuotaSold[i] + model.QuantitySuperQuotaSold[i] - (model.Yield[i] * model.DevotedAcreage[i]) <= 0.0

model.LimitAmountSold = Constraint(model.CROPS)

def EnforceQuotas_rule(model, i):
    return (0.0, model.QuantitySubQuotaSold[i], model.PriceQuota[i])

model.EnforceQuotas = Constraint(model.CROPS)

#
# Stage-specific cost computations
#

def ComputeFirstStageCost_rule(model):
    return model.FirstStageCost - summation(model.PlantingCostPerAcre, model.DevotedAcreage) == 0.0

model.ComputeFirstStageCost = Constraint()

def ComputeSecondStageCost_rule(model):
    expr = summation(model.PurchasePrice, model.QuantityPurchased)
    expr -= summation(model.SubQuotaSellingPrice, model.QuantitySubQuotaSold)
    expr -= summation(model.SuperQuotaSellingPrice, model.QuantitySuperQuotaSold)
    return (model.SecondStageCost - expr) == 0.0

model.ComputeSecondStageCost = Constraint()

#
# Objective
#

def Total_Cost_Objective_rule(model):
    return model.FirstStageCost + model.SecondStageCost

model.Total_Cost_Objective = Objective(sense=minimize)

12.2.2. ReferenceModel.dat

The data introduced here are in the file ReferenceModel.dat in the sub-directory examples/pysp/farmer/scenariodata that is distributed with Coopr.

set CROPS := WHEAT CORN SUGAR_BEETS ;

param TOTAL_ACREAGE := 500 ;

# no quotas on wheat or corn
param PriceQuota :=
            WHEAT 100000 CORN 100000 SUGAR_BEETS 6000 ;

param SubQuotaSellingPrice :=
            WHEAT 170 CORN 150 SUGAR_BEETS 36 ;

param SuperQuotaSellingPrice :=
            WHEAT 0 CORN 0 SUGAR_BEETS 10 ;

param CattleFeedRequirement :=
            WHEAT 200 CORN 240 SUGAR_BEETS 0 ;

# can't purchase beets (no need, as cattle don't eat them)
param PurchasePrice :=
            WHEAT 238 CORN 210 SUGAR_BEETS 100000 ;

param PlantingCostPerAcre :=
            WHEAT 150 CORN 230 SUGAR_BEETS 260 ;

param Yield := WHEAT 3.0 CORN 3.6 SUGAR_BEETS 24 ;

Any of these data could be modeled as uncertain, but we will consider only the possibility that the yield per acre could be higher or lower than expected. Assume that there is a probability of 1/3 that the yields will be the average values that were given (i.e., wheat 2.5; corn 3; and beets 20). Assume that there is a 1/3 probability that they will be lower (2, 2.4, 16) and 1/3 probability they will be higher (3, 3.6, 24). We refer to each full set of data as a scenario and collectively we call them a scenario tree. In this case the scenario tree is very simple: there is a root node and three leaf nodes: one corresponding to each scenario. The acreage-to-plant decisions are root node decisions because they must be made without knowing what the yield will be. The other variables are so-called second stage decisions, because they will depend on which scenario is realized.

12.2.3. ScenarioStructure.dat

PySP requires that users describe the scenario tree using specific constructs in a file named ScenarioStructure.dat; for the farmer problem, this file can be found in the coopr sub-directory examples/pysp/farmer/scenariodata that is distributed with Coopr.

# IMPORTANT - THE STAGES ARE ASSUMED TO BE IN TIME-ORDER.

set Stages := FirstStage SecondStage ;

set Nodes := RootNode
             BelowAverageNode
             AverageNode
             AboveAverageNode ;

param NodeStage := RootNode         FirstStage
                   BelowAverageNode SecondStage
                   AverageNode      SecondStage
                   AboveAverageNode SecondStage ;

set Children[RootNode] := BelowAverageNode
                          AverageNode
                          AboveAverageNode ;

param ConditionalProbability := RootNode          1.0
                                BelowAverageNode  0.33333333
                                AverageNode       0.33333334
                                AboveAverageNode  0.33333333 ;

set Scenarios := BelowAverageScenario
                 AverageScenario
                 AboveAverageScenario ;

param ScenarioLeafNode :=
                    BelowAverageScenario BelowAverageNode
                    AverageScenario      AverageNode
                    AboveAverageScenario AboveAverageNode ;

set StageVariables[FirstStage] :=  DevotedAcreage[*] ;
set StageVariables[SecondStage] := QuantitySubQuotaSold[*]
                                   QuantitySuperQuotaSold[*]
                                   QuantityPurchased[*] ;

param StageCostVariable := FirstStage  FirstStageCost
                           SecondStage SecondStageCost ;

This data file is verbose and somewhat redundant, but in most applications it is generated by software rather than by a person, so this is not an issue. Generally, the left-most part of each expression (e.g. “set Stages :=”) is required and uses reserved words (e.g., Stages) and the other names are supplied by the user (e.g., “FirstStage” could be any name). Every assignment is terminated with a semi-colon. We will now consider the assignments in this file one at a time.

The first assignments provides names for the stages and the words "set Stages" are required, as are the := symbols. Any names can be used. In this example, we used "FirstStage" and "SecondStage" but we could have used "EtapPrimero" and "ZweiteEtage" if we had wanted to. Whatever names are given here will continue to be used to refer to the stages in the rest of the file. The order of the names is important. A simple way to think of it is that generally, the names must be in time order (technically, they need to be in order of information discovery, but that is usually time-order). Stages refers to decision stages, which may, or may not, correspond directly with time stages. In the farmer example, decisions about how much to plant are made in the first stage and "decisions" (which are pretty obvious, but which are decision variables nonetheless) about how much to sell at each price and how much needs to be bought are second stage decisions because they are made after the yield is known.

set Stages := FirstStage SecondStage ;

Node names are constructed next. The words "set Nodes" are required, but any names may be assigned to the nodes. In two stage stochastic problems there is a root node, which we chose to name "RootNode" and then there is a node for each scenario.

set Nodes := RootNode
             BelowAverageNode
             AverageNode
             AboveAverageNode ;

Nodes are associated with time stages with an assignment beginning with the required words "param Nodestage." The assignments must make use of previously defined node and stage names. Every node must be assigned a stage.

param NodeStage := RootNode         FirstStage
                   BelowAverageNode SecondStage
                   AverageNode      SecondStage
                   AboveAverageNode SecondStage ;

The structure of the scenario tree is defined using assignment of children to each node that has them. Since this is a two stage problem, only the root node has children. The words "param Children" are required for every node that has children and the name of the node is in square brackets before the colon-equals assignment symbols. A list of children is assigned.

set Children[RootNode] := BelowAverageNode
                          AverageNode
                          AboveAverageNode ;

The probability for each node, conditional on observing the parent node is given in an assignment that begins with the required words "param ConditionalProbability." The root node always has a conditional probability of 1, but it must always be given anyway. In this example, the second stage nodes are equally likely.

param ConditionalProbability := RootNode          1.0
                                BelowAverageNode  0.33333333
                                AverageNode       0.33333334
                                AboveAverageNode  0.33333333 ;

Scenario names are given in an assignment that begins with the required words "set Scenarios" and provides a list of the names of the scenarios. Any names may be given. In many applications they are given unimaginative names generated by software such as "Scen1" and the like. In this example, there are three scenarios and the names reflect the relative values of the yields.

set Scenarios := BelowAverageScenario
                 AverageScenario
                 AboveAverageScenario ;

Leaf nodes, which are nodes with no children, are associated with scenarios. This assignment must be one-to-one and it is initiated with the words "param ScenarioLeafNode" followed by the colon-equals assignment characters.

param ScenarioLeafNode :=
                    BelowAverageScenario BelowAverageNode
                    AverageScenario      AverageNode
                    AboveAverageScenario AboveAverageNode ;

Variables are associated with stages using an assignment that begins with the required words "set StageVariables" and the name of a stage in square brackets followed by the colon-equals assignment characters. Variable names that have been defined in the file ReferenceModel.py can be assigned to stages. Any variables that are not assigned are assumed to be in the last stage. Variable indexes can be given explicitly and/or wildcards can be used. Note that the variable names appear without the prefix "model." In the farmer example, DevotedAcreage is the only first stage variable.

set StageVariables[FirstStage] :=  DevotedAcreage[*] ;
set StageVariables[SecondStage] := QuantitySubQuotaSold[*]
                                   QuantitySuperQuotaSold[*]
                                   QuantityPurchased[*] ;

For reporting purposes, it is useful to define auxiliary variables in ReferenceModel.py that will be assigned the cost associated with each stage. This variables do not impact algorithms, but the values are output by some software during execution as well as upon completion. The names of the variables are assigned to stages using the "param StageCostVariable" assignment. The stages are previously defined in ScenarioStructure.dat and the variables are previously defined in ReferenceModel.py. Note that the variable names appear without the prefix "model."

param StageCostVariable := FirstStage  FirstStageCost
                           SecondStage SecondStageCost ;

12.2.4. Scenario data specification

So far, we have given a model in the file named ReferenceModel.py, a set of deterministic data in the file named ReferenceModel.py, and a description of the stochastics in the file named ScenarioStructure.dat. All that remains is to give the data for each scenario. There are two ways to do that in PySP: scenario-based and node-based. The default is scenario-based so we will describe that first.

For scenario-based data, the full data for each scenario is given in a .dat file with the root name that is the name of the scenario. So, for example, the file named AverageScenario.dat must contain all the data for the model for the scenario named "AvererageScenario." It turns out that this file can be created by simply copying the file ReferenceModel.dat as shown above because it contains a full set of data for the "AverageScenario" scenario. The files BelowAverageScenario.dat and AboveAverageScenario.dat will differ from this file and from each other only in their last line, where the yield is specified. These three files are distributed with Coopr and are in the coopr sub-directory examples/pysp/farmer/scenariodata along with ScenarioStructure.dat and ReferenceModel.dat.

Scenario-based data wastes resources by specifying the same thing over and over again. In many cases, that does not matter and it is convenient to have full scenario data files available (for one thing, the scenarios can easily be run independently using the pyomo command). However, in many other settings, it is better to use a node-based specification where the data that is unique to each node is specified in a .dat file with a root name that matches the node name. In the farmer example, the file RootNode.dat will be the same as ReferenceModel.dat except that it will lack the last line that specifies the yield. The files BelowAverageNode.dat, AverageNode.dat, and AboveAverageNode.dat will contain only one line each to specify the yield. If node-based data is to be used, then the ScenarioStructure.dat file must contain the following line:

param ScenarioBasedData := False ;

An entire set of files for node-based data for the farmer problem are distributed with Coopr in the sub-directory examples/pysp/farmer/nodedata

12.3. Finding Solutions for Stochastic Models

PySP provides a variety of tools for finding solutions to stochastic programs.

12.3.1. runef

The runef command puts together the so-called extensive form version of the model. It creates a large model that has constraints to ensure that variables at a node have the same value. For example, in the farmer problem, all of the DevotedAcres variables must have the same value regardless of which scenario is ultimately realized. The objective can be the expected value of the objective function, or the CVaR, or a weighted combination of the two. Expected value is the default. A full set of options for runef can be obtained using the command:

runef --help

The coopr distribution contains the files need to run the farmer example in the sub-directories to the sub-directory examples/pysp/farmer so if this is the current directory and if CPLEX is installed, the following command will cause formation of the EF and its solution using CPLEX.

runef -m models -i nodedata --solver=cplex --solve

The option -m models has one dash and is short-hand for the option --model-directory=models and note that the full option uses two dashes. The -i is equivalent to --instance-directory= in the same fashion. The default solver is CPLEX, so the solver option is not really needed. With the --solve option, runef would simply write an .lp data file that could be passed to a solver.

12.3.2. runph

The runph command executes an implementation of Progressive Hedging (PH) that is intended to support scripting and extension.

The coopr distribution contains the files need to run the farmer example in the sub-directories to the sub-directory examples/pysp/farmer so if this is the current directory and if CPLEX is installed, the following command will cause PH to execute using the default sub-problem solver, which is CPLEX.

runph -m models -i nodedata

After about 33 iterations, the algorithm will achieve the default level of convergence and terminate. A lot of output is generated and among the output is the following solution information:

Variable=DevotedAcreage
        Index: [CORN]            (Scenarios:  BelowAverageScenario   AverageScenario   AboveAverageScenario   )
                Values:       79.9844      80.0000      79.9768     Max-Min=      0.0232     Avg=     79.9871
        Index: [SUGAR_BEETS]             (Scenarios:  BelowAverageScenario   AverageScenario   AboveAverageScenario   )
                Values:      249.9848     249.9770     250.0000     Max-Min=      0.0230     Avg=    249.9873
        Index: [WHEAT]           (Scenarios:  BelowAverageScenario   AverageScenario   AboveAverageScenario   )
                Values:      170.0308     170.0230     170.0232     Max-Min=      0.0078     Avg=    170.0256
Cost Variable=FirstStageCost
        Tree Node=RootNode               (Scenarios:  BelowAverageScenario   AverageScenario   AboveAverageScenario   )
        Values:   108897.0836  108897.4725  108898.1476     Max-Min=      1.0640     Avg= 108897.5679

For problems with no, or few, integer variables, the default level of convergence leaves root-node variables almost converged. Since the acreage to be planted cannot depend on the scenario that will be realized in the future, the average, which is labeled "Avg" in this output, would be used. A farmer would probably interpret acreages of 79.9871, 249.9873, and 170.0256 to be 80, 250, and 170. In real-world applications, PH is embedded in scripts that produce output in a format desired by a decision maker.

But in real-world applications, the default settings for PH seldom work well enough. In addition to post-processing the output, a number of parameters need to be adjusted and sometimes scripting to extend or augment the algorithm is needed to improve convergence rates. A full set of options can be obtained with the command.

runph --help

Note that there are two dashes before help.

By default, PH uses quadratic objective functions after iteration zero; in some settings it may be desirable to linearize the quadratic terms. This is required to use a solver such as glpk for MIPs because it does not support quadratic MIPs. The directive --linearize-nonbinary-penalty-terms=n causes linearization of the penalty terms using n pieces. For example, to use glpk on the farmer, assuming glpk is installed and the command is given when the current directory is the examples/pysp/farmer, the following command will use default settings for most parameters and four pieces to approximate quadratic terms in sub-problems:

runph -i nodedata -m models --solver=glpk --linearize-nonbinary-penalty-terms=4

Use of the linearize-nonbinary-penalty-terms option requires that all variables not in the final stage have bounds.

13. Scripts

13.1. Iterative Example

14. Coopr Solver Interfaces

This chapter describes how Coopr interfaces with different solvers.

15. Using Black-Box Optimizers with Coopr.Opt

Many optimization software packages contain black-box optimizers, which perform optimization without using detailed knowledge of the structure of an optimization problem. Thus, black-box optimizers require a generic interface for optimization problems that defines key features of problems, like objectives and constraints.

The coopr.opt package contains the coopr.opt.blackbox subpackage, which provides facilities for (a) integrating Coopr solvers with blackbox optimization applications and (b) wrapping Pyomo models for use by external blackbox optimizers. We illustrate these capabilities in this chapter with simple examples that illustrate the use of coopr.opt.blackbox.

15.1. Defining and Optimizing Simple Black-Box Applications

Many black-box optimizers interact with an optimization problem by executing a separate process that computes properties of the optimization problem. This process typically reads an input file that defines the requested properties and writes an output file that contains the computed values. Unfortunately, no standards have emerged for black-box optimizers that interact with problems in this manner. Thus, different file formats are used by different optimizer software packages.

15.1.1. Defining an Optimization Problem

The coopr.opt.blackbox package provides several Python classes for optimization problems that coordinates file I/O for the user and simplifies the definition of simple black-box problems. The RealOptProblem class provides a generic interface for continuous optimization problems (i.e. with real variables). The following example defines a simple continuous optimization problem:

This problem is equivalent to the following problem definition:

$\begin{array}{lll} \min & x_0 - x_1 + (x_2 - 1.5)^2 + (x_3+2)^4 & \\ \mathrm{s.t.} & 0 \leq x_0 & \\ & -1 \leq x_1 \leq 0 & \\ & 0 \leq x_2 \leq 2 & \\ & x_3 \leq -1 & \\ \end{array}$

Note that the problem class does not specify the sense of the optimization problem. These problem classes are not a complete specification of an optimization problem. Rather, an instance of a problem class can compute information about the problem that is used during optimization.

Similarly, the MixedIntOptProblem class provides a generic interface for mixed-integer optimization problems, which may contain real variables, integer variables and binary variables. The following example defines a simple mixed-integer optimization problem:

This problem is equivalent to the following problem definition:

$\begin{array}{lll} \min & \sum_{i=1}^4 (x_i-1)^2 + \sum_{i=1}^3 (y_i+1)^2 + \sum_{i=1}^2 z_i & \\ \mathrm{s.t.} & 0 \leq x_i \leq 2 & \\ & -2 \leq y_i \leq 0 & \\ & z_i \in \{0,1\} & \\ \end{array}$

15.1.2. Optimizating with Coliny Solvers

The Coliny software library supports interfaces to a variety of black-box optimizers <Coliny>. The coliny executable reads an XML specification of the optimization problem and solver, as well as a specification of how the optimizer is applied. Consider the following XML specification:

This XML specification defines a MINLP0 problem, which indicated that this is a mixed-integer problem that supports zero-order derivatives (i.e. no derivatives). This problem has four real variables with lower and upper bounds specified. The problem values are computed with the RealProblem1.py command-line, which defines and uses the RealProblem1 class defined above:

Note that this command is a Python script that includes the shebang character sequence on the first line. On Linux and Unix systems, this line indicates that this is a script that is executed using the python command that is found in the user environment. Thus, this example assumes that the python command has coopr.opt installed. Since multiple versions of Python can be installed on a single computer, the XML Command element may need to be defined with an explicitly Python version. For example, if Python 2.6 is installed in /usr/local with coopr.opt, then the Command element would look like:

Additionally, the duplication of bounds information between RealProblem1.py and RealProblem1.xml is not strictly necessary in this example. The bounds information in RealProblem1.py is used in the validate method to verify that the point being evaluated is consistent with the bounds information. We can generally assume that the Coliny solver will only evaluate feasible points, so a simpler problem definition can be used:

The last two lines of RealProblem1.py create a problem instance and then call the main method to parse the command-line arguments. This script has the following command-line syntax:

The first argument is the name of an XML input file, and the second argument is the name of an XML output file. The optimization problem class manages the parsing of the input and generation of the output file. For example, consider the following input file:

The RealProblem1.py script creates the following output file:

15.1.3. Optimizating with Dakota Solvers

TODO: The only difference is the format used in the script.

15.2. Diving Deeper

The previous section provided an overview of the how the coopr.opt.blackbox package supports the definition of optimization problems that are solved with black-box optimizers. In this section we provide more detail about how the Python problem class can be customized, as well as details about the XML file format used to communicate with Coliny optimizers. The Dakota User Manual <Dakota> provides documentation of the file format of the input and output files used with Dakota optimizers.

Table [Table-OptProblemMethods] summarizes the methods of the OptProblem class that a user is likely to either use or redefine when declaring a subclass. The MixedIntOptProblem class is a convenient base class for the problems solved by most black-box optimizers, and this class provides the definition of the main, create_point and validate methods. However, any of the remaining methods may need to be defined, depending on the problem.

Table 1. Methods in the `OptProblem` class.
Method	Description
_init_	The constructor, which may be redefined to specify problem properties.
main	Method that processes command-line options to create a results file from an input file.
create_point	Create an instances of the class that defines a point in the search domain.
function_value	Returns the value of the objective function.
function_values	Returns a list of objective function values.
gradient	Returns a list that represents the gradient vector at the given point.
hessian	Returns a Hessian matrix.
nonlinear_constraint_values	Returns a list of values for the constraint functions.
jacobian	Returns a Jacobian matrix.
validate	Returns `True` if the given point is feasible, and `False` otherwise.

The following detailed example illustrates the use of all of these methods in a simple application:

The response_types attribute defined in the constructor specifies the type of information that this class can compute. For example, consider the following input XML file:

This input file requests that the class compute all of the response values, and thus the following output is generated:

Note that the values for Jacobian and Hessian matrices are represented in a sparse manner. Currently, these are represented with a list of tuple values, though a sparse matrix representation might be supported in the future.

16. Coopr Reference Documentation

16.1. Components

TODO

16.2. Utility Functions