1"""
2An example of the lake problem using the ema workbench.
3
4The model itself is adapted from the Rhodium example by Dave Hadka,
5see https://gist.github.com/dhadka/a8d7095c98130d8f73bc
6
7"""
8
9import math
10
11import numpy as np
12from scipy.optimize import brentq
13
14from ema_workbench import (
15 Model,
16 RealParameter,
17 ScalarOutcome,
18 ema_logging,
19 MultiprocessingEvaluator,
20 Constraint,
21)
22from ema_workbench.em_framework.optimization import HyperVolume, EpsilonProgress
23
24
25def lake_problem(
26 b=0.42,
27 q=2.0,
28 mean=0.02,
29 stdev=0.0017,
30 delta=0.98,
31 alpha=0.4,
32 nsamples=100,
33 l0=0,
34 l1=0,
35 l2=0,
36 l3=0,
37 l4=0,
38 l5=0,
39 l6=0,
40 l7=0,
41 l8=0,
42 l9=0,
43 l10=0,
44 l11=0,
45 l12=0,
46 l13=0,
47 l14=0,
48 l15=0,
49 l16=0,
50 l17=0,
51 l18=0,
52 l19=0,
53 l20=0,
54 l21=0,
55 l22=0,
56 l23=0,
57 l24=0,
58 l25=0,
59 l26=0,
60 l27=0,
61 l28=0,
62 l29=0,
63 l30=0,
64 l31=0,
65 l32=0,
66 l33=0,
67 l34=0,
68 l35=0,
69 l36=0,
70 l37=0,
71 l38=0,
72 l39=0,
73 l40=0,
74 l41=0,
75 l42=0,
76 l43=0,
77 l44=0,
78 l45=0,
79 l46=0,
80 l47=0,
81 l48=0,
82 l49=0,
83 l50=0,
84 l51=0,
85 l52=0,
86 l53=0,
87 l54=0,
88 l55=0,
89 l56=0,
90 l57=0,
91 l58=0,
92 l59=0,
93 l60=0,
94 l61=0,
95 l62=0,
96 l63=0,
97 l64=0,
98 l65=0,
99 l66=0,
100 l67=0,
101 l68=0,
102 l69=0,
103 l70=0,
104 l71=0,
105 l72=0,
106 l73=0,
107 l74=0,
108 l75=0,
109 l76=0,
110 l77=0,
111 l78=0,
112 l79=0,
113 l80=0,
114 l81=0,
115 l82=0,
116 l83=0,
117 l84=0,
118 l85=0,
119 l86=0,
120 l87=0,
121 l88=0,
122 l89=0,
123 l90=0,
124 l91=0,
125 l92=0,
126 l93=0,
127 l94=0,
128 l95=0,
129 l96=0,
130 l97=0,
131 l98=0,
132 l99=0,
133):
134 decisions = np.array(
135 [
136 l0,
137 l1,
138 l2,
139 l3,
140 l4,
141 l5,
142 l6,
143 l7,
144 l8,
145 l9,
146 l10,
147 l11,
148 l12,
149 l13,
150 l14,
151 l15,
152 l16,
153 l17,
154 l18,
155 l19,
156 l20,
157 l21,
158 l22,
159 l23,
160 l24,
161 l25,
162 l26,
163 l27,
164 l28,
165 l29,
166 l30,
167 l31,
168 l32,
169 l33,
170 l34,
171 l35,
172 l36,
173 l37,
174 l38,
175 l39,
176 l40,
177 l41,
178 l42,
179 l43,
180 l44,
181 l45,
182 l46,
183 l47,
184 l48,
185 l49,
186 l50,
187 l51,
188 l52,
189 l53,
190 l54,
191 l55,
192 l56,
193 l57,
194 l58,
195 l59,
196 l60,
197 l61,
198 l62,
199 l63,
200 l64,
201 l65,
202 l66,
203 l67,
204 l68,
205 l69,
206 l70,
207 l71,
208 l72,
209 l73,
210 l74,
211 l75,
212 l76,
213 l77,
214 l78,
215 l79,
216 l80,
217 l81,
218 l82,
219 l83,
220 l84,
221 l85,
222 l86,
223 l87,
224 l88,
225 l89,
226 l90,
227 l91,
228 l92,
229 l93,
230 l94,
231 l95,
232 l96,
233 l97,
234 l98,
235 l99,
236 ]
237 )
238
239 Pcrit = brentq(lambda x: x**q / (1 + x**q) - b * x, 0.01, 1.5)
240 nvars = len(decisions)
241 X = np.zeros((nvars,))
242 average_daily_P = np.zeros((nvars,))
243 decisions = np.array(decisions)
244 reliability = 0.0
245
246 for _ in range(nsamples):
247 X[0] = 0.0
248
249 natural_inflows = np.random.lognormal(
250 math.log(mean**2 / math.sqrt(stdev**2 + mean**2)),
251 math.sqrt(math.log(1.0 + stdev**2 / mean**2)),
252 size=nvars,
253 )
254
255 for t in range(1, nvars):
256 X[t] = (
257 (1 - b) * X[t - 1]
258 + X[t - 1] ** q / (1 + X[t - 1] ** q)
259 + decisions[t - 1]
260 + natural_inflows[t - 1]
261 )
262 average_daily_P[t] += X[t] / float(nsamples)
263
264 reliability += np.sum(X < Pcrit) / float(nsamples * nvars)
265
266 max_P = np.max(average_daily_P)
267 utility = np.sum(alpha * decisions * np.power(delta, np.arange(nvars)))
268 inertia = np.sum(np.abs(np.diff(decisions)) > 0.02) / float(nvars - 1)
269
270 return max_P, utility, inertia, reliability
271
272
273if __name__ == "__main__":
274 ema_logging.log_to_stderr(ema_logging.INFO)
275
276 # instantiate the model
277 lake_model = Model("lakeproblem", function=lake_problem)
278 lake_model.time_horizon = 100 # used to specify the number of timesteps
279
280 # specify uncertainties
281 lake_model.uncertainties = [
282 RealParameter("mean", 0.01, 0.05),
283 RealParameter("stdev", 0.001, 0.005),
284 RealParameter("b", 0.1, 0.45),
285 RealParameter("q", 2.0, 4.5),
286 RealParameter("delta", 0.93, 0.99),
287 ]
288
289 # set levers, one for each time step
290 lake_model.levers = [RealParameter(f"l{i}", 0, 0.1) for i in range(lake_model.time_horizon)]
291
292 # specify outcomes
293 # specify outcomes
294 lake_model.outcomes = [
295 ScalarOutcome("max_P", kind=ScalarOutcome.MINIMIZE, expected_range=(0, 5)),
296 ScalarOutcome("utility", kind=ScalarOutcome.MAXIMIZE, expected_range=(0, 2)),
297 ScalarOutcome("inertia", kind=ScalarOutcome.MAXIMIZE, expected_range=(0, 1)),
298 ScalarOutcome("reliability", kind=ScalarOutcome.MAXIMIZE, expected_range=(0, 1)),
299 ]
300
301 convergence_metrics = [EpsilonProgress()]
302
303 constraints = [
304 Constraint("max pollution", outcome_names="max_P", function=lambda x: max(0, x - 5))
305 ]
306
307 with MultiprocessingEvaluator(lake_model) as evaluator:
308 results, convergence = evaluator.optimize(
309 nfe=5000,
310 searchover="levers",
311 epsilons=[0.125, 0.05, 0.01, 0.01],
312 convergence=convergence_metrics,
313 constraints=constraints,
314 )