Cobre

Open infrastructure for power system computation.

Cobre is an ecosystem of Rust crates for power system analysis and optimization. The name comes from the Portuguese word for copper — the metal that conducts electricity. The project is licensed under Apache-2.0.

This site is the methodology reference

This site documents the theory, design specifications, and roadmap behind Cobre. It is written for researchers, algorithm designers, and contributors who need to understand why the system works the way it does.

What you will find here:

• Theory — The mathematical foundations: SDDP and Benders decomposition, stochastic inflow modeling with PAR(p), cut management strategies, and risk measures including CVaR.
• Specifications archive — The design-phase artifacts that governed each implementation phase: data model schemas, architectural trait contracts, HPC communication patterns, and testing contracts.
• Roadmap — Planned features and deferred variants, with rationale for what was prioritized in the minimal viable solver.
• Reference — Glossary of domain terms, notation conventions, and bibliography.

Software Book

For installation instructions, tutorials, user guides, CLI reference, and per-crate API documentation, see the Software Book.

The software book describes what the software does today. This site describes the theory and design decisions behind it.

Quick links

Go to the Software Book >>

Go to the Software Book >>

SDDP Theory

Stochastic Dual Dynamic Programming (SDDP) is the core optimization algorithm behind Cobre’s hydrothermal dispatch solver. This page provides a practical overview of how the algorithm works. For the complete formal specification, see SDDP Algorithm (spec).

What problem does SDDP solve?

Power system operators must decide how much water to release from reservoirs and how much thermal generation to dispatch at each time step over a multi-month or multi-year horizon. The challenge is uncertainty: future river inflows are unknown, and using too much water now risks expensive thermal generation later, while hoarding water risks spilling it.

SDDP finds the least-cost generation schedule that accounts for this uncertainty across all stages simultaneously. The decision at each stage depends on the current reservoir levels, recent inflow history, and a probabilistic model of future inflows.

How SDDP works: forward and backward passes

SDDP is an iterative algorithm. Each iteration has two phases:

Forward pass. The algorithm simulates the system forward in time, from stage 1 to stage $T$. At each stage it samples a random inflow scenario, solves a linear program (LP) to determine optimal generation and reservoir operations, and records the resulting reservoir levels. These recorded states are called trial points. Multiple independent trajectories are simulated in parallel.

Backward pass. Walking backward from stage $T$ to stage 1, the algorithm re-solves each stage LP at the trial points collected during the forward pass – but now for every inflow scenario in the scenario set. From each solve it extracts dual variables (shadow prices) that measure how sensitive the future cost is to changes in reservoir levels. These sensitivities are assembled into Benders cuts: linear constraints that approximate the future cost function.

Each iteration adds new cuts, progressively tightening the approximation of future costs. The forward pass uses the latest cuts to make better decisions; the backward pass generates cuts at the states actually visited.

Lower and upper bounds

SDDP produces two bounds that bracket the true optimal cost:

Lower bound. The first-stage LP objective (including the cut-based future cost approximation) provides a deterministic lower bound ${\underline{z}}^k$. Because cuts can only tighten the approximation, this bound increases monotonically across iterations.

Upper bound. The average total cost across all forward pass trajectories provides a statistical upper bound ${\bar{z}}^k$. This estimate has sampling noise, so it is typically reported with a confidence interval.

Convergence

The algorithm monitors the optimality gap between the bounds:

$${\text{gap}}^k=\frac{{\bar{z}}^k-{\underline{z}}^k}{\max (1,|{\bar{z}}^k|)}$$

When the gap falls below a threshold, or the lower bound stalls, the algorithm terminates. The resulting set of cuts defines the operational policy: at any future state, solving the LP with these cuts yields the recommended dispatch.

For detailed stopping criteria, see Convergence and Stopping Rules (spec).

Policy graphs

SDDP organizes stages into a policy graph – a directed graph where nodes are decision stages and arcs represent transitions.

Finite horizon. The standard setup is a linear chain of $T$ stages (e.g., 60 monthly stages for a 5-year study). The terminal stage has zero future cost.

Infinite horizon. For long-term planning, the last stage can loop back to the first, creating a cyclic graph. A discount factor $d<1$ on the cycle arc ensures convergence. Cuts at equivalent positions in the cycle are shared.

See SDDP Algorithm (spec) – Policy Graph Structure for details.

State variables

The state at each stage captures everything the future cost depends on. In Cobre, state variables are:

• Reservoir storage volumes ($v_h$): one per hydro plant, representing end-of-stage water levels
• Autoregressive inflow lags ($a_{h,\ell }$): past inflow values needed by the PAR(p) stochastic model to generate current-stage inflows

Together, these can reach approximately 2000 dimensions for a system with 160+ hydro plants and multi-lag inflow models.

The Markov property – future costs depend only on the current state, not on the path that led to it – is what makes cut generation mathematically valid. The AR lags are included as state variables precisely to preserve this property.

Single-cut formulation

Cobre uses the single-cut formulation by default: at each iteration, the per-scenario cut coefficients are averaged into one aggregated cut per stage. This keeps the LP small and solve times fast, at the cost of potentially needing more iterations. A multi-cut variant (one cut per scenario) is planned for future releases.

Further reading

• SDDP Algorithm (spec) – Complete formal specification with all formulas, edge cases, and design rationale
• Benders Decomposition – The decomposition technique underlying cut generation
• Forward and Backward Passes – Detailed pass mechanics
• Convergence – Bound computation and stopping criteria
• Cut Management – Cut selection, deletion, and storage strategies
• Single-Cut vs Multi-Cut – Trade-offs between aggregation strategies

See also: For implementation details and usage, see the SDDP crate documentation in the software book.

Benders Decomposition

Overview

Stochastic Dual Dynamic Programming (SDDP) solves large multistage stochastic optimization problems by decomposing them into smaller single-stage subproblems. This decomposition is rooted in Benders decomposition (also known as the L-shaped method), which splits a monolithic problem into a master problem and subproblems connected through linear inequalities called cuts.

In the hydrothermal dispatch context, a direct formulation over all stages, scenarios, and equipment would produce an LP with millions of variables and constraints. Benders decomposition breaks this into $T$ stage subproblems, each solved independently for a given incoming state and scenario realization.

The Stage Subproblem

Each stage $t$ solves an LP of the form:

$$\underset{x_t,{\theta }_t}{\min }{c}_t^{\top }{x}_t+{\theta }_t$$

$$\text{s.t.}{A}_t{x}_t=b_t-E_t{x}_{t-1},x_t\in {\mathcal{X}}_t,{\theta }_t\ge {\underline{V}}_{t+1}(x_t)$$

where $x_{t-1}$ is the incoming state (reservoir volumes, AR inflow lags) fixed from the previous stage, $c_t^{\top }{x}_t$ captures the immediate cost (thermal generation, penalties, regularization), and ${\theta }_t$ approximates the expected future cost from stages $t+1$ through $T$.

The full LP structure — objective terms, constraint families, slack variables, and variable bounds — is described in LP Formulation.

How Cuts Approximate Future Cost

The true future cost function $V_{t+1}(x_t)$ is convex but unknown. SDDP builds a piecewise-linear lower approximation ${\underline{V}}_{t+1}$ through iterative cut generation:

$${\underline{V}}_{t+1}(x_t)=\underset{i\in \mathcal{K}}{\max }\left\{{\alpha }_i+{\pi }_i^{\top }{x}_t\right\}$$

Each cut $({\alpha }_i,{\pi }_i)$ is a supporting hyperplane derived from dual variables of the backward pass. When the stage $t+1$ subproblem is solved at a trial state ${\hat{x}}_t$ under scenario $\omega$, the optimal dual multipliers ${\pi }^{\ast }$ of the state-linking constraints yield:

• Slope: $\pi =E_{t+1}^{\top }{\pi }^{\ast }$ (sensitivity of future cost to incoming state)
• Intercept: $\alpha =V_{t+1}({\hat{x}}_t,\omega )-{\pi }^{\top }{\hat{x}}_t$

Notation note: The matrix form $E_{t+1}^{\top }{\pi }^{\ast }$ above is the standard Benders decomposition notation from the stochastic programming literature. In Cobre’s formal specs, the cut coefficients are expressed directly in terms of LP dual variables: ${\pi }_h^v={\pi }_h^{fix}$ (the storage fixing constraint dual) and ${\pi }_{h,\ell }^{lag}$ (the lag-fixing constraint dual), without explicit technology matrices. Both conventions are mathematically equivalent — the matrix form is compact for general derivations, while the direct-dual form is more practical for implementation. See Cut Coefficient Derivation below.

In the hydrothermal problem, the state variables are reservoir storage volumes $v_h$ and AR inflow lags $a_{h,\ell }$, so the cut coefficients come from the storage fixing dual ${\pi }_h^{fix}$ and the lag-fixing constraint dual ${\pi }_{h,\ell }^{lag}$.

Cut Coefficient Derivation

The key insight is how dual variables map to cut coefficients. Cobre introduces an explicit incoming-state LP variable for each state dimension, fixed to the trial value via an equality constraint:

$$v_h^{in}={\hat{v}}_h,a_{h,\ell }={\hat{a}}_{h,\ell }$$

By the LP envelope theorem, the dual of each fixing constraint equals the total derivative of the optimal objective with respect to the fixed value, accounting for all downstream effects:

For storage, the dual ${\pi }_h^{fix}$ captures contributions from the water balance, FPHA hyperplane constraints, and any generic constraints that reference incoming storage — all aggregated automatically by the LP solver through the single incoming-storage variable. No manual dual combination is required.

The cut intercept ensures the cut passes exactly through the trial point:

$${\alpha }_t=Q_t({\hat{x}}_{t-1},\omega )-\sum _h{\pi }_{t,h}^v\cdot {\hat{v}}_h-\sum _{h,\ell }{\pi }_{t,h,\ell }^{lag}\cdot {\hat{a}}_{h,\ell }$$

where $Q_t({\hat{x}}_{t-1},\omega )$ is the optimal stage $t$ objective value.

Single-Cut Aggregation

In the backward pass, each trial state is evaluated against all $|\Omega |$ scenario realizations. The resulting per-scenario cuts are aggregated into one cut per trial point using probability-weighted expectation:

$${\bar{\alpha }}_{t-1}=\sum _{\omega \in \Omega }p(\omega )\cdot {\alpha }_t(\omega )$$

$${\bar{\pi }}_{t-1,h}^v=\sum _{\omega \in \Omega }p(\omega )\cdot {\pi }_{t,h}^v(\omega )$$

$${\bar{\pi }}_{t-1,h,\ell }^{lag}=\sum _{\omega \in \Omega }p(\omega )\cdot {\pi }_{t,h,\ell }^{lag}(\omega )$$

where $p(\omega )$ is the probability of scenario $\omega$ (uniform for the fixed opening tree: $p(\omega )=1/|\Omega |$). The aggregated cut $(\bar{\alpha },{\bar{\pi }}^v,{\bar{\pi }}^{lag})$ is added to stage $t-1$’s cut pool. This single-cut formulation keeps the LP compact: after $K$ iterations the LP has $K$ cuts, regardless of the number of scenarios per iteration.

The Forward-Backward Iteration

SDDP alternates between two passes:

1. Forward pass: Simulate scenario paths from stage 1 to $T$, solving each stage LP with the current cut approximation. This produces trial states ${\hat{x}}_t$ and a statistical upper bound on optimal cost.

2. Backward pass: Traverse from stage $T$ back to stage 1. At each stage, solve subproblems for all scenario realizations at each trial state, extract dual multipliers, aggregate them into new cuts, and add each cut to the stage’s cut pool. The first-stage objective (including ${\theta }_1$) provides a deterministic lower bound.

Cuts accumulate over iterations, progressively tightening ${\underline{V}}_{t+1}$ toward the true $V_{t+1}$. The algorithm converges when the gap between lower and upper bounds falls below a tolerance.

Why the LP Formulation Matters

The quality of Benders cuts depends directly on the LP formulation:

• Recourse slacks (deficit, excess) guarantee feasibility for every scenario, ensuring duals always exist. Without relatively complete recourse, the backward pass would encounter infeasible subproblems.
• Penalty magnitudes shape the value function. Correct priority ordering ensures that cuts propagate economically meaningful signals across stages.
• State variable identification determines which dual multipliers contribute to cut coefficients. Every constraint linking the current stage to the incoming state produces a component of $\pi$.

Related Topics

• LP Formulation — complete stage subproblem: objective, constraints, penalties, and Benders cut interface
• Cut Management — cut selection, deletion, and storage strategies
• SDDP Theory — the iterative algorithm structure
• Forward-Backward Pass — detailed mechanics of the forward simulation and backward cut generation passes

Further Reading

