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
Discount factor on graph edges and padding stages in SDDP.jl (NEWAVE-style)
2

Hi @MagnoOcarlos welcome to the forum :slight_smile:

I’m away for the weekend so I’ll reply properly on Monday, but your approach looks good.

Instead of setting the inflow value in simulation you need to set the noise term.

3

@odow Thank you!

4

Hi @MagnoOcarlos,

Is this the correct approach for conditioning on lagged states in SDDP.jl

Yes. I need to write a tutorial that explains this in the documentation.

Are there numerical concerns with this pattern

Potentially, but not really. Assuming that your stochastic model is a reasonable it probably has some nice properties like mean reversion, etc.

Is there a better way to handle this, perhaps using SDDP.MarkovianGraph

In theory, yes, you could use a Markovian graph of some sort. But the way you’re doing it is probably better.

1 Like
5

Thank you a lot, @odow ! I am facing some challenges creating this complex model, and have some other questions about the discount factor in the cost-to-go function. I will open another post about it if you don’t mind :sweat_smile:

1 Like
6

Hi @MagnoOcarlos,

There is no easy way of passing the exact 2000 scenarios that you simulated in Python directly to SDDP while modeling the inflows as a PAR(p).

When you model the inflows as a PAR(p) your uncertainty in SDDP is not the inflow itself anymore but the noises that you add in the PAR(p) equations. This means that at every stage you will have a different set of noises that you can sample from. When you sample different noises at each stage there is no guarantee that you will recover the exact scenario that you simulated in Python.

A second possible problem is that when you keep sampling noises if you reach a sequence of particularly negative noises the resulting inflow can be less than 0 meaning that you are creating water in your model. Some implementation add penalties for negative inflows, other change the noises in a way such that the accumulation cannot reach a negative inflow and there is novel research on how to deal with this pattern.

1 Like
7

@guilhermebodin , thank you for your detailed answer! And yes, I am dealing with the possibility of negative noise by the lognormal transformation with the shift variable Gamma to avoid the negative inflow. Even with that, as you said, there is no guarantee to recover the exact scenario of ENA and could result in negative values. I saw some possibilities to define a maximum gamma shift, which takes the highest value of gamma per stage, but I think this one could generate unreal scenarios. So I think that the approach of the penalties would be more interesting.

8

Hi @odow, following your confirmation that the approach was correct, here is my full implementation. I would like to check if it makes sense end-to-end.

State variables

@variable(sp, 0 <= ear[r in codigos_ree], SDDP.State,
          initial_value = ear_initial_value[r][1])

# Lagged standardized inflows — PAR(p) memory
@variable(sp, eaf_lag[r in codigos_ree, τ in 1:p_max], SDDP.State,
          initial_value = init_std_last[m, r])

The stored energy ear is the main hydro state. The lagged standardized inflows eaf_lag[r, τ] carry the PAR(p) memory (up to p_max lags) for each reservoir equivalence area r.

Inflow variables

@variable(sp, ω_ena[r in codigos_ree])    # full inflow (MWmonth, original scale)
@variable(sp, ω_std[r in codigos_ree])    # inflow in standardized space
@variable(sp, ω_noise[r in codigos_ree])  # noise term (fixed by parameterize)

Stochastic parameterization — noise residuals

Instead of parameterizing the full inflow, I pre-compute the PAR(p) residuals \xi_t = EAF_t - \sum_k \phi_{m,k} \cdot EAF_{t-k} in Python and pass only the noise term:

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

PAR(p) reconstruction and denormalization

The full inflow is reconstructed inside the model in two steps.

First, the standardized inflow is rebuilt from the lagged states and the noise:

\tilde{\omega}_r = \sum_{k=1}^{p_{m,r}} \phi_{m,k,r} \cdot \widetilde{EAF}_{r,k}^{\text{in}} + \xi_r 
@constraint(sp, par_reconstruction[r in codigos_ree],
    ω_std[r] == sum(
        get_phi(m, k, r) * eaf_lag[r, k].in
        for k in 1:get_pm(m, r)
    ) + ω_noise[r]
)

Then it is denormalized back to the original scale:

\omega_r = \mu_{m,r} + \sigma_{m,r} \cdot \tilde{\omega}_r 
@constraint(sp, denorm[r in codigos_ree],
    ω_ena[r] == ω_std[r] * std_m_par[m, r] + mean_m_par[m, r]
)

Questions

  1. Does this structure look correct? In particular, is it valid to have ω_std and ω_ena as free (non-state) variables reconstructed each stage via constraints, or is there a reason to make them states?

  2. The lagged shift pattern (updating eaf_lag[r, τ].out) is done elsewhere in the subproblem. The key point is that ω_std[r] depends on eaf_lag[r, k].in, so the Benders cuts should correctly capture the dual value of past inflows — is that reasoning right?

Thanks again!

9

Seems good. I’m writing a tutorial for the SDDP docs that should help you.

Edit: How to implement a PAR model · SDDP.jl

2 Likes
10

Thanks!

1 Like