This example combines PyMC probabilistic forecasting with the FICO® Xpress Optimization capabilities for power system planning under demand uncertainty, using real Ontario electricity demand data (2024–2025) sourced from the IESO Data Directory.
The study is presented in the FICO Xpress Community Blog post Getting Started with Optimization Under Uncertainty Using PyMC and FICO Xpress.
The full end-to-end pipeline is implemented in a single script: stochastic_energy_planning.py
A Bayesian time-series model is fitted on historical Ontario daily demand (IESO, 2024–2025) and used to generate an out-of-sample posterior predictive forecast for January–June 2026.
Model structure:
mu(t) = intercept + trend(t) + fourier_biyearly(t) + fourier_weekly(t)
y(t) ~ LogNormal(mu(t), sigma)
Fourier seasonality uses sin/cos pairs at bi-yearly and weekly frequencies. The posterior is sampled with NUTS (4 chains × 1,000 draws).
Part 1 outputs used by Part 2:
ontario-demand-forecast-2026-jan-july.nc— a NetCDF file generated by the script containing the full posterior predictive distribution: 4,000 daily demand scenarios (MWh/day) for 181 days. This file is not included in the repository; it is created on the first run and cached locally to avoid re-sampling.
From these 4,000 scenarios, Part 2 derives the 95th-percentile demand series used in Approach 1, and samples a 1,000-scenario subset for the CVaR model in Approach 2.
Two robust dispatch approaches are compared over a multi-source generation fleet (nuclear, hydro, gas, wind, solar), both using the scenario distribution produced by Part 1:
Replaces stochastic demand with its 95th-percentile value per day (derived from the 4,000 PyMC scenarios), guaranteeing feasibility in at least 95% of scenarios by construction. Runs on the Xpress community licence.
Minimizes generation cost plus a weighted CVaR term over 1,000 sampled scenarios from the PyMC posterior, using the Rockafellar-Uryasev LP reformulation. Controls tail risk explicitly via a risk-aversion weight.
Licence note: Approach 2 uses 1,000 demand scenarios, which exceeds the Xpress community edition limit. A free 30-day trial licence is available at fico.com/en/fico-xpress-trial-and-licensing-options and provides full solver capability. Approach 1 runs without restriction on the community licence.
A Pareto analysis sweeps CVaR risk-aversion weights across two clean energy policy levels (80% and 90%) to trace cost-vs-reliability frontiers.
| Symbol | Description |
|---|---|
| S | Generation sources: {nuclear, hydro, gas, wind, solar} |
| C ⊆ S | Clean (non-emitting) sources: {nuclear, hydro, wind, solar} |
| T | Planning horizon: days t = 1, ..., 181 |
| Ω | Full set of demand scenarios from PyMC posterior: |Ω| = 4,000 (4 chains × 1,000 draws) |
| Ω' ⊆ Ω | Sampled subset used in the CVaR model: |Ω'| = 1,000 |
| Symbol | Description |
|---|---|
| cap_s | Installed capacity of source s (MW) |
| cf_s | Capacity factor of source s (fraction) |
| avail_s = cap_s × cf_s × 24 | Available energy from source s per day (MWh/day) |
| cost_s | Marginal cost of source s ($/MWh) |
| clean_min | Minimum clean energy fraction (e.g., 0.80) |
| demand^β_t | β-th percentile of demand on day t, computed from all scenarios in Ω |
| demand_{ω,t} | Demand in scenario ω on day t (MWh), ω ∈ Ω' |
| tail_prob | CVaR tail probability (e.g., 0.05) |
| penalty | Shortage penalty ($/MWh) |
| risk_weight | CVaR weight (risk aversion parameter) |
| Symbol | Description |
|---|---|
| gen_{s,t} ≥ 0 | Energy generated by source s on day t (MWh), ∀s ∈ S, ∀t ∈ T |
| # | Constraint | Description |
|---|---|---|
| (1) | Σ_{s∈S} gen_{s,t} ≥ demand^{0.95}_t, ∀t∈T | Meet 95th percentile demand each day |
| (2) | gen_{s,t} ≤ avail_s, ∀s∈S, ∀t∈T | Available capacity limit |
| (3) | Σ_{s∈C} gen_{s,t} ≥ clean_min · Σ_{s∈S} gen_{s,t}, ∀t∈T | Clean energy policy (per day) |
| (4) | gen_{s,t} ≥ 0, ∀s∈S, ∀t∈T | Non-negativity |
demand^{0.95}_t is computed from all |Ω| = 4,000 PyMC scenarios, so (1) is a deterministic linear constraint. The model is an LP.
| Symbol | Description |
|---|---|
| shortage_{ω,t} ≥ 0 | Unmet demand in scenario ω on day t (MWh), ω ∈ Ω', t ∈ T |
| var_threshold ∈ ℝ | Value at Risk (VaR) threshold |
| cvar_aux_ω ≥ 0 | CVaR auxiliary variable for scenario ω, ω ∈ Ω' |
| # | Constraint | Description |
|---|---|---|
| (1) | gen_{s,t} ≤ avail_s, ∀s∈S, ∀t∈T | Available capacity limit |
| (2) | Σ_{s∈C} gen_{s,t} ≥ clean_min · Σ_{s∈S} gen_{s,t}, ∀t∈T | Clean energy policy (per day) |
| (3) | shortage_{ω,t} ≥ demand_{ω,t} − Σ_{s∈S} gen_{s,t}, ∀ω∈Ω', ∀t∈T | Shortage linearization |
| (4) | cvar_aux_ω ≥ penalty · Σ_{t∈T} shortage_{ω,t} − var_threshold, ∀ω∈Ω' | CVaR linearization |
| (5) | gen_{s,t} ≥ 0, shortage_{ω,t} ≥ 0, cvar_aux_ω ≥ 0, var_threshold ∈ ℝ | Variable bounds |
Constraints (3) and (4) together linearize max(0, shortage_cost_ω − var_threshold). The model is an LP.
| Chance-Constrained | CVaR | |
|---|---|---|
| Demand constraint | Hard: must meet demand^{0.95}_t | Soft: penalize shortages via tail cost |
| Decision variables | gen_{s,t} | gen_{s,t}, shortage_{ω,t}, var_threshold, cvar_aux_ω |
| Scenario data used | 95th percentile from Ω (4,000 scenarios) | All scenarios in Ω' (1,000 scenarios) |
| Risk control | Fixed confidence level β | risk_weight |
| Problem type | LP | LP |
| File | Type | Description |
|---|---|---|
ontario-demand-data-2024-2025-daily.csv |
Real | Historical daily electricity demand (MWh) for Ontario, 2024–2025. Used to train the PyMC model. Sourced from the IESO Data Directory. |
The generation fleet parameters (capacity, cost, capacity factors) hardcoded in the script are illustrative values, chosen to represent a realistic Ontario-like mix for educational purposes.
- Python 3.10+
- FICO Xpress Python API (
xpress) - PyMC 6+ and ArviZ 1+ (only required for the forecasting step; skipped if the cached
.ncfile exists) - NumPy, pandas, matplotlib, xarray, scipy
Install exact pinned versions (recommended):
pip install -r requirements.txtOr install without pinning:
pip install xpress xpresslibs pymc arviz numpy pandas matplotlib xarray scipyNote for Windows users: PyMC compiles the model graph to C++ via PyTensor. Without a C++ compiler, NUTS sampling falls back to pure Python and becomes significantly slower. See the PyMC installation guide for details.
python stochastic_energy_planning.pyThe first run performs PyMC sampling (~5–10 minutes with a C++ compiler) and caches the
forecast to ontario-demand-forecast-2026-jan-july.nc. Subsequent runs skip sampling and
load from cache.
See legal and license requirements in the main repository README.