5. Exploratory Modeling

Mapping policy performance across scenarios

In the DMDU literature, exploratory modeling means systematically running a model across a wide range of plausible assumptions to understand where and why a policy succeeds or fails (Bankes 1993). Rather than searching for the single “optimal” policy, exploratory modeling asks: under what conditions does a given policy lead to unacceptable outcomes?

This section uses explore() to run all policy×scenario combinations and then analyzes the results to build a scenario map—a visualization of which regions of the uncertainty space are favorable or unfavorable for a given strategy.

Setup

using SimOptDecisions
using Distributions
using Random
using CairoMakie
using Statistics

# Physics
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)

# Types
Base.@kwdef struct HouseConfig{T<:AbstractFloat} <: AbstractConfig
    horizon::Int = 70
    house_height::T = 4.0
    house_area_ft2::T = 1500.0
    house_value::T = 200_000.0
end

@scenariodef HouseScenario begin
    @timeseries water_levels
    @continuous dd_threshold
    @continuous dd_saturation
    @continuous discount_rate
end

function sample_scenario(rng::AbstractRNG, horizon::Int)
    # Draw GEV parameters from prior distributions
    μ = rand(rng, Normal(2.8, 0.3))
    σ = rand(rng, truncated(Normal(1.0, 0.15); lower=0.3))
    ξ = rand(rng, truncated(Normal(0.15, 0.05); lower=-0.2, upper=0.5))

    # Generate water level time series from GEV
    surge_dist = GeneralizedExtremeValue(μ, σ, ξ)
    water_levels = [rand(rng, surge_dist) for _ in 1:horizon]

    HouseScenario(
        water_levels = water_levels,
        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)),
    )
end

struct HouseState{T<:AbstractFloat} <: AbstractState
    elevation_ft::T
end

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

@policydef ElevationPolicy begin
    @continuous elevation_ft 0.0 14.0
end

# Functions
SimOptDecisions.initialize(::HouseConfig, ::HouseScenario, ::AbstractRNG) = HouseState(0.0)
SimOptDecisions.time_axis(config::HouseConfig, ::HouseScenario) = 1:(config.horizon)
SimOptDecisions.get_action(policy::ElevationPolicy, ::HouseState, ::TimeStep, ::HouseScenario) = ElevationAction(value(policy.elevation_ft))

function SimOptDecisions.run_timestep(
    state::HouseState, action::ElevationAction, t::TimeStep,
    config::HouseConfig, scenario::HouseScenario, rng::AbstractRNG
)
    if is_first(t)
        new_elevation = state.elevation_ft + action.elevation_ft
        construction_cost = elevation_cost(action.elevation_ft, config.house_area_ft2, config.house_value)
    else
        new_elevation = state.elevation_ft
        construction_cost = 0.0
    end

    water_level = scenario.water_levels[t]
    floor_level = config.house_height + new_elevation
    flood_depth = water_level - floor_level
    damage = depth_damage(flood_depth, value(scenario.dd_threshold), value(scenario.dd_saturation))

    return (HouseState(new_elevation), (construction_cost=construction_cost, damage_fraction=damage))
end

# Outcome type: @outcomedef wraps fields in parameter types for explore()
@outcomedef HouseOutcome begin
    @continuous construction_cost
    @continuous npv_damages
    @continuous total_cost
end

function SimOptDecisions.compute_outcome(
    step_records::Vector, ::HouseConfig, scenario::HouseScenario
)
    construction_cost = step_records[1].construction_cost
    npv_damages = sum(
        step_records[t].damage_fraction * discount_factor(value(scenario.discount_rate), t)
        for t in eachindex(step_records)
    )
    return HouseOutcome(
        construction_cost = construction_cost,
        npv_damages = npv_damages,
        total_cost = construction_cost + npv_damages,
    )
end

The explore() Function

The explore() function runs simulate() for every (policy, scenario) combination and returns a YAXArray Dataset—a labeled multi-dimensional array.

NoteOutcome Types for explore()

In the previous section, compute_outcome returned a plain NamedTuple. For explore(), outcome fields must be parameter types so the framework can build labeled arrays. We use @outcomedef to define the outcome type—the macro auto-wraps plain values into ContinuousParameter, just like @scenariodef and @policydef do for their fields.

config = HouseConfig()
rng = Random.Xoshiro(42)

