From 22edecc24357ddff6269b16b8fc05d7e5196e8ab Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Thu, 31 Aug 2023 21:03:37 +1200 Subject: [PATCH] [docs] simplify milk producer example (#660) --- docs/src/tutorial/example_milk_producer.jl | 119 +++++++++++---------- 1 file changed, 63 insertions(+), 56 deletions(-) diff --git a/docs/src/tutorial/example_milk_producer.jl b/docs/src/tutorial/example_milk_producer.jl index f73488048..ea305ca99 100644 --- a/docs/src/tutorial/example_milk_producer.jl +++ b/docs/src/tutorial/example_milk_producer.jl @@ -18,22 +18,20 @@ import Plots # A company produces milk for sale on a spot market each month. The quantity of # milk they produce is uncertain, and so too is the price on the spot market. - -# The company can store the milk, but, over time, some milk spoils and must be -# discarded. Eventually the milk expires and is worthless. +# The company can store unsold milk in a stockpile of dried milk powder. # The spot price is determined by an auction system, and so varies from month to # month, but demonstrates serial correlation. In each auction, there is # sufficient demand that the milk producer finds a buyer for all their -# widgets, regardless of the quantity they supply. Furthermore, the spot price +# milk, regardless of the quantity they supply. Furthermore, the spot price # is independent of the milk producer (they are a small player in the market). # The spot price is highly volatile, and is the result of a process that is out # of the control of the company. To counteract their price risk, the company # engages in a forward contracting programme. -# The forward contracting programme is a deal for physical milk at a future -# date in time, up to four months in the future. +# The forward contracting programme is a deal for physical milk four months in +# the future. # The futures price is the current spot price, plus some forward contango (the # buyers gain certainty that they will receive the milk in the future). @@ -91,7 +89,7 @@ plot = Plots.plot( # We can't incorporate this price process directly into SDDP.jl, but we can fit # a [`SDDP.MarkovianGraph`](@ref) directly from the simulator: -graph = SDDP.MarkovianGraph(simulator; budget = 60, scenarios = 10_000); +graph = SDDP.MarkovianGraph(simulator; budget = 30, scenarios = 10_000); nothing # hide # Here `budget` is the number of nodes in the policy graph, and `scenarios` is @@ -125,55 +123,58 @@ plot model = SDDP.PolicyGraph( graph; sense = :Max, - upper_bound = 1e4, + upper_bound = 1e2, optimizer = HiGHS.Optimizer, ) do sp, node - t, price = node - c_contango = [1.0, 1.025, 1.05, 1.075] - c_transaction, c_perish_factor, c_buy_premium = 0.01, 0.95, 1.5 - F, P = length(c_contango), 5 - @variable(sp, 0 <= x_forward[1:F], SDDP.State, initial_value = 0) - @variable(sp, 0 <= x_stock[p = 1:P], SDDP.State, initial_value = 0) - @variable(sp, 0 <= u_spot_sell[1:P]) - @variable(sp, 0 <= u_spot_buy) - @variable(sp, 0 <= u_forward_deliver[1:P]) - @variable(sp, 0 <= u_forward_sell[1:F] <= 10) - @variable(sp, 0 <= u_production) - for f in 1:F - if t + f >= 12 - fix(u_forward_sell[f], 0.0; force = true) - end - end - @constraint( - sp, - [i = 1:F-1], - x_forward[i].out == x_forward[i+1].in + u_forward_sell[i], - ) - @constraint(sp, x_forward[F].out == u_forward_sell[F]) - @constraint(sp, x_forward[1].in == sum(u_forward_deliver)) + ## Decompose the node into the month (::Int) and spot price (::Float64) + t, price = node::Tuple{Int,Float64} + ## Transactions on the futures market cost 0.01 + c_transaction = 0.01 + ## It costs the company +50% to buy milk on the spot market and deliver to + ## their customers + c_buy_premium = 1.5 + ## Buyer is willing to pay +5% for certainty + c_contango = 1.05 + ## Distribution of production + Ω_production = range(0.1, 0.2; length = 5) + c_max_production = 12 * maximum(Ω_production) + ## x_stock: quantity of milk in stock pile + @variable(sp, 0 <= x_stock, SDDP.State, initial_value = 0) + ## x_forward[i]: quantity of milk for delivery in i months + @variable(sp, 0 <= x_forward[1:4], SDDP.State, initial_value = 0) + ## u_spot_sell: quantity of milk to sell on spot market + @variable(sp, 0 <= u_spot_sell <= c_max_production) + ## u_spot_buy: quantity of milk to buy on spot market + @variable(sp, 0 <= u_spot_buy <= c_max_production) + ## u_spot_buy: quantity of milk to sell on futures market + c_max_futures = t <= 8 ? c_max_production : 0.0 + @variable(sp, 0 <= u_forward_sell <= c_max_futures) + ## ω_production: production random variable + @variable(sp, ω_production) + ## Forward contracting constraints: + @constraint(sp, [i in 1:3], x_forward[i].out == x_forward[i+1].in) + @constraint(sp, x_forward[4].out == u_forward_sell) + ## Stockpile balance constraint @constraint( sp, - x_stock[1].out == - u_production - u_forward_deliver[1] - u_spot_sell[1] + u_spot_buy + x_stock.out == + x_stock.in + ω_production + u_spot_buy - x_forward[1].in - u_spot_sell ) - @constraint( - sp, - [i = 2:P], - x_stock[i].out == - c_perish_factor * x_stock[i-1].in - u_forward_deliver[i] - - u_spot_sell[i] - ) - Ω = [(price, (production = p,)) for p in range(0.1, 0.2; length = 5)] - SDDP.parameterize(sp, Ω) do (price, ω) - set_upper_bound(u_production, ω.production) + ## The random variables. `price` comes from the Markov node + Ω = [(price, p) for p in Ω_production] + SDDP.parameterize(sp, Ω) do ω::Tuple{Float64,Float64} + ## Fix the ω_production variable + fix(ω_production, ω[2]) @stageobjective( sp, - price * sum(u_spot_sell) + - price * sum(c_contango[f] * u_forward_sell[f] for f in 1:F) - - price * c_buy_premium * u_spot_buy - - c_transaction * sum(u_forward_sell) + ## Sales on spot market + ω[1] * (u_spot_sell - c_buy_premium * u_spot_buy) + + ## Sales on futures smarket + (ω[1] * c_contango - c_transaction) * u_forward_sell ) + return end + return end # ## Training a policy @@ -188,7 +189,7 @@ end SDDP.train( model; time_limit = 20, - risk_measure = SDDP.EAVaR(; lambda = 0.5, beta = 0.25), + risk_measure = 0.5 * SDDP.Expectation() + 0.5 * SDDP.AVaR(0.25), sampling_scheme = SDDP.SimulatorSamplingScheme(simulator), ) @@ -205,7 +206,7 @@ SDDP.train( simulations = SDDP.simulate( model, 200, - Symbol[:x_forward, :x_stock, :u_spot_sell, :u_spot_buy]; + Symbol[:x_stock, :u_forward_sell, :u_spot_sell, :u_spot_buy]; sampling_scheme = SDDP.SimulatorSamplingScheme(simulator), ); nothing # hide @@ -226,18 +227,24 @@ simulations[1][12][:noise_term] # terminated the training early, so we should run the re-train the policy for # more iterations before making too many judgements). -Plots.plot( - SDDP.publication_plot(simulations; title = "sum(x_stock.out)") do data - return sum(y.out for y in data[:x_stock]) +plot = Plots.plot( + SDDP.publication_plot(simulations; title = "x_stock.out") do data + return data[:x_stock].out end, - SDDP.publication_plot(simulations; title = "sum(x_forward.out)") do data - return sum(y.out for y in data[:x_forward]) + SDDP.publication_plot(simulations; title = "u_forward_sell") do data + return data[:u_forward_sell] end, SDDP.publication_plot(simulations; title = "u_spot_buy") do data return data[:u_spot_buy] end, - SDDP.publication_plot(simulations; title = "sum(u_spot_sell)") do data - return sum(data[:u_spot_sell]) + SDDP.publication_plot(simulations; title = "u_spot_sell") do data + return data[:u_spot_sell] end; layout = (2, 2), ) + +# ## Next steps + +# * Train the policy for longer. What do you observe? +# * Try creating different Markovian graphs. What happens if you add more nodes? +# * Try different risk measures