• Pereira, M. V. F. and Pinto, L. M. V. G. (1991). “Multi-stage stochastic optimization applied to energy planning.” Mathematical Programming, 52(1-3), 359-375.
• Benders, J. F. (1962). “Partitioning procedures for solving mixed-variables programming problems.” Numerische Mathematik, 4(1), 238-252.
• Birge, J. R. and Louveaux, F. V. (2011). Introduction to Stochastic Programming. Springer, 2nd edition.
• SDDP.jl documentation: https://sddp.dev/stable/
• Cut Management (spec): ../specs/math/cut-management.md

Forward and Backward Passes

SDDP iterates between two complementary phases: a forward pass that simulates the system under the current policy, and a backward pass that improves the policy by generating new Benders cuts. This page summarizes how the two passes work together. For the complete specification, see SDDP Algorithm (spec).

Forward pass

The forward pass simulates the system from stage 1 to stage $T$, producing trial points – the reservoir states that the backward pass will use to build cuts.

For each of $M$ independent trajectories:

1. Start from the known initial state $x_0$ (reservoir volumes and inflow history).
2. At each stage $t$, sample an inflow scenario ${\omega }_t$ from the opening tree and solve the stage LP with the incoming state ${\hat{x}}_{t-1}$ and the sampled scenario. The LP includes all Benders cuts accumulated so far, so the future cost variable $\theta$ reflects the best available approximation of costs from stage $t+1$ onward.
3. Record the optimal outgoing state ${\hat{x}}_t$ (end-of-stage storage volumes and updated autoregressive lags) as the trial point for stage $t$.
4. Pass ${\hat{x}}_t$ as the incoming state to stage $t+1$ and repeat.

The forward pass also records each trajectory’s total cost, which feeds the statistical upper bound estimate (see Convergence).

Backward pass

The backward pass walks stages in reverse order, from $T$ down to 1, generating new cuts that tighten the future cost approximation.

At each stage $t$, for each trial point ${\hat{x}}_{t-1}$ collected during the forward pass:

1. Solve the stage $t$ LP for every scenario $\omega$ in the opening tree, using ${\hat{x}}_{t-1}$ as the incoming state. This is called backward branching: unlike the forward pass which samples one scenario per stage, the backward pass evaluates all of them.
2. From each LP solution, extract the optimal objective value and the dual multipliers of the state-linking constraints (water balance, autoregressive lag fixing). These duals measure how sensitive the future cost is to small changes in the incoming state.
3. Compute per-scenario cut coefficients $(\alpha (\omega ),\pi (\omega ))$ from the duals and trial point.
4. Aggregate the per-scenario coefficients into a single cut via probability-weighted averaging (the single-cut formulation).
5. Add the new cut to stage $t-1$’s problem.

Each iteration adds one cut per stage per trial point, progressively tightening the piecewise-linear lower approximation of the true cost-to-go function.

How the passes work together

The two passes form a feedback loop:

• The forward pass uses the current cut set to make decisions. Better cuts lead to better trial points.
• The backward pass generates cuts at the states actually visited. More representative trial points lead to more useful cuts.

Over iterations, the cuts converge toward the true value function, and the forward pass policy converges toward the optimal dispatch schedule.

Parallelization

Both passes offer natural parallelism:

• Forward pass: The $M$ trajectories are independent. Cobre distributes them across MPI ranks and OpenMP threads.
• Backward pass: Within each stage, the branching solves across scenarios can run in parallel. However, there is a synchronization barrier at each stage boundary – all threads must finish generating cuts at stage $t$ before any thread moves to stage $t-1$, because the newly generated cuts must be available when solving stage $t-1$.

The forward pass LP solution at each stage also provides a warm-start basis for the corresponding backward pass solves, significantly reducing solver time.

For the full specification of the forward and backward pass mechanics, including scenario sampling schemes, feasibility guarantees, and the thread-trajectory affinity model, see SDDP Algorithm (spec).

Related topics

• SDDP Theory – High-level overview of the algorithm
• Benders Decomposition – The decomposition technique behind cut generation
• Convergence – How the bounds produced by these passes are used to monitor convergence
• Scenario Generation – How scenarios are sampled for forward and backward passes
• Solver Workspaces (spec) – Solver state management, basis persistence, and warm-starting across forward and backward passes

See also: For implementation details and usage, see the SDDP crate documentation in the software book.

Convergence

SDDP is an iterative algorithm, so it needs well-defined criteria for when to stop. This page explains how convergence is monitored through bounds, what the five available stopping rules are and how their formulas work, and how to obtain a deterministic upper bound independent of the risk measure.

Lower Bound

The first-stage LP objective — immediate cost plus the cut-based approximation of future cost — provides a deterministic lower bound ${\underline{z}}^k$ at iteration $k$:

$${\underline{z}}^k=c_1^{\top }{\hat{x}}_1^k+{\theta }_1^k$$

Because each iteration adds new cuts that can only tighten the approximation, the lower bound increases monotonically. It never decreases, making it a reliable progress indicator.

Statistical Upper Bound

The average total cost across all forward pass trajectories provides a statistical upper bound ${\bar{z}}^k$:

$${\bar{z}}^k=\frac{1}{M}\sum _{m=1}^M\sum _{t=1}^T{c}_t^{\top }{\hat{x}}_t^{k,m}$$

This is a Monte Carlo estimate of the true policy cost, so it carries sampling noise. It is typically reported with a confidence interval based on the sample standard deviation across trajectories.

For risk-averse problems (e.g., CVaR objectives), the statistical upper bound is not valid in the classical sense — the sample average is not a valid upper bound on a risk measure. In such cases, the deterministic upper bound described below is necessary.

Deterministic Upper Bound via Inner Approximation

For a convergence certificate valid under any risk measure, Cobre supports a deterministic upper bound through inner approximation (SIDP). The idea is to build a concave piecewise-linear overestimate of the value function from the inside, complementing the convex piecewise-linear underestimate built by cuts from the outside.

The inner approximation stores vertices — visited state-value pairs $(x^{(i)},{\bar{v}}^{(i)})$ collected during forward passes. At any query state $x$, the upper bound is interpolated via the Lipschitz condition:

$${\bar{V}}_t(x)=\underset{i\in {\mathcal{V}}_t}{\min }\left\{{\bar{v}}^{(i)}+L_t\cdot \Vert x-x^{(i)}{\Vert }_1\right\}$$

where $L_t$ is the Lipschitz constant for stage $t$, computed backward from the maximum penalty coefficient:

$$L_T=c_{max}^{penalty},L_t=L_{t+1}+c_{max}^{penalty,t}$$

This formula uses the largest penalty cost (the dominant term in the objective) as a bound on how fast the value function can change per unit change in state. The upper bound at stage 1 then gives:

$${\bar{z}}^k=c_1({\hat{x}}_1)+{\bar{V}}_2({\hat{x}}_1)$$

As iterations proceed, more vertices are collected and the inner approximation tightens. Both bounds converge to the same value, and the gap goes to zero.

Optimality Gap

The gap between the bounds measures how far the current policy may be from optimal:

$${\text{gap}}^k=\frac{{\bar{z}}^k-{\underline{z}}^k}{\max (1,|{\bar{z}}^k|)}$$

A small gap means the cut approximation is close to the true value function and the policy is near-optimal. In practice, the gap is monitored over iterations and used by the stopping rules described below.

Stopping Rules

Cobre supports five stopping criteria, evaluated at the end of each iteration. Multiple rules can be combined with OR logic (stop when any rule triggers) or AND logic (stop when all rules trigger simultaneously).

Iteration Limit

A hard cap on the number of iterations. This rule must always be present as a safety bound:

$$\text{STOP}⟺k\ge k_{max}$$

Time Limit

A wall-clock time budget. Checked at the end of each iteration:

$$\text{STOP}⟺t_{elapsed}\ge t_{max}$$

Bound Stalling

Detects when the lower bound has plateaued. The rule computes the relative improvement of ${\underline{z}}^k$ over a sliding window of $\tau$ iterations:

$${\Delta }_k=\frac{{\underline{z}}^k-{\underline{z}}^{k-\tau }}{\max (1,|{\underline{z}}^k|)}$$

$$\text{STOP}⟺|{\Delta }_k|<{\epsilon }_{stall}$$

If $|{\Delta }_k|$ falls below the tolerance, further iterations are unlikely to improve the policy significantly.

Gap Convergence

Stops when the optimality gap falls below a threshold:

$$\text{STOP}⟺{\text{gap}}^k=\frac{{\bar{z}}^k-{\underline{z}}^k}{\max (1,|{\bar{z}}^k|)}<{\epsilon }_{gap}$$

Requires an upper bound to be available (statistical or deterministic).

Simulation-Based Stopping (Recommended)

The most robust criterion combines bound stability with policy stability. Periodically (e.g., every 20 iterations), the rule:

1. Checks bound stability:

   $$\text{Bound stable}⟺\left|{\underline{z}}^k-{\underline{z}}^{k-w}\right|<{\epsilon }_{bound}\times \max (1,|{\underline{z}}^k|)$$

2. If stable, runs a batch of Monte Carlo simulations and computes the normalized distance between the per-stage cost profiles of consecutive simulation batches:

   $$d=\sqrt{\sum _t{\left(\frac{c_t^{new}-c_t^{old}}{\max (1,|c_t^{old}|)}\right)}^2}$$

3. Stops when both converge:

   $$\text{STOP}⟺\text{Bound stable}\wedge d<{\epsilon }_{policy}$$

This avoids premature termination from statistical noise and ensures the policy is genuinely stable, not just the bound.

Termination Output

When the algorithm stops, the output reports which rule triggered, the final iteration count, the lower and upper bounds, and the optimality gap. This information supports post-hoc analysis of whether the policy is sufficiently converged or whether additional iterations would be beneficial.

Related Topics

• SDDP Theory — Overview of the algorithm and its bound structure
• Forward and Backward Passes — The iteration mechanics that produce bounds
• Stopping Rules — Complete stopping rule formulas and configuration guidance
• Stopping Rules (spec) — Full specification with configuration schemas and formulas
• Upper Bound Evaluation (spec) — Deterministic upper bounds via inner approximation

Stopping Rules

Overview

SDDP is an iterative algorithm with no fixed termination point: it can always run another iteration and tighten the future cost approximation. Stopping rules provide the criteria for deciding when the current policy is good enough to stop.

Good stopping criteria balance two risks. Stopping too early leaves a suboptimal policy; stopping too late wastes computation on marginal improvements. Because the lower bound and the upper bound (statistical or deterministic) approach the true optimum from opposite sides, monitoring both gives a robust termination signal.

The Five Available Rules

Cobre provides five stopping criteria. Multiple rules can be combined: "any" mode (default) stops when the first rule triggers; "all" mode requires every rule to trigger simultaneously.

1. Iteration Limit (Mandatory)

A hard cap on the number of iterations:

$$\text{STOP}⟺k\ge k_{max}$$

where $k$ is the current iteration count and $k_{max}$ is the configured limit.

This rule must always be present. It acts as a safety bound preventing runaway computation if other rules fail to trigger. For production studies, $k_{max}$ is typically set high enough that convergence rules trigger first, but not so high that a non-converging run consumes unbounded resources.

2. Time Limit

A wall-clock budget:

$$\text{STOP}⟺t_{elapsed}\ge t_{max}$$

Checked at the end of each iteration. Useful in operational contexts where results must be available within a fixed window, regardless of convergence status.

3. Bound Stalling

Detects when the lower bound has plateaued — the outer approximation is no longer improving despite additional iterations. The rule computes the relative improvement of the deterministic lower bound ${\underline{z}}^k$ over a sliding window of $\tau$ iterations:

$${\Delta }_k=\frac{{\underline{z}}^k-{\underline{z}}^{k-\tau }}{\max (1,|{\underline{z}}^k|)}$$

$$\text{STOP}⟺|{\Delta }_k|<{\epsilon }_{stall}$$

where ${\epsilon }_{stall}$ is the relative tolerance (typically 0.01%). The window size $\tau$ controls sensitivity: a small window detects stalling faster but is more susceptible to transient plateaus; a larger window is more robust.

The intuition: if the LP approximation of future costs has not improved appreciably over the last $\tau$ iterations, the current cuts already provide a tight lower bound and further iteration is unlikely to change the dispatch policy.

4. Gap Convergence

When both a lower bound ${\underline{z}}^k$ and an upper bound ${\bar{z}}^k$ are available, their gap directly measures how far the current policy may be from optimal:

$${\text{gap}}^k=\frac{{\bar{z}}^k-{\underline{z}}^k}{\max (1,|{\bar{z}}^k|)}$$

$$\text{STOP}⟺{\text{gap}}^k<{\epsilon }_{gap}$$

The upper bound ${\bar{z}}^k$ can be the Monte Carlo average of forward pass costs (statistical, valid for risk-neutral problems) or the inner approximation via Lipschitz interpolation (deterministic, valid for all risk measures). A tight gap provides a verifiable certificate of near-optimality.

5. Simulation-Based Stopping (Recommended)