scenarios = [sample_scenario(rng, config.horizon) for _ in 1:200]
policies = [ElevationPolicy(elevation_ft=Float64(e)) for e in 0:2:14]

result = explore(config, scenarios, policies; progress=false)

println("Result type: $(typeof(result))")
println("Cubes: $(keys(result.cubes))")
println("Shape: $(size(result[:total_cost]))")  # (n_policies, n_scenarios)

Result type: YAXArrays.Datasets.Dataset
Cubes: [Symbol("scenario_water_levels[32]"), Symbol("scenario_water_levels[51]"), Symbol("scenario_water_levels[62]"), Symbol("scenario_water_levels[30]"), Symbol("scenario_water_levels[40]"), Symbol("scenario_water_levels[55]"), Symbol("scenario_water_levels[10]"), Symbol("scenario_water_levels[19]"), Symbol("scenario_water_levels[12]"), Symbol("scenario_water_levels[65]"), Symbol("scenario_water_levels[29]"), Symbol("scenario_water_levels[13]"), :scenario_dd_threshold, Symbol("scenario_water_levels[27]"), Symbol("scenario_water_levels[45]"), Symbol("scenario_water_levels[35]"), :construction_cost, Symbol("scenario_water_levels[24]"), Symbol("scenario_water_levels[15]"), Symbol("scenario_water_levels[39]"), Symbol("scenario_water_levels[18]"), Symbol("scenario_water_levels[52]"), Symbol("scenario_water_levels[1]"), Symbol("scenario_water_levels[33]"), Symbol("scenario_water_levels[42]"), Symbol("scenario_water_levels[61]"), :npv_damages, Symbol("scenario_water_levels[69]"), Symbol("scenario_water_levels[56]"), Symbol("scenario_water_levels[7]"), Symbol("scenario_water_levels[58]"), Symbol("scenario_water_levels[22]"), Symbol("scenario_water_levels[48]"), Symbol("scenario_water_levels[50]"), Symbol("scenario_water_levels[2]"), Symbol("scenario_water_levels[67]"), Symbol("scenario_water_levels[44]"), Symbol("scenario_water_levels[47]"), Symbol("scenario_water_levels[41]"), Symbol("scenario_water_levels[37]"), Symbol("scenario_water_levels[59]"), Symbol("scenario_water_levels[9]"), Symbol("scenario_water_levels[54]"), :scenario_dd_saturation, Symbol("scenario_water_levels[36]"), :policy_elevation_ft, Symbol("scenario_water_levels[31]"), Symbol("scenario_water_levels[43]"), Symbol("scenario_water_levels[16]"), Symbol("scenario_water_levels[4]"), Symbol("scenario_water_levels[14]"), Symbol("scenario_water_levels[70]"), Symbol("scenario_water_levels[38]"), Symbol("scenario_water_levels[57]"), Symbol("scenario_water_levels[63]"), Symbol("scenario_water_levels[11]"), Symbol("scenario_water_levels[6]"), Symbol("scenario_water_levels[23]"), Symbol("scenario_water_levels[49]"), Symbol("scenario_water_levels[3]"), Symbol("scenario_water_levels[28]"), Symbol("scenario_water_levels[17]"), Symbol("scenario_water_levels[20]"), Symbol("scenario_water_levels[5]"), Symbol("scenario_water_levels[21]"), Symbol("scenario_water_levels[25]"), Symbol("scenario_water_levels[68]"), Symbol("scenario_water_levels[46]"), :total_cost, Symbol("scenario_water_levels[60]"), :scenario_discount_rate, Symbol("scenario_water_levels[34]"), Symbol("scenario_water_levels[26]"), Symbol("scenario_water_levels[53]"), Symbol("scenario_water_levels[64]"), Symbol("scenario_water_levels[8]"), Symbol("scenario_water_levels[66]")]
Shape: (8, 200)

Accessing Results

# Get specific value
first_val = first(result[:total_cost][policy=1, scenario=1])
println("Policy 1, Scenario 1: $(round(first_val, digits=3))")

# Get all outcomes for one policy
policy_3_costs = result[:total_cost].data[3, :]
println("Policy 3 mean: $(round(mean(policy_3_costs), digits=3))")

Policy 1, Scenario 1: 0.712
Policy 3 mean: 0.916

Visualizing Results

