The farmer's problem

This tutorial was generated using Literate.jl. Download the source as a .jl file. Download the source as a .ipynb file.

This problem is taken from Section 1.1 of the book Birge, J. R., & Louveaux, F. (2011). Introduction to Stochastic Programming. New York, NY: Springer New York. Paragraphs in quotes are taken verbatim.

Problem description

Consider a European farmer who specializes in raising wheat, corn, and sugar beets on his 500 acres of land. During the winter, [they want] to decide how much land to devote to each crop.
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 [they expect] 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 [their] land is roughly 2.5 T, 3 T, and 20 T per acre for wheat, corn, and sugar beets, respectively.
[To introduce uncertainty,] assume some correlation among the yields of the different crops. A very simplified representation of this 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. To fix these ideas, above and below average indicate a yield 20% above or below the mean yield.

Problem data

The area of the farm.

MAX_AREA = 500.0

500.0

There are three crops:

CROPS = [:wheat, :corn, :sugar_beet]

3-element Vector{Symbol}:
 :wheat
 :corn
 :sugar_beet

Each of the crops has a different planting cost ($/acre).

PLANTING_COST = Dict(:wheat => 150.0, :corn => 230.0, :sugar_beet => 260.0)

Dict{Symbol, Float64} with 3 entries:
  :wheat      => 150.0
  :sugar_beet => 260.0
  :corn       => 230.0

The farmer requires a minimum quantity of wheat and corn, but not of sugar beet (tonnes).

MIN_QUANTITIES = Dict(:wheat => 200.0, :corn => 240.0, :sugar_beet => 0.0)

Dict{Symbol, Float64} with 3 entries:
  :wheat      => 200.0
  :sugar_beet => 0.0
  :corn       => 240.0

In Europe, there is a quota system for producing crops. The farmer owns the following quota for each crop (tonnes):

QUOTA_MAX = Dict(:wheat => Inf, :corn => Inf, :sugar_beet => 6_000.0)

Dict{Symbol, Float64} with 3 entries:
  :wheat      => Inf
  :sugar_beet => 6000.0
  :corn       => Inf

The farmer can sell crops produced under the quota for the following amounts ($/tonne):

SELL_IN_QUOTA = Dict(:wheat => 170.0, :corn => 150.0, :sugar_beet => 36.0)

Dict{Symbol, Float64} with 3 entries:
  :wheat      => 170.0
  :sugar_beet => 36.0
  :corn       => 150.0

If they sell more than their allotted quota, the farmer earns the following on each tonne of crop above the quota ($/tonne):

SELL_NO_QUOTA = Dict(:wheat => 0.0, :corn => 0.0, :sugar_beet => 10.0)

Dict{Symbol, Float64} with 3 entries:
  :wheat      => 0.0
  :sugar_beet => 10.0
  :corn       => 0.0

The purchase prices for wheat and corn are 40% more than their sales price. However, the description does not address the purchase price of sugar beet. Therefore, we use a large value of $1,000/tonne.

BUY_PRICE = Dict(:wheat => 238.0, :corn => 210.0, :sugar_beet => 1_000.0)

Dict{Symbol, Float64} with 3 entries:
  :wheat      => 238.0
  :sugar_beet => 1000.0
  :corn       => 210.0

On average, each crop has the following yield in tonnes/acre:

MEAN_YIELD = Dict(:wheat => 2.5, :corn => 3.0, :sugar_beet => 20.0)

Dict{Symbol, Float64} with 3 entries:
  :wheat      => 2.5
  :sugar_beet => 20.0
  :corn       => 3.0

However, the yield is random. In good years, the yield is +20% above average, and in bad years, the yield is -20% below average.

YIELD_MULTIPLIER = Dict(:good => 1.2, :fair => 1.0, :bad => 0.8)

Dict{Symbol, Float64} with 3 entries:
  :bad  => 0.8
  :good => 1.2
  :fair => 1.0

Mathematical formulation

SDDP.jl code

Note

In what follows, we make heavy use of the fact that you can look up variables by their symbol name in a JuMP model as follows:

@variable(model, x)
model[:x]

Read the JuMP documentation if this isn't familiar to you.

First up, load SDDP.jl and a solver. For this example, we use HiGHS.jl.

using SDDP, HiGHS

State variables

State variables are the information that flows between stages. In our example, the state variables are the areas of land devoted to growing each crop.

function add_state_variables(subproblem)
    @variable(subproblem, area[c = CROPS] >= 0, SDDP.State, initial_value = 0)
end

add_state_variables (generic function with 1 method)

First stage problem

We can only plant a maximum of 500 acres, and we want to minimize the planting cost

function create_first_stage_problem(subproblem)
    @constraint(
        subproblem,
        sum(subproblem[:area][c].out for c in CROPS) <= MAX_AREA
    )
    @stageobjective(
        subproblem,
        -sum(PLANTING_COST[c] * subproblem[:area][c].out for c in CROPS)
    )
end

create_first_stage_problem (generic function with 1 method)

Second stage problem

Now let's consider the second stage problem. This is more complicated than the first stage, so we've broken it down into four sections:

control variables
constraints
the objective
the uncertainty

First, let's add the second stage control variables.

Variables

