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
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 Policy,
24)
25from ema_workbench.em_framework import get_SALib_problem
26from ema_workbench.em_framework.evaluators import Samplers
27
28
29def lake_problem(
30 b=0.42, # decay rate for P in lake (0.42 = irreversible)
31 q=2.0, # recycling exponent
32 mean=0.02, # mean of natural inflows
33 stdev=0.001, # future utility discount rate
34 delta=0.98, # standard deviation of natural inflows
35 alpha=0.4, # utility from pollution
36 nsamples=100, # Monte Carlo sampling of natural inflows
37 **kwargs,
38):
39 try:
40 decisions = [kwargs[str(i)] for i in range(100)]
41 except KeyError:
42 decisions = [0] * 100
43
44 Pcrit = brentq(lambda x: x**q / (1 + x**q) - b * x, 0.01, 1.5)
45 nvars = len(decisions)
46 X = np.zeros((nvars,))
47 average_daily_P = np.zeros((nvars,))
48 decisions = np.array(decisions)
49 reliability = 0.0
50
51 for _ in range(nsamples):
52 X[0] = 0.0
53
54 natural_inflows = np.random.lognormal(
55 math.log(mean**2 / math.sqrt(stdev**2 + mean**2)),
56 math.sqrt(math.log(1.0 + stdev**2 / mean**2)),
57 size=nvars,
58 )
59
60 for t in range(1, nvars):
61 X[t] = (
62 (1 - b) * X[t - 1]
63 + X[t - 1] ** q / (1 + X[t - 1] ** q)
64 + decisions[t - 1]
65 + natural_inflows[t - 1]
66 )
67 average_daily_P[t] += X[t] / float(nsamples)
68
69 reliability += np.sum(X < Pcrit) / float(nsamples * nvars)
70
71 max_P = np.max(average_daily_P)
72 utility = np.sum(alpha * decisions * np.power(delta, np.arange(nvars)))
73 inertia = np.sum(np.absolute(np.diff(decisions)) < 0.02) / float(nvars - 1)
74
75 return max_P, utility, inertia, reliability
76
77
78def analyze(results, ooi):
79 """analyze results using SALib sobol, returns a dataframe"""
80
81 _, outcomes = results
82
83 problem = get_SALib_problem(lake_model.uncertainties)
84 y = outcomes[ooi]
85 sobol_indices = sobol.analyze(problem, y)
86 sobol_stats = {key: sobol_indices[key] for key in ["ST", "ST_conf", "S1", "S1_conf"]}
87 sobol_stats = pd.DataFrame(sobol_stats, index=problem["names"])
88 sobol_stats.sort_values(by="ST", ascending=False)
89 s2 = pd.DataFrame(sobol_indices["S2"], index=problem["names"], columns=problem["names"])
90 s2_conf = pd.DataFrame(
91 sobol_indices["S2_conf"], index=problem["names"], columns=problem["names"]
92 )
93
94 return sobol_stats, s2, s2_conf
95
96
97if __name__ == "__main__":
98 ema_logging.log_to_stderr(ema_logging.INFO)
99
100 # instantiate the model
101 lake_model = Model("lakeproblem", function=lake_problem)
102 lake_model.time_horizon = 100
103
104 # specify uncertainties
105 lake_model.uncertainties = [
106 RealParameter("b", 0.1, 0.45),
107 RealParameter("q", 2.0, 4.5),
108 RealParameter("mean", 0.01, 0.05),
109 RealParameter("stdev", 0.001, 0.005),
110 RealParameter("delta", 0.93, 0.99),
111 ]
112
113 # set levers, one for each time step
114 lake_model.levers = [RealParameter(str(i), 0, 0.1) for i in range(lake_model.time_horizon)]
115
116 # specify outcomes
117 lake_model.outcomes = [
118 ScalarOutcome("max_P"),
119 ScalarOutcome("utility"),
120 ScalarOutcome("inertia"),
121 ScalarOutcome("reliability"),
122 ]
123
124 # override some of the defaults of the model
125 lake_model.constants = [Constant("alpha", 0.41), Constant("nsamples", 150)]
126
127 # generate sa single default no release policy
128 policy = Policy("no release", **{str(i): 0.1 for i in range(100)})
129
130 n_scenarios = 1000
131
132 with MultiprocessingEvaluator(lake_model) as evaluator:
133 results = evaluator.perform_experiments(
134 n_scenarios, policy, uncertainty_sampling=Samplers.SOBOL
135 )
136
137 sobol_stats, s2, s2_conf = analyze(results, "max_P")
138 print(sobol_stats)
139 print(s2)
140 print(s2_conf)