Conditioning stochastic inflows on lagged state variables in SDDP.jl (PAR(p) model)

Specific Domains Optimization (Mathematical)
1

Conditioning stochastic inflows on lagged state variables in SDDP.jl (PAR(p) model)

Hi all,

I’m implementing a hydrothermal dispatch model (similar to the Brazilian NEWAVE) using SDDP.jl. The inflow process follows a Periodic Autoregressive model PAR(p):

EAF_t(k) = φ₁·EAF_{t-1}(k) + φ₂·EAF_{t-2}(k) + ... + φₚ·EAF_{t-p}(k) + ξ_t

where ξ_t is the innovation (noise) term.

In the NEWAVE formulation, the Benders cuts explicitly depend on past inflows as state variables (not just stored energy):

α_{t+1} - Σ π^EA · EA_{t+1} - Σ π^EAF1 · EAF_t - ... - Σ π^EAFp · EAF_{t-p+1} ≥ δ

I currently carry the past inflows as state variables using the shift pattern from the docs
[Access a first-stage decision in a future stage]

@variable(sp, eaf_lag[r in REEs, τ in 1:max_lag], SDDP.State, initial_value = eaf_init[r])
@constraint(sp, [r in REEs], eaf_lag[r, 1].out == ω_ena[r])
@constraint(sp, [r in REEs, τ in 2:max_lag], eaf_lag[r, τ].out == eaf_lag[r, τ-1].in)

However, I pre-generate 2000 full EAF scenarios in Python and pass them via SDDP.parameterize with fixed samples — so the realized inflow ω_ena[r] does NOT depend on eaf_lag[r,τ].in. As a result, the duals π^EAF are essentially zero and the lag states are dead weight (66 extra dimensions that slow convergence without improving the policy).

My question: What is the recommended way to make the inflow realization depend on the lagged state variables inside SDDP.jl?

My current thinking is:

  1. Pre-compute the PAR(p) residuals ξ_t = EAF_t - Σ φ·EAF_{t-τ} in Python
  2. Pass ξ_t (the noise) as the parameterized random variable instead of the full EAF_t
  3. Reconstruct the full inflow inside the model as a constraint:
@variable(sp, ω_noise[r in REEs])  # parameterized noise
@variable(sp, ω_ena[r in REEs])    # reconstructed full inflow

@constraint(sp, par_reconstruction[r in REEs],
    ω_ena[r] == sum(φ[m,r,τ] * eaf_lag[r,τ].in for τ in 1:p_order[m,r]) + ω_noise[r]
)

SDDP.parameterize(sp, noise_samples[t], probs[t]) do ξ
    for r in REEs
        fix(ω_noise[r], ξ[idx[r]]; force=true)
    end
end

This way the water balance and demand constraints would see ω_ena[r] which depends on eaf_lag, giving non-zero π^EAF duals in the Benders cuts.

Questions:

  • Is this the correct approach for conditioning on lagged states in SDDP.jl?
  • Are there numerical concerns with this pattern (e.g., the PAR coefficients φ creating large coefficient ranges in the cuts)?
  • Is there a better way to handle this, perhaps using SDDP.MarkovianGraph?

Thanks!

1 Like