Examples

If you installed directly from GitHub, the top-level examples directory contains sub-directories with examples. Each example can be run via generic_cylinders.py (see generic_cylinders.py).

If you did not get the code from GitHub (e.g., if you installed with pip), you will need to get the examples from: https://github.com/Pyomo/mpi-sppy/tree/master/examples

Note

Some examples also have historical *_cylinders.py driver scripts. These are retained for reference but generic_cylinders.py is the recommended way to run all examples.

Tutorial: Farmer Example

In this section, we step through a simple example—the farmer example of [birge2011] (Section 1.1). This model can be expressed as the following linear program (LP):

\begin{array}{rl} \min & (150x_1 + 230x_2 + 260x_3) + (238y_1+210y_2) - (170w_1 + 150w_2 + 36w_3 + 10 w_4) \\ \mathrm{s.t.} & x_1 + x_2 + x_3 \leq 500 \\ & (2.5)x_1 + y_1 - w_1 \geq 200 \\ & (3)x_2 + y_2 - w_2 \geq 240 \\ & (20)x_3 - w_3 - w_4 \geq 0 \\ & w_3 \leq 6000 \\ & x,y,w\geq0 \end{array}

The decision variables are as follows:

\(x_i\) = number of acres to devote to crop i (1=wheat, 2=corn, 3=sugar beets)
\(y_i\) = tons of crop i to purchase from a wholesaler (i=1,2 but not 3)
\(w_i\) = tons of crop i sold (i=1,2)
\(w_i\) = tons of beets sold at favorable (i=3) or unfavorable (i=4) price

The coefficients of the \(x_i\) variables in the second, third and fourth constraints are the number of tons per acre that each crop will yield (2.5 for wheat, 3 for corn, and 20 for sugar beets).

The following code in examples/farmer/archive/farmer.py (with similar code in examples/farmer/farmer.py) creates an instance of the farmer’s model:

import pyomo.environ as pyo

def build_model(yields):
    model = pyo.ConcreteModel()

    # Variables
    model.X = pyo.Var(["WHEAT", "CORN", "BEETS"], within=pyo.NonNegativeReals)
    model.Y = pyo.Var(["WHEAT", "CORN"], within=pyo.NonNegativeReals)
    model.W = pyo.Var(
        ["WHEAT", "CORN", "BEETS_FAVORABLE", "BEETS_UNFAVORABLE"],
        within=pyo.NonNegativeReals,
    )

    # Objective function
    model.PLANTING_COST = 150 * model.X["WHEAT"] + 230 * model.X["CORN"] + 260 * model.X["BEETS"]
    model.PURCHASE_COST = 238 * model.Y["WHEAT"] + 210 * model.Y["CORN"]
    model.SALES_REVENUE = (
        170 * model.W["WHEAT"] + 150 * model.W["CORN"]
        + 36 * model.W["BEETS_FAVORABLE"] + 10 * model.W["BEETS_UNFAVORABLE"]
    )
    model.OBJ = pyo.Objective(
        expr=model.PLANTING_COST + model.PURCHASE_COST - model.SALES_REVENUE,
        sense=pyo.minimize
    )

    # Constraints
    model.CONSTR= pyo.ConstraintList()

    model.CONSTR.add(pyo.summation(model.X) <= 500)
    model.CONSTR.add(
        yields[0] * model.X["WHEAT"] + model.Y["WHEAT"] - model.W["WHEAT"] >= 200
    )
    model.CONSTR.add(
        yields[1] * model.X["CORN"] + model.Y["CORN"] - model.W["CORN"] >= 240
    )
    model.CONSTR.add(
        yields[2] * model.X["BEETS"] - model.W["BEETS_FAVORABLE"] - model.W["BEETS_UNFAVORABLE"] >= 0
    )
    model.W["BEETS_FAVORABLE"].setub(6000)

    return model

Note that the build_model function takes a list of values, containing the yields for each crop. We can solve the model:

yields = [2.5, 3, 20]
model = build_model(yields)
solver = pyo.SolverFactory("cplex_direct")
solver.solve(model)

# Display the objective value to one decimal place
print(f"{pyo.value(model.OBJ):.1f}")

The optimal objective value is:

  -118600.0

In practice, the farmer does not know the number of tons that each crop will yield per acre planted–the yield depends on the weather, the quality of the seeds, and other stochastic factors. Consequently, we replace the deterministic model above with the stochastic LP:

\begin{array}{rl} \min & (150x_1 + 230x_2 + 260x_3) \\ & \quad+\sum_{\omega\in\Omega}Pr[\omega]\big[(238y_1^\omega+210y_2^\omega) - (170w_1^\omega + 150w_2^\omega + 36w_3^\omega + 10 w_4^\omega)\big] \\ \mathrm{s.t.} & x_1 + x_2 + x_3 \leq 500 \\ & \xi^\omega_1 x_1 + y^\omega_1 - w^\omega_1 \geq 200\;\forall\;\omega\in\Omega\\ & \xi^\omega_2 x_2 + y^\omega_2 - w^\omega_2 \geq 240\;\forall\;\omega\in\Omega\\ & \xi^\omega_3 x_3 - w^\omega_3 - w^\omega_4 \geq 0\;\forall\;\omega\in\Omega\\ & w^\omega_3 \leq 6000 \\ & x,y^\omega,w^\omega\geq0\;\forall\;\omega\in\Omega \end{array}

The variables \(y_i\) and \(w_i\) have been replaced with copies \(y_i^\omega\) and \(w_i^\omega\), corresponding to the values of each variable chosen under scenario \(\omega\in\Omega\), where \(\Omega\) is a finite set of scenarios. The parameter \(\xi^\omega_i\) is the number of tons of crop \(i\) yielded per acre under scenario \(\omega\).

We assume that there are three scenarios: “good”, “bad”, and “average”. We assume that each scenario is equally likely to occur. The yield values (\(\xi^\omega_i\)) are given here:

Table 1 Crop yields under each scenario (tons/acre)
	Wheat	Corn	Sugar Beets
Good	3	3.6	24
Average	2.5	3	20
Bad	2	2.4	16

In order to transform the code for the deterministic model above into a stochastic model which can be manipulated by MPI-SPPy, we need only incorporate a few extra elements (see scenario_creator function for full details). The scenario_creator function is told the name of the scenario to build, and builds a Pyomo model for that scenario appropriately:

import mpisppy.utils.sputils as sputils

def scenario_creator(scenario_name):
    if scenario_name == "good":
        yields = [3, 3.6, 24]
    elif scenario_name == "average":
        yields = [2.5, 3, 20]
    elif scenario_name == "bad":
        yields = [2, 2.4, 16]
    else:
        raise ValueError("Unrecognized scenario name")

    model = build_model(yields)
    sputils.attach_root_node(model, model.PLANTING_COST, [model.X])
    model._mpisppy_probability = 1.0 / 3
    return model

The scenario_creator accomplishes two important tasks

It calls the attach_root_node function. We tell this function which part of the objective function (model.PLANTING_COST) and which set of variables (model.X) belong to the first stage. In this case, the problem is only two stages, so we need only specify the root node and the first-stage information–MPI-SPPy assumes the remainder of the model belongs to the second stage.
It attaches an attribute called _mpisppy_probability to the model object. This is the probability that the specified scenario occurs. If this probability is not specified, MPI-SPPy will assume that all scenarios are equally likely.

Now that we have specified a scenario creator, we can use MPI-SPPy to solve the farmer’s stochastic program.

Solving the Extensive Form

The simplest approach is to solve the extensive form of the model directly. Assuming you are in the directory examples/farmer the following unix command will work.

python ../../mpisppy/generic_cylinders.py --module-name farmer --num-scens 3 --EF --EF-solver-name gurobi

We can extract the optimal solution itself using the --solution-base-name option:

python ../../mpisppy/generic_cylinders.py --module-name farmer --num-scens 3 --EF --EF-solver-name gurobi --solution-base-name farmersol

This command writes solution data for nonanticipative variables to two files with the base name farmersol and full scenario solutions to a directory named farmersol_soldir.

Note

Most command line options relevant to the EF start with –EF. Most other command line options will be silently ignored if --EF is specified (one exception is --solution-base-name).

Solving Using Progressive Hedging (PH)

Here is a simple command that uses PH as the hub algorithm and computes lower bounds using a Lagrangian spoke (--lagrangian) with upper bounds computed by randomly trying scenario solutions to fix the nonanticipative variables (--xhatshuffle).

mpiexec -np 3 python -m mpi4py ../../mpisppy/generic_cylinders.py --module-name farmer --num-scens 3 --solver-name gurobi_persistent --max-iterations 10 --max-solver-threads 4 --default-rho 1 --lagrangian --xhatshuffle --rel-gap 0.01

Solving Using Benders’ Decomposition

Finally, we can solve our example using Benders’ decomposition, known as the L-shaped method in stochastic programming. The setup code is similar to the previous methods:

from mpisppy.opt.lshaped import LShapedMethod

all_scenario_names = ["good", "average", "bad"]
bounds = {name: -432000 for name in all_scenario_names}
options = {
    "root_solver": "cplex_persistent",
    "sp_solver": "cplex_persistent",
    "sp_solver_options" : {"threads" : 1},
    "valid_eta_lb": bounds,
    "max_iter": 10,
}

ls = LShapedMethod(options, all_scenario_names, scenario_creator)
result = ls.lshaped_algorithm()

variables = ls.gather_var_values_to_rank0()
for ((scen_name, var_name), var_value) in variables.items():
    print(scen_name, var_name, var_value)