The most robust criterion combines two convergence indicators: stability of the outer approximation (bound) and stability of the policy itself (simulated cost profile). The rule executes in two steps, checked every period iterations:

Step 1 — Bound stability check:

$$\text{Bound stable}⟺\left|{\underline{z}}^k-{\underline{z}}^{k-w}\right|<{\epsilon }_{bound}\times \max (1,|{\underline{z}}^k|)$$

where $w$ is the bound window (number of past iterations to compare against).

Step 2 — Policy stability check (only if bound is stable): Run replications Monte Carlo forward simulations under the current policy. Compute the mean per-stage cost $c_t$ across replications and compare to the previous simulation batch:

$$d=\sqrt{\sum _{t=1}^T{\left(\frac{c_t^{new}-c_t^{old}}{\max (1,|c_t^{old}|)}\right)}^2}$$

Stopping condition:

$$\text{STOP}⟺\text{Bound stable}\wedge d<{\epsilon }_{policy}$$

This criterion requires both the bound and the simulated policy costs to have converged. The two-step design avoids premature termination: a temporarily stable bound might coincide with a policy that still changes significantly across iterations, and vice versa.

Why this rule is recommended: The bound monitors convergence of the mathematical approximation; the simulation monitors convergence of the economic dispatch decisions. Requiring both to stabilize ensures the solution is practically meaningful, not just numerically tight.

Combining Rules

A typical conservative setup uses an iteration limit as a hard cap with simulation-based stopping as the primary convergence criterion:

• Iteration limit: 500 (hard safety cap)
• Simulation-based stopping: every 20 iterations, 100 replications, 1% policy tolerance
• Mode: "any" (stop when either triggers)

The iteration limit ensures termination even if the simulation criterion never triggers (e.g., due to high variance). In practice, well-tuned simulation criteria typically trigger well before the iteration limit.

Termination Output

When any stopping rule triggers, the solver records which rule terminated the run, the final iteration count, the lower bound, the upper bound (if evaluated), and the optimality gap. This information is important for post-hoc analysis: a run terminated by the iteration limit before the simulation criterion converged may warrant rerunning with a higher limit.

Further Reading

• Convergence — how bounds are computed and the gap is interpreted
• SDDP Theory — the iteration structure that stopping rules monitor
• Stopping Rules (spec) — complete formal specification with configuration schemas
• Upper Bound Evaluation (spec) — deterministic upper bounds via inner approximation

Stage Subproblem LP Formulation

Overview

At each stage of the SDDP algorithm, the dispatcher solves a linear program (LP) that models the power system for that time period. This LP is the fundamental computational unit of the solver: the forward pass executes it once per stage to simulate the policy, and the backward pass executes it once per stage per scenario to generate Benders cuts.

Understanding the LP structure is important because it directly shapes the quality of the cuts and the economic signals propagated across stages. A well-structured LP ensures that cuts carry meaningful marginal water values and that every scenario is feasible.

Cost and Penalty Taxonomy

The objective function combines four types of costs. The distinction matters because each type serves a different algorithmic role.

Resource Costs

Resource costs represent actual operating expenditures with economic meaning:

Category 1: Recourse Slacks

These penalty terms guarantee that every LP is feasible regardless of the inflow scenario. Without them, the backward pass could encounter infeasible subproblems, making cut generation impossible.

Category 2: Constraint Violation Penalties

These provide soft enforcement of physical constraints that may be impossible to satisfy under extreme conditions. They must be high enough to propagate penalty signals into the future cost function and influence reservoir decisions multiple stages ahead.

Category 3: Regularization Costs

Very small costs (2–3 orders of magnitude below economic costs) that break degeneracy and guide the solver toward physically preferred solutions when the LP is otherwise indifferent:

Priority Ordering

The ordering must be maintained to ensure correct economic signals:

$$c^{fill}>c^{sv-}>c^{def}>c^{tv-},c^{ov\pm },c^{gv-}>c^{th},c^{ctr}\gg c^{spill},c^{fpha},c^{exch}$$

Filling a dead reservoir (dam safety) has the highest priority; load shedding ranks above generation costs; regularization terms are negligible compared to any economic quantity.

The Stage Subproblem LP

The complete objective minimizes total cost over all load blocks $k\in \mathcal{K}$:

$$\min \sum _{k\in \mathcal{K}}{\tau }_k[\sum _j{c}_j^{th}{g}_{j,k}+\sum _b\sum _s{c}_{b,s}^{def}{\delta }_{b,k,s}+\sum _b{c}_b^{exc}{ϵ}_{b,k}+\sum _h{c}_h^{spill}{s}_{h,k}+\text{other violation terms}]$$

$$+\sum _h\left[c_h^{sv-}{\sigma }_h^{v-}+c_h^{fill}{\sigma }_h^{fill}\right]+\theta$$

Storage violation penalties (${\sigma }_h^{v-}$, ${\sigma }_h^{fill}$) appear outside the block sum because they apply to end-of-stage storage (hm³), not per-block flow rates.

Key Constraint Families

Load Balance

For each bus $b$ and block $k$, supply must meet demand. Hydro generation, thermal dispatch, imports, and transmission flows contribute to supply; deficit and excess slack variables ensure feasibility:

$$\sum _{h\in {\mathcal{H}}_b}{g}_{h,k}+\sum _{j\in {\mathcal{T}}_b}{g}_{j,k}+\text{(flows in)}-\text{(flows out)}+\sum _s{\delta }_{b,k,s}-{ϵ}_{b,k}=D_{b,k}$$

Hydro Water Balance

End-of-stage storage equals incoming storage plus net inflows and minus all outflows, scaled by the time conversion factor $\zeta =0.0036\times {\sum }_k{\tau }_k$:

$$v_h=v_h^{in}+\zeta \left[a_h+\sum _k{w}_k\left(\sum _{i\in {\mathcal{U}}_h}(q_{i,k}+s_{i,k})-q_{h,k}-s_{h,k}-e_{h,k}\right)\right]$$

where $a_h$ is the incremental inflow from the PAR(p) model and $v_h^{in}$ is the incoming storage.

State-Fixing Constraints and the Benders Interface

A critical design choice is the introduction of explicit incoming-state LP variables $v_h^{in}$ and $a_{h,\ell }$, each fixed to the trial value via an equality constraint:

$$v_h^{in}={\hat{v}}_h,a_{h,\ell }={\hat{a}}_{h,\ell }$$

The dual of each fixing constraint is the cut coefficient for that state dimension. This “fishing constraint” technique collects all sensitivity information into a single dual value per state variable — the LP solver automatically propagates contributions from all downstream constraints (water balance, FPHA, generic constraints) through the incoming-state variable. The result appears in the stage LP as:

$$\theta \ge {\alpha }_i+\sum _h{\pi }_{i,h}^v\cdot v_h+\sum _{h,\ell }{\pi }_{i,h,\ell }^{lag}\cdot a_{h,\ell }$$

where ${\pi }_{i,h}^v$ is the dual of the storage fixing constraint and ${\pi }_{i,h,\ell }^{lag}$ is the dual of the lag fixing constraint. Each new Benders cut is derived from these duals after each backward-pass solve.

Hydro Production Models

Two production models are available during training:

Constant productivity: A linear equality linking generation to turbined flow:

$$g_{h,k}={\rho }_h\cdot q_{h,k}$$

where ${\rho }_h$ (MW per m³/s) depends on turbine efficiency and reference head. Simple and fast.

FPHA (piecewise-linear): A set of $M$ linear hyperplanes forming a concave upper bound on generation as a function of storage and flow:

$$g_{h,k}\le {\gamma }_0^m+{\gamma }_v^m\cdot v_h^{avg}+{\gamma }_q^m\cdot q_{h,k}+{\gamma }_s^m\cdot s_{h,k}$$

where $v_h^{avg}=(v_h^{in}+v_h)/2$ is average stage storage. FPHA captures head variation and the effect of spillage on tailrace level, at the cost of more constraints per hydro.

Further Reading

• Benders Decomposition — how cuts are derived from this LP structure
• Hydro Production — full derivation of constant productivity and FPHA models
• PAR(p) Autoregressive Models — the inflow model that produces $a_h$
• LP Formulation (spec) — complete formal specification with all constraint families, variable bounds, and dual variable conventions
• Cut Management (spec) — dual extraction and cut coefficient derivation

Hydro Production Models

Overview

The amount of electrical power a hydroelectric plant produces depends on two physical quantities: how much water flows through the turbines and the pressure difference (net head) driving that flow. Net head is the difference between the upstream reservoir level and the downstream tailrace level, minus hydraulic losses in the penstock.

This nonlinear relationship creates a modeling challenge for LP-based optimization: the exact production function $g=f(v,q)$ is a bilinear function of storage and flow, incompatible with linear programming. Cobre uses two linearization strategies during training (policy construction) that preserve LP tractability while capturing different levels of head-variation accuracy.

Constant Productivity Model

The simplest model assumes productivity is constant — independent of how full the reservoir is. Generation is a linear function of turbined flow:

$$g_{h,k}={\rho }_h\cdot q_{h,k}$$

where ${\rho }_h$ is the productivity coefficient (MW per m³/s), computed from the turbine efficiency and a reference net head:

$${\rho }_h=\frac{9.81\times {\eta }_h\times H_h^{ref}}{1000}$$

with ${\eta }_h$ the turbine-generator efficiency (typically 0.85–0.92) and $H_h^{ref}$ the net head at a representative operating point (typically 65% full storage).

This is an equality constraint: for each hydro $h$ and each load block $k$, generation is fully determined by turbined flow. The LP does not treat generation as a free variable — it is substituted out. The simplicity makes this model fast and numerically well-conditioned, which is valuable for long planning horizons where head variation is secondary to reservoir management.

When to use: Run-of-river plants, far-future stages, and initial algorithm development where computational speed matters more than fine-grained head variation.

FPHA: Piecewise-Linear Head Approximation

The exact hydroelectric production function is:

$$\varphi (v,q,q_{out})=\frac{9.81\times \eta \times q\times h_{net}(v,q,q_{out})}{1000}$$

where the net head captures the full hydraulic chain:

$$h_{net}(v,q,q_{out})=h_{fore}(v)-h_{tail}(q_{out})-h_{loss}(q)$$

• $h_{fore}(v)$: forebay (reservoir surface) level, a nonlinear function of storage
• $h_{tail}(q_{out})$: tailrace level, a function of total outflow $q_{out}=q+s$
• $h_{loss}(q)$: hydraulic losses in the penstock and turbines

Because $\varphi$ is a product of $q$ and $h_{net}(v,q,q_{out})$, it is nonlinear in the LP variables. FPHA (Função de Produção Hidrelétrica Aproximada) replaces this nonlinear surface with a concave piecewise-linear envelope: a set of $M$ hyperplanes that together bound generation from above everywhere in the operating region:

$$g_{h,k}\le {\gamma }_0^m+{\gamma }_v^m\cdot v_h^{avg}+{\gamma }_q^m\cdot q_{h,k}+{\gamma }_s^m\cdot s_{h,k},m=1,\dots ,M$$

where $v_h^{avg}=(v_h^{in}+v_h)/2$ is the average storage over the stage.

Physical Meaning of Coefficients

Each hyperplane coefficient has a clear hydraulic interpretation:

The negative sign on ${\gamma }_s^m$ reflects a real hydraulic effect: water released as spillage raises the downstream level, which partially cancels the head driving power generation. In systems with strong tailrace backwater effects, this is an important operational consideration.

Concave Envelope and Conservative Bias

The LP optimizer maximizes hydro generation subject to the FPHA constraints (since hydro has zero fuel cost). At the optimum, generation lies on one of the hyperplane faces. Because the envelope is a concave upper bound on the exact nonlinear surface, the FPHA model may overestimate generation compared to the physical plant.

To prevent this, a correction factor $\kappa \le 1$ scales the hyperplane intercepts:

$${\tilde{\gamma }}_0^m=\kappa \times {\gamma }_0^m$$

The default approach sets $\kappa$ as the worst-case ratio between the exact production function and the FPHA approximation across the operating grid, guaranteeing that FPHA never overestimates. For plants with modest head variation, typical values are $\kappa \approx 0.97$–$1.00$.

Impact on Water Values

Because the average storage $v_h^{avg}=(v_h^{in}+v_h)/2$ appears in the FPHA constraint, the incoming storage variable $v_h^{in}$ is linked to the generation upper bound. The dual of the incoming storage fixing constraint therefore captures both the water balance sensitivity and the FPHA generation sensitivity — all through a single LP variable. This is what makes the “fishing constraint” technique (see LP Formulation) particularly powerful: the FPHA contribution to cut coefficients is automatic.

FPHA vs. Constant Productivity: Trade-offs

In practice, many systems use FPHA for near-term stages (where head variation significantly affects near-term dispatch) and constant productivity for far-future stages (where approximation accuracy matters less than speed).

Further Reading