let
    elevations = 0:2:14
    n_policies = length(elevations)

    fig = Figure(; size=(800, 500))
    ax = Axis(fig[1, 1];
        xlabel="Elevation (ft)",
        ylabel="Total cost (fraction of house value)",
        title="Cost Distribution by Elevation ($(length(scenarios)) scenarios)",
        xticks=(1:n_policies, string.(elevations)))

    for (i, elev) in enumerate(elevations)
        costs = result[:total_cost].data[i, :]
        boxplot!(ax, fill(i, length(costs)), costs; width=0.6)
    end

    fig
end

Figure 1: Distribution of total costs by elevation

Executors

By default, explore() runs simulations sequentially. For large experiments (many policies × many scenarios), executors control how simulations are parallelized.

Executor	How it runs	Best for
SequentialExecutor()	Single-threaded loop	Debugging, small experiments
ThreadedExecutor()	Threads.@threads across Julia threads	Most production use
DistributedExecutor()	asyncmap across Distributed workers	Multi-node clusters

All three accept the same keyword arguments:

# Sequential (default) — single-threaded, good for debugging
result = explore(config, scenarios, policies;
    executor=SequentialExecutor(; crn=true, seed=1234))

# Multi-threaded — uses all available Julia threads
result = explore(config, scenarios, policies;
    executor=ThreadedExecutor(; crn=true, seed=1234))

# Distributed — across workers added with addprocs()
result = explore(config, scenarios, policies;
    executor=DistributedExecutor(; crn=true, seed=1234))

ThreadedExecutor also accepts ntasks to control the number of parallel tasks (defaults to Threads.nthreads()).

Common Random Numbers (CRN)

All executors support Common Random Numbers (CRN), enabled by default. CRN assigns each scenario index a deterministic random seed, so the rng passed to initialize and run_timestep produces the same stream regardless of which policy is being evaluated or which thread runs the simulation.

This matters because it reduces variance when comparing policies: any difference in outcomes is due to the policy, not to random noise. Without CRN, two policies evaluated on the “same” scenario might see different random draws, making comparisons noisier.

# CRN enabled (default) — same random stream per scenario across all policies
result = explore(config, scenarios, policies;
    executor=ThreadedExecutor(; crn=true, seed=42))

# CRN disabled — each simulation gets an independent random stream
result = explore(config, scenarios, policies;
    executor=ThreadedExecutor(; crn=false))

Note

CRN only affects the rng argument passed to callbacks. In this tutorial, our scenarios contain pre-generated water level time series, so run_timestep doesn’t use rng at all—the random draws happen during scenario generation. CRN becomes important when your model samples random values during simulation (e.g., stochastic demand, random failures).

Storage Backends

For very large experiments, stream results to disk:

# Save to Zarr format
result = explore(config, scenarios, policies;
    backend=ZarrBackend("results.zarr"))

# Load later
result = load_zarr_results("results.zarr")

Scenario Discovery

The core question in exploratory modeling is not “which policy is best?” but “under what conditions does a given policy fail?”

This reframing is powerful: instead of searching for an optimal decision, we map out the vulnerabilities of each candidate decision. This is sometimes called scenario discovery in the DMDU literature.

When does a moderate elevation fail?

Let’s take a concrete policy—elevating 4 feet—and ask: across our 200 scenarios, when does this policy lead to high total costs?

We’ll map the results against two key uncertain parameters: the discount rate (which controls how much we weight future damages) and the average storm severity (summarized as the mean of each scenario’s water level time series).

let
    # Extract scenario parameters
    discount_rates = [value(s.discount_rate) for s in scenarios]
    mean_water_levels = [mean(s.water_levels.values) for s in scenarios]

    # Get total costs for the 4ft elevation policy (index 3: 0,2,4,...)
    policy_idx = 3  # 4ft elevation (0:2:14 → index 3)
    costs = result[:total_cost].data[policy_idx, :]

    # Threshold for "unacceptable" outcome
    cost_threshold = 0.5

    fig = Figure(; size=(700, 550))
    ax = Axis(fig[1, 1];
        xlabel="Discount rate",
        ylabel="Mean annual water level (ft)",
        title="Scenario Map: 4ft Elevation Policy")

    sc = scatter!(ax, discount_rates, mean_water_levels;
        color=costs, colormap=Reverse(:RdYlGn), markersize=10, alpha=0.8)

    Colorbar(fig[1, 2], sc; label="Total cost (fraction of house value)")

    # Mark the threshold
    hlines!(ax, [mean(mean_water_levels)]; color=(:gray, 0.3), linestyle=:dot)
    vlines!(ax, [mean(discount_rates)]; color=(:gray, 0.3), linestyle=:dot)

    fig
