3. Running a Simulation

Physics, callbacks, and your first simulation

Now we’ll implement the physics and simulation logic, then run our first simulation.

Setup

First, we need all the type definitions from the previous section:

using SimOptDecisions
using Distributions
using Random
using CairoMakie

# Config
Base.@kwdef struct HouseElevationConfig{T<:AbstractFloat} <: AbstractConfig
    horizon::Int = 70
    gauge_height_above_ref::T = 0.0
    house_height_above_ref::T = 4.0
    house_area_ft2::T = 1500.0
    house_value::T = 200_000.0
end

# SOW
struct HouseElevationSOW{T<:AbstractFloat} <: AbstractSOW
    gev_μ::T
    gev_σ::T
    gev_ξ::T
    dd_threshold::T
    dd_saturation::T
    discount_rate::T
end

function sample_sow(rng::AbstractRNG)
    gev_μ = rand(rng, Normal(2.8, 0.3))
    gev_σ = rand(rng, truncated(Normal(1.0, 0.15); lower=0.3))
    gev_ξ = rand(rng, truncated(Normal(0.15, 0.05); lower=-0.2, upper=0.5))
    dd_threshold = rand(rng, Normal(0.0, 0.25))
    dd_saturation = rand(rng, Normal(8.0, 0.5))
    discount_rate = rand(rng, truncated(Normal(0.03, 0.015); lower=0.01, upper=0.07))
    return HouseElevationSOW(gev_μ, gev_σ, gev_ξ, dd_threshold, dd_saturation, discount_rate)
end

# State, Action, Policy
struct SeaLevelState{T<:AbstractFloat} <: AbstractState
    msl::T
end

struct ElevationAction{T<:AbstractFloat} <: AbstractAction
    elevation_ft::T
end

struct ElevationPolicy{T<:AbstractFloat} <: AbstractPolicy
    elevation_ft::T
end

function SimOptDecisions.get_action(
    policy::ElevationPolicy, state::SeaLevelState, sow::HouseElevationSOW, t::TimeStep
)
    return ElevationAction(policy.elevation_ft)
end

Domain Physics

Depth-Damage Function

The depth-damage function converts flood depth (at the house) to fractional damage. We use a logistic curve with a threshold and saturation point.

function depth_damage(depth::T, threshold::T, saturation::T) where {T<:AbstractFloat}
    depth <= threshold && return zero(T)
    depth >= saturation && return one(T)
    midpoint = (threshold + saturation) / 2
    steepness = T(6) / (saturation - threshold)
    return one(T) / (one(T) + exp(-steepness * (depth - midpoint)))
end
Code 
let
    depths = range(-1, 10; length=200)
    fig = Figure(; size=(600, 400))
    ax = Axis(fig[1, 1];
        xlabel="Flood depth at house (ft)",
        ylabel="Damage (fraction of house value)",
        title="Depth-Damage Function")

    damages = [depth_damage(d, 0.0, 8.0) for d in depths]
    lines!(ax, depths, damages; linewidth=2, label="Nominal (threshold=0, sat=8)")

    for (thresh, sat) in [(-0.5, 7.5), (0.5, 8.5)]
        damages_alt = [depth_damage(d, thresh, sat) for d in depths]
        lines!(ax, depths, damages_alt; linewidth=1, linestyle=:dash, alpha=0.5)
    end

    axislegend(ax; position=:rb)
    fig
end
Figure 1: Depth-damage function: fractional damage as a function of flood depth

Elevation Cost Function

The cost to elevate increases with height. We use the cost model from Doss-Gollin et al. (2022).

function elevation_cost(Δh::Real, area_ft2::Real, house_value::Real)
    Δh <= 0 && return 0.0
    Δh > 14 && error("Cannot elevate more than 14 ft")

    base_cost = 10_000 + 300 + 470 + 4_300 + 2_175 + 3_500  # = $20,745
    thresholds = [0.0, 5.0, 8.5, 12.0, 14.0]
    rates = [80.36, 82.5, 86.25, 103.75, 113.75]

    rate = rates[1]
    for i in 1:(length(thresholds) - 1)
        if Δh <= thresholds[i + 1]
            t = (Δh - thresholds[i]) / (thresholds[i + 1] - thresholds[i])
            rate = rates[i] + t * (rates[i + 1] - rates[i])
            break
        end
        rate = rates[i + 1]
    end

    cost_dollars = base_cost + area_ft2 * rate
    return cost_dollars / house_value
end

elevation_cost(Δh::Real) = elevation_cost(Δh, 1500.0, 200_000.0)
Code 
let
    heights = 0:14
    costs = [elevation_cost(h) for h in heights]

    fig = Figure(; size=(600, 400))
    ax = Axis(fig[1, 1];
        xlabel="Elevation height (ft)",
        ylabel="Cost (fraction of house value)",
        title="Elevation Cost Function (1500 ft², \$200k house)")
    scatterlines!(ax, heights, costs; markersize=10)
    fig
end
Figure 2: Construction cost as a function of elevation height

The Five Callbacks

SimOptDecisions uses five callbacks to run simulations:

Callback Purpose Returns
initialize Create initial state State
time_axis Define time points Iterable
get_action Map policy + state to action Action
run_timestep Execute one step (new_state, step_record)
finalize Aggregate into outcome Anything

Initialize

Create the initial state:

function SimOptDecisions.initialize(
    ::HouseElevationConfig, ::HouseElevationSOW, ::AbstractRNG
)
    return SeaLevelState(0.0)  # MSL at reference datum
end

Time Axis

Define the simulation time points:

function SimOptDecisions.time_axis(config::HouseElevationConfig, ::HouseElevationSOW)
    return 1:(config.horizon)
end

Run Timestep

Execute one year: sample a storm surge, compute damage:

function SimOptDecisions.run_timestep(
    state::SeaLevelState,
    action::ElevationAction,
    sow::HouseElevationSOW,
    config::HouseElevationConfig,
    t::TimeStep,
    rng::AbstractRNG,
)
    # Construction cost at t=1 only
    construction_cost = if SimOptDecisions.Utils.is_first(t)
        elevation_cost(action.elevation_ft, config.house_area_ft2, config.house_value)
    else
        0.0
    end

    # Sample annual maximum surge
    surge_dist = GeneralizedExtremeValue(sow.gev_μ, sow.gev_σ, sow.gev_ξ)
    surge_at_gauge = rand(rng, surge_dist)

    # Convert to water level relative to reference datum
    water_level = config.gauge_height_above_ref + surge_at_gauge + state.msl

    # House floor level with elevation
    house_floor_level = config.house_height_above_ref + action.elevation_ft

    # Flood depth at house
    flood_depth_at_house = water_level - house_floor_level

    # Compute damage
    damage_fraction = depth_damage(flood_depth_at_house, sow.dd_threshold, sow.dd_saturation)

    step_record = (construction_cost=construction_cost, damage_fraction=damage_fraction)
    return (state, step_record)
end

Finalize

Aggregate annual damages into NPV:

function SimOptDecisions.finalize(
    final_state::SeaLevelState,
    step_records::Vector,
    config::HouseElevationConfig,
    sow::HouseElevationSOW,
)
    construction_cost = step_records[1].construction_cost

    npv_damages = sum(
        step_records[t].damage_fraction * SimOptDecisions.Utils.discount_factor(sow.discount_rate, t)
        for t in eachindex(step_records)
    )

    return (
        construction_cost=construction_cost,
        npv_damages=npv_damages,
        annual_damages=[o.damage_fraction for o in step_records],
    )
end

Running Your First Simulation

Now we can run a simulation! Let’s try one policy with one SOW:

# Create instances
config = HouseElevationConfig()
rng = Random.Xoshiro(42)
sow = sample_sow(rng)
policy = ElevationPolicy(5.0)  # Elevate 5 feet

# Run simulation
outcome = simulate(config, sow, policy)

# Inspect results
println("Construction cost: $(round(outcome.construction_cost, digits=3)) (fraction of house value)")
println("NPV damages: $(round(outcome.npv_damages, digits=3)) (fraction of house value)")
println("Years with damage: $(count(d -> d > 0, outcome.annual_damages))")
Construction cost: 0.722 (fraction of house value)
NPV damages: 0.036 (fraction of house value)
Years with damage: 1

Visualizing Annual Damages

Let’s see when damage occurred over the 70-year horizon:

Code 
let
    fig = Figure(; size=(700, 400))
    ax = Axis(fig[1, 1];
        xlabel="Year",
        ylabel="Damage (fraction of house value)",
        title="Annual Damages (5 ft elevation, one SOW)")

    barplot!(ax, 1:length(outcome.annual_damages), outcome.annual_damages; color=:steelblue)
    fig
end
Figure 3: Annual flood damages over the 70-year simulation

Comparing Elevation Levels

Let’s compare different elevation choices on the same SOW:

Code 
let
    elevations = 0:14
    total_costs = Float64[]

    for elev in elevations
        policy = ElevationPolicy(Float64(elev))
        # Use same RNG seed for each to make comparison fair
        outcome = simulate(config, sow, policy, Random.Xoshiro(42))
        push!(total_costs, outcome.construction_cost + outcome.npv_damages)
    end

    fig = Figure(; size=(700, 400))
    ax = Axis(fig[1, 1];
        xlabel="Elevation (ft)",
        ylabel="Total cost (fraction of house value)",
        title="Total Cost by Elevation (same SOW)")

    barplot!(ax, collect(elevations), total_costs; color=:steelblue)

    # Mark minimum
    min_idx = argmin(total_costs)
    scatter!(ax, [elevations[min_idx]], [total_costs[min_idx]];
        color=:red, markersize=20, marker=:star5)

    fig
end
Figure 4: Total cost (construction + NPV damages) for different elevations

Summary

We’ve now:

  1. Implemented domain physics (depth-damage, elevation cost)
  2. Written the five callbacks (initialize, time_axis, get_action, run_timestep, finalize)
  3. Run our first simulation with simulate(config, sow, policy)
  4. Inspected the outcome and visualized results

But this was just one SOW (one possible future). The “optimal” elevation depends on which future actually occurs.

Next Steps

In the next section, we’ll evaluate a single policy across many SOWs to understand how it performs under uncertainty.