• LP Formulation — how production constraints fit into the assembled stage LP
• Benders Decomposition — how FPHA affects cut coefficient computation
• Hydro Production Models (spec) — complete formal specification including topology function interpolation, correction factor computation methods, and the simulation-only linearized head model

Stochastic Modeling

Why Inflows Are Uncertain

Hydrothermal dispatch optimizes the operation of interconnected hydroelectric and thermal power plants over a multi-year horizon. The dominant source of uncertainty is river inflow – the volume of water arriving at each reservoir in each time period. Inflows depend on rainfall, snowmelt, and upstream conditions that cannot be predicted with precision months or years ahead. Because water stored today displaces expensive thermal generation tomorrow (and vice versa), the dispatch decision at every stage depends critically on the distribution of future inflows, not just a single forecast.

Role of Autoregressive Models

SDDP requires a probabilistic model of inflows that captures two key properties:

1. Seasonal patterns: Inflows are strongly periodic – wet and dry seasons drive large swings in mean and variance across the year.
2. Temporal correlation: A wet month is more likely to be followed by another wet month than by a dry one. Ignoring this correlation would underestimate hydrological risk.

The Periodic Autoregressive model of order p (PAR(p)) addresses both properties. It fits separate autoregressive parameters for each season (month or week), allowing the mean, variance, and lag structure to vary across the hydrological cycle. At each stage, the model expresses the current inflow as a linear function of recent past inflows plus a random innovation term:

$$a_{h,t}={\mu }_{m(t)}+\sum _{\ell =1}^p{\psi }_{m(t),\ell }\left(a_{h,t-\ell }-{\mu }_{m(t-\ell )}\right)+{\sigma }_{m(t)}\cdot {\epsilon }_t$$

The lagged inflows $a_{h,t-1},\dots ,a_{h,t-p}$ become state variables in the SDDP subproblem, linking stages through the Bellman recursion.

Spatial Correlation

In systems with many river basins, inflows at different plants are correlated – a drought typically affects an entire region, not a single plant. Cobre handles spatial correlation by sampling the innovation vector ${\epsilon }_t$ from a multivariate distribution that preserves the observed cross-plant correlation structure. The PAR(p) model captures temporal dependence per plant, while the joint sampling of innovations captures spatial dependence across plants. See Spatial Correlation for details on the correlation modeling approach.

Integration with SDDP

During the forward pass, inflow scenarios are sampled from the PAR(p) model (or from external scenario files). During the backward pass, the opening tree of inflow realizations at each stage is generated from the same model, ensuring consistency between the policy being evaluated and the uncertainty it was trained against. When external scenarios are provided for training, a PAR(p) model is fitted to those scenarios so that the backward pass opening tree remains well-defined.

Further Reading

• PAR(p) Autoregressive Models – accessible explanation of the PAR(p) model, its parameters, and how lags become state variables
• PAR Inflow Model Specification – full mathematical specification including fitting procedure, model order selection, and validation invariants
• Scenario Generation – how inflow scenarios are produced for forward and backward passes

See also: For implementation details and usage, see the Stochastic Modeling guide in the software book.

PAR(p) Autoregressive Models

What Is a PAR(p) Model?

A Periodic Autoregressive model of order p (PAR(p)) is a time series model designed for data with strong seasonal patterns. It extends the classical autoregressive (AR) model by allowing every parameter to vary by season — the coefficients that govern January inflows are different from those that govern July inflows.

The “order p” indicates how many past time steps the model looks back. A PAR(3) model for a given month predicts the current inflow using the inflows from the previous three months. The order can differ by season: January might need only one lag while April might need four, reflecting different hydrological dynamics across the year.

The PAR(p) Equation

For hydro plant $h$ at stage $t$ falling in season $m(t)$, the PAR(p) model is:

$$a_{h,t}={\mu }_{m(t)}+\sum _{\ell =1}^p{\psi }_{m(t),\ell }\left(a_{h,t-\ell }-{\mu }_{m(t-\ell )}\right)+{\sigma }_{m(t)}\cdot {\epsilon }_t$$

In words: the inflow at stage $t$ equals the seasonal mean, plus a weighted combination of how much recent inflows deviated from their seasonal means, plus a random shock.

Parameters by Season

Each season $m$ has its own set of parameters:

The seasonal mean ${\mu }_m$ and sample standard deviation $s_m$ are estimated from historical data. The AR coefficients ${\psi }_{m,\ell }$ are fitted using the Yule-Walker equations (see below). The residual standard deviation ${\sigma }_m$ is derived at runtime from the other parameters (it is not stored independently).

How Lags Become State Variables

In the SDDP framework, decisions at each stage depend on a set of state variables that summarize everything the optimizer needs to know from the past. For the PAR(p) model, the state variables are the lagged inflows:

$$\text{State at stage }t:(v_{h,t},a_{h,t-1},a_{h,t-2},\dots ,a_{h,t-p_{\max }})$$

where $v_{h,t}$ is the reservoir volume and $a_{h,t-\ell }$ are the lagged inflows needed by the autoregressive equation. Each lag adds one state variable per hydro plant to the SDDP subproblem.

This is significant for problem size: a system with 150 hydro plants and a maximum PAR order of 6 adds up to $150\times 6=900$ state variables beyond the reservoir volumes. The LP formulation includes constraints that “shift” lagged inflows forward from one stage to the next, ensuring the autoregressive structure is respected across the Bellman recursion.

Stored vs. Computed Quantities

Cobre stores the natural outputs of the fitting process:

• Stored: seasonal means (${\mu }_m$), seasonal sample standard deviations ($s_m$), AR order ($p_m$), and AR coefficients in original units (${\psi }_{m,\ell }$)
• Computed at runtime: the residual standard deviation ${\sigma }_m$, derived from the stored quantities to guarantee consistency

This design avoids redundancy — ${\sigma }_m$ is fully determined by the other parameters and recomputing it is inexpensive.

Yule-Walker Fitting Procedure

When fitting PAR(p) parameters from historical inflow data, the AR coefficients are estimated by solving the Yule-Walker equations — a linear system that relates the autocorrelations of the data to the model coefficients. The procedure has five steps.

Step 1 — Seasonal Statistics

For each season $m$, compute the sample mean and standard deviation from historical observations $\{a_{h,t}:m(t)=m\}$:

$${\hat{\mu }}_m=\frac{1}{N_m}\sum _{t:m(t)=m}{a}_{h,t}$$

$${\hat{s}}_m=\sqrt{\frac{1}{N_m-1}\sum _{t:m(t)=m}{\left(a_{h,t}-{\hat{\mu }}_m\right)}^2}$$

where $N_m$ is the number of historical observations for season $m$.

Step 2 — Seasonal Autocorrelations

Compute the cross-seasonal autocorrelation at lag $\ell$ for season $m$. The cross-seasonal structure arises because lag $\ell$ at season $m$ reaches back to season $m-\ell$ (cyclically):

$${\hat{\gamma }}_m(\ell )=\frac{1}{N_m-1}\sum _{t:m(t)=m}\left(a_{h,t}-{\hat{\mu }}_m\right)\left(a_{h,t-\ell }-{\hat{\mu }}_{m-\ell }\right)$$

$${\hat{\rho }}_m(\ell )=\frac{{\hat{\gamma }}_m(\ell )}{{\hat{s}}_m\cdot {\hat{s}}_{m-\ell }}$$

Note that ${\hat{s}}_{m-\ell }$ is the standard deviation of season $m-\ell$, not of season $m$. This is the defining feature of a periodic (as opposed to stationary) autoregressive model.

Step 3 — Yule-Walker System

For each season $m$, the coefficients in standardized form ${\psi }_{m,1}^{\ast },\dots ,{\psi }_{m,p}^{\ast }$ satisfy:

$${\mathbf{R}}_m{\mathbf{\psi }}_m^{\ast }={\mathbf{r}}_m$$

where:

$${\mathbf{R}}_m=\left(\begin{array}{cccc}1& {\hat{\rho }}_{m-1}(1)& \cdots & {\hat{\rho }}_{m-p+1}(p-1)\\ {\hat{\rho }}_m(1)& 1& \cdots & {\hat{\rho }}_{m-p+2}(p-2)\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{\rho }}_m(p-1)& {\hat{\rho }}_{m-1}(p-2)& \cdots & 1\end{array}\right),{\mathbf{r}}_m=\left(\begin{array}{c}{\hat{\rho }}_m(1)\\ {\hat{\rho }}_m(2)\\ ⋮\\ {\hat{\rho }}_m(p)\end{array}\right)$$

The solution is:

$${\hat{\mathbf{\psi }}}_m^{\ast }={\mathbf{R}}_m^{-1}{\mathbf{r}}_m$$

The matrix ${\mathbf{R}}_m$ is not a standard Toeplitz matrix (because consecutive rows use different seasons’ correlations), but it has a similar structure. The correlation matrix must be positive definite for the solution to exist; if not, the historical record may be too short for the requested order.

Step 4 — Residual Standard Deviation

After solving the Yule-Walker system, the residual standard deviation for season $m$ is:

$${\hat{\sigma }}_m={\hat{s}}_m\sqrt{1-{\mathbf{\psi }}_m^{\ast \top }{\mathbf{r}}_m}$$

This equals ${\hat{s}}_m$ times the square root of the unexplained variance fraction. If ${\hat{\sigma }}_m^2\le 0$, the model overfits — it explains all historical variance, leaving no room for the noise term.

Step 5 — Convert to Original Units

The Yule-Walker solution yields coefficients in standardized form ${\psi }_{m,\ell }^{\ast }$ (dimensionless, relating standardized deviations). The LP requires original-unit coefficients:

$${\psi }_{m,\ell }={\psi }_{m,\ell }^{\ast }\cdot \frac{s_m}{s_{m-\ell }}$$

These are computed once at initialization and used directly as LP constraint matrix entries.

Key Properties

• Periodicity: All parameters vary by season, matching the strong seasonality of hydrological data.
• Parsimony: The model order $p$ is selected per season using AIC, BIC, or coefficient significance tests, avoiding unnecessary lags.
• Stationarity: Fitted models are validated to ensure the AR process does not diverge — the characteristic polynomial roots must lie outside the unit circle.
• Positive residual variance: After fitting, ${\sigma }_m^2>0$ must hold for all seasons. A zero or negative residual variance indicates overfitting.

Further Reading

• Stochastic Modeling — overview of why inflow uncertainty matters and how PAR(p) fits into the SDDP workflow
• Inflow Non-Negativity — how negative PAR(p) realizations are handled
• PAR Inflow Model Specification — complete mathematical specification with unit conversion formulas, runtime preprocessing steps, and all validation invariants

See also: For implementation details and usage, see the Stochastic Modeling guide in the software book.

Inflow Non-Negativity

Overview

The PAR(p) autoregressive model generates inflow realizations by adding a random innovation to a weighted combination of past inflows. Because the innovation term can be arbitrarily negative in the Gaussian model, the total inflow can sometimes be negative — which is physically impossible. Negative inflows would cause the water balance to add water to the reservoir rather than drain it, producing absurd dispatch decisions.

Handling this situation correctly is important for both LP feasibility and for maintaining the statistical properties of the inflow model. Three strategies are available, with different trade-offs between model fidelity, LP complexity, and computational cost.

How Negative Inflows Arise

The PAR(p) model at stage $t$ for hydro $h$ generates:

$$a_h=\underset{\text{deterministic base}}{\underbrace{{\mu }_m-\sum _{\ell =1}^p{\psi }_{\ell }{\mu }_{m-\ell }}}+\underset{\text{lag contribution}}{\underbrace{\sum _{\ell =1}^p{\psi }_{\ell }\cdot {\hat{a}}_{h,\ell }}}+\underset{\text{noise term}}{\underbrace{{\sigma }_m\cdot \eta }}$$

When the noise realization $\eta$ is sufficiently negative (e.g., $\eta <-2$), the total can become negative. For systems with low seasonal mean flows and high variability, this can occur with non-negligible probability in dry seasons.

If $a_h<0$ enters the water balance constraint uncorrected, the LP solver treats it as a negative inflow (water added to the river system), which has no physical meaning and will cause the optimizer to mismanage reservoir storage.

The Penalty Method (Recommended)

The recommended approach introduces a non-negative slack variable ${\sigma }_h^{inf}$ that absorbs any negative inflow, keeping the effective inflow non-negative while preserving the LP formulation and the AR model structure for positive realizations.

Modified AR constraint:

$$a_h+{\sigma }_h^{inf}=\text{deterministic base}+\sum _{\ell =1}^p{\psi }_{\ell }\cdot {\hat{a}}_{h,\ell }+{\sigma }_m\cdot \eta$$

with ${\sigma }_h^{inf}\ge 0$. When $\eta$ is mild, the optimizer sets ${\sigma }_h^{inf}=0$ and $a_h$ takes its natural value. When $\eta$ is extreme, ${\sigma }_h^{inf}$ absorbs the excess, so the effective inflow:

