Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Stochastic Programming

This notebook was prepared by Jialu Wang and revised by Maddie Watson at the University of Notre Dame.

# Imports
import sys

if "google.colab" in sys.modules:
    !wget "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    import helper

    helper.easy_install()
else:
    sys.path.insert(0, "../")
    import helper
helper.set_plotting_style()

Farmers Example

Here is the handout for lecture

Consider a European farmer who specializes in raising wheat, corn, and sugar beets on his 500 acres of land. During the winter, he wants to decide how much land to devote to each crop. (We refer to the farmer as “he” for convenience and not to imply anything about the gender of European farmers)

The farmer knows that at least 200 tons (T) of wheat and 240 T of corn are needed for cattle feed. These amounts can be raised on the farm or bought from a wholesaler. Any production in excess of the feeding requirement would be sold. Over the last decade, mean selling prices have been $170 and $150 per ton of wheat and corn, respectively. The purchase prices are 40% more than this due to the wholesaler’s margin and transportation costs.

Another profitable crop is sugar beet, which he expects to sell at $36/T; however, the European Commission imposes a quota on sugar beet production. Any amount in excess of the quota can be sold only at $10/T. The farmer’s quota for next year is 6000 T.

Based on past experience, the farmer knows that the mean yield on his land is roughly 2.5 T, 3 T, and 20 T per acre for wheat, corn, and sugar beets, respectively. Table 1 summarizes these data and the planting costs for these crops.

WheatCornSugar Beets
Yield (T/acre)2.5320
Planting cost ($/acre)150230260
Selling price ($/T)17015036 under 6000 T, 10 above 6000 T
Purchase price ($/T)238210–
Minimum requirement (T)200240–

Total available land: 500 acres

To help the farmer make up his mind, we can set up the following model. Let

  • x1x_1 = acres of land devoted to wheat,

  • x2x_2 = acres of land devoted to corn,

  • x3x_3 = acres of land devoted to sugar beets,

  • w1w_1 = tons of wheat sold,

  • y1y_1 = tons of wheat purchased,

  • w2w_2 = tons of corn sold,

  • y2y_2 = tons of corn purchased,

  • w3w_3 = tons of sugar beets sold at the favorable price,

  • w4w_4 = tons of sugar beets sold at the unfavorable price.

Farmer's problem: acres allocated to each crop, with required animal feed and quantities sold to or purchased from the market.

Problem Formulation In Words

Minimize total cost, subject to the following constraints:

  • Plant up to 500 acres.

  • Need at least 200 tons of wheat (for animals).

  • Need at least 240 tons of corn (for animals).

  • Sugar beet sales must be less than or equal to the yield from the farm.

  • All variables are positive.

  • Up to 6000 tons of sugar beets can be sold at a favorable price.

Perfect information

With the ‘perfect’ information shown in the table above, the optimization problem is formed as:

min⁡150x1+230x2+260x3⏟planting costs+238y1−170w1⏟wheat purchases less sales+210y2−150w2⏟corn purchases less sales −36w3−10w4⏟sugar beet saless.t.x1+x2+x3≤500,(plant up to 500 acres)2.5x1+y1−w1≥200,(satisfy wheat demand)3x2+y2−w2≥240,(satisfy corn demand)w3+w4≤20x3,(beet sales cannot exceed production)w3≤6000,(cap beet sales at favorable price)x1,x2,x3,y1,y2,w1,w2,w3,w4≥0.\begin{align*} \min \quad & \underbrace{150x_1 + 230x_2 + 260x_3}_{\text{planting costs}} \underbrace{+ 238y_1 - 170w_1}_{\text{wheat purchases less sales}} \\ & \underbrace{+ 210y_2 - 150w_2}_\text{corn purchases less sales} ~ \underbrace{- 36w_3 - 10w_4}_{\text{sugar beet sales}} \\ \text{s.t.} \quad & x_1 + x_2 + x_3 \leq 500, \quad \text{(plant up to 500 acres)} \\ & 2.5x_1 + y_1 - w_1 \geq 200, \quad \text{(satisfy wheat demand)}\\ & 3x_2 + y_2 - w_2 \geq 240, \quad \text{(satisfy corn demand)}\\ & w_3 + w_4 \leq 20x_3, \quad \text{(beet sales cannot exceed production)}\\ & w_3 \leq 6000, \quad \text{(cap beet sales at favorable price)}\\ & x_1, x_2, x_3, y_1, y_2, w_1, w_2, w_3, w_4 \geq 0. \end{align*}
(1)

