-
Notifications
You must be signed in to change notification settings - Fork 59
/
example_milk_producer.jl
255 lines (209 loc) · 9.14 KB
/
example_milk_producer.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
# Copyright (c) 2017-23, Oscar Dowson and SDDP.jl contributors. #src
# This Source Code Form is subject to the terms of the Mozilla Public #src
# License, v. 2.0. If a copy of the MPL was not distributed with this #src
# file, You can obtain one at http://mozilla.org/MPL/2.0/. #src
# # Example: the milk producer
# The purpose of this tutorial is to demonstrate how to fit a Markovian policy
# graph to a univariate stochastic process.
# This tutorial uses the following packages:
using SDDP
import HiGHS
import Plots
# ## Background
# 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 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
# 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 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).
# In general, the milk company should forward contract (since they reduce
# their price risk), however they also have production risk. Therefore, it may
# be the case that they forward contract a fixed amount, but find that they do
# not produce enough milk to meet the fixed demand. They are then forced to
# buy additional milk on the spot market.
# The goal of the milk company is to choose the extent to which they forward
# contract in order to maximise (risk-adjusted) revenues, whilst managing their
# production risk.
# ## A stochastic process for price
# It is outside the scope of this tutorial, but assume that we have gone away
# and analysed historical data to fit a stochastic process to the sequence of
# monthly auction spot prices.
# One plausible model is a multiplicative auto-regressive model of order one,
# where the white noise term is modeled by a finite distribution of empirical
# residuals. We can simulate this stochastic process as follows:
function simulator()
residuals = [0.0987, 0.199, 0.303, 0.412, 0.530, 0.661, 0.814, 1.010, 1.290]
residuals = 0.1 * vcat(-residuals, 0.0, residuals)
scenario = zeros(12)
y, μ, α = 4.5, 6.0, 0.05
for t in 1:12
y = exp((1 - α) * log(y) + α * log(μ) + rand(residuals))
scenario[t] = clamp(y, 3.0, 9.0)
end
return scenario
end
simulator()
# It may be helpful to visualize a number of simulations of the price process:
plot = Plots.plot(
[simulator() for _ in 1:500];
color = "gray",
opacity = 0.2,
legend = false,
xlabel = "Month",
ylabel = "Price [\$/kg]",
xlims = (1, 12),
ylims = (3, 9),
)
# The prices gradually revert to the mean of \$6/kg, and there is high
# volatility.
# 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 = 30, scenarios = 10_000);
nothing # hide
# Here `budget` is the number of nodes in the policy graph, and `scenarios` is
# the number of simulations to use when estimating the transition probabilities.
# The graph contains too many nodes to be show, but we can plot it:
for ((t, price), edges) in graph.nodes
for ((t′, price′), probability) in edges
Plots.plot!(
plot,
[t, t′],
[price, price′];
color = "red",
width = 3 * probability,
)
end
end
plot
# That looks okay. Try changing `budget` and `scenarios` to see how different
# Markovian policy graphs can be created.
# ## Model
# Now that we have a Markovian graph, we can build the model. See if you can
# work out how we arrived at this formulation by reading the background
# description. Do all the variables and constraints make sense?
model = SDDP.PolicyGraph(
graph;
sense = :Max,
upper_bound = 1e2,
optimizer = HiGHS.Optimizer,
) do sp, node
## 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.out ==
x_stock.in + ω_production + u_spot_buy - x_forward[1].in - u_spot_sell
)
## The random variables. `price` comes from the Markov node
##
## !!! warning
## The elements in Ω MUST be a tuple with 1 or 2 values, where the first
## value is `price` and the second value is the random variable for the
## current node. If the node is deterministic, use Ω = [(price,)].
Ω = [(price, p) for p in Ω_production]
SDDP.parameterize(sp, Ω) do ω
## Fix the ω_production variable
fix(ω_production, ω[2])
@stageobjective(
sp,
## 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
# Now we have a model, we train a policy. The [`SDDP.SimulatorSamplingScheme`](@ref)
# is used in the forward pass. It generates an out-of-sample sequence of prices
# using `simulator` and traverses the closest sequence of nodes in the policy
# graph. When calling [`SDDP.parameterize`](@ref) for each subproblem, it uses
# the new out-of-sample price instead of the price associated with the Markov
# node.
SDDP.train(
model;
time_limit = 20,
risk_measure = 0.5 * SDDP.Expectation() + 0.5 * SDDP.AVaR(0.25),
sampling_scheme = SDDP.SimulatorSamplingScheme(simulator),
)
# !!! warning
# We're intentionally terminating the training early so that the
# documentation doesn't take too long to build. If you run this example
# locally, increase the time limit.
# ## Simulating the policy
# When simulating the policy, we can also use the
# [`SDDP.SimulatorSamplingScheme`](@ref).
simulations = SDDP.simulate(
model,
200,
Symbol[:x_stock, :u_forward_sell, :u_spot_sell, :u_spot_buy];
sampling_scheme = SDDP.SimulatorSamplingScheme(simulator),
);
nothing # hide
# To show how the sampling scheme uses the new out-of-sample price instead of
# the price associated with the Markov node, compare the index of the Markov
# state visited in stage 12 of the first simulation:
simulations[1][12][:node_index]
# to the realization of the noise `(price, ω)` passed to [`SDDP.parameterize`](@ref):
simulations[1][12][:noise_term]
# ## Visualizing the policy
# Finally, we can plot the policy to gain insight (although note that we
# terminated the training early, so we should run the re-train the policy for
# more iterations before making too many judgements).
plot = Plots.plot(
SDDP.publication_plot(simulations; title = "x_stock.out") do data
return data[:x_stock].out
end,
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 = "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