GitHub

Walkthrough

This walkthrough builds the Bolivian hydrothermal planning problem, constructs the feasibility-guaranteeing policy that strict TS-DDR needs, and takes a few gradient steps — on CPU, in a few minutes, on a horizon short enough to watch.

It is deliberately not the full-scale run, which trains for hours on a GPU over 126 stages. The point here is that every piece of it is visible in a script you can execute, and that the pieces are the ones the full run depends on. Where this simplifies, it says so.

The mathematics is in The long-term hydrothermal planning problem; how to train and evaluate at full scale is Reproducing the study.

Setup

Run from examples/HydroPowerModels, whose Project.toml carries everything used below.

using DecisionRules
using JuMP, DiffOpt, Ipopt
using Flux
using Random
using Statistics

1. The system

The case is the Bolivian national grid operated over weekly stages: a transmission network with thermal units, and a set of reservoirs linked into cascades, driven by historical inflow scenarios. Three files describe it — the network, the hydro topology, and the inflows.

CASE_DIR = joinpath(HYDRO_DIR, "bolivia")
"/home/runner/work/DecisionRules.jl/DecisionRules.jl/examples/HydroPowerModels/bolivia"

2. Build the stage problems

build_hydropowermodels reads one serialized stage model per stage — produced from the case by export_subproblem_mof.jl through HydroPowerModels — and re-parameterizes it: the incoming reservoir state becomes a parameter, the inflow becomes a parameter, and in strict mode the outgoing state is bound to a target parameter by a hard equality.

A short horizon keeps this runnable; the published run uses 126.

include(joinpath(HYDRO_DIR, "load_hydropowermodels.jl"))
include(joinpath(HYDRO_DIR, "hydro_reachable_policy.jl"))

NUM_STAGES = 3

diff_optimizer = () -> DiffOpt.diff_optimizer(
    optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0),
)

subproblems, state_params_in, state_params_out, uncertainty_samples,
    initial_volumes, max_volume, hydro_meta = build_hydropowermodels(
    CASE_DIR, "ACPPowerModel.mof.json";
    num_stages = NUM_STAGES, optimizer = diff_optimizer, strict = true,
)

(stages = length(subproblems), reservoirs = hydro_meta.nHyd,
 inflow_scenarios = length(uncertainty_samples[1]), K = hydro_meta.K)
(stages = 3, reservoirs = 11, inflow_scenarios = 15, K = 0.6048)

In strict mode there are no slack variables on the target and no penalty term. The dual of reservoir_out == target is therefore the clean marginal value of water, $\partial Q_t / \partial \hat v_t$, with no penalty noise mixed into it. That is the entire reason strict mode exists.

It is only well posed if every target the policy emits is reachable in one stage — otherwise the equality makes the stage infeasible. Hence the policy.

3. The reachable policy

hydro_reachable_policy maps an unconstrained network output into the one-stage reachable interval of each reservoir,

\[\hat v_r \;=\; \ell_r(v, w) + \bigl(u_r(v, w) - \ell_r(v, w)\bigr)\,\sigma(z_r),\]

and then applies a cascade clamp, so a downstream target can never assume more water than the upstream unit actually released.

Two details carry the published result:

  • the activation is a stretched sigmoid onto [0, 1 - 1e-3], not a plain one. A plain sigmoid cannot attain the ends of the interval, and the good policy on this case puts a substantial share of its targets exactly at a feasibility extreme;
  • the bounds $\ell_r, u_r$ depend on the incoming state, and that dependence is differentiated. Treating it as constant still trains and still lowers the loss, while descending a direction 48 degrees off the true gradient.
Random.seed!(42)
policy = hydro_reachable_policy(hydro_meta, [128, 128]; combiner_layers = [256, 256])

The encoder is an LSTM over the inflow sequence only; the reservoir state enters through the state-conditioned head, not through the recurrence.

4. One rollout

A rollout threads the realized state: the policy sees the state it actually reached, emits a target, the stage problem projects that target onto the feasible set, and the realized outgoing state becomes the next stage's input.

Written as a function rather than a bare loop: at top level (and inside a documentation @example block) a for introduces its own scope, so a loop-carried total_cost += ... would fail with UndefVarError. This is a recurring Julia trap in exactly this kind of script.

function rollout(policy, scenario)
    state = Float64.(initial_volumes)
    total = 0.0
    for t in 1:NUM_STAGES
        for (j, param) in enumerate(state_params_in[t])
            set_parameter_value(param, state[j])
        end
        for (param, val) in scenario[t]
            set_parameter_value(param, val)
        end

        w = Float32.([val for (_, val) in scenario[t]])
        target = policy(vcat(w, Float32.(state)))
        for j in 1:hydro_meta.nHyd
            set_parameter_value(state_params_out[t][j][1], Float64(target[j]))
        end

        optimize!(subproblems[t])
        @assert termination_status(subproblems[t]) in
                (MOI.LOCALLY_SOLVED, MOI.OPTIMAL, MOI.ALMOST_LOCALLY_SOLVED)
        total += objective_value(subproblems[t])

        for j in 1:hydro_meta.nHyd
            state[j] = value(state_params_out[t][j][2])
        end
    end
    return total, state
end

scenario = [uncertainty_samples[t][1] for t in 1:NUM_STAGES]
total_cost, final_state = rollout(policy, scenario)
total_cost
5013.770547840853

