generated from CEVE-421-521/lab-06
-
Notifications
You must be signed in to change notification settings - Fork 0
/
template.qmd
330 lines (273 loc) · 11.6 KB
/
template.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
---
title: "Lab 6: Policy Search"
author: "Lucia Roemro-Alston (lmr12)"
jupyter: julia-1.10
date: 2024-03-01
week: 7
categories: [Lab]
format:
html: default
# YOU DO NOT NEED BOTH PDF AND DOCX.
# COMMENT OR DELETE THE ONE YOU DON'T WANT TO USE.
pdf:
documentclass: article
fontsize: 11pt
geometry:
- margin=1in
number-sections: true
code-line-numbers: true
#docx:
# toc: true
# fig-format: png
# number-sections: true
# code-line-numbers: true
date-format: "ddd., MMM. D"
execute:
cache: true
freeze: auto
bibliography: references.bib
---
```{julia}
using Revise
using HouseElevation
using CSV
using DataFrames
using DataFramesMeta
using Distributions
using LaTeXStrings
using Metaheuristics
using Plots
using Random
using Unitful
Plots.default(; margin=5Plots.mm)
```
```{julia}
function objective_function(a::AbstractFloat)
return true # PLACEHOLDER
end
```
## Objective Function
We want to keep the same objective function that we have been using in order to find the net present value taking in the cost of heightening the house and the discounted expected costs of future flood damage. We also consider and define the states of the world over which we will analyze.
The building I am studying is Fisherman's Wharf at 2200 Harborside Drive Galveston, TX.
I got the information for the building area from the cvent website for event space. The value of the property results from searches on Zillow for average building prices in the area. Finally, the depth-damage curve which I chose is a result of location, building type, and what I am looking to analyze for damages.
```{julia}
slr_scenarios = let
df = CSV.read("data/slr_oddo.csv", DataFrame)
[Oddo17SLR(a, b, c, tstar, cstar) for (a, b, c, tstar, cstar) in eachrow(df)]
end
house = let
haz_fl_dept = CSV.read("data/haz_fl_dept.csv", DataFrame) # read in the file
desc = "Cafeteria Restaurant, structure"
row = @rsubset(haz_fl_dept, :Description == desc)[1, :] # select the row I want
area = 4004u"ft^2"
height_above_gauge = 4*u"ft"
House(row; area=area, height_above_gauge=height_above_gauge, value_usd=400_000)
end
p = ModelParams(; house=house, years=2024:2100)
function draw_surge_distribution()
μ = rand(Normal(5, 1))
σ = rand(Exponential(1.5))
ξ = rand(Normal(0.1, 0.05))
return GeneralizedExtremeValue(μ, σ, ξ)
end
function draw_discount_rate()
return 0.0
end
N_SOW = 100_000
sows = [
SOW(rand(slr_scenarios), draw_surge_distribution(), draw_discount_rate()) for
_ in 1:N_SOW
] # for 100,000 SOWs
N_SOW_opt = 10
sows_opt = collect(Iterators.take(sows, 10))
```
## Validation
We can validate our answer by comparing it to a simple Monte Carlo estimate using 25,000 samples and timing how long this takes us.
```{julia}
a = Action(3.0u"ft")
sow = first(sows)
_ = run_sim(a, sow, p)
@time run_sim(a, sow, p)
```
```{julia}
_ = HouseElevation.run_sim_old(a, sow, p)
@time HouseElevation.run_sim_old(a, sow, p)
```
## Optomization Problem Example
We are looking to solve a non-linear, non-convex optomization problem using the Metaheuristics.jl package.
The following plots the function we are looking to minimize for D=2.
```{julia}
f(x) = 10length(x) + sum(x .^ 2 - 10cos.(2π * x))
let
# Generate a grid of points for the surface plot
x = range(-5; stop=5, length=1000)
y = range(-5; stop=5, length=1000)
z = [f([i, j]) for i in x, j in y]
# Create the surface plot
surface(
x, y, z; xlabel="x1", ylabel="x2", zlabel="f(x)", title=L"Minimize $f(x)$ for $D=2$"
)
end
```
Now, if we want to minimize with D=10, we have to define bounds that constrians the search space to the decision variables.
```{julia}
D = 10
bounds = boxconstraints(; lb=-5ones(D), ub=5ones(D))
```
This can be optomized using the optomize function.
```{julia}
result = optimize(f, bounds)
```
The minimum of the objective function:
```{julia}
minimum(result)
```
The value of the decision variable that achieves the minimum:
```{julia}
minimizer(result)
```
It can be helpful to set constraints and specifications to the optomizer such as time limits, especially for problems that may take very long to run.
```{julia}
options = Options(; time_limit=10.0)
```
We will use the ECA algorithm which is suggested as default.
```{julia}
algorithm = ECA(; options=options)
```
We will set a random seed to make our results more reproducable and to understand the sensitivity of our results with respect to the random seed.
```{julia}
Random.seed!(918)
result = optimize(f, bounds, algorithm)
```
And check for a different iteration:
```{julia}
Random.seed!(952)
result = optimize(f, bounds, algorithm)
```
## Varying D
We can play around with some of the conditions of the optomization above in order to gain some knowledge on the effects of certian values on the optomization. For example, we can vary D to develop our understanding of this parameter.
```{julia}
function optomization()
trials = 1:2:19
optomization_res = []
for D in trials
bounds = boxconstraints(; lb=-5ones(D), ub=5ones(D))
# Optomization:
result = optimize(f, bounds)
# Minimum:
minimum(result)
# Value of the decision variable that achieves the minimum:
minimizer(result)
# Specifications:
options = Options(; time_limit=10.0)
algorithm = ECA(; options=options)
# Random Seed:
Random.seed!(918)
result2 = optimize(f, bounds, algorithm)
push!(optomization_res, result2)
end
return optomization_res
end
x =optomization()
```
From this analysis with D as the varying parameter, one general trend is an increase in the amount of iterations as the value of D increaes. There are some deviations to this trend, but for the most part this correlation is strong. Another interesting observation is the increase in the number of decision variables that result in the minimum value as the value for D increases. This means that there are more options for an ideal scenario that cuases minimization. The number of function calls is another result that increases along with the value of D in thsi analysis. Interestingly, the value of the minimum has no consistent trend, and the highest values of it occur around D values between 13 and 15.
## Optimization
In order to optimize our problem we must define an objective function that will include the objective and the SOWs over which we will analyse the problem.
```{julia}
# Sea Level Rise Scenarios:
slr_scenarios = let
df = CSV.read("data/slr_oddo.csv", DataFrame)
[Oddo17SLR(a, b, c, tstar, cstar) for (a, b, c, tstar, cstar) in eachrow(df)]
end
# Flood Depth and Home-Specific Information:
house = let
haz_fl_dept = CSV.read("data/haz_fl_dept.csv", DataFrame) # read in the file
desc = "Cafeteria Restaurant, structure"
row = @rsubset(haz_fl_dept, :Description == desc)[1, :] # select the row I want
area = 4004u"ft^2"
height_above_gauge = 4*u"ft"
House(row; area=area, height_above_gauge=height_above_gauge, value_usd=400_000)
end
p = ModelParams(; house=house, years=2024:2100)
# Surge Distribiton:
function draw_surge_distribution()
μ = rand(Normal(5, 1))
σ = rand(Exponential(1.5))
ξ = rand(Normal(0.1, 0.05))
return GeneralizedExtremeValue(μ, σ, ξ)
end
# Discount Rate:
function draw_discount_rate()
return 0.0
end
# Combination of information above to create different states of the world:
N_SOW = 100_000
sows = [
SOW(rand(slr_scenarios), draw_surge_distribution(), draw_discount_rate()) for
_ in 1:N_SOW
] # for 100,000 SOWs
# taking the first 10 SOWs drom those generated above
N_SOW_opt = 10
sows_opt = collect(Iterators.take(sows, 10))
```
This is the objective function which will take in the SOW information from the cell above and use it to analyze a sat of actions, or home elevations over all of these SOWs. This will provide information on how a given home elvation will act across "all" scenarios. We then take the negative sum of these values to have a comparison over which to consider different actions.
```{julia}
function objective(a::Vector{Float64})
action = Action(a[1]*u"ft")
# finding the npv of each N_SOW_opt
npv_each = [run_sim(action,s,p) for s in sows_opt]
# the negative sum of the expected values
min = -1.0 * mean(npv_each)
return min
end
```
Testing this function for different heights:
∆h = 1.0ft
```{julia}
a = [1.0]
npv_min = objective(a)
```
∆h = 7.0ft
```{julia}
a = [7.0]
npv_min = objective(a)
```
∆h = 13.0ft
```{julia}
a = [13.0]
npv_min = objective(a)
```
We can run an optimization of the objective function to find what elevation it recommends.
```{julia}
bounds = boxconstraints(; lb=0.0, ub=14.0)
# Optomization:
result = optimize(objective, bounds)
# Minimum:
minimum(result)
# Value of the decision variable that achieves the minimum:
minimizer(result)
# Specifications:
options = Options(; time_limit=10.0)
algorithm = ECA(; options=options)
# Random Seed:
Random.seed!(2024)
result2 = optimize(objective, bounds, algorithm)
```
Based on the optimization of our objective function, we are finding that the minimizing action is elevating our building to ∆h = 14.0ft.
```{julia}
NPVs = []
possible_elevations = [0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,9.5,10.0,10.5,11.0,11.5,12.0,12.5,13.0,13.5,14.0]
a =[[getindex(possible_elevations,i)] for i in 1:length(possible_elevations)]
npv_a = objective.(a)
push!(NPVs, npv_a)
return NPVs
scatter(possible_elevations, NPVs, legend=false)
```
As is clear in this graph, the recommended elevation, 14.0ft, is the minimum, validating our results.
## Reflection
How are we framing this problem? What are the decision variables, the objective function, and the states of the world over which we optimize?
In order to optimize this problem, we are setting it up as a means to figure out which one of our actions are policies has the best results when the simulation is run on them and parameters are varied. The decision variable is how high to elevate our home, and the objective function gives us data about our home and the flooding and surge distributions surrounding my building. The states of the world over which we optimize are different situations composed of real life considerations that make up the environment for a trial of our solutions. The state of the world can include future economy (discount rate), surge situations, and sea level rise.
Diggning deeper, we are averaging the objective function computed over a finite number of states of the world. This assumes that they are all drawn from a distribution representing the “true” distribution of states of the world. Is this a good assumption?
This is a good assumption for our current modeling capabilities and with limited computer efficiency. If we had the ability to do so, we would ideally run infinate simulations so that the function converges onto one value. Because this is not possible, we can get a good unferstanding of what would happen if this were the case and develop further knowledge on a topic that is difficult to predict.
What’s not being considered in this analysis that might be important?
One thing that is not being considered in this analysis that can have a significant imppact is fwater contamination and quality in our home's area. This can affect the way that the flooding intereacts with the home and as a result, the damages produced, making home elevation more or less important. There are also more personal and individual influences to whether or not someone will elevate their home like work promotions or job losses that can impact a persons ability to raise.