[    0.00] Start SPBase.__init__
Current Iteration: 1 Time Elapsed:    0.00 Current Objective: -Inf
Current Iteration: 2 Time Elapsed:    0.01 Time Spent on Last Master: 0.00 Time Spent Generating Last Cut Set:    0.01 Current Objective: -1296000.00
Current Iteration: 3 Time Elapsed:    0.02 Time Spent on Last Master: 0.00 Time Spent Generating Last Cut Set:    0.01 Current Objective: -160000.00
Current Iteration: 4 Time Elapsed:    0.02 Time Spent on Last Master: 0.00 Time Spent Generating Last Cut Set:    0.00 Current Objective: -113750.00
Converged in 4 iterations.
Total Time Elapsed:    0.03 Time Spent on Last Master:    0.00 Time spent verifying second stage:    0.00 Final Objective: -108390.00
good X[BEETS] 250.0
good X[CORN] 80.0
good X[WHEAT] 170.0
average X[BEETS] 250.0
average X[CORN] 80.0
average X[WHEAT] 170.0
bad X[BEETS] 250.0
bad X[CORN] 80.0
bad X[WHEAT] 170.0

We see that, for this toy example, the L-shaped method has converged to the optimal solution within just 10 iterations.

aircond

This is fairly complicated example because it is multi-stage and the model itself offers a lot of flexibility. The aircond example is unusual in that the model file, aircond.py, lives in mpisppy.tests.examples directory. Scripts and bash files that use it live in examples.aircond. Because aircond.py provides scenario_creator, inparser_adder, kw_creator, scenario_names_creator and sample_tree_scen_creator, runs can go through mpisppy/generic_cylinders.py by pointing --module-name at mpisppy/tests/examples/aircond.

A good place to start reading is the aircond_cylinders.py file, which is retained as a worked example of the hand-rolled cylinders pattern (and, together with bundle_pickler.py, as the canonical proper-bundles demo). It starts with some functions that support the main program, and the main program makes use of the Config object called cfg that creates a parser and gets arguments. The configuration data obtained by the parser are passed directly to the vanilla hub and spoke creator which knows how to use the arguments from a Config object. The arguments unique to aircond are processed by the create_kwargs function in the reference model file.

A simple example that uses a few of the options is shown in aircond_zhat.bash, which also calls the xhat4xhat program to estimate confidence intervals for the solution obtained.

aircondMulti

This is a multi-product version of aircond. The model module lives in mpisppy/tests/examples/aircondMulti.py; runs go through mpisppy/generic_cylinders.py with --module-name ../../mpisppy/tests/examples/aircondMulti. See examples/aircondMulti/README for a pointer and examples/run_all.py for the regression-test command line.

Note: sample_tree_scen_creator is not implemented for aircondMulti, so MMW / sequential-sampling confidence intervals are not available for this example.

hydro

Hydro is a three stage example that was originally coded in PySP and we make extensive use of the PySP files. Unlike farmer and aircond where the scenario data are created from distributions, for this problem the scenario data are provided in files.

The model module hydro.py provides the scenario_creator and the other hooks expected by mpisppy/generic_cylinders.py, so most runs are launched through the generic driver (see examples/run_all.py or examples/generic_tester.py for exact invocations). The file hydro_cylinders.py is retained alongside it as an illustration of the hand-rolled cylinder pattern for multistage problems; it reads hydro.py and builds the hub/spoke dictionaries directly.

An older example using PySPModel is preserved under examples/hydro/archive/hydro_cylinders_pysp.py for reference; new code should skip PySPModel and use generic_cylinders.py.

netdes

This is a very challenging network design problem, which has many instances each defined by a data file. netdes.py provides the model-module hooks so most runs go through mpisppy/generic_cylinders.py. For this problem, cross scenario cuts are helpful; the use of that spoke (together with --slammax) is illustrated in netdes_cylinders.py, which is retained as a worked example of the hand-rolled cylinders pattern.

sslp

This is a classic problem from Ntaimo and Sen with data in PySP format. sslp.py provides the model-module hooks so most runs go through mpisppy/generic_cylinders.py. sslp_cylinders.py is retained as a worked example of the hand-rolled cylinders driver; it is similar to the hydro example except sslp is simpler because it is just two stages.

UC

This example uses the egret package for the underlying unit commitment model and reads PySP format data using the pyomo dataportal. Data files for a variety of numbers of scenarios are provided.

sizes

The sizes example (Jorjani et al, IJPR, 1999) is a two-stage problem with general integers in each stage. sizes.py provides the model-module hooks so most runs go through mpisppy/generic_cylinders.py. The file sizes_cylinders.py is retained as a worked example of the hand-rolled cylinders driver pattern, and sizes_demo.py provides an example of serial execution (no cylinders) used by the extensions documentation.