Create a function to build a pyomo model for the farmer’s problem with crop yields as an input

Units. Land is in acres, crops in metric tons, money in dollars, and yields in tons per acre. The literal coefficients carry their units inline --- 150 * USD_PER_ACRE rather than a bare 150 --- and each builder ends with assert_units_consistent, which raises if any of the eleven constraints or the objective does not balance. A units declaration nobody checks only looks verified.

import pyomo.environ as pyo
from pyomo.environ import (
    Block,
    ConcreteModel,
    Constraint,
    ConstraintList,
    Var,
    NonNegativeReals,
    Objective,
    minimize,
    summation,
    SolverFactory,
    value,
)
from pyomo.environ import units as u
from pyomo.util.check_units import assert_units_consistent

# Pyomo's unit library has acres and metric tons but no money, so declare one.
u.load_definitions_from_strings(["USD = [currency]"])

# Shorthands for the three compound units this problem uses
USD_PER_ACRE = u.USD / u.acre
USD_PER_TON = u.USD / u.metric_ton
TON_PER_ACRE = u.metric_ton / u.acre

### Create a function to build a pyomo model for the farmer's problem with crop yields as an input


def build_model(yields):
    """
    Code adapted from https://mpi-sppy.readthedocs.io/en/latest/examples.html#examples

    Arguments:
        yields: Yield information as a list, following the rank [wheat, corn, beets]

    Return:
        model: farmer problem model
    """
    model = ConcreteModel()

    # Define sets
    all_crops = ["WHEAT", "CORN", "BEETS"]
    purchase_crops = ["WHEAT", "CORN"]
    sell_crops = ["WHEAT", "CORN", "BEETS_FAVORABLE", "BEETS_UNFAVORABLE"]

    # Crops field allocation [acre]
    model.X = Var(all_crops, domain=NonNegativeReals, units=u.acre)
    # How many tons of crops to purchase [t]
    model.Y = Var(purchase_crops, domain=NonNegativeReals, units=u.metric_ton)
    # How many tons of crops to sell [t]
    model.W = Var(sell_crops, domain=NonNegativeReals, units=u.metric_ton)

    # Objective function [USD]
    model.PLANTING_COST = (
        150 * USD_PER_ACRE * model.X["WHEAT"]
        + 230 * USD_PER_ACRE * model.X["CORN"]
        + 260 * USD_PER_ACRE * model.X["BEETS"]
    )
    model.PURCHASE_COST = (
        238 * USD_PER_TON * model.Y["WHEAT"] + 210 * USD_PER_TON * model.Y["CORN"]
    )
    model.SALES_REVENUE = (
        170 * USD_PER_TON * model.W["WHEAT"]
        + 150 * USD_PER_TON * model.W["CORN"]
        + 36 * USD_PER_TON * model.W["BEETS_FAVORABLE"]
        + 10 * USD_PER_TON * model.W["BEETS_UNFAVORABLE"]
    )
    # Maximize the Obj is to minimize the negative of the Obj
    model.OBJ = Objective(
        expr=model.PLANTING_COST + model.PURCHASE_COST - model.SALES_REVENUE,
        sense=minimize,
    )

    # Constraints
    model.CONSTR = ConstraintList()

    # Plant at most 500 acres [acre]
    model.CONSTR.add(summation(model.X) <= 500 * u.acre)
    # Animal feed requirements and beet sales limits [t]
    model.CONSTR.add(
        yields[0] * TON_PER_ACRE * model.X["WHEAT"]
        + model.Y["WHEAT"]
        - model.W["WHEAT"]
        >= 200 * u.metric_ton
    )
    model.CONSTR.add(
        yields[1] * TON_PER_ACRE * model.X["CORN"] + model.Y["CORN"] - model.W["CORN"]
        >= 240 * u.metric_ton
    )
    model.CONSTR.add(
        yields[2] * TON_PER_ACRE * model.X["BEETS"]
        - model.W["BEETS_FAVORABLE"]
        - model.W["BEETS_UNFAVORABLE"]
        >= 0 * u.metric_ton
    )
    model.W["BEETS_FAVORABLE"].setub(6000)

    # Raises UnitsError if any constraint or the objective is inconsistent
    assert_units_consistent(model)

    return model
