Capacity Expansion (Master Problem)¶
The Master Problem (master_problem.jl) determines optimal investment and retirement decisions over a multi-year planning horizon by minimizing the net present value (NPV) of total system costs. Representative days approximate operational costs and validate that investment plans are operationally feasible. Notation follows the conventions in the Formulation Overview.
1. Overview¶
The Master Problem is a two-stage linear program [24], [25]:
- First stage (investment): Capacity additions for generators, batteries, transmission, and primary energy infrastructure across all planning years.
- Second stage (operational validation): Representative-day dispatch subproblems embedded within the master model to link investment decisions to operational feasibility.
Retirement is handled through age-based expiry -- units automatically leave service when their age exceeds their configured lifetime. No binary retirement decisions are needed, keeping the formulation as a pure LP.
Entry point: create_master_problem(input; use_representative_days=true)
2. Sets¶
| Symbol | Description | Source |
|---|---|---|
| \(\mathcal{Y} = \{1, \ldots, Y\}\) | Planning years | input.years |
| \(\mathcal{G} = \{1, \ldots, G\}\) | Generator types | input.generators |
| \(\mathcal{G}_{RE} \subseteq \mathcal{G}\) | Renewable generators | gen.type == "Renewable" |
| \(\mathcal{B} = \{1, \ldots, B\}\) | Battery types | input.batteries |
| \(\mathcal{N} = \{1, \ldots, N\}\) | Nodes (buses) | input.network.num_buses |
| \(\mathcal{T}_d = \{1, \ldots, H\}\) | Hours within a representative day | 24 / temporal_resolution_hours |
| \(\mathcal{D}_y\) | Representative days for year \(y\) | select_representative_days() |
| \(\mathcal{F}\) | Fuel types (primary energy) | input.pe_configs |
| \(\mathcal{S}\) | Demand sectors | input.sectoral_demand |
3. Decision Variables¶
3.1 Investment Variables¶
| Variable | Domain | Units | Description | Julia name |
|---|---|---|---|---|
| \(I^{gen}_{y,g,n}\) | \(\mathbb{R}_+\) | MW | Generator capacity investment | gen_investment[y][g][n] |
| \(I^{bat,P}_{y,b,n}\) | \(\mathbb{R}_+\) | MW | Battery power investment | bat_power_investment[y][b][n] |
| \(I^{bat,E}_{y,b,n}\) | \(\mathbb{R}_+\) | MWh | Battery energy investment | bat_capacity_investment[y][b][n] |
| \(I^{tr}_{y,i,j}\) | \(\mathbb{R}_+\) | MW | Transmission expansion | transfer_investment[y][(i,j)] |
| \(I^{fs}_{y,f,n}\) | \(\mathbb{R}_+\) | units | Fuel storage investment | fuel_storage_investment[f][y][n] |
| \(I^{ft}_{y,f,i,j}\) | \(\mathbb{R}_+\) | units/day | Fuel transport investment | fuel_transport_investment[f][y][(i,j)] |
3.2 RE Tracking and Slack Variables¶
| Variable | Domain | Units | Description | Julia name |
|---|---|---|---|---|
| \(\rho_y\) | \([0, 1]\) | -- | RE penetration ratio | re_penetration_ratio[y] |
| \(s^{re}_y\) | \(\mathbb{R}_+\) | -- | RE target slack | slack_re_target[y] |
| \(s^{cap}_{y,n}\) | \(\mathbb{R}_+\) | MW | Capacity adequacy slack | slack_capacity[(y,n)] |
| \(s^{bud}_y\) | \(\mathbb{R}_+\) | $ | Budget slack | slack_budget[y] |
3.3 Representative Day Operational Variables¶
For each year \(y\) and representative day \(d\), a set of dispatch variables is created (struct RepresentativeDayVariables):
| Variable | Domain | Units | Description | Julia name |
|---|---|---|---|---|
| \(P^{day}_{g,n,t}\) | \(\mathbb{R}_+\) | MW | Generator output | gen_output[g,n,t] |
| \(C^{day}_{g,n,t}\) | \(\mathbb{R}_+\) | MW | Curtailment per generator | curtailment[g,n,t] |
| \(P^{ch,day}_{b,n,t}\) | \(\mathbb{R}_+\) | MW | Battery charge | bat_charge[b,n,t] |
| \(P^{dis,day}_{b,n,t}\) | \(\mathbb{R}_+\) | MW | Battery discharge | bat_discharge[b,n,t] |
| \(E^{day}_{b,n,t}\) | \(\mathbb{R}_+\) | MWh | Battery SOC | bat_soc[b,n,t] |
| \(L^{day}_{n,t}\) | \(\mathbb{R}_+\) | MW | Loss of load | loss_load[n,t] |
| \(\phi^{day}_{i,j,t}\) | \(\mathbb{R}\) | MW | Transfer flow | transfer[(i,j)][t] |
| \(\lambda^{fre}_{n,t}\) | \(\mathbb{R}_+\) | MW | FRE penetration loss slack | fre_penetration_loss[n,t] |
| \(L^{sec,day}_{s,n,t}\) | \(\mathbb{R}_+\) | MW | Sectoral loss of load | loss_of_load_sectoral[s][n,t] |
4. Parameters¶
4.1 Generator Parameters¶
| Symbol | Units | Description | Julia field |
|---|---|---|---|
| \(\bar{P}_{g,n}\) | MW | Rated power | gen.rated_power[n] |
| \(c^{inv}_{g,n}\) | $/MW | Investment cost | gen.invest_cost[n] |
| \(\bar{I}_{g,n}\) | MW | Maximum investment | gen.invest_max[n] |
| \(\tau_{g,n}\) | years | Lifetime | gen.life_time[n] |
| \(a^0_{g,n}\) | years | Initial age | gen.initial_age[n] |
| \(\delta_{g,n}\) | [0,1] | Annual degradation rate | gen.degradation_rate[n] |
| \(\alpha_{g,t,n}\) | [0,1] | Availability factor | gen.availability[t,n] |
| \(c^{fuel}_{g,n}\) | $/MWh | Fuel cost | gen.fuel_cost[n] |
| \(c^{fix}_{g,n}\) | $/MWh | Fixed O&M cost | gen.fixed_cost[n] |
| \(c^{maint}_{g,n}\) | $/MWh | Maintenance cost | gen.maintenance_cost[n] |
4.2 Battery Parameters¶
| Symbol | Units | Description | Julia field |
|---|---|---|---|
| \(\bar{P}^{bat}_{b,n}\) | MW | Max discharge power | bat.max_discharge_power[n] |
| \(\bar{E}_{b,n}\) | MWh | Energy capacity | bat.capacity[n] |
| \(\eta^{ch}_{b,n}\) | [0,1] | Charge efficiency | bat.charge_efficiency[n] |
| \(\eta^{dis}_{b,n}\) | [0,1] | Discharge efficiency | bat.discharge_efficiency[n] |
| \(E^0_{b,n}\) | [0,1] | Initial SOC fraction | bat.soc_initial[n] |
| \(h^{min}_{b}\) | hours | Minimum duration | bat.min_duration_hours |
| \(h^{max}_{b}\) | hours | Maximum duration | bat.max_duration_hours |
| \(c^{\text{thr}}_{b,n}\) | $/MWh | Throughput degradation cost | bat.throughput_degradation_cost[n] |
4.3 System Parameters¶
| Symbol | Units | Description | Julia field |
|---|---|---|---|
| \(r\) | [0,1] | Discount rate | input.discount_rate |
| \(\bar{B}_y\) | $ | Max annual investment budget | input.max_annual_investment |
| \(\rho^{target}\) | [0,1] | Target RE penetration (final year) | input.target_re_penetration |
| \(\rho^{init}\) | [0,1] | Initial RE penetration | input.initial_re_penetration |
| \(\Delta\rho^{min}\) | [0,1] | Minimum annual RE increment | input.min_re_increment |
| \(\Delta\rho^{max}\) | [0,1] | Maximum annual RE increment | input.max_re_increment |
| \(M_{res}\) | -- | Reserve margin multiplier | input.reserve_margin (default 1.15) |
| \(\gamma\) | -- | Annual demand growth rate | input.demand_growth |
| \(c^{VOLL}\) | $/MW | Value of lost load penalty | input.loss_of_load_penalty |
| \(c^{fre}\) | $/MWh | FRE penetration loss penalty | input.fre_penetration_loss_penalty |
| \(c^{slack}\) | $ | Generic slack penalty | input.slack_penalty |
| \(\bar{\kappa}\) | [0,1] | Maximum curtailment ratio | input.max_curtailment_ratio (default 0.05) |
| \(D_r\) | -- | Representative days per year | input.representative_days_per_year |
5. Objective Function¶
Julia function: build_master_objective!(model, vars, input)
The objective minimizes the NPV of total system costs over the planning horizon:
5.1 Investment Cost¶
where \(d_{ij}\) is the distance between nodes \(i\) and \(j\), and \(c^{ft}_f\) is the transport investment cost per unit per km for fuel \(f\).
5.2 Operational Cost¶
Operational costs come from the representative day subproblems, scaled to annual:
where the cost for each representative day is (from calculate_day_operational_cost):
where \(\kappa_s\) is the criticality weight for sector \(s\) and \(\Delta t\) is temporal_resolution_hours.
PWL cost curves in the master problem
When generators or technologies have a fuel_cost_curve configured, the representative day operational cost (OBJ-5) uses the same piecewise-linear (PWL) fuel cost decomposition as the operational dispatch formulation. Generator output is decomposed into segments with non-decreasing marginal costs (see Operational Dispatch -- PWL Fuel Cost). The same applies to batteries with discharge_cost_curve. This ensures that investment decisions in the master problem account for the actual shape of generator cost curves rather than assuming flat marginal costs.
5.3 Slack Penalties¶
6. Cumulative Capacity Expressions¶
Julia function: build_cumulative_capacity_expressions(vars, input, year_idx)
A key building block for all constraints. For each technology and node, the cumulative available capacity in year \(y\) combines existing (degraded) capacity with invested capacity, both subject to age-based retirement.
6.1 Generator Cumulative Capacity¶
6.2 Age Calculations¶
| Unit Type | Age Formula | Indicator Variable | Active Condition |
|---|---|---|---|
| Existing | \(a^{exist}_{y,g,n} = a^0_{g,n} + (y - 1)\) | -- | \(a^{exist}_{y,g,n} < \tau_{g,n}\) |
| Investment (year \(y'\)) | \(a^{inv}_{y,y'} = y - y'\) | -- | \(a^{inv}_{y,y'} < \tau_{g,n}\) |
The indicator function \(\mathbb{1}[\cdot]\) is evaluated at model construction time (not as a binary variable), which keeps the formulation as a pure LP.
6.3 Battery Cumulative Capacity¶
Battery cumulative power and energy follow the same structure but without degradation:
6.4 Transmission Cumulative Capacity¶
7. Constraint Families¶
7.1 Investment Limits (INV)¶
Julia function: add_investment_constraints!(model, vars, input)
INV-1: Cumulative investment limit per generator. Total investment across all years must not exceed the maximum allowed:
INV-2: Cumulative battery power investment limit.
INV-3: Cumulative battery energy investment limit.
INV-4: Battery duration constraints. Applied per year to each year's investment independently (not cumulative):
7.2 Budget Constraint (BUD)¶
Julia function: add_budget_constraints!(model, vars, input)
The annual investment expenditure must not exceed the budget (with slack):
7.3 Transmission Symmetry (TXN)¶
Julia function: add_transmission_symmetry_constraints!(model, vars, input)
Bidirectional transmission investment is symmetric:
7.4 Capacity Adequacy (CAP)¶
Julia function: add_capacity_adequacy_constraints!(model, vars, input)
Total generation and storage capacity must meet peak demand with a reserve margin at each node:
where the peak demand at node \(n\) in year \(y\) is:
The full expression for total capacity in CAP-1 sums across all generators and batteries using the cumulative expressions from equations (CUM-1) and (CUM-2).
7.5 RE Target Constraints (RE)¶
Julia functions: calculate_target_ratios(input), add_re_target_constraints!(), add_re_increment_constraints!()
RE-1: Target ratio calculation. Linear interpolation from initial to target penetration:
with annual increment clamping:
RE-2: RE penetration equality. The ratio variable is forced to match the target exactly:
RE-3: RE increment bounds. Annual change in the ratio variable is bounded:
7.6 Age-Based Retirement (RET)¶
Julia function: add_retirement_cascade_constraints!() (no-op; retirement handled implicitly)
Retirement is enforced structurally through the cumulative capacity expressions (Section 6). A unit contributes capacity only if its age is below its lifetime:
RET-1: Existing unit retirement.
RET-2: Investment retirement. An investment made in year \(y'\) retires when:
These conditions are evaluated at model construction time. When the condition is false, the corresponding term is simply omitted from the capacity expression -- no binary variables are needed.
8. Representative Day Validation¶
Representative day subproblems link strategic investment decisions to operational feasibility. Without them, the Master Problem could select investments that are infeasible in practice.
8.1 Representative Day Selection¶
Julia function: select_representative_days(demand, year_idx, num_days, min_separation, timesteps_per_day, timesteps_per_year)
The algorithm selects high-demand days with temporal diversity:
- Compute daily peak demand for each day in the year: \(D^{peak}_d = \max_t \sum_n D_{d,t,n}\)
- Select the global peak day first.
- Divide the remaining year into \(D_r\) segments for seasonal diversity.
- Within each segment, select the day with the highest peak that satisfies a minimum separation constraint (\(\geq\)
min_day_separationdays from any already-selected day).
8.2 Operational Variables¶
Julia function: create_day_operational_vars!(model, input, year_idx, day_idx, hours)
For each year \(y\) and representative day \(d\), a full set of hourly dispatch variables is created as described in Section 3.3. These variables are embedded within the same JuMP model as the investment variables, allowing the solver to jointly optimize investment and dispatch.
8.3 Operational Constraints¶
Julia function: add_day_operational_constraints!(model, day_vars, vars, input, year_idx, day_idx, demand, start_hour)
8.3.1 Power Balance¶
where the demand at bus \(n\), hour \(t\) in year \(y\) is:
Note that curtailment does not appear in the power balance. It is handled separately per generator for renewables.
8.3.2 Generator Capacity Linked to Investment¶
For renewable generators:
where \(t' = \text{mod1}(\text{start\_hour} + t - 1,\; 8760)\) maps the representative day hour to the annual availability profile.
For conventional generators:
The cumulative capacity \(\bar{P}^{cum}_{y,g,n}\) is the expression from Eq. (CUM-1), which is a linear expression in the investment variables. This is the key coupling between investment decisions and dispatch feasibility.
8.3.3 Battery Constraints Linked to Investment¶
Power limits:
SOC capacity limit:
SOC dynamics:
Cyclic SOC constraint (prevents batteries from acting as infinite energy sources):
8.3.4 Transmission Constraints Linked to Investment¶
PWL Transmission Losses in the Master Problem
When DC power flow is enabled, the master problem uses the same piecewise linear (PWL) loss model as operational dispatch, but with fewer segments (default 2 vs. 3) for computational performance. Losses are split 50/50 between the two bus endpoints of each line (half-loss split), ensuring symmetric loss allocation regardless of flow direction. The number of segments is configurable via dc_power_flow.pwl_loss_segments_master. See DC Power Flow -- Transmission Losses for the full PWL formulation.
8.3.5 RE Penetration per Day¶
Total renewable generation must meet the target ratio for each representative day:
where \(D^{total}_{y,d} = \sum_{n,t} D^{day}_{n,t}\) is the total demand for this day (with growth). The slack variable \(\lambda^{fre}\) allows soft enforcement, penalized in the objective.
Additionally, the ratio variable is linked to actual generation:
8.3.6 Curtailment Limit per Day¶
When max_curtailment_ratio < 1.0, curtailment is limited to a fraction of renewable generation:
This forces investment in storage rather than relying on curtailment to manage surplus RE.
8.3.7 Sectoral Loss-of-Load Bounds¶
Each sector's load shedding cannot exceed the sector's demand at that node and hour:
8.4 Cost Scaling¶
Each representative day's operational cost is scaled to approximate annual costs:
Julia function: add_representative_days_validation!(model, vars, input, targets) orchestrates the creation of variables, constraints, and cost expressions for all representative days across all years.
8.5 Time-Series Aggregation Method (TSAM)¶
When use_tsam = true, the peak-demand-based day selection (Section 8.1) is replaced by data-driven clustering with inter-period SOC linking.
Python function: compute_tsam_periods() in models/tsam.py
Julia function: add_tsam_periods_validation!(model, vars, input, targets)
8.5.1 Clustering Algorithm¶
The annual demand series is reshaped into daily blocks and clustered using k-medoids (default) or k-means:
- Build feature matrix \(F \in \mathbb{R}^{365 \times (H \cdot N)}\) from daily demand blocks
- Optionally concatenate normalized availability profiles as additional features
- Standardize features (zero mean, unit variance)
- Apply k-medoids: select \(K\) cluster medoids as representative periods
- Assign each original day to the nearest cluster
8.5.2 Period Weights¶
Each representative period \(p\) has a weight equal to its cluster size:
The operational cost scaling becomes:
This replaces the uniform \(365/N\) scaling (SCALE-1).
8.5.3 Inter-Period SOC Linking¶
To enable seasonal storage representation, inter-period SOC boundary variables are introduced:
Variables: \(\text{SOC}^{bnd}_{y,b,n,p}\) for \(p = 0, 1, \ldots, K\), where \(p=0\) represents the year start.
SOC chain (periods ordered chronologically):
Year-cyclic constraint:
Capacity bound:
When tsam_inter_period_linking = false, the standard cyclic SOC constraint (SOC-CYC) is used per period instead, and only the weighted cost scaling is active.
9. NPV-Based Iterative Retirement¶
Julia function: solve_with_npv_iteration(input; max_iterations=5, npv_threshold=0.0)
After the initial solve, this procedure identifies units that are economically unviable and forces their retirement through iterative re-optimization.
9.1 Algorithm¶
1. Solve initial Master Problem
2. For iter = 1 to max_iterations:
a. For each unit (generator g at node n, battery b at node n):
- Calculate NPV over remaining lifetime
b. Identify units with NPV < threshold
c. If none found: CONVERGED, return
d. Force retirement of negative-NPV units
e. Re-solve Master Problem
3. Return result (converged or max iterations reached)
9.2 Unit NPV Calculation¶
Julia function: calculate_unit_npv(result, input, unit_type, unit_idx, node)
For a generator with remaining life \(L_{rem} = \max(0, \tau_{g,n} - a^0_{g,n})\):
where:
- Revenue estimate: \(R_y = \bar{P}_{g,n} \cdot c^{inv}_{g,n} \cdot r_{return}\) (capacity value at configurable return rate, default 10% for generators)
- Cost: \(C_y = \bar{P}_{g,n} \cdot (c^{fix}_{g,n} + c^{maint}_{g,n}) \cdot 8760\)
Units with \(\text{NPV}_g < 0\) are flagged for retirement.
9.3 Result¶
The function returns an NPVIterationResult containing:
- Number of iterations performed
- Whether convergence was achieved
- Final
MasterProblemResult - List of all forced retirements (
UnitNPVrecords) - NPV history across iterations
10. Multi-System Master Problem¶
Julia function: create_multi_system_master_problem(input::MultiSystemMasterInput)
For planning across multiple interconnected power systems (e.g., countries or regions).
10.1 Structure¶
Each system \(s\) has its own set of investment variables, operational sub-models, and RE targets. Systems are coupled through inter-system transmission links.
10.2 Inter-System Investment Variables¶
For each link \(l\) connecting system \(s_1\) (node \(n_1\)) to system \(s_2\) (node \(n_2\)):
Cumulative limit:
Symmetry for bidirectional links (where link \(l'\) is the reverse of \(l\)):
10.3 Inter-System Operational Coupling¶
Julia function: add_inter_system_operational_coupling!(model, ext_vars, input)
For each representative day, forward and reverse flow variables are created for each inter-system link:
where the cumulative inter-system capacity is:
Border node injection:
At the source node:
At the destination node:
where \(\ell_l\) is the loss factor for link \(l\).
10.4 Multi-System Objective¶
Julia function: build_multi_system_objective!(model, ext_vars, input)
where \(Z_s\) is the per-system cost from Eq. (OBJ-1), \(c^{inter}_l\) is the inter-system investment cost per MW, and \(c^{flow}_l \cdot d^{km}_l\) is the distance-dependent operational flow cost.
11. Stochastic Extension¶
Julia function: create_stochastic_master_problem(input::StochasticMasterInput)
The stochastic formulation extends the deterministic Master Problem to handle uncertainty through scenario-based optimization.
11.1 Scenario Definition¶
Each scenario \(\omega\) has:
- A probability weight \(\pi_\omega\) where \(\sum_\omega \pi_\omega = 1\)
- A set of multipliers (
ScenarioMultipliers) that scale costs and parameters:
| Multiplier | Description |
|---|---|
invest_cost_renewables |
Scales RE investment costs |
invest_cost_conventional |
Scales conventional investment costs |
fuel_cost |
Scales fuel costs |
maintenance_cost |
Scales maintenance costs |
invest_cost_storage |
Scales storage investment costs |
invest_cost_transmission |
Scales transmission investment costs |
discount_rate |
Scenario-specific discount rate adjustment |
demand_growth |
Scenario-specific demand growth adjustment |
11.2 Two-Stage Structure¶
First stage (here-and-now): Investment decisions \(I^{gen}_{y,g,n}\), \(I^{bat,P}_{y,b,n}\), \(I^{bat,E}_{y,b,n}\), \(I^{tr}_{y,i,j}\) are common across all scenarios. They use base (unmodified) investment costs.
Second stage (wait-and-see): Operational costs are scenario-dependent. Modified generators and batteries are created by applying ScenarioMultipliers via apply_scenario_multipliers().
11.3 Stochastic Objective¶
Julia function: build_stochastic_objective!(model, vars, input, scenarios)
where \(r_\omega = r \cdot m^{discount}_\omega\) is the scenario-adjusted discount rate and \(C^{op}_{y,\omega}\) are the operational costs under scenario \(\omega\).
12. Solution Extraction¶
Julia function: extract_master_solution(model, vars, input)
After optimization, the following quantities are extracted:
| Output | Description | Computation |
|---|---|---|
gen_investment[y][g] |
Generator investment per year and node (MW) | Direct variable values |
bat_power_investment[y][b] |
Battery power investment (MW) | Direct variable values |
bat_capacity_investment[y][b] |
Battery energy investment (MWh) | Direct variable values |
transfer_investment[y][(i,j)] |
Transmission investment (MW) | Direct variable values |
cumulative_gen_capacity[y][g] |
Cumulative gen capacity with degradation | Eq. (CUM-1) evaluated |
cumulative_bat_capacity[y][b] |
Cumulative battery energy capacity | Eq. (CUM-3) evaluated |
re_penetration_by_year[y] |
Achieved RE penetration ratio | value(re_penetration_ratio[y]) |
total_investment_by_year[y] |
Total investment cost per year ($) | Sum of all investment costs |
total_operational_cost_by_year[y] |
Total operational cost per year ($) | From representative days |
gen_life_extension[y][g] |
Retirement status (1.0 = active, 0.0 = retired) | Age check against lifetime |
13. Infeasibility Diagnostics¶
Julia function: diagnose_infeasibility(model, vars, input)
If the model solves successfully but with non-zero slack, the diagnostics report which constraints required relaxation:
| Slack Variable | Indicates |
|---|---|
| \(s^{re}_y > 0\) | RE target cannot be met in year \(y\) |
| \(s^{bud}_y > 0\) | Investment budget exceeded in year \(y\) |
| \(s^{cap}_{y,n} > 0\) | Capacity inadequacy at node \(n\) in year \(y\) |
14. Implementation Summary¶
Function Call Graph¶
create_master_problem(input)
|-- calculate_target_ratios(input) # RE-1
|-- build_master_variables!(model, input) # All variables (Sec. 3)
|-- add_investment_constraints!() # INV-1 to INV-3
|-- add_budget_constraints!() # BUD-1
|-- add_retirement_cascade_constraints!() # No-op (RET handled in CUM)
|-- add_capacity_adequacy_constraints!() # CAP-1
|-- add_transmission_symmetry_constraints!() # TXN-1
|-- add_representative_days_validation!() # Sec. 8
| |-- select_representative_days() # Sec. 8.1
| |-- create_day_operational_vars!() # Sec. 8.2
| |-- add_day_operational_constraints!() # Sec. 8.3
| +-- calculate_day_operational_cost() # OBJ-5
|-- add_re_target_constraints!() # RE-2
|-- add_re_increment_constraints!() # RE-3
+-- build_master_objective!() # OBJ-1
solve_with_npv_iteration(input)
|-- create_master_problem()
|-- optimize!()
+-- loop:
|-- get_units_with_negative_npv() # Sec. 9.2
|-- force_unit_retirements!()
+-- optimize!()
15. MGA/SPORES: Near-Optimal Alternative Exploration¶
Julia file: mga.jl
MGA (Modeling to Generate Alternatives) [8] and SPORES (Spatially-explicit Practically Optimal REsultS, Lombardi et al. 2020 [7]) explore the space of near-optimal investment plans. In energy planning, multiple investment configurations can have similar total costs but very different compositions (e.g., more wind vs. more solar, distributed vs. centralized storage). These diverse alternatives inform policy decisions.
ESFEX implements both methods as distinct paths through the same machinery:
- MGA (\(\S 15.2\)–\(\S 15.8\)): the classical Hop-Skip-Jump (HSJ) loop. One diversity objective is applied \(K\) times, each iteration penalising investment variables seen in the previous solutions.
- SPORES (\(\S 15.9\)–\(\S 15.15\)): a menu of distinct objectives. Each declared objective produces one alternative under the shared cost-slack constraint; the alternative count equals
len(objectives).
Both share the near-optimal constraint \(Z \leq (1+\varepsilon) C^*\) and the per-alternative cost recovery; only the objective family differs.
15.1 Algorithm¶
1. Solve cost-optimal Master Problem → C*, x*₀
2. Add near-optimal constraint: total_cost ≤ (1 + ε) × C*
3. For k = 1, ..., K:
a. Compute frequency scores from {x*₀, ..., x*_{k-1}}
b. Set diversity-maximizing LP objective
c. Solve → x*_k
Key property: The diversity objective is LP-compatible -- no binary variables are introduced.
15.2 Near-Optimal Constraint¶
After obtaining the cost-optimal solution with cost \(C^*\), a single constraint is added:
where \(\varepsilon\) is the slack_fraction (e.g., 0.05 for 5% cost increase).
15.3 Frequency Scoring¶
Julia function: compute_frequency_scores(alternatives, input; investment_threshold=0.1)
For each investment variable indexed by technology \(g\), node \(n\), and year \(y\), the frequency of investment across all previous alternatives is:
where \(\tau\) is the investment_threshold (MW) and \(K_{prev}\) is the number of alternatives found so far. The diversity score is:
| Score range | Meaning | Effect |
|---|---|---|
| \(\sigma \approx +1\) | Rarely invested | Encourage |
| \(\sigma \approx 0\) | Invested ~50% of the time | Neutral |
| \(\sigma \approx -1\) | Always invested | Discourage |
Variables that have never appeared receive a default score of +1 (maximum encouragement).
15.4 Diversity Objective¶
Julia function: set_spores_objective!(model, vars, input, frequency_scores)
The diversity-maximizing objective replaces the cost objective:
Each investment variable is weighted by its diversity score and normalized by maximum investment capacity \(\bar{I}\). This normalization ensures technologies with different scales (e.g., 10 MW solar vs. 1000 MW wind) are treated comparably.
15.5 Cost Evaluation¶
After solving with the diversity objective, the actual system cost is recovered by evaluating the original cost expression:
Note: objective_value(model) returns the diversity objective, not the cost. The cost is obtained via value(total_cost_expr) where total_cost_expr was saved before modifying the objective.
15.6 Computational Cost¶
For \(K\) alternatives:
- Solves: \(K + 1\) (1 cost-optimal + \(K\) diversity)
- Time: Approximately \((K+1) \times T_{\text{master}}\)
- Memory: The JuMP model is reused; only the objective changes and one constraint is added
- LP size: Unchanged from the base Master Problem
15.7 Function Call Graph¶
run_mga_spores(input)
|-- create_master_problem(input) # Build base model
|-- optimize!(model) # Step 0: cost-optimal
|-- objective_function(model) # Save cost expression
|-- extract_master_solution(...) # Alternative 0
|-- @constraint(cost ≤ (1+ε)×C*) # Near-optimal bound
+-- loop k = 1..K:
|-- compute_frequency_scores() # MGA-2, MGA-3
|-- set_spores_objective!() # MGA-4 (replaces obj)
|-- optimize!(model) # Diversity solve
|-- value(total_cost_expr) # MGA-5 (actual cost)
+-- extract_master_solution(...) # Alternative k
15.8 Configuration¶
| Parameter | Default | Description |
|---|---|---|
mga.enabled |
false |
Enable MGA / SPORES |
mga.method |
"mga" |
Generation method: "mga" (HSJ loop, \(\S 15.2\)–\(\S 15.7\)) or "spores" (per-objective sweep, \(\S 15.9\)–\(\S 15.15\)) |
mga.num_alternatives |
10 | Number of diversity alternatives \(K\). Used only when method = "mga" — ignored under SPORES (the count equals len(objectives)) |
mga.slack_fraction |
0.05 | Near-optimal slack \(\varepsilon\) — shared by both methods |
mga.investment_threshold |
0.1 MW | Threshold \(\tau\) for frequency counting. Used by HSJ and by the SPORES hsj_diversity objective; ignored by the other SPORES objectives |
mga.objectives |
[] |
SPORES objective menu (list of SporesObjective). Required when method = "spores" — ignored under MGA |
15. (continued) SPORES: per-objective alternative generation¶
15.9 Overview¶
Julia entry point: run_spores(input; objectives, slack_fraction, …)
SPORES replaces the HSJ frequency loop with a menu of distinct LP objectives. Each entry in objectives produces one alternative under the same cost-slack constraint \(Z \leq (1+\varepsilon) C^*\). ESFEX ships five canonical objectives, the four "classical" SPORES objectives from Lombardi et al. (2020) plus the HSJ score retained as a special case:
| Symbol | Sense | Section |
|---|---|---|
:hsj_diversity |
\(\max\) | § 15.4 (same as MGA — reusable inside a SPORES sweep) |
:min_total_build |
\(\min\) | § 15.10 |
:max_tech_equity |
\(\min\) (min-max) | § 15.11 |
:max_regional_equity |
\(\min\) (min-max) | § 15.12 |
:evolutionary_dist |
\(\max\) | § 15.13 |
All objectives are linear; auxiliary variables introduced by the min-max equity objectives and by the L1-linearised distance objective are tracked in model[:_spores_objective_aux] and deleted between sweep iterations (see § 15.14) so the JuMP model never accumulates dead variables across the sweep.
15.10 Minimum total build objective¶
Julia function: set_min_build_objective!(model, vars, input)
Selects the smallest near-optimal portfolio. Useful when the cost slack admits a deployment-light plan that meets the same demand and renewable-energy targets:
where the investment-period gating from § 16 still applies (the sums only run over the period-start years \(y_{\mathrm{idx}} \in \{1, 1+y_{pp}, 1+2 y_{pp}, \ldots\}\)). Battery energy investments are deliberately excluded from the MW sum because the units (MWh vs MW) do not combine cleanly; the energy variable is implicitly constrained through the per-tech duration limit.
15.11 Technology equity objective¶
Julia function: set_tech_equity_objective!(model, vars, input)
Minimises the largest per-technology share of the build (a min-max equity formulation, sometimes called Gini-min in the SPORES literature):
[ \min \; M \tag{SPORES-2a} ] [ \text{s.t.}\quad \sum_{y, n} \frac{I^{tech}{y,t,n}}{\bar{I} ]}} \;\leq\; M \qquad \forall\, t \in \mathcal{T} \tag{SPORES-2b
Investments are normalised by the per-tech-per-node investment cap \(\bar{I}_{t,n}\) so technologies with different scales (10 MW vs 1000 MW) compete on the same footing. The auxiliary scalar \(M \geq 0\) and the \(|\mathcal{T}|\) per-tech constraints are anonymous (no base_name) and stashed in the aux registry. Technologies whose investment cap is zero everywhere contribute an empty sum that reduces to the trivial constraint \(0 \leq M\).
15.12 Regional equity objective¶
Julia function: set_regional_equity_objective!(model, vars, input)
The spatially-explicit twin of § 15.11 — minimises the largest per-node share of the build. This is the canonical SPORES objective from which the name spatially-explicit derives:
[ \min \; M \tag{SPORES-3a} ] [ \text{s.t.}\quad \sum_{y, t} \frac{I^{tech}{y,t,n}}{\bar{I} \;+\; \sum_{y, b} \frac{I^{bat,P}}{y,b,n}}{\bar{I}^{P} \;\leq\; M \qquad \forall\, n \in \mathcal{N} \tag{SPORES-3b} ]}
Battery power investments are included because they materially affect spatial siting; battery energy is omitted for the same unit-mixing reason as in (SPORES-1).
15.13 Evolutionary distance objective¶
Julia function: set_evolutionary_distance_objective!(model, vars, input, reference_solution)
Maximises the L1 distance (normalised) from a reference solution \(x_{\mathrm{ref}}\) — typically the cost-optimal plan \(x_0^*\). Used to surface the maximally different feasible plan when several visually similar near-optima exist:
The L1 norm is linearised via positive / negative deviation auxiliaries:
so the objective becomes \(\max \sum (d^+ + d^-)/\bar{I}\). The Euclidean L2 distance would force the model out of LP into QP and is therefore not used.
15.14 Dispatcher and aux-variable cleanup¶
Julia function: apply_spores_objective!(model, vars, input, objective::Symbol; frequency_scores = nothing, reference_solution = nothing)
Routes a SPORES objective symbol to the matching set_*_objective! function and validates required kwargs:
:hsj_diversityrequiresfrequency_scores(aDict{String, Float64}).:evolutionary_distrequiresreference_solution(aMasterProblemResult).- The other three need neither.
Each call begins with _clear_spores_aux!(model) which deletes every variable / constraint installed by the previous objective. Aux references live in model[:_spores_objective_aux] as a Vector{Any}; JuMP.delete is wrapped in try / catch to stay idempotent across repeated cleanups. This is what allows the sweep loop (§ 15.15) to re-use a single model across all objectives without unbounded growth.
15.15 SPORES sweep loop¶
Julia function: run_spores(input; objectives, slack_fraction, use_representative_days, investment_threshold)
run_spores(input, objectives = [:min_total_build, :max_tech_equity, …])
|-- create_master_problem(input) # Build base model
|-- optimize!(model) # Step 0: cost-optimal
|-- objective_function(model) # Save cost expression
|-- extract_master_solution(...) # Reference / alt 0
|-- @constraint(cost ≤ (1+ε)×C*) # Near-optimal envelope
+-- for k = 1..len(objectives):
|-- _clear_spores_aux!(model) # Drop previous aux
|-- apply_spores_objective!(...) # Install objective k
|-- optimize!(model) # Solve under cost cap
|-- value(total_cost_expr) # MGA-5 (actual cost)
+-- extract_master_solution(...) # Alternative k
The result is an MGAResult whose alternatives[1] is the cost-optimal seed, alternatives[2:end] are the SPORES solutions, and objective_labels[k] carries the symbol that produced alternatives[k+1].
15.16 When to use MGA vs SPORES¶
| Use case | Use MGA | Use SPORES |
|---|---|---|
| "Show me K alternatives, surprise me" | ✅ | — |
| "Show me the smallest near-optimal plan" | — | ✅ (min_total_build) |
| "How spatially flexible is the optimum?" | — | ✅ (max_regional_equity) |
| "Quantify the technology-substitution envelope" | indirect | ✅ (max_tech_equity) |
| Need named alternatives for policy discussion | — | ✅ |
| Need a large set (K = 20+) for statistical robustness | ✅ | ❌ (limited by objective count) |
| Cost cap is the only constraint | ✅ | ✅ |
In practice, SPORES is preferable when each alternative needs to answer a specific question; MGA is preferable when the goal is to map the breadth of the near-optimal space with a larger sample.
Solver Configuration¶
The Master Problem uses HiGHS with the following settings:
| Parameter | Default | Julia field |
|---|---|---|
| Threads | 4 | input.threads |
| Time limit | 3600 s | input.time_limit |
| MIP gap | 0.01 | input.gap |
| Verbose | false | input.verbose |
Since the formulation is a pure LP (all continuous variables, age-based retirement evaluated at construction time), the MIP gap setting is effectively unused. The solver employs the simplex or interior point method for LP.
16. Per-Technology Investment Model¶
Investment decisions can be modeled at the technology level rather than per-generator. This aggregates investment capacity across all generators of the same type, enabling more realistic technology-level constraints.
16.1 Technology Configuration¶
Technologies are defined in the YAML configuration:
technologies:
- name: "Solar PV"
type: "Renewable"
generators: ["Solar_Farm_1", "Solar_Farm_2"]
invest_min_mw: 0.0
invest_max_mw: 500.0
investment_cost_per_mw: 800000.0
fixed_om_per_mw: 12000.0
lifetime_years: 25
battery_technologies:
- name: "Li-Ion Storage"
batteries: ["Battery_1", "Battery_2"]
invest_min_power_mw: 0.0
invest_max_power_mw: 200.0
invest_max_energy_mwh: 800.0
power_cost_per_mw: 150000.0
energy_cost_per_mwh: 200000.0
lifetime_years: 15
16.2 Technology-Level Variables¶
| Variable | Domain | Units | Description |
|---|---|---|---|
| \(I^{tech}_{y,k}\) | \(\mathbb{R}_+\) | MW | Total technology \(k\) investment in year \(y\) |
| \(I^{bat,tech,P}_{y,k}\) | \(\mathbb{R}_+\) | MW | Battery technology \(k\) power investment |
| \(I^{bat,tech,E}_{y,k}\) | \(\mathbb{R}_+\) | MWh | Battery technology \(k\) energy investment |
16.3 Cumulative Capacity¶
The cumulative installed capacity for technology \(k\) at year \(y\) includes both existing capacity and all investments up to year \(y\), minus retirements:
where \(R^{tech}_{y',k}\) is retirement computed by age-based expiry.
16.4 Virtual Generators¶
Technology-level investments create virtual generators for operational dispatch. After the master problem solves, the runner:
- Computes cumulative investment per technology per year
- Creates synthetic generator entries with
rated_power = cumulative_investment - Passes these virtual generators to the operational dispatch adapter
- Virtual generators use the same availability profiles and cost parameters as the technology definition
This bridge between strategic (technology-level) and operational (generator-level) modeling ensures investment decisions are reflected in detailed dispatch simulation.
17. Piecewise-Linear Cost Curves in Investment¶
Investment costs can vary with scale using piecewise-linear (PWL) cost curves. This captures economies of scale (decreasing marginal cost for larger projects) or diseconomies (increasing cost for difficult sites).
17.1 Cost Curve Structure¶
Each technology can define cost curve blocks:
technologies:
- name: "Wind Onshore"
cost_curve:
blocks:
- power_mw: 100.0
cost_per_mw: 1200000.0 # First 100 MW at $1.2M/MW
- power_mw: 200.0
cost_per_mw: 1000000.0 # Next 200 MW at $1.0M/MW
- power_mw: 200.0
cost_per_mw: 900000.0 # Next 200 MW at $0.9M/MW
17.2 Formulation¶
The total investment cost for technology \(k\) in year \(y\) is:
where: - \(c^{inv}_{k,s}\) is the cost per MW in segment \(s\) - \(\delta^{inv}_{k,y,s}\) is the investment in segment \(s\) (bounded by segment width)
The total investment equals the sum of segments:
When costs are non-increasing (\(c^{inv}_{k,1} \geq c^{inv}_{k,2} \geq \ldots\)), the formulation is convex and the LP solver naturally fills cheaper segments first without requiring binary variables.
References¶
The capacity expansion formulation draws on the general framework for generation expansion planning reviewed by Koltsaklis and Dagoumas [24]. The importance of embedding operational flexibility constraints within planning models is analyzed by Palmintier and Webster [25] and Poncelet et al. [26]. The interaction between unit commitment constraints and generation expansion is studied by Schwele et al. [35], which motivates ESFEX's two-stage decomposition approach. The role of operational detail in planning models with high variable RE shares is reviewed by Helistö et al. [36]. Representative day selection and time series aggregation methods follow Kotzur et al. [22] and Nahmmacher et al. [37]; for a broader review see Hoffmann et al. [23]. The MGA/SPORES methodology for near-optimal space exploration follows Lombardi et al. [7] and DeCarolis [8]; see also Neumann and Brown [38]. The optimization is formulated via JuMP [20] and solved with HiGHS [21] by default.
See the full bibliography for complete citation details.