end

Figure 2: Scenario map for a 4ft elevation policy. Each point is one scenario, colored by total cost. The dashed line marks a cost threshold of 0.5 (half the house value). Scenarios above and to the left tend to produce worse outcomes: low discount rates weight future damages heavily, and high water levels cause more flooding.

The scenario map reveals the structure of vulnerability:

• Low discount rates (left side) make future damages count more, pushing total costs up
• Higher mean water levels (top) mean more frequent flooding, increasing damage regardless of discounting
• The combination of low discounting and high water levels produces the worst outcomes

This is more informative than asking “what is the optimal elevation for each scenario?”—it tells us where in the uncertainty space our chosen policy is at risk.

Comparing vulnerability across policies

We can repeat this analysis for multiple elevation levels to see how vulnerability shifts:

let
    discount_rates = [value(s.discount_rate) for s in scenarios]
    median_dr = median(discount_rates)
    elevations = 0:2:14
    cost_threshold = 0.5

    low_dr = discount_rates .< median_dr
    high_dr = .!low_dr

    fig = Figure(; size=(700, 450))
    ax = Axis(fig[1, 1];
        xlabel="Elevation (ft)",
        ylabel="Fraction of scenarios with cost > $(cost_threshold)",
        title="Policy Vulnerability by Discount Rate Region",
        xticks=(1:length(elevations), string.(elevations)))

    for (label, mask, color) in [
        ("Low discount rate", low_dr, :firebrick),
        ("High discount rate", high_dr, :steelblue)
    ]
        fracs = Float64[]
        for (i, _) in enumerate(elevations)
            costs_subset = result[:total_cost].data[i, mask]
            push!(fracs, mean(costs_subset .> cost_threshold))
        end
        barplot!(ax, 1:length(elevations), fracs;
            dodge=label == "Low discount rate" ? fill(1, length(elevations)) : fill(2, length(elevations)),
            color=color, label=label)
    end

    axislegend(ax; position=:rt)
    fig
end

Figure 3: Fraction of scenarios exceeding a cost threshold of 0.5, by elevation and scenario region. Low discount rate scenarios (below median) are more sensitive to elevation choice.

TipKey Insight: Scenario Discovery

The scenario map doesn’t tell us which policy is “best”—it tells us under what conditions each policy is vulnerable. A decision-maker who believes discount rates will be low should choose a higher elevation. One who expects high discount rates can accept lower elevation. This is the essence of exploratory modeling: understanding the relationship between assumptions and outcomes, rather than collapsing everything into a single recommendation.

Iterative Analysis

Exploratory modeling is not a one-shot process. A typical DMDU workflow iterates:

1. Explore: Run policies across scenarios and build scenario maps
2. Discover: Identify which uncertainties drive poor outcomes
3. Refine: Generate new scenarios that stress-test the vulnerable regions, or design new policies that address the discovered vulnerabilities
4. Repeat: Explore the refined set

For example, the scenario maps above might lead you to:

• Add more scenarios with low discount rates (to better resolve that region)
• Design an adaptive policy that performs differently depending on observed flood frequency
• Include additional uncertain parameters (e.g., sea-level rise) that were initially omitted

The explore() output provides the raw data for this iterative process. External tools (classification trees, PRIM, etc.) can be applied to the result matrix to automate the scenario discovery step.

Summary

The explore() function:

• explore(config, scenarios, policies) → YAXArray Dataset
• Results indexed by :policy and :scenario dimensions
• Executors for parallelization (Sequential, Threaded, Distributed)
• CRN for variance reduction when comparing policies
• Storage backends for large experiments (InMemory, Zarr)

Exploratory modeling concepts:

• A scenario map shows how outcomes vary across the uncertainty space for a given policy
• Scenario discovery identifies the conditions under which a policy leads to unacceptable outcomes
• The goal is to understand where a policy is vulnerable, not just whether it’s optimal on average

Exploratory modeling gives us a rich picture of how our decisions interact with uncertainty. In the next section, we’ll use optimize() to automatically search for Pareto-optimal policies—but the scenario maps from this section will help us interpret those results.

References

Bankes, Steve. 1993. “Exploratory Modeling for Policy Analysis.” Operations Research 41 (3): 435–49. https://doi.org/10/c7rgcr.