6. Policy Search

Multi-objective optimization to find the best trade-offs

Instead of manually searching over elevation levels, we can use optimization to automatically find the best policies. When objectives conflict (cost vs. damage), we get a Pareto front of optimal trade-offs.

Setup

using SimOptDecisions
using Distributions
using Random
using CairoMakie
using Statistics

# === All type definitions (self-contained) ===

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

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

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

# Required for optimization
ElevationPolicy(params::AbstractVector) = ElevationPolicy(params[1])
SimOptDecisions.params(p::ElevationPolicy) = [p.elevation_ft]
SimOptDecisions.param_bounds(::Type{<:ElevationPolicy}) = [(0.0, 14.0)]

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

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

function elevation_cost(Δh::Real, area_ft2::Real, house_value::Real)
    Δh <= 0 && return 0.0
    base_cost = 20_745.0
    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
    return (base_cost + area_ft2 * rate) / house_value
end
elevation_cost(Δh::Real) = elevation_cost(Δh, 1500.0, 200_000.0)

function SimOptDecisions.initialize(::HouseElevationConfig, ::HouseElevationSOW, ::AbstractRNG)
    return SeaLevelState(0.0)
end

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

function SimOptDecisions.run_timestep(
    state::SeaLevelState, action::ElevationAction, sow::HouseElevationSOW,
    config::HouseElevationConfig, t::TimeStep, rng::AbstractRNG
)
    construction_cost = SimOptDecisions.Utils.is_first(t) ?
        elevation_cost(action.elevation_ft, config.house_area_ft2, config.house_value) : 0.0

    surge_dist = GeneralizedExtremeValue(sow.gev_μ, sow.gev_σ, sow.gev_ξ)
    surge_at_gauge = rand(rng, surge_dist)
    water_level = config.gauge_height_above_ref + surge_at_gauge + state.msl
    house_floor_level = config.house_height_above_ref + action.elevation_ft
    flood_depth_at_house = water_level - house_floor_level
    damage_fraction = depth_damage(flood_depth_at_house, sow.dd_threshold, sow.dd_saturation)

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

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)
end

Outcomes vs Metrics

Before optimization, understand two key concepts:

  • Outcome: Result of a single simulation (one policy × one SOW)
  • Metric: Aggregated statistic across many SOWs for a given policy

The optimizer evaluates policies by computing metrics across the SOW ensemble.

"""
Compute metrics from outcomes across many SOWs.
Returns a NamedTuple that the optimizer can reference.
"""
function calculate_metrics(outcomes)
    damages = [o.npv_damages for o in outcomes]
    construction = outcomes[1].construction_cost

    return (
        construction_cost=construction,
        expected_damage=mean(damages),
        cvar_95=mean(sort(damages)[(end - max(1, length(damages) ÷ 20)):end]),
        worst_case=maximum(damages),
    )
end

Setting Up the Problem

config = HouseElevationConfig()
n_sows = 1000
rng = Random.Xoshiro(42)
sows = [sample_sow(rng) for _ in 1:n_sows]

println("Config: $(config.horizon)-year horizon, \$$(Int(config.house_value)) house")
println("SOWs: $(n_sows) sampled futures")
Config: 70-year horizon, $200000 house
SOWs: 1000 sampled futures

Multi-Objective Problem

We trade off construction cost against expected damage—conflicting objectives:

import Metaheuristics  # loads the optimization extension