$$a_h^{eff}=a_h+{\sigma }_h^{inf}\ge 0$$

The violation is penalized in the objective outside the block summation:

$$+\sum _{h\in \mathcal{H}}{c}^{inf}\cdot {\sigma }_h^{inf}\cdot T$$

where $T={\sum }_k{\tau }_k$ is the total stage duration in hours. The factor $T$ converts the slack rate (m³/s) to an equivalent energy dimension over the full stage, making the penalty dimensionally consistent with storage violation penalties that appear in the same part of the objective.

The penalty $c^{inf}$ is a Category 2 constraint violation penalty — it sits above generation costs in the penalty hierarchy, ensuring the optimizer “pays” for negative inflow events in the value function and propagates this signal back through the Benders cuts.

Properties

• LP always feasible: the slack prevents infeasibility for any noise realization
• AR dynamics preserved for positive realizations: when ${\sigma }_h^{inf}=0$, the model is unchanged
• Clear cost signal: the penalty magnitude tracks how frequently and severely the threshold is violated

Truncation (Simple Alternative)

The simplest method truncates negative inflows to zero before they enter the LP:

$$a_h=\max \left(0,\text{deterministic base}+\sum _{\ell =1}^p{\psi }_{\ell }\cdot {\hat{a}}_{h,\ell }+{\sigma }_m\cdot \eta \right)$$

No additional LP variables or constraints are needed. However, truncation has two statistical drawbacks: it shifts the mean inflow upward (positive bias) and disrupts temporal correlation at truncation events, because the extreme noise that triggered truncation is discarded rather than tracked.

For quick studies or debugging, truncation is acceptable. For production use, the bias can distort long-run cost estimates in drought scenarios.

Deferred: Truncation with Penalty

A more sophisticated hybrid approach (deferred to post-v0.1.0) combines truncation with explicit noise adjustment tracking. It introduces a dimensionless slack ${\xi }_h\ge 0$ that adjusts the noise term rather than the inflow directly:

$$a_h=\text{deterministic base}+\sum _{\ell =1}^p{\psi }_{\ell }\cdot {\hat{a}}_{h,\ell }+{\sigma }_m\cdot (\eta +{\xi }_h),a_h\ge 0$$

The penalty is proportional to the actual inflow correction ${\sigma }_m\cdot {\xi }_h$, which is larger in high-variability seasons. This approach preserves AR model structure better than pure truncation and provides a more precise signal of how much the noise realization was adjusted. It is deferred because the interaction with the noise generation pipeline adds complexity that is not needed for the minimal viable solver.

Summary

Further Reading

• PAR(p) Autoregressive Models — the inflow model that generates the realizations handled here
• Inflow Non-Negativity (spec) — complete formal specification including the full LP formulations for all four methods and penalty hierarchy placement
• LP Formulation (spec) — penalty taxonomy and priority ordering

Spatial Correlation

Hydro plants in the same river basin or geographic region tend to experience similar hydrological conditions – when one plant receives high inflows, nearby plants often do too. Capturing this spatial correlation is essential for realistic scenario generation. This page summarizes how Cobre handles it. For the full autoregressive model specification, see PAR Inflow Model (spec).

Temporal correlation: the PAR(p) model

Each hydro plant’s inflow series is modeled independently by a Periodic Autoregressive model of order $p$ (PAR(p)). This model captures temporal correlation – the tendency for high-inflow months to follow high-inflow months – through autoregressive coefficients that link the current inflow to past values:

$$a_{h,t}={\mu }_{m(t)}+\sum _{\ell =1}^p{\psi }_{m(t),\ell }(a_{h,t-\ell }-{\mu }_{m(t-\ell )})+{\sigma }_{m(t)}\cdot {\epsilon }_{h,t}$$

The parameters (seasonal means $\mu$, AR coefficients $\psi$, residual standard deviations $\sigma$) vary by season, reflecting the strong seasonality of hydrological processes.

Spatial correlation: correlated noise terms

The PAR(p) model for each plant produces a residual (innovation term) ${\epsilon }_{h,t}$ that represents the unexplained randomness after accounting for the autoregressive structure. Spatial correlation is introduced by making these residuals correlated across plants.

Instead of drawing independent noise for each plant, the scenario generator samples a vector of correlated noise terms. Plants in the same river basin share a large portion of their randomness, while plants in distant basins are less correlated. The correlation structure is typically estimated from the historical residuals of the fitted PAR(p) models.

Effect on scenario generation

Correlated noise means that the scenarios used in the forward and backward passes reflect realistic joint behavior across plants:

• In a wet scenario, most plants in a basin receive above-average inflows simultaneously, leading to system-wide surplus.
• In a dry scenario, most plants are deficient together, stressing the thermal fleet.

Without spatial correlation, the scenario generator would produce unrealistic combinations – some plants wet while their upstream neighbors are dry – leading to policies that underestimate system-level risk.

Implementation notes

The spatial correlation structure is encoded in the noise vectors stored in the opening tree. Each noise vector contains entries for all hydro plants, and the cross-plant correlations are baked into the vector generation process. The PAR(p) model then transforms these correlated noise vectors into inflow realizations that respect both the temporal dynamics (autoregressive lags) and the spatial structure (correlated residuals).

Related topics

• PAR Inflow Model (spec) – Full specification of the autoregressive model, fitting procedure, and validation
• Scenario Generation – How scenarios incorporating spatial correlation are used in the SDDP passes
• SDDP Theory – The broader algorithmic context in which inflow modeling operates

Scenario Generation

SDDP requires scenarios – discrete realizations of uncertain quantities (primarily inflows) – at multiple points in the algorithm. This page provides an overview of how scenarios are generated and used. For the underlying stochastic model, see PAR Inflow Model (spec). For the forward/backward pass mechanics that consume these scenarios, see SDDP Algorithm (spec).

The opening tree

Before the SDDP iteration loop begins, Cobre generates a fixed opening tree: a set of pre-sampled noise vectors for each stage. Each noise vector contains correlated random draws for all hydro plants, reflecting the spatial correlation structure estimated from historical data (see Spatial Correlation).

The opening tree is generated once and remains fixed throughout training. This ensures that every iteration’s backward pass evaluates the same set of realizations, which is important for cut validity – cuts are valid lower bounds on the expected cost-to-go computed over a specific discrete distribution.

The number of openings per stage is a key configuration parameter. More openings give a finer approximation of the continuous inflow distribution but increase the computational cost of each backward pass proportionally.

Forward pass sampling

During the forward pass, each trajectory needs one scenario realization per stage. The default scheme (InSample) draws uniformly at random from the opening tree – the same set used by the backward pass. This means the forward pass explores states that are consistent with the distribution the cuts were built for.

Alternative sampling schemes are available:

• External: Draws from user-provided scenario data, allowing the forward pass to explore states driven by historical or synthetically generated inflow sequences.
• Historical: Uses actual historical inflow records as forward pass scenarios, useful for validating the policy against observed conditions.

Regardless of the sampling scheme, the backward pass always uses the complete opening tree.

Backward pass branching

The backward pass evaluates all openings in the tree at each visited trial point. This complete enumeration is what makes the resulting Benders cut a valid expected-value inequality: the cut coefficients are computed as probability-weighted averages over the full discrete distribution.

For each trial point ${\hat{x}}_{t-1}$ at stage $t$:

1. Solve the stage LP once per opening $\omega \in {\Omega }_t$, using ${\hat{x}}_{t-1}$ as the incoming state and $\omega$ as the inflow realization.
2. Extract dual multipliers and objective values from each solve.
3. Aggregate into a single cut (probability-weighted average of per-scenario coefficients).

The uniform probability $p(\omega )=1/|{\Omega }_t|$ applies when all openings are equally likely, which is the default.

From noise to inflows

The opening tree stores noise vectors (standardized innovations), not inflows directly. At each stage, the PAR(p) model transforms these noise vectors into physical inflow values using the autoregressive equation:

$$a_{h,t}=\text{deterministic base}+\text{lag contribution}+{\sigma }_{m(t)}\cdot {\eta }_t$$

where ${\eta }_t$ is the noise draw from the opening tree. The deterministic base and lag contribution depend on the incoming state (autoregressive lags), so the same noise vector produces different inflows depending on the system’s history. This is how the Markov property is maintained – the state carries enough information to reconstruct the inflow dynamics.

Scenario count and computational cost

The number of openings directly affects the backward pass cost: each backward stage solve requires $|{\Omega }_t|$ LP solves per trial point. Typical configurations use 20 to 200 openings per stage, balancing approximation quality against runtime.

The forward pass cost is independent of the opening count – each trajectory samples exactly one realization per stage.

Related topics

• PAR Inflow Model (spec) – The autoregressive model that transforms noise into inflows
• Spatial Correlation – How cross-plant correlation is embedded in the noise vectors
• Forward and Backward Passes – How the passes consume scenarios
• SDDP Algorithm (spec) – Full specification including scenario sampling schemes and opening tree structure

Cut Management

Cut management controls the lifecycle of Benders cuts in the SDDP algorithm: how cuts are generated, aggregated, stored, and selectively pruned. Effective cut management is essential for solver performance – without it, the number of constraints in each stage LP grows unboundedly, slowing every forward and backward pass solve.

What is a Benders cut?

A Benders cut is a linear inequality added to a stage’s LP that approximates the cost of future decisions. Each cut encodes information learned during a backward pass solve: how sensitive the total future cost is to changes in reservoir storage volumes and inflow history at that stage.

Mathematically, a cut takes the form:

$$\theta \ge \alpha +\sum _h{\pi }_h^v\cdot v_h+\sum _{h,\ell }{\pi }_{h,\ell }^{lag}\cdot a_{h,\ell }$$

where $\theta$ is the future cost variable, $\alpha$ is the intercept, ${\pi }_h^v$ captures the marginal value of water in reservoir $h$, and ${\pi }_{h,\ell }^{lag}$ captures the marginal value of inflow history. The collection of all cuts at a stage forms a piecewise-linear lower approximation of the true future cost function.

How cuts are generated

During the backward pass, the algorithm solves each stage’s LP at trial states collected in the forward pass, for every scenario in the opening tree. The LP dual variables (shadow prices) from the water balance and AR lag-fixing constraints provide the cut coefficients. The intercept is computed so that the cut passes exactly through the trial point.

Because multiple scenarios are evaluated, the per-scenario coefficients are aggregated into a single cut per trial point using probability-weighted averaging. This is the single-cut formulation used by Cobre.

Why cut management matters

Each SDDP iteration adds new cuts, so after hundreds of iterations the cut pool can contain thousands of constraints per stage. Most older cuts become redundant as newer, tighter cuts supersede them. Without pruning:

• LP solve time increases linearly with the number of active cuts
• Memory consumption grows proportionally
• Numerical conditioning degrades as near-parallel cuts accumulate

Cut selection strategies identify and deactivate redundant cuts, keeping the LP compact while preserving the quality of the future cost approximation. The available strategies (Level-1, LML1, dominated cut detection) offer different trade-offs between aggressiveness and safety. See Cut Selection for a summary of each strategy.

Importantly, deactivated cuts are not deleted from memory – they are relaxed to $-\infty$ so their indices remain stable. This preserves reproducibility across checkpoint/resume cycles and allows reactivation if needed.

Further reading

• Cut Management (spec) – complete formal specification: cut definition, dual extraction, aggregation formulas, validity conditions, selection strategies, convergence guarantees, and configuration parameters
• Single-Cut vs Multi-Cut – comparison of cut aggregation strategies
• Cut Selection – overview of the three selection strategies
• Benders Decomposition – the decomposition framework underlying cut generation
• Forward and Backward Passes – the iterative structure that drives cut generation

Single-Cut vs Multi-Cut

In SDDP, the backward pass evaluates multiple inflow scenarios at each trial state and must combine the resulting information into cuts for the previous stage. There are two approaches: the single-cut formulation (used by Cobre) and the multi-cut formulation (deferred to a future release).

Single-cut formulation

The single-cut approach aggregates per-scenario cut coefficients into one cut per trial point by computing the probability-weighted average:

$$\bar{\alpha }=\sum _{\omega }p(\omega )\cdot \alpha (\omega ),\bar{\pi }=\sum _{\omega }p(\omega )\cdot \pi (\omega )$$

This produces a single constraint added to the stage LP:

$$\theta \ge \bar{\alpha }+{\bar{\pi }}^{\top }x$$

Advantages:

• Each iteration adds exactly one cut per stage per trial point, keeping the LP small
• LP solve times remain fast, especially important for systems with many stages
• Simpler implementation and lower memory footprint

Trade-off:

• Each cut carries averaged information, so more iterations may be needed to achieve the same approximation quality

Multi-cut formulation

The multi-cut approach keeps one cut per scenario, introducing scenario-specific future cost variables ${\theta }_{\omega }$:

$${\theta }_{\omega }\ge \alpha (\omega )+\pi (\omega {)}^{\top }x,\forall \omega \in \Omega$$