# Solve the Optimimization Problem with Perfect yields
yields_perfect = [2.5, 3, 20]
model = build_model(yields_perfect)
solver = SolverFactory("ipopt")
results = solver.solve(model)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

# Define a function for printing the optimal solution


def print_opt_sol(model, yields):
    """
    Arguments:
        model: solved farmer problem model
        yields: the yields [wheat, corn, beets] this model was built with, in
            T/acre. These are needed to report the production, and they differ
            from scenario to scenario.

    Return:
        Prints the optimal solution
    """
    print("===Optimal solutions based on perfect information===")

    print("Culture.         | ", "Wheat |", "Corn  |", "Sugar Beets |")
    print(
        "Surface (acres)  | ",
        f'{value(model.X["WHEAT"]):.1f}',
        "|",
        f'{value(model.X["CORN"]):.1f}',
        " |",
        f'{value(model.X["BEETS"]):.1f}',
        " |",
    )
    print(
        "Yield (T)        | ",
        f'{value(model.X["WHEAT"])*yields[0]:.1f}',
        "|",
        f'{value(model.X["CORN"])*yields[1]:.1f}',
        "|",
        f'{value(model.X["BEETS"])*yields[2]:.1f}',
        "|",
    )
    print(
        "Sales (T)        | ",
        f'{value(model.W["WHEAT"]):.1f}',
        "|",
        f'{value(model.W["CORN"]):.1f}',
        "  |",
        f'{value(model.W["BEETS_FAVORABLE"]) + value(model.W["BEETS_UNFAVORABLE"]):.1f}',
        "|",
    )
    print(
        "Purchases (T)    | ",
        f'{value(model.Y["WHEAT"]):.1f}',
        "  |",
        f'{value(model.Y["CORN"]):.1f}',
        "  |",
        "-",
        "     |",
    )

    profit = -value(model.OBJ)
    print("Overall profit: $", f"{profit:.1f}")

    return profit


profit_perfect = print_opt_sol(model, yields_perfect)
===Optimal solutions based on perfect information===
Culture.         |  Wheat | Corn  | Sugar Beets |
Surface (acres)  |  120.0 | 80.0  | 300.0  |
Yield (T)        |  300.0 | 240.0 | 6000.0 |
Sales (T)        |  100.0 | -0.0   | 6000.0 |
Purchases (T)    |  -0.0   | -0.0   | -      |
Overall profit: $ 118600.0

The optimal solution based on perfect information is:

ex1.2

This solution is easy to understand:

  • The farmer devotes enough land to sugar beets to reach the quota of 6000 T

  • Devote enough land to wheat and corn production to meet the feeding requirement

  • Plant wheat in the rest of the land

However, there are often some ‘real world’ constraints that break the perfect information heuristic:

  • Market prices change

  • Yield is uncertain

  • Planting cost materials, water, labor...

  • Crop rotation

A representation of the uncertainty would be to assume that years are good, fair, or bad for all crops, resulting in above average, average, or below average yields for all crops. Three scenarios are defined as:

  • Above average yield (+20%)

  • Average yield (base case)

  • Below average yield (-20%)

### Run Above average case
yields_above = [2.5 * 1.2, 3 * 1.2, 20 * 1.2]
model = build_model(yields_above)
solver = SolverFactory("ipopt")
results = solver.solve(model)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

profit_above = print_opt_sol(model, yields_above)
===Optimal solutions based on perfect information===
Culture.         |  Wheat | Corn  | Sugar Beets |
Surface (acres)  |  183.3 | 66.7  | 250.0  |
Yield (T)        |  550.0 | 240.0 | 6000.0 |
Sales (T)        |  350.0 | -0.0   | 6000.0 |
Purchases (T)    |  -0.0   | -0.0   | -      |
Overall profit: $ 167666.7

### Run Below average case
yields_below = [2.5 * 0.8, 3 * 0.8, 20 * 0.8]
model = build_model(yields_below)
solver = SolverFactory("ipopt")
results = solver.solve(model)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