# Optimize over [3, 14] ft (continuous range)
# We'll add 0 ft separately to handle the discrete gap
prob = OptimizationProblem(
    config,
    sows,
    ElevationPolicy,
    calculate_metrics,
    [minimize(:construction_cost), minimize(:expected_damage)];
    bounds=[(3.0, 14.0)],  # custom bounds for smooth optimization
)
OptimizationProblem{HouseElevationConfig{Float64}, HouseElevationSOW{Float64}, ElevationPolicy, typeof(calculate_metrics)}(HouseElevationConfig{Float64}(70, 0.0, 4.0, 1500.0, 200000.0), HouseElevationSOW{Float64}[HouseElevationSOW{Float64}(3.0365066804812875, 0.8680212106068401, 0.10631032588951869, -0.18331374110872317, 7.64890240062116, 0.02891127197931195), HouseElevationSOW{Float64}(2.5261229840980217, 1.094743124667513, 0.22193416378557068, -0.08854246084261151, 8.398063459639017, 0.012369433001736197), HouseElevationSOW{Float64}(2.62198148820797, 1.0644409547457927, 0.19761642315606365, 0.032284802915086036, 7.852613053219432, 0.024385973471723183), HouseElevationSOW{Float64}(3.1508535220741116, 1.0385272744163991, 0.16004625443194756, 0.2834748872756562, 8.144894197218777, 0.03327601966065405), HouseElevationSOW{Float64}(2.617265231083102, 1.0973421041907703, 0.2167934731792261, -0.16648013540557244, 7.497201004812148, 0.04212737631233758), HouseElevationSOW{Float64}(2.6857762671706715, 0.9690080328278211, 0.13446280634581328, -0.010118350062088726, 8.052385395894461, 0.015517642694704159), HouseElevationSOW{Float64}(3.0896982248669187, 1.1966526033600249, 0.170208770037958, 0.2994899174531351, 7.385455643459239, 0.05190114622109096), HouseElevationSOW{Float64}(2.7887796907768525, 0.8926094029958801, 0.10204855303721891, 0.23871132138863554, 8.189284201920376, 0.026139587209633768), HouseElevationSOW{Float64}(2.4488340058761127, 1.1396098524494624, 0.1763496994434102, -0.023277577703625953, 8.077624832158886, 0.02579488418677791), HouseElevationSOW{Float64}(2.5471607034554085, 0.9695640248226942, 0.1356734277226387, -0.05794016963802915, 8.527877033976434, 0.012784669743705207)  …  HouseElevationSOW{Float64}(3.154943776317375, 0.9942331861823888, 0.16098948268650382, 0.04762699225156271, 7.5302773642838, 0.04503504140041249), HouseElevationSOW{Float64}(2.9135749666117077, 1.0177890877980103, 0.18947726286777905, -0.4207742718598112, 8.213295395679053, 0.03310545691725598), HouseElevationSOW{Float64}(2.5449337058968684, 1.2630902096579346, 0.18558769596389765, 0.08282287259428643, 7.459876304894185, 0.04729096442977641), HouseElevationSOW{Float64}(3.0371365288802483, 0.9732238479028417, 0.1305851227184358, 0.2790027556754343, 7.581911813476502, 0.01887313719248586), HouseElevationSOW{Float64}(2.579426994648304, 0.7087510296964865, 0.16688708885715364, -0.09774641513196625, 7.71523460368626, 0.0465985994114465), HouseElevationSOW{Float64}(2.4908412719644164, 1.112873596622874, 0.09705221448106119, -0.06456692252890037, 7.716307728287777, 0.0413284357926222), HouseElevationSOW{Float64}(2.6990457438573445, 0.9446613535979722, 0.1308037201578582, -0.10060726588121222, 7.773259437218494, 0.04192163666589968), HouseElevationSOW{Float64}(2.8240845764361193, 0.7311811457344246, 0.20354460010299635, 0.12778582436684355, 7.422199429957107, 0.029279655192239042), HouseElevationSOW{Float64}(2.972506213918957, 0.9519997609075798, 0.14419626426405974, 0.004968072619444291, 9.031183591355582, 0.038376220748816914), HouseElevationSOW{Float64}(2.7309177703267533, 0.7497132417794324, 0.22088824549170816, -0.3315442613704665, 7.918125612398264, 0.02764433807571119)], ElevationPolicy, calculate_metrics, Objective[Objective(:construction_cost, Minimize), Objective(:expected_damage, Minimize)], FullBatch(), AbstractConstraint[], [(3.0, 14.0)])
TipHandling Discrete Choices

