Examples

Farmer problem

The following defines the well-known "Farmer problem", first outlined in Introduction to Stochastic Programming, in StochasticPrograms. The problem revolves around a farmer who needs to decide how to partition his land to sow three different crops. The uncertainty comes from not knowing what the future yield of each crop will be. Recourse decisions involve purchasing/selling crops at the market.

using StochasticPrograms
using GLPKMathProgInterface

We begin by introducing some variable indices:

Crops = [:wheat, :corn, :beets];
Purchased = [:wheat, :corn];
Sold = [:wheat,:corn,:beets_quota,:beets_extra];

The price of beets drops after a certain quantity (6000), so we introduce an extra variable to handle the excess beets. Using the variable indices, we define the deterministic problem parameters:

Cost = Dict(:wheat=>150, :corn=>230, :beets=>260);
Required = Dict(:wheat=>200, :corn=>240, :beets=>0);
PurchasePrice = Dict(:wheat=>238, :corn=>210);
SellPrice = Dict(:wheat=>170, :corn=>150, :beets_quota=>36, :beets_extra=>10);
Budget = 500;

In the first stage, the farmer needs to know what crops to plant, the cost of planting them, and the available land. Therefore, we introduce the first stage data:

first_stage_data = (Crops, Cost, Budget)

(Symbol[:wheat, :corn, :beets], Dict(:wheat=>150,:beets=>260,:corn=>230), 500)

In the second stage, the farmer needs to know the required quantity of each crop, the purchase price, and the sell price:

second_stage_data = (Required, PurchasePrice, SellPrice)

(Dict(:wheat=>200,:beets=>0,:corn=>240), Dict(:wheat=>238,:corn=>210), Dict(:wheat=>170,:beets_extra=>10,:corn=>150,:beets_quota=>36))

The uncertainty lies in the future yield of each crop. We define a scenario type to capture this:

@scenario Farmer = begin
    Yield::Dict{Symbol, Float64}

    @zero begin
        return FarmerScenario(Dict(:wheat=>0., :corn=>0., :beets=>0.))
    end

    @expectation begin
        return FarmerScenario(Dict(:wheat=>sum([probability(s)*s.Yield[:wheat] for s in scenarios]),
                                   :corn=>sum([probability(s)*s.Yield[:corn] for s in scenarios]),
                                   :beets=>sum([probability(s)*s.Yield[:beets] for s in scenarios])))
    end
end

We provide an implementation of expected since it can not be autogenerated for the internal Dict type. The three predicted outcomes can be defined through:

s₁ = FarmerScenario(Dict(:wheat=>3.0,:corn=>3.6,:beets=>24.0), probability = 1/3)
s₂ = FarmerScenario(Dict(:wheat=>2.5,:corn=>3.0,:beets=>20.0), probability = 1/3)
s₃ = FarmerScenario(Dict(:wheat=>2.0,:corn=>2.4,:beets=>16.0), probability = 1/3)

FarmerScenario with probability 0.3333333333333333

Now, we create a stochastic program with the defined data:

farmer = StochasticProgram(first_stage_data, second_stage_data, [s₁,s₂,s₃], solver=GLPKSolverLP())

Stochastic program with:
 * 3 scenarios of type FarmerScenario
 * undefined first stage
 * undefined second stage
Solver is GLPKInterfaceLP

Finally, we define the optimization models:

@first_stage farmer = begin
    (Crops,Cost,Budget) = stage
    @variable(model, x[c = Crops] >= 0)
    @objective(model, Min, sum(Cost[c]*x[c] for c in Crops))
    @constraint(model, sum(x[c] for c in Crops) <= Budget)
end
@second_stage farmer = begin
    @decision x
    (Required, PurchasePrice, SellPrice) = stage
    @variable(model, y[p = Purchased] >= 0)
    @variable(model, w[s = Sold] >= 0)
    @objective(model, Min, sum( PurchasePrice[p] * y[p] for p = Purchased) - sum( SellPrice[s] * w[s] for s in Sold))

    @constraint(model, const_minreq[p=Purchased],
                   scenario.Yield[p] * x[p] + y[p] - w[p] >= Required[p])
    @constraint(model, const_minreq_beets,
                   scenario.Yield[:beets] * x[:beets] - w[:beets_quota] - w[:beets_extra] >= Required[:beets])
    @constraint(model, const_aux, w[:beets_quota] <= 6000)
end

Stochastic program with:
 * 3 scenarios of type FarmerScenario
 * 3 decision variables
 * 6 recourse variables
Solver is GLPKInterfaceLP

We can now optimize the model:

optimize!(farmer)
x = optimal_decision(farmer, :x)
println("Wheat: $(x[:wheat])")
println("Corn: $(x[:corn])")
println("Beets: $(x[:beets])")
println("Profit: $(optimal_value(farmer))")

Wheat: 170.00000000000006
Corn: 80.00000000000001
Beets: 250.0
Profit: -108390.00000000001

Finally, we calculate the stochastic performance of the model:

println("EVPI: $(EVPI(farmer))")
println("VSS: $(VSS(farmer))")

EVPI: 7015.555555555533
VSS: 1150.0000000000437