profit_below = print_opt_sol(model, yields_below)
===Optimal solutions based on perfect information===
Culture.         |  Wheat | Corn  | Sugar Beets |
Surface (acres)  |  100.0 | 25.0  | 375.0  |
Yield (T)        |  200.0 | 60.0 | 6000.0 |
Sales (T)        |  -0.0 | -0.0   | 6000.0 |
Purchases (T)    |  -0.0   | 180.0   | -      |
Overall profit: $ 59950.0

Running the optimization problem based on above average and below average yields gives optimal solutions:

ex1.2

The solutions again seem natural. When yields are high, smaller surfaces are needed to raise the minimum requirements in wheat and corn and the sugar beet quota. The remaining land is devoted to wheat, whose extra production is sold. When yields are low, larger surfaces are needed to raise the minimum requirements and the sugar beet quota.

Unfortunately, weather conditions cannot be accurately predicted six months ahead. The farmer must make up his mind without perfect information on yields!

Stochastic Programming

Stochastic programming optimizes when some parameters are uncertain (i.e., crop yields), but defined with a probability distribution. It is opposed to deterministic programming where all parameters are known.

Ξ\Xi: random variable, results of an ‘experiment’

ξ\xi: a realization, i.e., outcome of a simple experiment

Each realization has an associated probability p(ξ)p(\xi). Θ\Theta is the set of all possible realizations, ξ∈Θ\xi \in \Theta.

Let xx be stage 1 decision variables (make now, i.e., land allocation), yy be stage 2 decision variables (wait-and-see, i.e., buy/sell crops), a deterministic optimization problem is formed as:

min⁡f(x,y,θ)\min f(x,y,\theta)

      g(x,y,θ)≤0\ \ \ \ \ \ g(x,y,\theta) \leq 0

      u(x,y,θ)=0\ \ \ \ \ \ u(x,y,\theta) = 0

When the parameters θ\theta are uncertain, the corresponding stochastic optimization problem is formed as:

min⁡E[f(x,y,Ξ)]\min E[f(x,y,\Xi)]

      g(x,y,Ξ)≤0\ \ \ \ \ \ g(x,y,\Xi) \leq 0

      u(x,y,Ξ)=0\ \ \ \ \ \ u(x,y,\Xi) = 0

Key Concepts

Probability Review

Example: roll a 6-sided die.

  • ξ \xi is the realization (outcome) of the experiment (e.g., roll the die once).

  • p(ξ) p(\xi) represents the probability for each realization.

  • Θ\Theta represents the set of all possible realizations/outcomes (e.g., Θ\Theta = {1,2,3,4,5,6}, numbers on a die).

  • The probability density p(ξ) p(\xi) encodes the probability for all realizations ξ∈Θ\xi \in \Theta

Key properties:

  • 0≤p(ξ)≤1∀ξ∈Θ 0 \leq p(\xi) \leq 1 \quad \forall \quad \xi \in \Theta

  • ∫ξ∈Θp(ξ)dξ=1\int_{\xi \in \Theta} p(\xi)d\xi = 1 for continuous probabilities

  • ∑ξ∈Θp(ξ)=1\sum_{\xi \in \Theta} p(\xi) = 1 for discrete probabilities.

Example: rolling a 6-sided die.

  • P(ξ)=1/6 P(\xi) = 1/6 for each side of the die.

Infinite Dimensional Formulation

Continuous random variables have the key property:

0≤p(ξ)≤10 \leq p(\xi) \leq 1 for all ξ∈Θ\xi \in \Theta

∫p(ξ)dξ=1\int p(\xi) d\xi = 1

The expectation is formed as:

EΞ=∫Ξf(ξ)PΞ(ξ)dξE_{\Xi} = \int_{\Xi} f(\xi) P_{\Xi}(\xi) d\xi

Discrete (Finite Dimensional) Approximations

Discrete random variables have the key property:

0≤p(ξ)≤10 \leq p(\xi) \leq 1 for all ξ∈Θ\xi \in \Theta

∑p(ξ)=1\sum p(\xi) = 1

The expectation is formed as:

EΞ=∑Ξf(ξ)w(ξ)E_{\Xi} = \sum_{\Xi} f(\xi) w(\xi)

General Two-Stage Stochastic Programming Formulation

min⁡∑ξw(ξ)f(x,yξ,ξ) \min \sum_{\xi} w(\xi) f(x, y_\xi, \xi)

