1"""
2This example takes the direct policy search formulation of the lake problem as
3found in Quinn et al (2017), but embeds in in a robust optimization.
4
5Quinn, J.D., Reed, P.M., Keller, K. (2017)
6Direct policy search for robust multi-objective management of deeply
7uncertain socio-ecological tipping points. Environmental Modelling &
8Software 92, 125-141.
9
10
11"""
12
13import math
14
15import numpy as np
16from scipy.optimize import brentq
17
18from ema_workbench import (
19 Model,
20 RealParameter,
21 ScalarOutcome,
22 Constant,
23 ema_logging,
24 MultiprocessingEvaluator,
25)
26from ema_workbench.em_framework.samplers import sample_uncertainties
27
28# Created on 1 Jun 2017
29#
30# .. codeauthor::jhkwakkel <j.h.kwakkel (at) tudelft (dot) nl>
31
32__all__ = []
33
34
35def get_antropogenic_release(xt, c1, c2, r1, r2, w1):
36 """
37
38 Parameters
39 ----------
40 xt : float
41 pollution in lake at time t
42 c1 : float
43 center rbf 1
44 c2 : float
45 center rbf 2
46 r1 : float
47 ratius rbf 1
48 r2 : float
49 ratius rbf 2
50 w1 : float
51 weight of rbf 1
52
53 note:: w2 = 1 - w1
54
55 """
56
57 rule = w1 * (abs(xt - c1 / r1)) ** 3 + (1 - w1) * (abs(xt - c2 / r2) ** 3)
58 at1 = max(rule, 0.01)
59 at = min(at1, 0.1)
60
61 return at
62
63
64def lake_problem(
65 b=0.42, # decay rate for P in lake (0.42 = irreversible)
66 q=2.0, # recycling exponent
67 mean=0.02, # mean of natural inflows
68 stdev=0.001, # future utility discount rate
69 delta=0.98, # standard deviation of natural inflows
70 alpha=0.4, # utility from pollution
71 nsamples=100, # Monte Carlo sampling of natural inflows
72 myears=1, # the runtime of the simulation model
73 c1=0.25,
74 c2=0.25,
75 r1=0.5,
76 r2=0.5,
77 w1=0.5,
78):
79 Pcrit = brentq(lambda x: x**q / (1 + x**q) - b * x, 0.01, 1.5)
80
81 X = np.zeros(myears)
82 average_daily_P = np.zeros(myears)
83 reliability = 0.0
84 inertia = 0
85 utility = 0
86
87 for _ in range(nsamples):
88 X[0] = 0.0
89 decision = 0.1
90
91 decisions = np.zeros(myears)
92 decisions[0] = decision
93
94 natural_inflows = np.random.lognormal(
95 math.log(mean**2 / math.sqrt(stdev**2 + mean**2)),
96 math.sqrt(math.log(1.0 + stdev**2 / mean**2)),
97 size=myears,
98 )
99
100 for t in range(1, myears):
101 # here we use the decision rule
102 decision = get_antropogenic_release(X[t - 1], c1, c2, r1, r2, w1)
103 decisions[t] = decision
104
105 X[t] = (
106 (1 - b) * X[t - 1]
107 + X[t - 1] ** q / (1 + X[t - 1] ** q)
108 + decision
109 + natural_inflows[t - 1]
110 )
111 average_daily_P[t] += X[t] / nsamples
112
113 reliability += np.sum(X < Pcrit) / (nsamples * myears)
114 inertia += np.sum(np.absolute(np.diff(decisions) < 0.02)) / (nsamples * myears)
115 utility += np.sum(alpha * decisions * np.power(delta, np.arange(myears))) / nsamples
116 max_P = np.max(average_daily_P)
117
118 return max_P, utility, inertia, reliability
119
120
121if __name__ == "__main__":
122 ema_logging.log_to_stderr(ema_logging.INFO)
123
124 # instantiate the model
125 lake_model = Model("lakeproblem", function=lake_problem)
126 # specify uncertainties
127 lake_model.uncertainties = [
128 RealParameter("b", 0.1, 0.45),
129 RealParameter("q", 2.0, 4.5),
130 RealParameter("mean", 0.01, 0.05),
131 RealParameter("stdev", 0.001, 0.005),
132 RealParameter("delta", 0.93, 0.99),
133 ]
134
135 # set levers
136 lake_model.levers = [
137 RealParameter("c1", -2, 2),
138 RealParameter("c2", -2, 2),
139 RealParameter("r1", 0, 2),
140 RealParameter("r2", 0, 2),
141 RealParameter("w1", 0, 1),
142 ]
143
144 # specify outcomes
145 lake_model.outcomes = [
146 ScalarOutcome("max_P"),
147 ScalarOutcome("utility"),
148 ScalarOutcome("inertia"),
149 ScalarOutcome("reliability"),
150 ]
151
152 # override some of the defaults of the model
153 lake_model.constants = [
154 Constant("alpha", 0.41),
155 Constant("nsamples", 100),
156 Constant("myears", 100),
157 ]
158
159 # setup and execute the robust optimization
160 def signal_to_noise(data):
161 mean = np.mean(data)
162 std = np.std(data)
163 sn = mean / std
164 return sn
165
166 MAXIMIZE = ScalarOutcome.MAXIMIZE # @UndefinedVariable
167 MINIMIZE = ScalarOutcome.MINIMIZE # @UndefinedVariable
168 robustnes_functions = [
169 ScalarOutcome("mean p", kind=MINIMIZE, variable_name="max_P", function=np.mean),
170 ScalarOutcome("std p", kind=MINIMIZE, variable_name="max_P", function=np.std),
171 ScalarOutcome(
172 "sn reliability", kind=MAXIMIZE, variable_name="reliability", function=signal_to_noise
173 ),
174 ]
175 n_scenarios = 10
176 scenarios = sample_uncertainties(lake_model, n_scenarios)
177 nfe = 1000
178
179 with MultiprocessingEvaluator(lake_model) as evaluator:
180 evaluator.robust_optimize(
181 robustnes_functions,
182 scenarios,
183 nfe=nfe,
184 epsilons=[0.1] * len(robustnes_functions),
185 population_size=5,
186 )