Real elevation choices are discrete: you either don’t elevate (0 ft) or elevate to at least 3 ft. Rather than encoding this discontinuity in the policy (which confuses the optimizer), we:

  1. Optimize over the continuous [3, 14] range
  2. Add the 0 ft baseline using merge_into_pareto!

This gives the optimizer a smooth landscape while including discrete alternatives.

Running Optimization

# NSGA-II is a popular multi-objective evolutionary algorithm
backend = MetaheuristicsBackend(;
    algorithm=:NSGA2,
    max_iterations=100,
    population_size=30,
    parallel=true
)
MetaheuristicsBackend(:NSGA2, 100, 30, true, Dict{Symbol, Any}())
opt_result = SimOptDecisions.optimize(prob, backend)

# Add the "no elevation" baseline to the Pareto front
merge_into_pareto!(opt_result, prob, ElevationPolicy(0.0))

println("Found $(length(opt_result.pareto_params)) Pareto-optimal solutions")
Found 31 Pareto-optimal solutions

Understanding the Pareto Front

A solution is Pareto-optimal if you can’t improve one objective without worsening another. The Pareto front shows the efficient frontier of trade-offs.

Code 
let
    fig = Figure(; size=(600, 600))
    ax = Axis(fig[1, 1];
        xlabel="Construction cost (fraction of house value)",
        ylabel="Expected NPV damages (fraction of house value)",
        title="Pareto Front: Cost vs Damage Trade-off",
        aspect=1)

    pareto_pts = collect(SimOptDecisions.pareto_front(opt_result))

    if !isempty(pareto_pts)
        construction = [obj[1] for (_, obj) in pareto_pts]
        expected = [obj[2] for (_, obj) in pareto_pts]
        elevations = [round(p[1]; digits=1) for (p, _) in pareto_pts]

        # Iso-cost lines (x + y = constant)
        all_costs = construction .+ expected
        cost_range = range(minimum(all_costs) * 0.8, maximum(all_costs) * 1.1; length=6)
        for total_cost in cost_range
            lines!(ax, [0.0, total_cost], [total_cost, 0.0];
                color=(:gray, 0.3), linewidth=0.5, linestyle=:dash)
        end

        scatter!(ax, construction, expected; markersize=12, color=:steelblue)

        # Label selected points
        labeled = Set{Float64}()
        for (i, elev) in enumerate(elevations)
            if all(abs(elev - l) > 1.5 for l in labeled)
                text!(ax, construction[i], expected[i];
                    text="$(round(Int, elev)) ft", align=(:left, :bottom), offset=(5, 5))
                push!(labeled, elev)
            end
        end
    end

    fig
end
Figure 1: Pareto front showing trade-off between construction cost and expected damage

Interpreting the Results

The diagonal lines show iso-cost contours (total cost = construction + damage). Moving along the Pareto front:

  • Lower left: Infeasible—you can’t get low cost and low damage
  • Upper left: Low construction, high damage (0 ft)
  • Lower right: High construction, low damage (14 ft)
  • The curve: Efficient trade-offs between extremes
println("Pareto-optimal solutions:")
for (params, objectives) in pareto_front(opt_result)
    elev = round(params[1], digits=1)
    cost = round(objectives[1], digits=3)
    damage = round(objectives[2], digits=3)
    total = round(objectives[1] + objectives[2], digits=3)
    println("  $(elev) ft → cost=$cost, damage=$damage (total=$total)")