Subject to: g(x,yξ,ξ)≤0,h(x,yξ,ξ)=0 g(x, y_\xi, \xi) \leq 0, \quad h(x, y_\xi, \xi) = 0

Key features:

  • Replace expectation with summation.

  • Constraints enforced only for ξ \xi realizations.

  • y y reacts to uncertainty in ξ \xi .

Include uncertainty in the Farmer’s Problem (two-stage stochastic program)

Now the farmer wants to assess the benefits and losses of each decision in each situation.

Decisions on land assignment (x1x_1, x2x_2, x3x_3) have to be taken now, but sales and purchases (wi,i=1,...,4w_i, i=1,...,4, yj,j=1,2y_j, j=1,2) depend on the yields. This forms the two-stage stochastic program:

  1. Stage 1 decisions: land assignments (x1x_1, x2x_2, x3x_3)

  2. Uncertainty is realized

  3. Stage 2 decisions: wait-and-see (sales and purchases)

It is useful to index those decisions by a scenario index s=1,2,3s=1,2,3 according to above average, average or below average yields, respectively. This creates a new set of variables wi,s,i=1,2,3,4,s=1,2,3w_{i,s}, i=1,2,3,4, s=1,2,3 and yj,s,j=1,2,s=1,2,3.y_{j,s}, j=1,2, s=1,2,3. For e.g., w3,2w_{3,2} represents the amount of sugar beets sold at the favorable price if yields are average.

If the three scenarios have an equal probability of 1/3, the farmer’s problem is formed as:

min⁡150x1+230x2+260x3−13(170w1,1−238y1,1+150w2,1−210y2,1+36w3,1+10w4,1)−13(170w1,2−238y1,2+150w2,2−210y2,2+36w3,2+10w4,2)−13(170w1,3−238y1,3+150w2,3−210y2,3+36w3,3+10w4,3)s.t.scenario 1:x1+x2+x3≤500,3x1+y1,1−w1,1≥200,3.6x2+y2,1−w2,1≥240,w3,1+w4,1≤24x3,w3,1≤6000,scenario 2:2.5x1+y1,2−w1,2≥200,3x2+y2,2−w2,2≥240,w3,2+w4,2≤20x3,w3,2≤6000,scenario 3:2x1+y1,3−w1,3≥200,2.4x2+y2,3−w2,3≥240,w3,3+w4,3≤16x3,w3,3≤6000,x,y,w≥0.\begin{align*} \min \quad & 150x_1 + 230x_2 + 260x_3\\ & -\frac{1}{3}(170w_{1,1} - 238y_{1,1} + 150w_{2,1} - 210y_{2,1} + 36w_{3,1} + 10w_{4,1}) \\ & -\frac{1}{3}(170w_{1,2} - 238y_{1,2} + 150w_{2,2} - 210y_{2,2} + 36w_{3,2} + 10w_{4,2}) \\ & -\frac{1}{3}(170w_{1,3} - 238y_{1,3} + 150w_{2,3} - 210y_{2,3} + 36w_{3,3} + 10w_{4,3}) \\ \text{s.t.} \quad & \text{scenario 1:} \\ & x_1 + x_2 + x_3 \leq 500, \\ & 3x_1 + y_{1,1} - w_{1,1} \geq 200, \\ & 3.6x_2 + y_{2,1} - w_{2,1} \geq 240, \\ & w_{3,1} + w_{4,1} \leq 24x_3, \\ & w_{3,1} \leq 6000, \\ & \text{scenario 2:} \\ & 2.5x_1 + y_{1,2} - w_{1,2} \geq 200, \\ & 3x_2 + y_{2,2} - w_{2,2} \geq 240, \\ & w_{3,2} + w_{4,2} \leq 20x_3, \\ & w_{3,2} \leq 6000, \\ & \text{scenario 3:} \\ & 2x_1 + y_{1,3} - w_{1,3} \geq 200, \\ & 2.4x_2 + y_{2,3} - w_{2,3} \geq 240, \\ & w_{3,3} + w_{4,3} \leq 16x_3, \\ & w_{3,3} \leq 6000, \\ & x, y, w \geq 0. \end{align*}
(2)

Blocks make the two-stage structure explicit. In the code below, the stage 1 decision X sits on the model itself, and each scenario gets its own Block holding that scenario’s stage 2 decisions Y and W together with the three constraints they appear in. This is the indexed-Block pattern from Chapter 8 of Pyomo -- Optimization Modeling in Python: one rule describes a single scenario, and Pyomo builds one copy of it per scenario.