with the linking constraint $\theta ={\sum }_{\omega }p(\omega )\cdot {\theta }_{\omega }$.

Advantages:

• Preserves more information per iteration – each scenario’s sensitivity is represented individually
• Can converge in fewer iterations for problems with high scenario variability

Trade-off:

• Each iteration adds $|\Omega |$ cuts per stage per trial point, leading to larger LPs
• LP solve time per iteration is higher

Which does Cobre use?

Cobre implements the single-cut formulation. For the system sizes typical of Brazilian hydrothermal dispatch (160+ hydro plants, hundreds of scenarios), the single-cut approach provides the best balance between iteration count and per-iteration solve time.

The multi-cut formulation is planned for a future release. See Deferred Features for status.

Further reading

• Cut Management (spec) – formal definition of single-cut aggregation (section 3) and the multi-cut deferral note
• Cut Management – overview of the full cut lifecycle
• Benders Decomposition – how cuts are derived from LP duals

Cut Selection

As SDDP iterates, the number of Benders cuts grows with each backward pass. Many older cuts become redundant once newer, tighter cuts are generated at nearby states. Cut selection strategies identify and deactivate these redundant cuts, keeping the stage LPs compact without sacrificing approximation quality.

Cobre supports three selection strategies, listed from least to most aggressive.

Cut activity

All three strategies rely on the concept of cut activity. A cut is active at a given state if it is binding (or nearly binding) at the LP optimum – meaning it is actually shaping the future cost approximation at that point. Inactive cuts sit below the current approximation surface and contribute nothing to the solution.

Activity is determined by checking whether the gap between the future cost variable ${\theta }^{\ast }$ and the cut’s value is below a configurable threshold $ϵ$:

$$\text{cut }k\text{ is active at }\hat{x}⟺{\theta }^{\ast }-({\alpha }_k+{\pi }_k^{\top }\hat{x})<ϵ$$

A threshold of zero means only strictly binding cuts count; a small positive value (e.g., 1e-6) accounts for numerical tolerance.

Level-1

The simplest strategy. A cut is retained if it was active at least once during the entire algorithm execution. Periodically, cuts that have never been active are deactivated.

• Retention policy: ever-active
• Aggressiveness: low – keeps cuts that were useful at any point in training
• Risk: may retain cuts that were active early but are now permanently dominated by newer cuts

Limited Memory Level-1 (LML1)

A refinement of Level-1 that adds a time horizon. Each cut is timestamped with the most recent iteration at which it was active. Cuts whose timestamp is older than a configurable memory window $W$ (in iterations) are deactivated.

• Retention policy: active within the last $W$ iterations
• Aggressiveness: moderate – controlled by the window size
• Risk: a very small window may prematurely remove cuts that are infrequently but meaningfully active

Dominated cut detection

The most aggressive strategy. A cut is dominated if, at every visited state, some other cut in the pool achieves an equal or higher value. Dominated cuts provide no value at any known operating point and can be safely deactivated.

Since verifying domination over the entire state space is intractable, Cobre checks domination only at the states visited during recent forward passes. As the set of visited states grows denser over iterations, the domination check becomes more reliable.

• Retention policy: not dominated at any visited state
• Aggressiveness: high – directly targets cuts that contribute nothing
• Cost: $\mathcal{O}(|\text{cuts}|\times |\text{visited states}|)$ per stage per check

Convergence guarantees

Level-1 and LML1 both preserve the finite convergence guarantee of SDDP: removing cuts that are never active at visited states does not affect the outer approximation quality at those states. This result is established in Bandarra and Guigues (2021).

Configuration

Cut selection is controlled by four parameters: the strategy (level1, lml1, or domination), the activity threshold, the frequency of selection checks (in iterations), and the memory window (for LML1). See Cut Management (spec) – section 9 for the parameter table.

Further reading

• Cut Management (spec) – formal definitions of cut activity, all three selection strategies, convergence theorem, and configuration parameters
• Cut Management – overview of the full cut lifecycle
• Convergence – how cut quality affects stopping criteria and bound computation

Risk Measures

Why risk matters in power system planning

Standard SDDP minimizes the expected total cost across all scenarios. This produces policies that perform well on average but can lead to extremely high costs in adverse scenarios – for example, a prolonged drought that depletes reservoirs and forces expensive thermal generation or load shedding.

Risk-averse SDDP addresses this by incorporating a risk measure into the optimization objective. Instead of minimizing only the average cost, the algorithm also penalizes outcomes in the tail of the cost distribution. The result is a policy that sacrifices a small amount of average-case performance for substantially better worst-case behavior.

Risk-neutral vs risk-averse

In risk-neutral SDDP, each backward pass aggregates per-scenario cut coefficients using the scenario probabilities $p(\omega )$. All scenarios contribute equally (weighted by probability) to the cut that approximates the future cost function.

In risk-averse SDDP, the scenario probabilities are replaced by risk-adjusted weights ${\mu }_{\omega }^{\ast }$ that place more emphasis on expensive (adverse) scenarios. The cut generation mechanics are otherwise identical – only the aggregation weights change.

The convex combination risk measure

Cobre uses a convex combination of expectation and Conditional Value-at-Risk (CVaR):

$${\rho }^{\lambda ,\alpha }[Z]=(1-\lambda )\mathbb{E}[Z]+\lambda \cdot {\text{CVaR}}_{\alpha }[Z]$$

This formulation has two parameters:

• $\lambda$ (risk aversion weight): Controls the blend between expected value and CVaR. Setting $\lambda =0$ recovers risk-neutral SDDP; setting $\lambda =1$ uses pure CVaR.
• $\alpha$ (confidence level): Controls how deep into the tail CVaR looks. Smaller $\alpha$ means focusing on more extreme adverse scenarios.

Both parameters can vary by stage, allowing operators to apply stronger risk aversion to near-term decisions (where consequences are immediate) and weaker risk aversion to distant stages.

Implications for convergence monitoring

A critical subtlety of risk-averse SDDP is that the first-stage LP objective is not a valid lower bound on the true risk-averse optimal cost. It remains a useful convergence indicator (it increases monotonically and plateaus), but it cannot be interpreted as a bound in the same way as in risk-neutral SDDP. Valid upper bounds require the inner approximation (SIDP) method.

Further reading

• CVaR – Intuitive explanation of Conditional Value-at-Risk
• Risk Measures (spec) – Complete formal specification: dual representations, subgradient theorem, risk-averse Bellman equation, cut generation algorithm, and lower bound validity warning
• Cut Management – How risk-adjusted weights integrate into cut aggregation
• Convergence – Bound computation and stopping criteria

Conditional Value at Risk (CVaR)

What is CVaR?

Conditional Value-at-Risk (CVaR) is a risk measure that captures the expected cost in the worst-case tail of a probability distribution. While the ordinary expected value averages over all outcomes equally, CVaR focuses on the most adverse fraction.

Formally, for a random cost $Z$ and a confidence level $\alpha \in (0,1]$:

$${\text{CVaR}}_{\alpha }(Z)=\underset{\eta \in \mathbb{R}}{\min }\left\{\eta +\frac{1}{\alpha }\mathbb{E}\left[(Z-\eta {)}^{+}\right]\right\}$$

The parameter $\alpha$ determines how deep into the tail the measure looks. Smaller $\alpha$ means a more pessimistic assessment.

Intuitive interpretation

Imagine sorting all possible cost outcomes from lowest to highest. CVaR${}_{\alpha }$ is the average cost among the worst $\alpha$-fraction of those outcomes.

• $\alpha =1$: The “worst 100%” is the entire distribution, so CVaR${}_1$ equals the ordinary expected value $\mathbb{E}[Z]$.
• $\alpha =0.5$: The average cost of the worst half of outcomes.
• $\alpha =0.05$: The average cost of the worst 5% of outcomes – a strongly pessimistic measure.

Relationship to VaR

Value-at-Risk (VaR) is the threshold below which the worst $\alpha$-fraction of outcomes falls. CVaR goes further: it averages all costs beyond that threshold. This makes CVaR strictly more conservative than VaR and, unlike VaR, it is a coherent risk measure (satisfying monotonicity, translation equivariance, positive homogeneity, and subadditivity).

Why CVaR for hydrothermal dispatch?

In hydrothermal systems, the primary source of uncertainty is river inflows. A sequence of dry years can deplete reservoirs and force the system to rely on expensive thermal generation or, in extreme cases, impose load curtailment. The expected cost may look acceptable, but the cost in drought scenarios can be catastrophic.

By incorporating CVaR into the SDDP objective, the optimizer builds policies that maintain higher reservoir levels as a hedge against drought. The resulting policy costs slightly more on average but avoids the worst tail outcomes – a trade-off that power system operators typically prefer.

How CVaR enters the SDDP algorithm

Cobre does not use pure CVaR in isolation. Instead, it uses a convex combination of expectation and CVaR:

$${\rho }^{\lambda ,\alpha }[Z]=(1-\lambda )\mathbb{E}[Z]+\lambda \cdot {\text{CVaR}}_{\alpha }[Z]$$

During the backward pass, this combination changes only the aggregation weights used to combine per-scenario cut coefficients into a single cut. Scenarios with higher costs receive proportionally more weight, while low-cost scenarios receive less. The LP formulation, dual extraction, and cut structure remain unchanged.

Further reading

• Risk Measures – Overview of risk-averse vs risk-neutral SDDP
• Risk Measures (spec) – Complete formal specification including the dual representation of CVaR, the sorting-based weight computation algorithm, and the critical warning about lower bound validity
• Cut Management (spec) – How CVaR-modified aggregation weights affect cut coefficient computation
• SDDP Theory – How the forward-backward iteration works

Roadmap

This section describes planned features for Cobre that are not yet implemented. The features listed here represent the team’s current thinking about where the solver is headed — they are directional plans, not committed timelines or version targets.

Each feature has been analyzed for feasibility and grouped by theme. The thematic pages below provide detail on motivation, prerequisites, and planned approach for each feature.

Note: Feature C.5 (Non-Controllable Sources) was originally deferred but is now fully implemented. It appears in the table below for historical reference with status “Implemented”. All other features remain planned.

Feature Index

Each feature page below provides details on motivation, prerequisites, and planned approach.

Equipment & Modeling

This page covers planned extensions to the set of power system equipment types and production models that Cobre can simulate. These features expand the physical scope of the solver to cover generation technologies and accuracy improvements not included in the initial release.

C.1 GNL Thermal Plants

Gas Natural Liquefeito (GNL) thermal plants operate under complex fuel supply contracts with minimum take-or-pay obligations, variable spot market fuel costs, and start-up and shutdown constraints. The existing thermal model covers simple dispatchable units with fixed cost; GNL requires modeling the fuel inventory as an additional state variable and enforcing contract fulfillment constraints over time.

Why it is deferred: GNL modeling requires integer variables for unit commitment decisions. Extending SDDP to handle integer subproblems requires SDDiP infrastructure (Lagrangian relaxation or other duality handlers) — a significant architectural change that is out of scope for the initial release.

Prerequisites:

• Duality handler infrastructure for MIP subproblems (Lagrangian relaxation or strengthened Benders cuts)
• MIP solver integration (Gurobi, CPLEX, or an open-source alternative such as Cbc)
• GNL-specific data model extensions (fuel inventory, contract parameters)

Estimated effort: Large (3-4 weeks). Requires SDDiP infrastructure before GNL-specific modeling can begin.

See also: Deferred Features §C.1

C.2 Battery Energy Storage Systems

Grid-scale batteries introduce a state-of-charge dimension alongside the existing reservoir storage state. Unlike hydro reservoirs, batteries have symmetric charge/discharge efficiency losses and a capacity that degrades over operational lifetime. The LP formulation is linear — no integer variables are required — making this a more tractable extension than GNL.

Why it is deferred: Batteries require a new state variable dimension (state of charge per battery), integration of charge and discharge flows into bus load balance constraints, and output schema additions for battery simulation results. The data model schemas for batteries are already fully defined, so the main effort is LP construction and testing.

Prerequisites:

• State variable infrastructure supports a dynamic number of storage dimensions
• Bus balance constraint generation handles battery charge and discharge contributions
• Output schema for battery simulation results implemented in cobre-io

Estimated effort: Medium (2-3 weeks). LP formulation is straightforward; the main effort is the data pipeline and realistic test cases.

See also: Deferred Features §C.2

C.5 Non-Controllable Sources (Wind/Solar)

Status: Implemented. Non-controllable sources (wind, solar, and other variable renewable generation) are fully specified and implemented. The LP formulation — including generation bounds, curtailment variables, and bus balance participation — is in Equipment Formulations §6. The entity schema, scenario pipeline, and curtailment penalty are also fully specified. This entry is retained for completeness.

See also: Equipment Formulations §6, Deferred Features §C.5

C.6 FPHA Enhancements

The core Four-Point Head Approximation (FPHA) model for hydro generation is implemented with a constant turbine efficiency. Three accuracy extensions are deferred:

Variable efficiency curves — Replace the constant efficiency coefficient with a flow-dependent hill chart approximation. The efficiency varies with turbine load as a polynomial function of the normalized flow. This captures the operating curve of real turbines and improves generation estimates away from the nominal operating point.

Pumped hydro production function — Extend the FPHA to model reversible hydro units that can operate in both generation and pumping mode. The pumping mode uses a reversed head approximation. Mutual exclusion between generation and pumping is enforced via a high-penalty relaxation (the alternative of binary SOS1 constraints is deferred along with general MIP support).

Dynamic FPHA recomputation — Periodically refit the FPHA hyperplanes based on the observed operating region during training, narrowing the volume window to match where the algorithm actually explores. This improves approximation quality in practice but requires careful management of cut validity across refitting events.

Why it is deferred: Core FPHA is operational and validated. These extensions are accuracy improvements with problem-dependent benefit. Variable efficiency requires hill chart data infrastructure; pumped hydro requires the unified gen/pump production function; dynamic recomputation requires online hyperplane fitting with cut validity tracking.

Prerequisites:

• Core FPHA operational and validated against reference implementations
• Hyperplane fitting infrastructure supports refitting with a new volume window
• Performance baseline established to measure accuracy improvements

Estimated effort: Medium-Large (2-4 weeks total across all three sub-features).

See also: Deferred Features §C.6, Hydro Production Models

Algorithm Variants

This page covers planned extensions to the core SDDP algorithm: alternative cut formulations, richer policy graph structures, and solver modes that address problem types not covered by the initial release.

C.3 Multi-Cut Formulation

The standard SDDP implementation uses a single-cut formulation: one aggregated Benders cut is added per backward pass iteration. The multi-cut formulation instead generates one cut per scenario in the opening tree. This provides more information per iteration at the cost of a larger LP at each stage.

The trade-off is well-studied: multi-cut converges in fewer iterations but each iteration is slower and uses more memory. Single-cut is better for large scenario counts and risk-averse configurations with CVaR; multi-cut is better for small scenario counts and risk-neutral settings.

Why it is deferred: Multi-cut requires the LP builder to handle a variable number of future cost variables (one per scenario instead of one aggregate), a cut pool indexed by scenario, and adapted cut selection logic. The interaction with CVaR risk measures requires careful analysis to ensure the cut validity guarantees are preserved.

Prerequisites:

• LP builder supports a configurable number of future cost variables per stage
• Cut pool is indexed by scenario identifier
• Cut selection strategy is adapted for per-scenario pools

Estimated effort: Medium (2-3 weeks). The core algorithm change has wide-reaching effects on cut management.

See also: Deferred Features §C.3, Cut Management

C.4 Markovian Policy Graphs

The current SDDP implementation assumes stage-wise independence: the same scenario tree structure is used at every stage regardless of which scenario was realized in the previous stage. Markovian policy graphs relax this assumption by introducing a Markov chain over “regimes” (for example, wet/dry/normal inflow regimes), where the transition to the next regime depends on the current one.

In this extension, the policy graph has nodes indexed by (stage, regime) rather than just stage. Cuts are regime-specific — they are only shared within nodes of the same regime. The state space grows by a factor equal to the number of Markov states.

Why it is deferred: Markovian policy graphs substantially increase algorithm complexity. The forward pass must track and sample regime transitions. The backward pass generates cuts for each (stage, markov_state) node. The cut pool, binary format, and cut selection all need to be extended with regime indexing. The data model is already prepared (the stages.json schema includes optional markov_states fields), but the algorithmic work is large.

Prerequisites:

• Policy graph supports multi-node-per-stage structure
• Cut pool indexed by (stage, markov_state) tuples
• Forward and backward pass logic handles Markov regime transitions
• Scenario generation respects regime-dependent inflow distributions

Estimated effort: Large (3-4 weeks). Fundamental change to the policy graph and forward/backward pass logic.

See also: Deferred Features §C.4, Input Scenarios

C.7 Temporal Scope Decoupling

The current SDDP formulation tightly couples three temporal scopes: the SDDP decomposition stage (where Benders cuts are generated), the physical time resolution of the constraints, and the stochastic process realization period. This forces a trade-off between fine temporal resolution (accurate physics, exponential cut growth) and coarse resolution (manageable cuts, poor short-term dynamics).

Temporal scope decoupling, inspired by the SPARHTACUS framework, allows each SDDP stage to contain multiple physical decision periods. Within a single stage, the LP chains multiple sub-period problems (period 1 feeds period 2, and so on). Cuts are still generated only at stage boundaries — they reference the state at the stage boundary, not intermediate periods — so the cut pool size does not increase even as temporal resolution within each stage improves.

The key advantage is that hourly or weekly resolution can be modeled within monthly SDDP stages without multiplying the cut pool size. The LP per stage grows proportionally to the number of sub-periods, but cut growth (the usual bottleneck) is unchanged.

Why it is deferred: This requires significant architectural changes to LP construction (multi-period sub-problem chaining within a stage), modified forward and backward pass logic, data model changes (a periods array within stages), and careful interaction with chronological block formulations.

Prerequisites:

• Core SDDP training loop validated and stable
• Chronological block formulation operational
• Performance profiling showing cut growth as the binding bottleneck

Estimated effort: Large (4-6 weeks). Fundamental change to LP construction, the data model, and the algorithm.

See also: Deferred Features §C.7

C.12 Complete Tree Solver

The complete tree solver replaces SDDP’s sample-based forward and backward passes with deterministic traversal of an explicit scenario tree. Instead of sampling scenarios, all nodes in the tree are enumerated and a subproblem is solved at each node. The result is an exact solution to the stochastic program (relative to the given tree discretization), matching the behavior of CEPEL’s DECOMP model.

The LP construction per node is identical to the SDDP stage subproblem — the complete tree mode reuses the existing LP builder. The difference is purely in the traversal logic (BFS/DFS enumeration of all tree nodes vs. random sampling) and cut aggregation (tree-weighted aggregation vs. sample-averaged aggregation).

Why it is deferred: Full tree traversal requires infrastructure for enumerating tree nodes, distributing them across MPI ranks, and storing results for potentially millions of nodes. The convergence criterion also changes: rather than a statistical lower bound, convergence means closing the Benders gap on the tree exactly.

Prerequisites:

• Core SDDP training loop operational
• Opening tree infrastructure supports variable per-stage branching counts
• External scenario integration validated
• Performance profiling establishes tractability bounds for target problem sizes

Estimated effort: Medium-Large (3-4 weeks). LP construction is reused; main effort is tree traversal, MPI node distribution, and result aggregation.

See also: Deferred Features §C.12, Scenario Generation §7

Stochastic Enhancements

This page covers planned improvements to the stochastic modeling and sampling layers of Cobre. These features extend the inflow model, allow finer temporal resolution, add flexibility to the noise generation pipeline, and introduce sampling variants for the backward pass and forward pass.

C.8 CEPEL PAR(p)-A Variant

The standard PAR(p) inflow model is implemented and validated. The CEPEL PAR(p)-A variant extends it with four refinements primarily relevant to basins with highly skewed inflow distributions:

• Maximum AR order fixed at 12 to enforce an annual cycle
• Stationarity enforcement via a coefficient constraint
• Lognormal transformation: the model fits on log(inflow) rather than the raw inflow, which prevents synthetic inflow values from going negative
• Regional correlation between hydro plants in the same river basin

The lognormal variant is mainly needed when inflow distributions are highly right-skewed and the standard model produces negative synthetic values with non-negligible probability. The inflow non-negativity handling strategies already mitigate this for the standard model in most basins.

Why it is deferred: The standard PAR(p) model covers most practical use cases. The PAR(p)-A variant is a refinement for specific basin characteristics. The current input format (inflow_seasonal_stats.parquet and inflow_ar_coefficients.parquet) is already compatible with PAR(p)-A — the lognormal transformation is applied during history preprocessing before computing the seasonal statistics.

Prerequisites:

• Standard PAR(p) fitting and validation operational
• Lognormal transformation infrastructure (log-space fitting, back-transformation)
• Validation suite comparing PAR(p) and PAR(p)-A on representative basins

Estimated effort: Small-Medium (1-2 weeks). The mathematical extensions are straightforward; the main effort is validation and testing.

See also: Deferred Features §C.8, PAR Inflow Model

C.10 Fine-Grained Temporal Resolution (Typical Days)

The current block formulation represents each SDDP stage with a small number of load blocks (for example, peak, shoulder, and off-peak) that aggregate monthly hours by load level. This is adequate for long-term planning but cannot model daily cycling patterns — daily storage behavior for pumped hydro and batteries, hourly renewable generation profiles, time-of-use pricing, or intra-day ramp constraints.

The planned extension introduces representative typical day types (for example, weekday and weekend) within each SDDP stage, each containing a sequence of chronological hourly blocks. This enables accurate daily sub-structure without increasing the cut pool size (cuts are still generated at stage boundaries, not within day types).

Key open design questions that require research before implementation:

• Day-type chaining: Should the end-of-day storage from one day type carry forward to the next day type within the same stage, or should each day type operate independently with a weighted average end-of-stage storage?
• Objective scaling: Block durations are actual hours (for water balance), but the objective contribution must be scaled by how many days each day type represents.
• Two-phase architecture: Training with aggregated blocks followed by simulation with full typical-day resolution. Cut compatibility conditions between the two phases need formal verification.

Why it is deferred: This feature requires research on the open design questions before implementation can begin. The LP size per stage grows substantially (from 3 blocks to 72 or more), and the interaction with existing block formulations is complex.

Prerequisites:

• Core parallel and chronological block modes operational and validated
• Simulation architecture supports variable block configurations
• Research on day-type chaining and two-phase cut compatibility completed

Estimated effort: Large (4-6 weeks). Significant input schema design, LP construction changes, and research required.

See also: Deferred Features §C.10

C.11 User-Supplied Noise Openings

The standard noise generation pipeline samples independent Gaussian noise vectors and applies Cholesky correlation. This pipeline covers most use cases, but some workflows require direct control over the noise values:

• Importing noise realizations from external stochastic models that use non-Gaussian distributions
• Reproducing exact noise sequences from legacy tools for validation
• Using domain-specific spatial correlation structures not captured by the Cholesky approach
• Research workflows where specific noise patterns are under study

The planned mechanism is a scenarios/noise_openings.parquet file. When this file is present, the scenario generator skips internal noise sampling and Cholesky correlation entirely, loading the user-supplied values directly into the opening tree.

Open design questions to resolve before implementation: the relationship to the existing external scenario mechanism (with noise inversion), what validation checks should be applied to user-supplied noise, whether load entity noise should be included alongside inflow noise, and whether separate noise sets are needed for forward and backward passes.

Prerequisites:

• Opening tree architecture implemented and validated
• Forward/backward noise separation operational
• Clear use case catalog that cannot be served by the existing external scenario mechanism

Estimated effort: Small (1 week). Input loading and validation are straightforward once the design questions are resolved.

See also: Deferred Features §C.11, Scenario Generation §2.3

C.14 Monte Carlo Backward Sampling

The standard backward pass evaluates all openings in the opening tree at each backward stage. When the opening count is large (500 or more, for high-fidelity tail representation), the backward pass dominates iteration time. Monte Carlo backward sampling replaces full enumeration with a sample of n openings drawn with replacement from the tree, reducing backward pass cost from O(N_openings) to O(n) LP solves per stage per trial point.

The resulting cut is an unbiased estimator of the full cut but with higher variance. The trade-off relative to the default full evaluation: fewer solves per iteration, but non-monotonic lower bound behavior (which complicates convergence monitoring) and slower per-iteration convergence.

Why it is deferred: The default complete evaluation is reliable and performant for typical production opening counts (50-200). Monte Carlo sampling introduces non-monotonic lower bounds, which require adaptations to convergence monitoring. The benefit depends on the chosen n relative to N_openings, which is problem-dependent.

Prerequisites:

• Core SDDP with complete backward sampling validated
• Convergence monitoring supports non-monotonic lower bounds
• Empirical study of the n vs. convergence rate trade-off

Estimated effort: Small (1 week). Sampling infrastructure is trivial; the main effort is convergence monitoring adaptation.

See also: Deferred Features §C.14

C.15 Risk-Adjusted Forward Sampling

Standard forward sampling draws scenarios uniformly from the opening tree. For strongly risk-averse policies (CVaR with weight lambda close to 1.0), uniform sampling under-explores the worst-case tail scenarios that matter most for the risk measure.

Risk-adjusted forward sampling biases the forward pass toward tail scenarios by evaluating the risk measure over candidate noise terms at each forward stage and re-weighting or re-sampling the next-stage noise accordingly. A complementary importance sampling variant weights forward trajectories by their likelihood ratio under the risk-adjusted distribution, enabling tighter upper bound estimates for risk-averse settings.