end
Pareto-optimal solutions:
  3.0 ft → cost=0.716, damage=0.261 (total=0.977)
  14.0 ft → cost=0.957, damage=0.007 (total=0.964)
  4.0 ft → cost=0.719, damage=0.165 (total=0.884)
  9.5 ft → cost=0.79, damage=0.027 (total=0.816)
  3.1 ft → cost=0.716, damage=0.254 (total=0.97)
  3.3 ft → cost=0.717, damage=0.229 (total=0.946)
  13.8 ft → cost=0.949, damage=0.008 (total=0.957)
  7.6 ft → cost=0.744, damage=0.046 (total=0.79)
  13.2 ft → cost=0.927, damage=0.009 (total=0.936)
  7.0 ft → cost=0.738, damage=0.057 (total=0.795)
  12.8 ft → cost=0.913, damage=0.01 (total=0.923)
  12.5 ft → cost=0.902, damage=0.011 (total=0.914)
  12.2 ft → cost=0.891, damage=0.012 (total=0.903)
  12.1 ft → cost=0.886, damage=0.013 (total=0.899)
  3.4 ft → cost=0.717, damage=0.212 (total=0.93)
  3.7 ft → cost=0.718, damage=0.192 (total=0.911)
  10.1 ft → cost=0.81, damage=0.023 (total=0.833)
  6.2 ft → cost=0.732, damage=0.072 (total=0.804)
  11.6 ft → cost=0.868, damage=0.015 (total=0.883)
  5.7 ft → cost=0.728, damage=0.084 (total=0.813)
  5.2 ft → cost=0.724, damage=0.104 (total=0.828)
  4.8 ft → cost=0.722, damage=0.119 (total=0.841)
  11.0 ft → cost=0.846, damage=0.017 (total=0.864)
  10.7 ft → cost=0.835, damage=0.019 (total=0.854)
  4.5 ft → cost=0.721, damage=0.135 (total=0.856)
  4.2 ft → cost=0.72, damage=0.151 (total=0.871)
  8.6 ft → cost=0.755, damage=0.035 (total=0.79)
  9.7 ft → cost=0.797, damage=0.025 (total=0.822)
  13.0 ft → cost=0.918, damage=0.01 (total=0.928)
  8.8 ft → cost=0.762, damage=0.033 (total=0.795)
  0.0 ft → cost=0.0, damage=1.444 (total=1.444)

Selecting from the Pareto Front

The optimizer finds the efficient frontier, but choosing a point is a decision, not an optimization. Different stakeholders might choose differently:

Stakeholder Preference Likely Choice
Risk-averse homeowner Minimize worst case Higher elevation
Budget-constrained Minimize upfront cost Lower elevation
Expected-value optimizer Minimize total expected cost Middle of curve 
# Find the solution with minimum total expected cost
pareto_pts = collect(pareto_front(opt_result))
total_costs = [obj[1] + obj[2] for (_, obj) in pareto_pts]
best_idx = argmin(total_costs)
best_params, best_obj = pareto_pts[best_idx]

println("Minimum expected total cost:")
println("  Elevation: $(round(best_params[1], digits=1)) ft")
println("  Construction: $(round(best_obj[1], digits=3))")
println("  Expected damage: $(round(best_obj[2], digits=3))")
println("  Total: $(round(best_obj[1] + best_obj[2], digits=3))")
Minimum expected total cost:
  Elevation: 7.6 ft
  Construction: 0.744
  Expected damage: 0.046
  Total: 0.79

Summary

We’ve covered the complete SimOptDecisions workflow:

  1. Define the problem with types (Config, SOW, State, Action, Policy)
  2. Implement physics (depth-damage, elevation cost)
  3. Write callbacks (initialize, time_axis, get_action, run_timestep, finalize)
  4. Run simulations with simulate()
  5. Evaluate policies across many SOWs
  6. Explore systematically with grid search or explore()
  7. Optimize with multi-objective algorithms to find Pareto-optimal solutions

Key Takeaways

  1. Structure your uncertainty in the SOW type
  2. Separate physics from framework in standalone functions
  3. Validate with exploration before optimizing
  4. Use multi-objective optimization when objectives conflict
  5. The Pareto front shows trade-offs; choosing a point requires judgment

Discussion

This model makes simplifying assumptions worth noting:

  • Stationary climate: Storm distribution doesn’t change over time
  • Constant sea level: No sea-level rise (though the state structure supports it)
  • Independent damages: No cumulative effects year-to-year
  • Single decision: Elevation chosen once, not adaptively

For real applications, you’d want to relax these assumptions. The framework supports time-varying SOWs, evolving state, and more complex policies.

Next Steps