BCNexus (CLEWs) Model Code

Overview

The BCNexus model is a CLEWs (Climate–Land–Energy–Water systems) implementation built on the OSeMOSYS (Open Source energy MOdelling SYStem) optimisation framework. It is a linear / mixed-integer cost-minimisation model written in GNU MathProg (GLPK's algebraic modelling language) and solved with glpsol. The objective is to find the least-cost portfolio of capacity investments and operating decisions that satisfies exogenous energy-service demand over a multi-year horizon, subject to physical, policy, and (in the CLEWs extension) land-use constraints.

Two model files are documented here. They are identical except for how intra-year storage chronology is represented; both carry the same CLEWs land-use extension:

The two files map directly to the models/model_Niet and models/model_Kotzur folders of the DeltaE/BC_Nexus repository. Each ships with an otoole YAML configuration that defines the CSV schema for every set, parameter, variable, and result.

OSeMOSYS vs CLEWs — what actually differs

OSeMOSYS at its core is a generic reference energy system optimiser: technologies convert fuels into other fuels and final demands, with costs, capacities, and emissions accounted across years and timeslices. CLEWs reuses this identical engine and the same mathematical primitives. The distinguishing feature in this codebase is a compact land-use accounting layer bolted onto the standard formulation.

Conceptually, the land system is represented through the energy-modelling abstraction itself: a parcel of land is a "technology", and the way that land is used (a crop, forest, pasture, conversion process) is a "mode of operation". By constraining activity per mode you constrain land allocation, land-use transitions, and their associated water and emission flows — without adding a separate land sub-model. This is why the user note is accurate: "the OSeMOSYS and CLEWs differences are simply in some land-use accounting inclusion in CLEWs." The land-use block is documented in full in Section 14.

Kotzur vs Niet — the storage chronology difference

Both files implement storage with the same investment, salvage, and limit logic. They differ only in the temporal bookkeeping of the state of charge, which matters when typical-day clustering is used to compress the timeslice structure. The relevant declarations and constraints differ as follows.

Practical implication: the Kotzur formulation lets a small set of clustered "typical days" still capture seasonal (inter-day) storage cycling — important for long-duration storage — at far lower variable count than tracking every timeslice chronologically across a full year.

Anatomy of the model file

A MathProg OSeMOSYS file is read top-to-bottom in five phases. Understanding this order makes the block-by-block documentation below easy to navigate.

# 1. SETS        — index domains (years, technologies, fuels, timeslices, ...)
# 2. PARAMETERS  — exogenous data, indexed over sets (costs, demands, efficiencies)
# 3. VARIABLES   — decision & accounting unknowns the solver chooses
# 4. minimize cost: ...   — the single objective function
# 5. s.t. <LABEL>: ...    — constraint blocks
# then:  solve;  +  table tout ... OUT "CSV" ...   — emit results

Lines beginning # are comments; large numbers of constraints in these files are commented out (#s.t. ...) because the BC team pre-computes some accounting in Python (preprocess_data.py / otoole) and keeps only the constraints needed for the solve. Where a block is mostly disabled, this is noted.

What is actually active in the base model used vs available

OSeMOSYS ships with far more sets, parameters and constraints than any one study uses. The BC_Nexus wiki's Model Structure page records which are wired up in the reference (base) configuration — this is why so many constraints appear commented out. The counts below are used / total.

1 · Sets index domains

Sets are the index ranges over which everything else is defined. The core OSeMOSYS sets describe the spatial, temporal, and technological dimensions of the reference energy system.

set YEAR;              # modelling horizon (e.g. 2020..2050)
set TECHNOLOGY;        # every converter incl. land parcels (CLEWs)
set TIMESLICE;         # intra-year time resolution (season × daytype × bracket)
set FUEL;              # energy carriers, plus water/material flows in CLEWs
set EMISSION;          # CO2, CH4, ... tracked pollutants
set MODE_OF_OPERATION; # operating modes; in CLEWs encode land/crop categories
set REGION;            # spatial nodes (BC sub-regions)
set SEASON; set DAYTYPE; set DAILYTIMEBRACKET;
set FLEXIBLEDEMANDTYPE;
set STORAGE;           # storage facilities
set DAYSCRO;           # KOTZUR ONLY: chronological order of typical days

# Derived membership sets that pre-filter mode×technology pairs by the fuel
# they produce/consume or the storage they charge/discharge — these make the
# energy-balance summations far cheaper:
set MODExTECHNOLOGYperFUELout{FUEL} within MODE_OF_OPERATION cross TECHNOLOGY;
set MODExTECHNOLOGYperFUELin{FUEL}  within MODE_OF_OPERATION cross TECHNOLOGY;
set MODExTECHNOLOGYperSTORAGEto{STORAGE}   ...;
set MODExTECHNOLOGYperSTORAGEfrom{STORAGE} ...;

2 · Parameters exogenous data

Parameters hold all input data, grouped by the source comment banners in the file. They are populated from CSVs by otoole before the solve.

Global & demand

Performance & cost

Land-use parameters CLEWs

#### Land use ####
param TechnologyActivityByModeUpperLimit{r in REGION, t in TECHNOLOGY, m in MODE_OF_OPERATION, y in YEAR};
param TechnologyActivityByModeLowerLimit{r in REGION, t in TECHNOLOGY, m in MODE_OF_OPERATION, y in YEAR};
param TechnologyActivityIncreaseByModeLimit{...};  # max year-on-year growth (fraction)
param TechnologyActivityDecreaseByModeLimit{...};  # max year-on-year decline (fraction)

These four parameters supply the right-hand-side bounds consumed by the LU1–LU4 constraints.

3 · Variables decision unknowns

Variables are what the solver chooses. They split into primary decisions (investment and operation) and accounting variables that merely tally derived quantities for reporting and for use inside constraints.

Many classic OSeMOSYS variables (RateOfProduction, Use, ProductionByTechnology, TotalDiscountedCost, …) are commented out here: the BC workflow pre-computes or omits them to shrink the problem.

4 · Objective function minimize cost

A single statement minimises the total discounted system cost over all regions, technologies, and years. Each term is discounted to the first model year.

minimize cost:
  sum{r in REGION, t in TECHNOLOGY, y in YEAR} (
      ( /* fixed O&M on installed capacity */
        ((sum{yy: y-yy < OperationalLife && y-yy>=0} NewCapacity[r,t,yy]) + ResidualCapacity[r,t,y]) * FixedCost[r,t,y]
        + /* variable O&M on activity */
        sum{m in MODEperTECHNOLOGY[t], l in TIMESLICE} RateOfActivity[r,l,t,m,y]*YearSplit[l,y]*VariableCost[r,t,m,y]
      ) / (1+DiscountRate[r])^(y-y0+0.5)                       # mid-year discounting of O&M
      + CapitalCost[r,t,y] * NewCapacity[r,t,y] / (1+DiscountRate[r])^(y-y0)   # investment
      + DiscountedTechnologyEmissionsPenalty[r,t,y]            # carbon penalty
      - DiscountedSalvageValue[r,t,y]                          # credit for residual life
  )
  + sum{s in STORAGE} (
      CapitalCostStorage[r,s,y]*NewStorageCapacity[r,s,y]/(1+DiscountRate[r])^(y-y0)
      - DiscountedSalvageValueStorage[r,s,y]
  );

Reading it term by term: (1) fixed O&M scales with total installed capacity (new + residual); (2) variable O&M scales with operating activity; (3) capital cost is paid on new capacity in its build year; (4) emission penalties add the carbon price; (5) salvage value credits capacity whose operational life extends beyond the horizon; (6) the storage block adds storage capital cost net of storage salvage. y0 = min(YEAR); the +0.5 on O&M represents mid-year cash flow.

5 · Capacity adequacy A & B capacity ↔ activity

These constraints connect installed capacity to permissible activity, ensuring the system is physically able to deliver what it operates.

s.t. CAa4_Constraint_Capacity{r,l,t,y}:
   sum{m in MODEperTECHNOLOGY[t]} RateOfActivity[r,l,t,m,y]
   <= ((sum{yy: y-yy<OperationalLife && y-yy>=0} NewCapacity[r,t,yy]) + ResidualCapacity[r,t,y])
      * CapacityFactor[r,t,l,y] * CapacityToActivityUnit[r,t];

6 · Energy balance A & B supply ≥ demand

The heart of the reference energy system: for every fuel, production must cover demand, downstream use, and net trade — at each timeslice (A) and over the year (B). Most fine-grained accounting equations are disabled; the two binding inequalities below are what actually enforce balance.

# Per-timeslice balance: production >= demand + use + exports
s.t. EBa11_EnergyBalanceEachTS5{r,l,f,y}:
   sum{(m,t) in MODExTECHNOLOGYperFUELout[f]} RateOfActivity[r,l,t,m,y]*OutputActivityRatio*YearSplit[l,y]
   >= SpecifiedAnnualDemand*SpecifiedDemandProfile
      + sum{(m,t) in MODExTECHNOLOGYperFUELin[f]} RateOfActivity*InputActivityRatio*YearSplit
      + sum{rr in REGION} Trade[r,rr,l,f,y]*TradeRoute[r,rr,f,y];

# Annual balance for accumulated (non-time-sliced) demand
s.t. EBb4_EnergyBalanceEachYear4{r,f,y}: production_annual >= use_annual + trade_annual + AccumulatedAnnualDemand;

Supporting equalities: EBa9 defines time-sliced Demand; EBa10 enforces trade antisymmetry (Trade[r,rr]=-Trade[rr,r]).

7 · Accounting: technology production / use bookkeeping

This short block tallies derived quantities. In these files most of it is commented out; the one active and important equation is Acc3, which feeds the land-use block.

s.t. Acc3_AverageAnnualRateOfActivity{r,t,m,y}:
   sum{l in TIMESLICE} RateOfActivity[r,l,t,m,y]*YearSplit[l,y]
   = TotalAnnualTechnologyActivityByMode[r,t,m,y];   # ← the per-mode annual activity used by LU1–LU4

8 · Storage equations, constraints & investment state of charge

The storage block tracks the state of charge, keeps it within capacity, prevents free energy creation, and handles storage investment and salvage. This is the only block that differs between the Niet and Kotzur files (see variants).

Niet (standard, timeslice chronology)

s.t. S1_StorageLevelYearStart{r,s,y}: level carried from previous year;
s.t. S2_StorageLevelTSStart{r,s,l,y}: # chains charge across timeslices within the year
   if l = first timeslice then StorageLevelYearStart else previous level + net charge;
s.t. SC1/SC2_Lower/UpperLimit: MinCharge*cap <= StorageLevelTSStart <= installed storage capacity;
s.t. SC8_StorageRefilling: net annual charge over all timeslices = 0;

Kotzur (typical-day chronology)

s.t. S2_StorageLevelDayTypeFinish{r,s,ld,y}: net charge accumulated within each day-type;
s.t. S1_StorageLevelYearStart{r,s,y}: carries the summed day-type finishes across years;
s.t. S3_StorageLevelChronoDayStart{r,s,ldc,y}: # chains day-types in chronological order DAYSCRO
   if ldc = first day then StorageLevelYearStart else prev chrono-day + that day's finish;
s.t. SC1/SC2: limits applied to StorageLevelChronoDayStart;
s.t. SC4_StorageRefilling: net charge over all day-types & years = 0;

Shared investment / salvage logic (SI6–SI9): storage gets zero salvage if its life ends within the horizon, otherwise a depreciation-method-dependent salvage value discounted to the start year.

9 · Capital, salvage & operating costs cost accounting

These blocks define the cost variables that the objective consumes. Several are commented out because the objective computes them inline; the active ones supply reported cost streams.

The "Total Discounted Costs" block (TDC1/TDC2) is entirely commented out — the objective already aggregates discounted costs, so the per-technology breakdown is skipped for speed.

10 · Capacity & activity limits policy bounds

Optional upper/lower bounds the modeller imposes on build-out and operation. A sentinel value of -1 means "no limit", so each constraint carries a guard condition (e.g. : TotalAnnualMaxCapacity[r,t,y] <> -1).

11 · Reserve margin reliability

Ensures enough tagged firm capacity exists to cover demand plus a reserve. Only the binding form RM3 is active; the accounting variants RM1/RM2 are disabled.

s.t. RM3_ReserveMargin_Constraint{r,l,y}:
   (tagged production) * ReserveMarginTagFuel * ReserveMargin
   <= sum{t} TotalCapacity * ReserveMarginTagTechnology[r,t,y] * CapacityToActivityUnit[r,t];

12 · Renewable-energy production target policy

Forces a minimum share of production from technologies/fuels tagged renewable. Active constraint: RE4_EnergyConstraint.

s.t. RE4_EnergyConstraint{r,y}:
   REMinProductionTarget[r,y] * (production of RE-tagged fuels)
   <= (production from RE-tagged technologies, year-weighted);

Companion accounting equations RE1–RE3, RE5 are commented out.

13 · Emissions accounting climate

Translates activity into emissions, prices them, and enforces annual and model-period caps.

In the BC base model the tracked pollutant is CO₂-equivalent, with emission factors sourced from the US EPA (2018) and BC's carbon-intensity records. The carbon price fed through EmissionsPenalty follows BC's legislated tax trajectory — roughly $45/t (2020) rising to $170/t (2030+) — which is also what makes the negative-cost forest CCS credit (see land-use) economically meaningful.

14 · Land-use constraints ★ the CLEWs extension

This is the block that makes BCNexus a CLEWs model rather than plain OSeMOSYS. Four constraints bound the per-mode annual activity (TotalAnnualTechnologyActivityByMode, defined by Acc3) in both level and rate of change. Because a mode encodes a land/crop/cover category, these are effectively land-allocation and land-transition limits.

######## Land use constraints ########

# LU1 — upper bound on land-use activity per mode (skipped where limit = 0)
s.t. LU1_TechnologyActivityByModeUL{r,t,m,y: TechnologyActivityByModeUpperLimit[r,t,m,y] <> 0}:
   TotalAnnualTechnologyActivityByMode[r,t,m,y] <= TechnologyActivityByModeUpperLimit[r,t,m,y];

# LU2 — lower bound on land-use activity per mode
s.t. LU2_TechnologyActivityByModeLL{r,t,m,y}:
   TotalAnnualTechnologyActivityByMode[r,t,m,y] >= TechnologyActivityByModeLowerLimit[r,t,m,y];

# LU3 — cap year-on-year INCREASE (e.g. how fast cropland can expand)
s.t. LU3_TechnologyActivityIncreaseByMode{r,t,m,y,yy: y-yy==1 && limit<>0}:
   activity[y] <= (1 + TechnologyActivityIncreaseByModeLimit[r,t,m,yy]) * activity[yy];

# LU4 — cap year-on-year DECREASE (e.g. how fast forest can be lost)
s.t. LU4_TechnologyActivityDecreaseByMode{r,t,m,y,yy: y-yy==1 && limit<>0}:
   activity[y] >= (1 - TechnologyActivityDecreaseByModeLimit[r,t,m,yy]) * activity[yy];

How BC's land is set up (base model)

The land data behind these modes comes from the GAEZ (Global Agro-Ecological Zoning) model using agglomerative hierarchical clustering of land with similar achievable yield. BC is divided into 7 clustered zones, and 9 crops covering ~90% of provincial production (alfalfa, barley, maize, oat, pea, potato, rapeseed, rye, wheat) are modelled explicitly, with the remainder grouped as "other". Land is measured in units of 1000 km² (BC ≈ 925 units), farmland is given a generic 15-year operational life, and each crop appears in mode/technology combinations of irrigated vs rain-fed × low/intermediate/high input intensity — visible in the commodity codes (e.g. LWHEIIBC1 = wheat, irrigated, intermediate).

15 · Solve & output tables results

After the constraints, a single solve; invokes GLPK. A series of table tout … OUT "CSV" … statements then writes each result variable to a CSV under res/csv/, which the otoole / BCNexus post-processing reads back. Exported results include capacity (NewCapacity, TotalCapacityAnnual), operation (RateOfActivity, ProductionByTechnologyAnnual), costs, emissions, trade, storage levels, and — for CLEWs — TotalAnnualTechnologyActivityByMode.

solve;
table tout {r,t,y} OUT "CSV" "…/NewCapacity.csv" : r, t, y, NewCapacity[r,t,y];
table tout {r,t,m,y} OUT "CSV" "…/TotalAnnualTechnologyActivityByMode.csv" : r, t, m, y, ...;  # land use output
... (≈30 result tables) ...
end;

References

• BCNexus source — github.com/DeltaE/BC_Nexus (models: model_Kotzur, model_Niet)
• OSeMOSYS documentation — osemosys.readthedocs.io
• OSeMOSYS reference implementation: Howells et al. (2011), Energy Policy 39(10):5850–5870, "OSeMOSYS: The Open Source Energy Modeling System."
• Kotzur et al. (2018) typical-period storage formulation — basis of the seasonal-storage variant.
• Storage with temporal aggregation — DeltaE/storage-in-OSeMOSYS (Cluster · Kotzur · Welsch · Niet methods; Borba, Islam & Niet, 2025).
• Part of the Combined Modelling BC suite — BCNexus + RESource + PyPSA-BC with a soft-linking tool.
• Source files documented: model_BCNexus_Kotzur.m (533 lines), model_BCNexus_Niet.m (522 lines), with their otoole_config_*.yaml.