Why it is deferred: Integration with the CVaR risk measure configuration is non-trivial, and importance weights must be handled carefully in upper bound estimation. Default uniform sampling is sufficient for most applications; the benefit is primarily for strongly risk-averse configurations.

Prerequisites:

• CVaR risk measure implemented and validated
• Forward sampling scheme abstraction supports per-stage re-weighting
• Upper bound evaluation handles importance-weighted trajectories

Estimated effort: Medium (2-3 weeks). The algorithm is well-documented in the literature; the main effort is integration with the risk measure and convergence analysis.

See also: Deferred Features §C.15, Risk Measures

Performance & Parallelism

This page covers planned optimizations to the SDDP training loop that reduce computation time, improve state space exploration, or exploit parallel hardware more effectively. These features are motivated by profiling observations rather than new algorithmic capabilities — they make the existing algorithm faster or more thorough, without changing its fundamental behavior.

C.13 Alternative Forward Pass

The default forward pass solves the training LP, which excludes simulation-only features (linearized head, unit commitment, bottom discharge) to preserve convexity for valid cut generation. Trial points from this simplified model may not represent states that are realistic under full operational modeling, potentially slowing convergence.

The alternative forward pass maintains two LP models per stage: the training LP (convex, used for backward pass cut generation) and an alternative LP (includes simulation-only features, used for the forward pass to generate more realistic trial points). After each backward pass, new cuts are added to both models.

Cuts remain valid because they are always generated from the convex training LP. Trial points are more realistic because they come from the richer alternative LP. The trade-off: memory approximately doubles (two LP models per stage), and convergence may be slower because forward trial points are not optimal for the training model.

Why it is deferred: Maintaining two parallel LP models per stage doubles memory and complicates solver workspace management. The benefit depends on how much simulation-only features change the trial points, which is problem-dependent and requires empirical investigation.

Prerequisites:

• Core SDDP training loop operational and validated
• Simulation-only LP construction (linearized head, unit commitment) implemented
• Solver workspace infrastructure supports multiple LP models per stage

Estimated effort: Medium (2-3 weeks). LP construction infrastructure is reused; the main effort is dual LP management, cut copying, and convergence testing.

See also: Deferred Features §C.13, Simulation Architecture

C.16 Revisiting Forward Pass

Standard forward passes always start from the initial state at stage 1. In problems with large state spaces, repeated iterations can visit similar states and generate cuts in the same regions, leaving parts of the state space poorly covered.

The revisiting forward pass maintains a buffer of previously visited states from past iterations. With some probability, a forward pass starts from a randomly selected historical state at a randomly selected stage instead of from stage 1. The resulting trial points generate cuts in regions that standard forward passes have not recently visited.

Why it is deferred: Starting a forward pass from an intermediate stage produces a partial trajectory that does not contribute a full-horizon cost estimate. This complicates upper bound estimation. The benefit is primarily for problems where standard forward passes provably concentrate in a narrow region of the state space.

Prerequisites:

• Core SDDP training loop validated
• State history buffer infrastructure
• Convergence monitoring handles partial-trajectory iterations

Estimated effort: Small-Medium (1-2 weeks). The state buffer is straightforward; the main effort is convergence analysis for partial trajectories.

See also: Deferred Features §C.16

C.17 Forward Pass State Deduplication

When multiple forward scenarios visit identical or near-identical states at a given stage, the backward pass redundantly evaluates the cost-to-go from duplicate trial points. Identical states produce identical cuts, so the duplication adds no information.

Two approaches are planned:

• Exact deduplication — Hash-based deduplication of state vectors with identical values. Straightforward to implement but limited in benefit, since exact duplicates are uncommon in problems with continuous state spaces.
• Approximate deduplication — Merge states within a tolerance using spatial indexing (for example, k-d tree). Greater reduction in backward solves, but introduces approximation in cut coefficients. Requires analysis of convergence implications before deployment.

Why it is deferred: Exact duplicates are uncommon, making exact deduplication a marginal optimization. Approximate deduplication requires careful convergence analysis. This feature should be implemented only when profiling shows backward pass cost is dominated by trial point count rather than per-solve cost.

Prerequisites:

• Core SDDP training loop validated and profiled
• Backward pass timing data identifying trial point count as the dominant cost factor

Estimated effort: Small (less than 1 week) for exact deduplication; Medium (2-3 weeks) for approximate deduplication with convergence analysis.

See also: Deferred Features §C.17, Training Loop §5.2

C.18 Pipelined Backward Pass

The standard sequential backward pass uses the cut approximation from the current iteration at each stage. This requires a synchronization barrier between stages: stage t cannot begin its backward solve until stage t+1 has finished and broadcast its cuts.

The pipelined backward pass uses the previous iteration’s approximation instead. This eliminates per-stage barriers and allows backward stages to proceed in a pipeline — stage t begins its solve while stage t+1 is still communicating cuts via non-blocking MPI collectives (MPI_Iallgatherv). The result is overlapped communication and computation.

The trade-off: pipelined cuts are looser (older information), so more iterations are needed to converge. Whether the shorter wall-clock time per iteration outweighs the slower per-iteration convergence depends on hardware-specific communication latency, which cannot be determined without profiling production workloads.

When to consider: Profiling shows backward pass per-stage barrier synchronization accounts for more than 20% of backward pass time. Most likely with very large stage counts (more than 100 stages) and high inter-node network latency.

Why it is deferred: The sequential mode is simpler, produces tighter cuts, and is the default in SDDP.jl and most production implementations. The pipelining benefit is hardware-specific.

Prerequisites:

• Core SDDP training loop validated and profiled
• Backward pass timing data showing barrier synchronization as a bottleneck
• Non-blocking MPI collective support in ferrompi (MPI_Iallgatherv)

Estimated effort: Medium (2-3 weeks). Requires async state management, non-blocking collective wrappers, and convergence regression testing.

See also: Deferred Features §C.18, Communication Patterns

C.19 Dynamic Forward-Passes Scheduler

The number of forward passes per iteration (forward_passes) is currently set statically in the configuration. A dynamic scheduler would adjust this count during training based on convergence metrics (for example, increasing forward passes when the lower bound is stagnant) or wall-clock budget (for example, maintaining a fixed seconds-per-iteration target).

Why it is deferred: The cut pool is pre-allocated using the formula warm_start_cuts + max_iterations × forward_passes with a deterministic slot assignment. A dynamic scheduler would require either a worst-case static allocation (wasteful for most runs) or dynamic reallocation (which invalidates the deterministic slot assignment and complicates reproducibility). The right approach requires profiling data from production workloads.

See also: Deferred Features §C.19, Cut Management Implementation §1.3

Tooling & Interfaces

This page covers planned additions to the tooling ecosystem around the Cobre solver: validation tooling for trained policies, and the three interface crates that expose Cobre to different consumer contexts (AI agents via MCP, Python workflows via PyO3, and interactive terminal monitoring via ratatui).

C.9 Policy Compatibility Validation

When the solver runs in simulation-only mode, warm-start mode, or resumes from a checkpoint, it must validate that the current input data is compatible with the previously trained policy. The policy (cut coefficients) encodes LP structural information: the number of state variables, the constraint topology, the block configuration, the penalty values. If any of these change between training and simulation, the cuts become invalid and produce silently incorrect results.

The planned mechanism has two components:

1. Policy metadata record — During training, the solver persists a metadata record alongside the cut data, capturing all input properties that affect LP structure. The record is versioned for forward compatibility.
2. Validation phase — At the start of simulation or warm-start, the current input is compared against the metadata. Any mismatch produces a hard error with a clear description of what changed.

Properties that must be validated include: block mode and block count per stage, number of hydro plants and their AR orders, cascade topology, bus and line count, thermal and contract counts, production model per hydro, and penalty configuration. Some properties may allow controlled relaxation in the future (for example, adding a new thermal plant at the end of the generator list might be compatible if the state variable dimensions are unchanged), but the initial implementation validates everything strictly.

Why it is deferred: This is a cross-cutting concern touching training, simulation, checkpointing, and input loading. All prerequisite specs are now approved, and C.9 is unblocked pending a dedicated specification that defines the complete property list, the metadata format, the validation algorithm, and the error reporting format.

Prerequisites (all now met):

• Training loop architecture finalized — Training Loop approved
• Checkpoint format defined — Checkpointing and Binary Formats §3-§4 approved
• Simulation architecture finalized — Simulation Architecture approved
• Input loading pipeline defined — Input Loading Pipeline approved

Estimated effort: Medium (2-3 weeks). Primarily specification and testing; implementation is straightforward once the property list is defined.

See also: Deferred Features §C.9, Binary Formats

MCP Server (cobre-mcp)

The MCP server exposes Cobre as a first-class tool for AI coding agents via the Model Context Protocol. It runs as a standalone binary (cobre-mcp) or as a cobre serve subcommand, in single-process mode without MPI. Computation is delegated to cobre-sddp; all tool responses use the structured output envelope defined in Structured Output.

The server exposes Cobre’s run, validate, and report capabilities as MCP tools, with streaming progress updates for long-running training jobs. Multi-process execution via TCP or shared memory backends is architecturally possible but deferred pending demonstrated demand from MCP-based agent workflows.

The full specification is in MCP Server.

Python Bindings (cobre-python)

The Python bindings expose Cobre’s solver to Python workflows as the cobre module via PyO3. The module compiles to a shared library (cdylib) that is installed and imported like any other Python package.

Key design points: the GIL is released during all Rust computation; NumPy and Arrow FFI provide zero-copy data paths; MPI is prohibited from the Python interface (multi-process execution uses TCP or shared memory backends instead). The Python exception hierarchy maps to Cobre’s structured error kind registry.

The full specification is in Python Bindings.

Terminal UI (cobre-tui)

The terminal UI provides real-time monitoring of SDDP training runs in the terminal. It is built on ratatui and renders convergence plots, iteration timing breakdowns, cut statistics, resource usage, and MPI communication metrics.

The TUI operates in two modes: co-hosted within the cobre binary via an in-process broadcast channel (sharing the same process as the solver), or standalone reading JSON-lines from stdin for monitoring remote or already-running jobs. Interactive features — pause/resume, cut inspection, scenario comparison, stopping rule adjustment — all operate at iteration boundaries to preserve training loop atomicity.

The full specification is in Terminal UI.

Overview

Note: These specifications were written before implementation to define contracts and interfaces. With all 8 phases complete, the source of truth for behavior is the source code and the software book. These specs remain as reference for the design rationale behind implementation decisions.

Internal development artifacts: This section is intended for Cobre contributors and maintainers. The specs document design decisions, trait contracts, and implementation constraints that were authored before the code was written. If you are here to understand the mathematics or algorithms, the Theory section is the recommended starting point.

This section contains the foundational specifications that underpin all other sections of the Cobre documentation. It establishes the design philosophy that guides every architectural and algorithmic decision, the notation conventions that ensure consistency across all mathematical and data-model specs, and the production-scale reference dimensions that anchor performance analysis and memory budget calculations throughout the HPC specs.

These three specs are prerequisites for reading any other specification section. They are referenced pervasively – design principles inform trade-off discussions in the architecture specs, notation conventions define the symbols used in every equation, and production-scale dimensions appear whenever a spec quantifies communication volume, memory footprint, or computational cost.

Reading Order

The specs are short and self-contained; reading them in order takes under fifteen minutes:

1. Design Principles – The 9 foundational principles guiding the Cobre solver design. These principles are cited throughout the architecture and HPC specs to justify specific design choices.
2. Notation Conventions – Index sets, symbols, and unit conventions used throughout all mathematical formulations. Required context before reading any spec in the Mathematical Formulations section.
3. Production Scale Reference – Representative problem dimensions for a national-scale hydrothermal system (number of plants, stages, scenarios, state dimension). Referenced by every spec that performs sizing analysis or communication volume estimation.

Spec Index

Navigation

The Specifications section contains 7 subsections: Overview (this section), Mathematical Formulations, Data Model, Architecture, High-Performance Computing, Configuration, and Deferred Features. Readers encountering the specs for the first time should complete the three overview specs before proceeding to any other subsection.

Cross-Reference Index

This page provides a centralized navigation index for the 79 specification files in the Cobre documentation. It is designed to help implementers determine which specs to read for a given crate, in what order, and how specs relate to each other through cross-references.

The index has five components:

1. Spec-to-Crate Mapping Table – which crate(s) each spec belongs to
2. Per-Crate Reading Lists – ordered reading lists for each of the 11 crates
3. Outgoing Cross-Reference Table – for each spec, which other specs it references
4. Incoming Cross-Reference Table – for each spec, which other specs reference it
5. Dependency Ordering – a global topological reading order for the entire corpus

Errata resolved: Seven HIGH-severity section-number errors (F-1, F-5, F-6, F-7, F-8, F-9, F-11) were identified in the Epic 03 cross-reference audit and corrected in Epic 06 remediation. All section-number references now point to the correct targets.

1. Spec-to-Crate Mapping Table