Every stage solved, from an untrained policy: the reachable map did its job and no emitted target was infeasible. That property is what makes hard target equalities usable at all.

5. The marginal value of water

The multiplier of the target equality is what TS-DDR differentiates, and it comes straight off the solved stage — no extra machinery:

lambda = [DecisionRules.pdual(state_params_out[NUM_STAGES][j][1])
          for j in 1:hydro_meta.nHyd]
round.(lambda; digits = 3)
11-element Vector{Float64}:
   0.0
   0.0
 209.346
   0.0
   0.0
   0.0
   0.0
  86.045
   0.0
   0.0
  57.79

A negative entry means holding one more unit in that reservoir lowers future cost — water has value there. The spread across reservoirs is the locational content that the production factors and the cascade topology create: a cubic metre of water is not a fungible commodity.

6. A few training steps

train_multistage assembles the loop: sample a scenario, roll out, collect the multipliers, backpropagate through the policy, step the optimizer. Three iterations here; the published run took 890 across three restarted stages.

uncertainty_samples is passed DIRECTLY: DecisionRules.sample has an overload for it that draws one inflow scenario per stage, which is exactly the per-stage joint sampling this case needs.

DecisionRules.train_multistage(
    policy, initial_volumes, subproblems,
    state_params_in, state_params_out, uncertainty_samples;
    num_train_per_batch = 2,
    num_batches = 3,
    optimizer = Flux.Adam(1e-3),
)
Main.HydroReachablePolicy{Flux.Chain{Tuple{Flux.LSTMCell{Matrix{Float32}, Matrix{Float32}, Vector{Float32}}, Flux.LSTMCell{Matrix{Float32}, Matrix{Float32}, Vector{Float32}}}}, Flux.Chain{Tuple{Flux.Dense{typeof(Main.stretchedsigmoid), Matrix{Float32}, Vector{Float32}}, Flux.Dense{typeof(Main.stretchedsigmoid), Matrix{Float32}, Vector{Float32}}, Flux.Dense{typeof(Main.stretchedsigmoid), Matrix{Float32}, Vector{Float32}}}}, Tuple{Tuple{Vector{Float32}, Vector{Float32}}, Tuple{Vector{Float32}, Vector{Float32}}}, Vector{Float32}, Nothing}(Chain(LSTMCell(11 => 128), LSTMCell(128 => 128)), Chain(Dense(139 => 256, stretchedsigmoid), Dense(256 => 256, stretchedsigmoid), Dense(256 => 11, stretchedsigmoid)), ((Float32[0.062251348, 0.19234937, 0.1234384, 0.11134807, -0.18531075, 0.2368241, 0.2202905, 0.045568105, -0.037735876, 0.18668492  …  0.027903153, -0.030833833, -0.021179313, -0.040772814, -0.014631278, -0.3843649, -0.3674907, -0.0060699694, -0.08303242, -0.04760731], Float32[0.2100807, 0.51544386, 0.32200336, 1.8758478, -0.86023045, 0.54381037, 1.5813458, 0.11377677, -0.14232308, 1.1232488  …  0.049564824, -0.15488139, -0.50934637, -0.88616645, -0.035394035, -1.4389149, -0.5914538, -0.009402217, -0.32100695, -0.47423267]), (Float32[0.078988984, -0.08059676, 0.021467313, -0.12199448, -0.18590732, 0.12625715, 0.05937218, -0.032897923, -0.04388441, 0.10670452  …  0.09982549, 0.09407575, -0.034589745, -0.03134629, -0.034946054, 0.03563281, 0.046024486, 0.030403953, 0.13255543, -0.0169413], Float32[0.1789464, -0.13874146, 0.045570377, -0.25449112, -0.3610112, 0.24651009, 0.12795624, -0.070383996, -0.085475296, 0.20946114  …  0.20102271, 0.18579787, -0.07020476, -0.05776476, -0.07909844, 0.066298306, 0.08091117, 0.059018314, 0.3270558, -0.034148563])), 0, 11, 11, Float32[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], Float32[0.1, 137.99, 0.042, 16.76, 10.35, 0.32, 9.72, 0.07, 0.04, 0.0, 2.93], Float32[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], Float32[10.2926, 10.5531, 0.7854, 3.63, 8.74, 17.01, 1.25, 7.0, 11.0, 0.84, 3.24], Float32[6.382515, 0.0, 0.0, 0.0, 0.0, 5.285952, 0.0, 0.0, 4.2336, 0.0, 0.0], nothing, 0.6048, Main.CascadeLink[Main.CascadeLink(1, 2, true, 6.382515f0), Main.CascadeLink(6, 5, false, 5.285952f0), Main.CascadeLink(9, 8, false, 4.2336f0)])

Three updates on three stages will not produce a good policy, and the number above is not meaningful on its own — it is one noisy sample of a stochastic objective, and the trap this case study keeps returning to: the training loss, the training objective and the fixed-panel evaluation are three different quantities, and only the last one selects a policy.

Where to go from here

The full-scale run is the same construction with a longer horizon, a training schedule of several phases, and a GPU. How each method arrives at a price for water is Valuing water: two approaches; how to train, evaluate and compare at full scale is Reproducing the study, and the example READMEs of DecisionRules.jl and DecisionRulesExa.jl carry the commands.