The alternative is to hang a scenario index off every stage 2 variable --- model.W[crop, scenario] --- which gives the same optimization problem, but nothing in the code then says which variables are decided now and which wait for the yield to be realized.

def build_sp_model(yields):
    """
    Code adapted from https://mpi-sppy.readthedocs.io/en/latest/examples.html#examples
    It specifies the extensive form of the two-stage stochastic programming

    Arguments:
        yields: Yield information as a list, following the rank [wheat, corn, beets]

    Return:
        model: farmer problem model
    """
    model = ConcreteModel()

    all_crops = ["WHEAT", "CORN", "BEETS"]
    purchase_crops = ["WHEAT", "CORN"]
    sell_crops = ["WHEAT", "CORN", "BEETS_FAVORABLE", "BEETS_UNFAVORABLE"]
    scenarios = ["ABOVE", "AVERAGE", "BELOW"]

    # Yield multiplier for each scenario, +/- 20% around the mean [dimensionless]
    yield_factor = {"ABOVE": 1.2, "AVERAGE": 1.0, "BELOW": 0.8}

    # Stage 1: fields allocation, chosen before the yields are known [acre]
    model.X = Var(all_crops, domain=NonNegativeReals, units=u.acre)

    # Stage 1 cost [USD]
    model.PLANTING_COST = (
        150 * USD_PER_ACRE * model.X["WHEAT"]
        + 230 * USD_PER_ACRE * model.X["CORN"]
        + 260 * USD_PER_ACRE * model.X["BEETS"]
    )

    # Stage 1 constraint: plant at most 500 acres [acre]
    model.total_land_constraint = Constraint(expr=summation(model.X) <= 500 * u.acre)

    def scenario_block_rule(b, scenario):
        """Populate the block for one scenario.

        The first argument of a block rule is the block being built, not the
        model, so the stage 1 variable is reached through `model` from the
        enclosing scope.
        """
        # Stage 2: how many tons of crops to purchase in this scenario [t]
        b.Y = Var(purchase_crops, domain=NonNegativeReals, units=u.metric_ton)
        # Stage 2: how many tons of crops to sell in this scenario [t]
        b.W = Var(sell_crops, domain=NonNegativeReals, units=u.metric_ton)

        # Purchase cost in this scenario [USD]
        b.PURCHASE_COST = (
            238 * USD_PER_TON * b.Y["WHEAT"] + 210 * USD_PER_TON * b.Y["CORN"]
        )
        # Sales revenue in this scenario [USD]
        b.SALES_REVENUE = (
            170 * USD_PER_TON * b.W["WHEAT"]
            + 150 * USD_PER_TON * b.W["CORN"]
            + 36 * USD_PER_TON * b.W["BEETS_FAVORABLE"]
            + 10 * USD_PER_TON * b.W["BEETS_UNFAVORABLE"]
        )

        # Yields realized in this scenario [t/acre]
        crop_yield = [y * yield_factor[scenario] for y in yields]

        # Satisfy the wheat feed requirement [t]
        b.wheat_constraint = Constraint(
            expr=crop_yield[0] * TON_PER_ACRE * model.X["WHEAT"]
            + b.Y["WHEAT"]
            - b.W["WHEAT"]
            >= 200 * u.metric_ton
        )
        # Satisfy the corn feed requirement [t]
        b.corn_constraint = Constraint(
            expr=crop_yield[1] * TON_PER_ACRE * model.X["CORN"]
            + b.Y["CORN"]
            - b.W["CORN"]
            >= 240 * u.metric_ton
        )
        # Beet sales cannot exceed beet production [t]
        b.beets_constraint = Constraint(
            expr=crop_yield[2] * TON_PER_ACRE * model.X["BEETS"]
            - b.W["BEETS_FAVORABLE"]
            - b.W["BEETS_UNFAVORABLE"]
            >= 0 * u.metric_ton
        )

        # Cap sales at the favorable price [t]
        b.W["BEETS_FAVORABLE"].setub(6000)

    # One block per scenario: identical structure, different yield data
    model.scenarios = Block(scenarios, rule=scenario_block_rule)

    # Objective function [USD]
    model.OBJ = Objective(
        expr=model.PLANTING_COST
        + 1 / 3 * sum(model.scenarios[s].PURCHASE_COST for s in scenarios)
        - 1 / 3 * sum(model.scenarios[s].SALES_REVENUE for s in scenarios),
        sense=minimize,
    )

    # Raises UnitsError if any constraint or the objective is inconsistent
    assert_units_consistent(model)

    return model
