1"""
2
3"""
4
5# Created on 12 Mar 2020
6#
7# .. codeauthor:: jhkwakkel
8
9import bisect
10import functools
11import math
12import operator
13
14import numpy as np
15import scipy as sp
16
17from ema_workbench import Model, RealParameter, ScalarOutcome, ema_logging, MultiprocessingEvaluator
18
19##==============================================================================
20## Implement the model described by Eijgenraam et al. (2012)
21# code taken from Rhodium Eijgenraam example
22##------------------------------------------------------------------------------
23
24# Parameters pulled from the paper describing each dike ring
25params = ("c", "b", "lam", "alpha", "eta", "zeta", "V0", "P0", "max_Pf")
26raw_data = {
27 10: (16.6939, 0.6258, 0.0014, 0.033027, 0.320, 0.003774, 1564.9, 0.00044, 1 / 2000),
28 11: (42.6200, 1.7068, 0.0000, 0.032000, 0.320, 0.003469, 1700.1, 0.00117, 1 / 2000),
29 15: (125.6422, 1.1268, 0.0098, 0.050200, 0.760, 0.003764, 11810.4, 0.00137, 1 / 2000),
30 16: (324.6287, 2.1304, 0.0100, 0.057400, 0.760, 0.002032, 22656.5, 0.00110, 1 / 2000),
31 22: (154.4388, 0.9325, 0.0066, 0.070000, 0.620, 0.002893, 9641.1, 0.00055, 1 / 2000),
32 23: (26.4653, 0.5250, 0.0034, 0.053400, 0.800, 0.002031, 61.6, 0.00137, 1 / 2000),
33 24: (71.6923, 1.0750, 0.0059, 0.043900, 1.060, 0.003733, 2706.4, 0.00188, 1 / 2000),
34 35: (49.7384, 0.6888, 0.0088, 0.036000, 1.060, 0.004105, 4534.7, 0.00196, 1 / 2000),
35 38: (24.3404, 0.7000, 0.0040, 0.025321, 0.412, 0.004153, 3062.6, 0.00171, 1 / 1250),
36 41: (58.8110, 0.9250, 0.0033, 0.025321, 0.422, 0.002749, 10013.1, 0.00171, 1 / 1250),
37 42: (21.8254, 0.4625, 0.0019, 0.026194, 0.442, 0.001241, 1090.8, 0.00171, 1 / 1250),
38 43: (340.5081, 4.2975, 0.0043, 0.025321, 0.448, 0.002043, 19767.6, 0.00171, 1 / 1250),
39 44: (24.0977, 0.7300, 0.0054, 0.031651, 0.316, 0.003485, 37596.3, 0.00033, 1 / 1250),
40 45: (3.4375, 0.1375, 0.0069, 0.033027, 0.320, 0.002397, 10421.2, 0.00016, 1 / 1250),
41 47: (8.7813, 0.3513, 0.0026, 0.029000, 0.358, 0.003257, 1369.0, 0.00171, 1 / 1250),
42 48: (35.6250, 1.4250, 0.0063, 0.023019, 0.496, 0.003076, 7046.4, 0.00171, 1 / 1250),
43 49: (20.0000, 0.8000, 0.0046, 0.034529, 0.304, 0.003744, 823.3, 0.00171, 1 / 1250),
44 50: (8.1250, 0.3250, 0.0000, 0.033027, 0.320, 0.004033, 2118.5, 0.00171, 1 / 1250),
45 51: (15.0000, 0.6000, 0.0071, 0.036173, 0.294, 0.004315, 570.4, 0.00171, 1 / 1250),
46 52: (49.2200, 1.6075, 0.0047, 0.036173, 0.304, 0.001716, 4025.6, 0.00171, 1 / 1250),
47 53: (69.4565, 1.1625, 0.0028, 0.031651, 0.336, 0.002700, 9819.5, 0.00171, 1 / 1250),
48}
49data = {i: dict(zip(params, raw_data[i])) for i in raw_data}
50
51# Set the ring we are analyzing
52ring = 15
53max_failure_probability = data[ring]["max_Pf"]
54
55
56# Compute the investment cost to increase the dike height
57def exponential_investment_cost(
58 u, # increase in dike height
59 h0, # original height of the dike
60 c, # constant from Table 1
61 b, # constant from Table 1
62 lam,
63): # constant from Table 1
64 if u == 0:
65 return 0
66 else:
67 return (c + b * u) * math.exp(lam * (h0 + u))
68
69
70def eijgenraam_model(
71 X1,
72 X2,
73 X3,
74 X4,
75 X5,
76 X6,
77 T1,
78 T2,
79 T3,
80 T4,
81 T5,
82 T6,
83 T=300,
84 P0=data[ring]["P0"],
85 V0=data[ring]["V0"],
86 alpha=data[ring]["alpha"],
87 delta=0.04,
88 eta=data[ring]["eta"],
89 gamma=0.035,
90 rho=0.015,
91 zeta=data[ring]["zeta"],
92 c=data[ring]["c"],
93 b=data[ring]["b"],
94 lam=data[ring]["lam"],
95):
96 """Python implementation of the Eijgenraam model
97
98 Params
99 ------
100 Xs : list
101 list of dike heightenings
102 Ts : list
103 time of dike heightenings
104 T : int, optional
105 planning horizon
106 P0 : <>, optional
107 constant from Table 1
108 V0 : <>, optional
109 constant from Table 1
110 alpha : <>, optional
111 constant from Table 1
112 delta : float, optional
113 discount rate, mentioned in Section 2.2
114 eta : <>, optional
115 constant from Table 1
116 gamma : float, optional
117 paper says this is taken from government report, but no indication
118 of actual value
119 rho : float, optional
120 risk-free rate, mentioned in Section 2.2
121 zeta : <>, optional
122 constant from Table 1
123 c : <>, optional
124 constant from Table 1
125 b : <>, optional
126 constant from Table 1
127 lam : <>, optional
128 constant from Table 1
129
130 """
131 Ts = [T1, T2, T3, T4, T5, T6]
132 Xs = [X1, X2, X3, X4, X5, X6]
133
134 Ts = [int(Ts[i] + sum(Ts[:i])) for i in range(len(Ts)) if Ts[i] + sum(Ts[:i]) < T]
135 Xs = Xs[: len(Ts)]
136
137 if len(Ts) == 0:
138 Ts = [0]
139 Xs = [0]
140
141 if Ts[0] > 0:
142 Ts.insert(0, 0)
143 Xs.insert(0, 0)
144
145 S0 = P0 * V0
146 beta = alpha * eta + gamma - rho
147 theta = alpha - zeta
148
149 # calculate investment
150 investment = 0
151
152 for i in range(len(Xs)):
153 step_cost = exponential_investment_cost(Xs[i], 0 if i == 0 else sum(Xs[:i]), c, b, lam)
154 step_discount = math.exp(-delta * Ts[i])
155 investment += step_cost * step_discount
156
157 # calculate expected losses
158 losses = 0
159
160 for i in range(len(Xs) - 1):
161 losses += math.exp(-theta * sum(Xs[: (i + 1)])) * (
162 math.exp((beta - delta) * Ts[i + 1]) - math.exp((beta - delta) * Ts[i])
163 )
164
165 if Ts[-1] < T:
166 losses += math.exp(-theta * sum(Xs)) * (
167 math.exp((beta - delta) * T) - math.exp((beta - delta) * Ts[-1])
168 )
169
170 losses = losses * S0 / (beta - delta)
171
172 # salvage term
173 losses += S0 * math.exp(beta * T) * math.exp(-theta * sum(Xs)) * math.exp(-delta * T) / delta
174
175 def find_height(t):
176 if t < Ts[0]:
177 return 0
178 elif t > Ts[-1]:
179 return sum(Xs)
180 else:
181 return sum(Xs[: bisect.bisect_right(Ts, t)])
182
183 failure_probability = [
184 P0 * np.exp(alpha * eta * t) * np.exp(-alpha * find_height(t)) for t in range(T + 1)
185 ]
186 total_failure = 1 - functools.reduce(operator.mul, [1 - p for p in failure_probability], 1)
187 mean_failure = sum(failure_probability) / (T + 1)
188 max_failure = max(failure_probability)
189
190 return (investment, losses, investment + losses, total_failure, mean_failure, max_failure)
191
192
193if __name__ == "__main__":
194 model = Model("eijgenraam", eijgenraam_model)
195
196 model.responses = [
197 ScalarOutcome("TotalInvestment", ScalarOutcome.INFO),
198 ScalarOutcome("TotalLoss", ScalarOutcome.INFO),
199 ScalarOutcome("TotalCost", ScalarOutcome.MINIMIZE),
200 ScalarOutcome("TotalFailureProb", ScalarOutcome.INFO),
201 ScalarOutcome("AvgFailureProb", ScalarOutcome.MINIMIZE),
202 ScalarOutcome("MaxFailureProb", ScalarOutcome.MINIMIZE),
203 ]
204
205 # Set uncertainties
206 model.uncertainties = [
207 RealParameter.from_dist("P0", sp.stats.lognorm(scale=0.00137, s=0.25)),
208 # @UndefinedVariable
209 RealParameter.from_dist("alpha", sp.stats.norm(loc=0.0502, scale=0.01)),
210 # @UndefinedVariable
211 RealParameter.from_dist("eta", sp.stats.lognorm(scale=0.76, s=0.1)),
212 ] # @UndefinedVariable
213
214 # having a list like parameter were values are automagically wrappen
215 # into a list can be quite useful.....
216 model.levers = [RealParameter(f"X{i}", 0, 500) for i in range(1, 7)] + [
217 RealParameter(f"T{i}", 0, 300) for i in range(1, 7)
218 ]
219
220 ema_logging.log_to_stderr(ema_logging.INFO)
221
222 with MultiprocessingEvaluator(model, n_processes=-1) as evaluator:
223 results = evaluator.perform_experiments(1000, 4)