We add four types of control variables. Technically, the yield isn't a control variable. However, we add it as a dummy "helper" variable because it will be used when we add uncertainty.

function second_stage_variables(subproblem)
    @variables(subproblem, begin
        0 <= yield[c = CROPS]                          # tonnes/acre
        0 <= buy[c = CROPS]                            # tonnes
        0 <= sell_in_quota[c = CROPS] <= QUOTA_MAX[c]  # tonnes
        0 <= sell_no_quota[c = CROPS]                  # tonnes
    end)
end

second_stage_variables (generic function with 1 method)

Constraints

We need to define is the minimum quantity constraint. This ensures that MIN_QUANTITIES[c] of each crop is produced.

function second_stage_constraint_min_quantity(subproblem)
    @constraint(
        subproblem,
        [c = CROPS],
        subproblem[:yield][c] + subproblem[:buy][c] -
        subproblem[:sell_in_quota][c] - subproblem[:sell_no_quota][c] >=
        MIN_QUANTITIES[c]
    )
end

second_stage_constraint_min_quantity (generic function with 1 method)

Objective

The objective of the second stage is to maximise revenue from selling crops, less the cost of buying corn and wheat if necessary to meet the minimum quantity constraint.

function second_stage_objective(subproblem)
    @stageobjective(
        subproblem,
        sum(
            SELL_IN_QUOTA[c] * subproblem[:sell_in_quota][c] +
            SELL_NO_QUOTA[c] * subproblem[:sell_no_quota][c] -
            BUY_PRICE[c] * subproblem[:buy][c] for c in CROPS
        )
    )
end

second_stage_objective (generic function with 1 method)

Random variables

Then, in the SDDP.parameterize function, we set the coefficient using JuMP.set_normalized_coefficient.

function second_stage_uncertainty(subproblem)
    @constraint(
        subproblem,
        uncertainty[c = CROPS],
        1.0 * subproblem[:area][c].in == subproblem[:yield][c]
    )
    SDDP.parameterize(subproblem, [:good, :fair, :bad]) do ω
        for c in CROPS
            JuMP.set_normalized_coefficient(
                uncertainty[c],
                subproblem[:area][c].in,
                MEAN_YIELD[c] * YIELD_MULTIPLIER[ω],
            )
        end
    end
end

second_stage_uncertainty (generic function with 1 method)

Putting it all together

Now we're ready to build the multistage stochastic programming model. In addition to the things already discussed, we need a few extra pieces of information.

First, we are maximizing, so we set sense = :Max. Second, we need to provide a valid upper bound. (See Choosing an initial bound for more on this.) We know from Birge and Louveaux that the optimal solution is $108,390. So, let's choose $500,000 just to be safe.

Here is the full model.

model = SDDP.LinearPolicyGraph(;
    stages = 2,
    sense = :Max,
    upper_bound = 500_000.0,
    optimizer = HiGHS.Optimizer,
) do subproblem, stage
    add_state_variables(subproblem)
    if stage == 1
        create_first_stage_problem(subproblem)
    else
        second_stage_variables(subproblem)
        second_stage_constraint_min_quantity(subproblem)
        second_stage_uncertainty(subproblem)
        second_stage_objective(subproblem)
    end
end

A policy graph with 2 nodes.
 Node indices: 1, 2

Training a policy

Now that we've built a model, we need to train it using SDDP.train. The keyword iteration_limit stops the training after 40 iterations. See Choose a stopping rule for other ways to stop the training.

SDDP.train(model; iteration_limit = 40)

-------------------------------------------------------------------
         SDDP.jl (c) Oscar Dowson and contributors, 2017-24
-------------------------------------------------------------------
problem
  nodes           : 2
  state variables : 3
  scenarios       : 3.00000e+00
  existing cuts   : false
options
  solver          : serial mode
  risk measure    : SDDP.Expectation()
  sampling scheme : SDDP.InSampleMonteCarlo
subproblem structure
  VariableRef                             : [7, 19]
  AffExpr in MOI.EqualTo{Float64}         : [3, 3]
  AffExpr in MOI.GreaterThan{Float64}     : [3, 3]
  AffExpr in MOI.LessThan{Float64}        : [1, 1]
  VariableRef in MOI.GreaterThan{Float64} : [3, 16]
  VariableRef in MOI.LessThan{Float64}    : [1, 2]
numerical stability report
  matrix range     [1e+00, 2e+01]
  objective range  [1e+00, 1e+03]
  bounds range     [6e+03, 5e+05]
  rhs range        [2e+02, 5e+02]
-------------------------------------------------------------------
 iteration    simulation      bound        time (s)     solves  pid
-------------------------------------------------------------------
         1  -9.800000e+04  4.922260e+05  8.943510e-02         6   1
        40   1.093500e+05  1.083900e+05  1.190240e-01       240   1
-------------------------------------------------------------------
status         : iteration_limit
total time (s) : 1.190240e-01
total solves   : 240
best bound     :  1.083900e+05
simulation ci  :  8.444504e+04 ± 1.819355e+04
numeric issues : 0
-------------------------------------------------------------------

Checking the policy

Birge and Louveaux report that the optimal objective value is $108,390. Check that we got the correct solution using SDDP.calculate_bound:

@assert isapprox(SDDP.calculate_bound(model), 108_390.0, atol = 0.1)