### calculate two-stage stochastic problem
yields_perfect = [2.5, 3, 20]
model = build_sp_model(yields_perfect)
solver = SolverFactory("ipopt")
results = solver.solve(model)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)

profit_2stage = -value(model.OBJ)

print("===Optimal solutions of two-stage stochastic problem===")
print("Culture.         | ", "Wheat |", "Corn  |", "Sugar Beets |")
print(
    "Surface (acres)  | ",
    f'{value(model.X["WHEAT"]):.1f}',
    "|",
    f'{value(model.X["CORN"]):.1f}',
    " |",
    f'{value(model.X["BEETS"]):.1f}',
    " |",
)

# One reporting block per scenario, read off the corresponding Pyomo block
for label, name, factor in [
    ("s=1 (Above average)", "ABOVE", 1.2),
    ("s=2 (Average average)", "AVERAGE", 1.0),
    ("s=3 (Below average)", "BELOW", 0.8),
]:
    b = model.scenarios[name]
    print("First stage:", label)
    print("Culture.         | ", "Wheat |", "Corn  |", "Sugar Beets |")
    print(
        "Yield (T)        | ",
        f'{value(model.X["WHEAT"])*yields_perfect[0]*factor:.1f}',
        "|",
        f'{value(model.X["CORN"])*yields_perfect[1]*factor:.1f}',
        "|",
        f'{value(model.X["BEETS"])*yields_perfect[2]*factor:.1f}',
        "|",
    )
    print(
        "Sales (T)        | ",
        f'{value(b.W["WHEAT"]):.1f}',
        "|",
        f'{value(b.W["CORN"]):.1f}',
        "  |",
        f'{value(b.W["BEETS_FAVORABLE"]) + value(b.W["BEETS_UNFAVORABLE"]):.1f}',
        "|",
    )
    print(
        "Purchases (T)    | ",
        f'{value(b.Y["WHEAT"]):.1f}',
        "  |",
        f'{value(b.Y["CORN"]):.1f}',
        "  |",
        "-",
        "     |",
    )

print("Overall profit: $", f"{profit_2stage:.1f}")
===Optimal solutions of two-stage stochastic problem===
Culture.         |  Wheat | Corn  | Sugar Beets |
Surface (acres)  |  170.0 | 80.0  | 250.0  |
First stage: s=1 (Above average)
Culture.         |  Wheat | Corn  | Sugar Beets |
Yield (T)        |  510.0 | 288.0 | 6000.0 |
Sales (T)        |  310.0 | 48.0   | 6000.0 |
Purchases (T)    |  -0.0   | -0.0   | -      |
First stage: s=2 (Average average)
Culture.         |  Wheat | Corn  | Sugar Beets |
Yield (T)        |  425.0 | 240.0 | 5000.0 |
Sales (T)        |  225.0 | -0.0   | 5000.0 |
Purchases (T)    |  -0.0   | -0.0   | -      |
First stage: s=3 (Below average)
Culture.         |  Wheat | Corn  | Sugar Beets |
Yield (T)        |  340.0 | 192.0 | 4000.0 |
Sales (T)        |  140.0 | -0.0   | 4000.0 |
Purchases (T)    |  -0.0   | 48.0   | -      |
Overall profit: $ 108390.0
ex1.2

Such a model of a stochastic decision program is known as the extensiveextensive formform of the stochastic program because it explicitly describes the second-stage decision variables for all scenarios.

This solution illustrates that it is impossible to find a solution that is ideal under all circumstances under uncertainty.

Comparing Perfect Information and the Stochastic Solutions

Perfect Information

YieldLow YieldAverage YieldHigh Yield
Profit$59,950$118,600$167,667

Average: $115,406 \$115,406 , representing the WS (wait-and-see problem)

Stochastic Programming

WS=Eξ[min⁡Xf(X,ξ)] WS = E_{\xi} \left[ \min_{X} f(X, \xi) \right]

ProfitLowAverageHigh
Stochastic$48,820$109,350$167,000

