1"""
2An example of the lake problem using the ema workbench. This example
3illustrated how you can control more finely how samples are being generated.
4In this particular case, we want to apply Sobol analysis over both the
5uncertainties and levers at the same time.
6
7"""
8
9import math
10
11import numpy as np
12import pandas as pd
13from SALib.analyze import sobol
14from scipy.optimize import brentq
15
16from ema_workbench import (
17 Model,
18 RealParameter,
19 ScalarOutcome,
20 Constant,
21 ema_logging,
22 MultiprocessingEvaluator,
23 Scenario,
24)
25from ema_workbench.em_framework import get_SALib_problem, sample_parameters
26from ema_workbench.em_framework import SobolSampler
27
28# from ema_workbench.em_framework.evaluators import Samplers
29
30
31def get_antropogenic_release(xt, c1, c2, r1, r2, w1):
32 """
33
34 Parameters
35 ----------
36 xt : float
37 pollution in lake at time t
38 c1 : float
39 center rbf 1
40 c2 : float
41 center rbf 2
42 r1 : float
43 ratius rbf 1
44 r2 : float
45 ratius rbf 2
46 w1 : float
47 weight of rbf 1
48
49 note:: w2 = 1 - w1
50
51 """
52
53 rule = w1 * (abs(xt - c1) / r1) ** 3 + (1 - w1) * (abs(xt - c2) / r2) ** 3
54 at1 = max(rule, 0.01)
55 at = min(at1, 0.1)
56
57 return at
58
59
60def lake_problem(
61 b=0.42, # decay rate for P in lake (0.42 = irreversible)
62 q=2.0, # recycling exponent
63 mean=0.02, # mean of natural inflows
64 stdev=0.001, # future utility discount rate
65 delta=0.98, # standard deviation of natural inflows
66 alpha=0.4, # utility from pollution
67 nsamples=100, # Monte Carlo sampling of natural inflows
68 myears=1, # the runtime of the simulation model
69 c1=0.25,
70 c2=0.25,
71 r1=0.5,
72 r2=0.5,
73 w1=0.5,
74):
75 Pcrit = brentq(lambda x: x**q / (1 + x**q) - b * x, 0.01, 1.5)
76
77 X = np.zeros((myears,))
78 average_daily_P = np.zeros((myears,))
79 reliability = 0.0
80 inertia = 0
81 utility = 0
82
83 for _ in range(nsamples):
84 X[0] = 0.0
85 decision = 0.1
86
87 decisions = np.zeros(myears)
88 decisions[0] = decision
89
90 natural_inflows = np.random.lognormal(
91 math.log(mean**2 / math.sqrt(stdev**2 + mean**2)),
92 math.sqrt(math.log(1.0 + stdev**2 / mean**2)),
93 size=myears,
94 )
95
96 for t in range(1, myears):
97 # here we use the decision rule
98 decision = get_antropogenic_release(X[t - 1], c1, c2, r1, r2, w1)
99 decisions[t] = decision
100
101 X[t] = (
102 (1 - b) * X[t - 1]
103 + X[t - 1] ** q / (1 + X[t - 1] ** q)
104 + decision
105 + natural_inflows[t - 1]
106 )
107 average_daily_P[t] += X[t] / nsamples
108
109 reliability += np.sum(X < Pcrit) / (nsamples * myears)
110 inertia += np.sum(np.absolute(np.diff(decisions) < 0.02)) / (nsamples * myears)
111 utility += np.sum(alpha * decisions * np.power(delta, np.arange(myears))) / nsamples
112 max_P = np.max(average_daily_P)
113
114 return max_P, utility, inertia, reliability
115
116
117def analyze(results, ooi):
118 """analyze results using SALib sobol, returns a dataframe"""
119
120 _, outcomes = results
121
122 parameters = lake_model.uncertainties.copy() + lake_model.levers.copy()
123
124 problem = get_SALib_problem(parameters)
125 y = outcomes[ooi]
126 sobol_indices = sobol.analyze(problem, y)
127 sobol_stats = {key: sobol_indices[key] for key in ["ST", "ST_conf", "S1", "S1_conf"]}
128 sobol_stats = pd.DataFrame(sobol_stats, index=problem["names"])
129 sobol_stats.sort_values(by="ST", ascending=False)
130 s2 = pd.DataFrame(sobol_indices["S2"], index=problem["names"], columns=problem["names"])
131 s2_conf = pd.DataFrame(
132 sobol_indices["S2_conf"], index=problem["names"], columns=problem["names"]
133 )
134
135 return sobol_stats, s2, s2_conf
136
137
138if __name__ == "__main__":
139 ema_logging.log_to_stderr(ema_logging.INFO)
140
141 # instantiate the model
142 lake_model = Model("lakeproblem", function=lake_problem)
143 # specify uncertainties
144 lake_model.uncertainties = [
145 RealParameter("b", 0.1, 0.45),
146 RealParameter("q", 2.0, 4.5),
147 RealParameter("mean", 0.01, 0.05),
148 RealParameter("stdev", 0.001, 0.005),
149 RealParameter("delta", 0.93, 0.99),
150 ]
151
152 # set levers
153 lake_model.levers = [
154 RealParameter("c1", -2, 2),
155 RealParameter("c2", -2, 2),
156 RealParameter("r1", 0, 2),
157 RealParameter("r2", 0, 2),
158 RealParameter("w1", 0, 1),
159 ]
160 # specify outcomes
161 lake_model.outcomes = [
162 ScalarOutcome("max_P", kind=ScalarOutcome.MINIMIZE),
163 # @UndefinedVariable
164 ScalarOutcome("utility", kind=ScalarOutcome.MAXIMIZE),
165 # @UndefinedVariable
166 ScalarOutcome("inertia", kind=ScalarOutcome.MAXIMIZE),
167 # @UndefinedVariable
168 ScalarOutcome("reliability", kind=ScalarOutcome.MAXIMIZE),
169 ] # @UndefinedVariable
170
171 # override some of the defaults of the model
172 lake_model.constants = [
173 Constant("alpha", 0.41),
174 Constant("nsamples", 100),
175 Constant("myears", 100),
176 ]
177
178 # combine parameters and uncertainties prior to sampling
179 n_scenarios = 1000
180 parameters = lake_model.uncertainties + lake_model.levers
181 scenarios = sample_parameters(parameters, n_scenarios, SobolSampler(), Scenario)
182
183 with MultiprocessingEvaluator(lake_model) as evaluator:
184 results = evaluator.perform_experiments(scenarios)
185
186 sobol_stats, s2, s2_conf = analyze(results, "max_P")
187 print(sobol_stats)
188 print(s2)
189 print(s2_conf)