Average: $108,390 \$108,390 , representing the RP (recourse problem).

RP=min⁡XEξ[f(X,ξ)] RP = \min_{X} E_{\xi} \left[ f(X, \xi) \right]

Expected Value of Perfect Information (EVPI)

Suppose yields vary over years but are cyclical. A year with above average yields is always followed by a year with average yields and then a year with below average yields. The farmer would take optimal solutions as given in perfect information chapter respectively. The mean profit in the long run will be the mean of the three figures, namely $115,406 per year.

Now assume again the yields vary over years but on a random basis. The farmer does not get prior information on the yields. So, the best he can do in the long run is to take the solution as given in the two-stage stochastic program. The difference between this figure and the value in the case of perfect information is thethe expectedexpected valuevalue ofof perfectperfect informationinformation (EVPI). It represents the loss of profit due to the presence of uncertainty.

Note: In the tables above, we reported profit. But our objective was actually to minimize negative profit. We added the negative sign in the calculation below. This is because the formulas are defined for minimization problems.

EVPI=RP−WS=−$108,390−(−$115,406)=$7,016 EVPI = RP - WS = -\$108,390 - (-\$115,406) = \$7,016

This represents how much the farmer is willing to pay for a perfect forecast.

Another approach is to assume expected yields and always to allocate the optimal planting surface according to these yields, which represents the expected values solution. The loss by not considering the random variation is the difference between this and the stochastic model profit, which is called the value of the stochastic solution (VSS).

Expected Value Solution (EV)

How good is the stochastic solution compared to not considering uncertainty?

  • Expected Value Solution: xˉ=arg⁡min⁡Xf(X,ξˉ) \bar{x} = \arg \min_{X} f(X, \bar{\xi})

  • Expected result using EV solution: EEV=Eξ[f(Xˉ,ξ)] EEV = E_{\xi} \left[ f(\bar{X}, \xi) \right]

Farmer’s: EEV=−$107,240 EEV = -\$107,240

We get this value by recomputing the objective considering the scenario data but using the deterministic decision.

Value of Stochastic Solution (VSS)

What is the cost of ignoring uncertainty?

VSS = EEV - RP = (-$107,240) - (-$108,390) = $ 1,150

# calculated EVPI
EVPI = (profit_perfect + profit_above + profit_below) / 3 - profit_2stage

# calculate expectation value
expected = build_sp_model(yields_perfect)
# fix variables with solutions
expected.X["WHEAT"].fix(120)
expected.X["CORN"].fix(80)
expected.X["BEETS"].fix(300)
# solve the model
solver = SolverFactory("ipopt")
results = solver.solve(expected)
assert pyo.check_optimal_termination(results), (
    f"Solve failed: status={results.solver.status}, "
    f"termination={results.solver.termination_condition}"
)
# calculate expected value
profit_expect = -value(expected.OBJ)
print("Expectation:", round(profit_expect, 2), "USD")

VSS = profit_2stage - profit_expect

print("EVPI:", round(EVPI, 2), "USD")
print("VSS:", round(VSS, 2), "USD")
Expectation: 107240.0 USD
EVPI: 7015.56 USD
VSS: 1150.0 USD

EVPIEVPI measures the value of knowing the future with certainty, while VSSVSS assesses the value of knowing and using distributions on future outcomes.

Reference

Biegler, L.T., 2010. Nonlinear programming: concepts, algorithms, and applications to chemical processes. Society for Industrial and Applied Mathematics.

Birge, J.R. and Louveaux, F., 2011. Introduction to stochastic programming. Springer Science & Business Media.

Code partly adapted from:

https://mpi-sppy.readthedocs.io/en/latest/quick_start.html

https://mpi-sppy.readthedocs.io/en/latest/examples.html#examples

References
  1. Bynum, M. L., Hackebeil, G. A., Hart, W. E., Laird, C. D., Nicholson, B. L., Siirola, J. D., Watson, J.-P., & Woodruff, D. L. (2021). Structured Modeling with Blocks. In Pyomo — Optimization Modeling in Python (pp. 111–122). Springer International Publishing. 10.1007/978-3-030-68928-5_8
Optimization Modeling in Pyomo
Optimization Under Uncertainty
Optimization Modeling in Pyomo
Blocks and Other Pyomo Best Practices