1 · Introduction

This document presents a linear programming (LP) model to evaluate whether Australia can achieve net-zero greenhouse gas emissions from its electricity sector by 2050. The model is built using the ompr algebraic modelling interface and solved with the open-source GLPK solver.

1.1 Modelling Philosophy

The model is deliberately kept transparent and modular:

  • All parameters are declared in a single section and are easy to swap for scenario analysis.
  • The mathematical structure closely mirrors the problem description, making the code straightforward to audit.
  • Outputs are presented as tidy data frames, feeding directly into ggplot2 visualisations.

1.2 Energy Sources Modelled

Source Type Notes
Coal Fossil No new builds; existing fleet depreciates
Gas Fossil Flexible backup; phased down
Solar Renewable High build potential
Wind Renewable High build potential
Hydro Renewable Limited expansion potential
Nuclear Low-carbon Long lead times; included as option

2 · Dependencies

library(ompr)            # Algebraic LP/MIP modelling interface
library(ompr.roi)        # Bridge between ompr and ROI solvers
library(ROI.plugin.glpk) # GLPK solver plugin
library(ggplot2)         # Visualisation
library(dplyr)           # Data wrangling
library(tidyr)           # Data reshaping
library(scales)          # Axis formatting helpers
library(knitr)           # Table rendering
library(kableExtra)      # Enhanced table styling

3 · Sets & Indices

# Energy source labels
sources <- c("coal", "gas", "solar", "wind", "hydro", "nuclear")

# Planning horizon
years <- 2025:2050

# Integer indices used inside ompr expressions
n_sources <- length(sources)   # I = 6
n_years   <- length(years)     # T = 26
src_idx   <- seq_len(n_sources)
yr_idx    <- seq_len(n_years)

cat("Sources:", paste(sources, collapse = ", "), "\n")
#> Sources: coal, gas, solar, wind, hydro, nuclear
cat("Years:  ", min(years), "–", max(years), "(", n_years, "periods )\n")
#> Years:   2025 – 2050 ( 26 periods )

4 · Parameters

Scenario analysis: All parameters live in this section. To run an alternative scenario, change the values here and re-knit the document. The model structure remains unchanged.

4.1 Demand \(D_t\) (TWh / year)

# Demand grows linearly from 300 TWh (2025) to ~425 TWh (2050)
# reflecting electrification of transport, industry, and heating.
D <- 300 + 5 * (yr_idx - 1)

data.frame(year = years, demand_TWh = D) |>
  filter(year %in% seq(2025, 2050, 5)) |>
  kable(caption = "Demand projection (selected years)") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Demand projection (selected years)
year demand_TWh
2025 300
2030 325
2035 350
2040 375
2045 400
2050 425

4.2 Cost Parameters

# Variable (production) cost  c[i]  — $M per TWh
c_prod <- c(coal = 40, gas = 60, solar = 15, wind = 20, hydro = 10, nuclear = 35)

# Capital (build) cost  k[i]  — $M per GW installed
k_cap  <- c(coal = 3000, gas = 1000, solar = 1200, wind = 1500, hydro = 4000, nuclear = 8000)

data.frame(
  Source          = sources,
  `Production cost ($/MWh)`  = c_prod,
  `Capital cost ($M/GW)`     = k_cap,
  check.names = FALSE
) |>
  kable(caption = "Cost parameters") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Cost parameters
Source Production cost (\(/MWh) </th> <th style="text-align:right;"> Capital cost (\)M/GW)
coal coal 40 3000
gas gas 60 1000
solar solar 15 1200
wind wind 20 1500
hydro hydro 10 4000
nuclear nuclear 35 8000

4.3 Emission Intensities \(e_i\) (Mt CO₂ / TWh)

# Renewables and nuclear treated as zero lifecycle emissions (simplified)
e_int <- c(coal = 0.88, gas = 0.45, solar = 0.00, wind = 0.00, hydro = 0.00, nuclear = 0.00)

data.frame(Source = sources, `Emission intensity (Mt CO2/TWh)` = e_int,
           check.names = FALSE) |>
  kable(caption = "Emission intensities") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Emission intensities
Source Emission intensity (Mt CO2/TWh)
coal coal 0.88
gas gas 0.45
solar solar 0.00
wind wind 0.00
hydro hydro 0.00
nuclear nuclear 0.00

4.4 Technical Parameters

# Annual depreciation rate  delta[i]  (fraction of capacity retired per year)
delta <- c(coal = 0.02, gas = 0.03, solar = 0.04, wind = 0.04, hydro = 0.01, nuclear = 0.02)

# Maximum new capacity built per year  B_max[i]  (GW/yr) — supply-chain constraint
B_max <- c(coal = 0.0, gas = 2.0, solar = 10.0, wind = 8.0, hydro = 0.5, nuclear = 1.0)

# Initial installed capacity  K0[i]  in 2025  (GW)
K0    <- c(coal = 20.0, gas = 12.0, solar = 20.0, wind = 12.0, hydro = 8.0, nuclear = 0.0)

# Capacity factor  cf[i]  — fraction of hours at full output
cf    <- c(coal = 0.70, gas = 0.60, solar = 0.22, wind = 0.35, hydro = 0.45, nuclear = 0.90)

# Maximum annual production per GW of installed capacity  (TWh / GW / year)
TWh_per_GW <- 8.76 * cf

data.frame(
  Source             = sources,
  `K0 (GW)`          = K0,
  `Depreciation`     = delta,
  `Build limit (GW/yr)` = B_max,
  `Capacity factor`  = cf,
  `TWh/GW/yr`        = round(TWh_per_GW, 2),
  check.names = FALSE
) |>
  kable(caption = "Technical parameters") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Technical parameters
Source K0 (GW) Depreciation Build limit (GW/yr) Capacity factor TWh/GW/yr
coal coal 20 0.02 0.0 0.70 6.13
gas gas 12 0.03 2.0 0.60 5.26
solar solar 20 0.04 10.0 0.22 1.93
wind wind 12 0.04 8.0 0.35 3.07
hydro hydro 8 0.01 0.5 0.45 3.94
nuclear nuclear 0 0.02 1.0 0.90 7.88

5 · Mathematical Formulation

5.1 Decision Variables

Variable Domain Description
\(x_{i,t}\) \(\mathbb{R}_+\) Energy produced by source \(i\) in year \(t\) (TWh)
\(K_{i,t}\) \(\mathbb{R}_+\) Installed capacity of source \(i\) in year \(t\) (GW)
\(B_{i,t}\) \(\mathbb{R}_+\) New capacity built for source \(i\) in year \(t\) (GW)

5.2 Objective Function

\[\min \sum_{i} \sum_{t} c_i \cdot x_{i,t} + \sum_{i} \sum_{t} k_i \cdot B_{i,t}\]

5.3 Constraints

Demand — supply must meet or exceed demand in every year: \[\sum_{i} x_{i,t} \geq D_t \quad \forall\, t\]

Capacity — production bounded by available capacity: \[x_{i,t} \leq K_{i,t} \cdot \text{TWh/GW}_i \quad \forall\, i, t\]

Capacity evolution — stock accounting with depreciation: \[K_{i,t+1} = K_{i,t} + B_{i,t} - \delta_i \cdot K_{i,t} \quad \forall\, i,\; t < T\]

Initial capacity: \[K_{i,1} = K^0_i \quad \forall\, i\]

Build limits — supply-chain constraint: \[B_{i,t} \leq B^{\max}_i \quad \forall\, i, t\]

No new coal — policy constraint: \[B_{\text{coal},\, t} = 0 \quad \forall\, t\]

Net-zero target — hard emissions ceiling in 2050: \[\sum_{i} e_i \cdot x_{i,T} \leq 0\]

6 · Build & Solve the Model

# Helper functions to extract scalar parameters inside ompr lambda expressions
# (ompr passes i, t as plain integer indices)
get_c  <- function(i) c_prod[i]
get_k  <- function(i) k_cap[i]
get_e  <- function(i) e_int[i]
get_d  <- function(t) D[t]
get_Bm <- function(i) B_max[i]
get_K0 <- function(i) K0[i]
get_tw <- function(i) TWh_per_GW[i]
get_dl <- function(i) delta[i]

model <- MIPModel() |>

  # ── Decision Variables ──────────────────────────────────────────────────────
  add_variable(x[i, t], i = src_idx, t = yr_idx, type = "continuous", lb = 0) |>
  add_variable(K[i, t], i = src_idx, t = yr_idx, type = "continuous", lb = 0) |>
  add_variable(B[i, t], i = src_idx, t = yr_idx, type = "continuous", lb = 0) |>

  # ── Objective ───────────────────────────────────────────────────────────────
  set_objective(
    sum_expr(get_c(i) * x[i, t], i = src_idx, t = yr_idx) +
    sum_expr(get_k(i) * B[i, t], i = src_idx, t = yr_idx),
    sense = "min"
  ) |>

  # ── C1: Demand ──────────────────────────────────────────────────────────────
  add_constraint(
    sum_expr(x[i, t], i = src_idx) >= get_d(t),
    t = yr_idx
  ) |>

  # ── C2: Capacity bound ──────────────────────────────────────────────────────
  add_constraint(
    x[i, t] <= get_tw(i) * K[i, t],
    i = src_idx, t = yr_idx
  ) |>

  # ── C3a: Initial capacity ───────────────────────────────────────────────────
  add_constraint(K[i, 1] == get_K0(i), i = src_idx) |>

  # ── C3b: Capacity evolution ─────────────────────────────────────────────────
  add_constraint(
    K[i, t + 1] == K[i, t] + B[i, t] - get_dl(i) * K[i, t],
    i = src_idx, t = yr_idx[-n_years]
  ) |>

  # ── C4: Build limits ────────────────────────────────────────────────────────
  add_constraint(B[i, t] <= get_Bm(i), i = src_idx, t = yr_idx) |>

  # ── C5: No new coal ─────────────────────────────────────────────────────────
  add_constraint(B[1, t] == 0, t = yr_idx) |>

  # ── C6: Net-zero in 2050 ────────────────────────────────────────────────────
  add_constraint(
    sum_expr(get_e(i) * x[i, n_years], i = src_idx) <= 0
  )

E0 <- sum(K0 * TWh_per_GW * e_int)

model <- model |>
  add_constraint(
    sum_expr(e_int[i] * x[i, t], i = src_idx) <= 
      E0 * (1 - (t - 1)/(n_years - 1)),
    t = yr_idx
  )
result <- solve_model(model, with_ROI(solver = "glpk", verbose = FALSE))

cat("Solver status      :", result$status, "\n")
#> Solver status      : error
cat("Total cost ($M)    :", format(round(result$objective_value), big.mark = ","), "\n")
#> Total cost ($M)    : 566,037

7 · Results

7.1 Extract Solution

# Extract ALL solutions first
x_sol <- get_solution(result, x[i, t])
K_sol <- get_solution(result, K[i, t])
B_sol <- get_solution(result, B[i, t])

# Map indices back to labels
prod_df <- x_sol |>
  mutate(
    source = sources[i],
    year   = years[t],
    production = value
  ) |>
  select(source, year, production)

cap_df <- K_sol |>
  mutate(
    source = sources[i],
    year   = years[t],
    capacity = value
  ) |>
  select(source, year, capacity)

build_df <- B_sol |>
  mutate(
    source = sources[i],
    year   = years[t],
    build = value
  ) |>
  select(source, year, build)

# Emissions
emissions_df <- prod_df |>
  mutate(ef = e_int[source]) |>
  group_by(year) |>
  summarise(emissions_Mt = sum(production * ef), .groups = "drop")

7.2 Energy Production by Source (TWh)

prod_df |>
  filter(year %in% c(2025, 2030, 2035, 2040, 2045, 2050)) |>
  mutate(production = round(production, 1)) |>
  pivot_wider(names_from = year, values_from = production) |>
  kable(caption = "Optimal energy production by source (TWh)") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Optimal energy production by source (TWh)
source 2025 2030 2035 2040 2045 2050
coal 122.6 110.9 92.9 62.0 31.0 0.0
gas 63.1 25.5 0.0 0.0 0.0 0.0
solar 38.5 120.4 187.1 241.5 285.9 302.8
wind 36.8 68.2 70.0 71.5 83.1 94.5
hydro 31.5 0.0 0.0 0.0 0.0 24.5
nuclear 0.0 0.0 0.0 0.0 0.0 3.2

7.3 Installed Capacity by Source (GW)

cap_df |>
  filter(year %in% c(2025, 2030, 2035, 2040, 2045, 2050)) |>
  mutate(capacity = round(capacity, 1)) |>
  pivot_wider(names_from = year, values_from = capacity) |>
  kable(caption = "Installed capacity by source (GW)") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Installed capacity by source (GW)
source 2025 2030 2035 2040 2045 2050
coal 20 18.1 16.3 14.8 13.4 12.1
gas 12 10.3 8.8 7.6 6.5 5.6
solar 20 62.5 97.1 125.3 148.3 157.1
wind 12 22.3 22.8 23.3 27.1 30.8
hydro 8 7.6 7.2 6.9 6.5 6.2
nuclear 0 0.0 0.0 0.0 0.0 0.4

7.4 Emissions Trajectory (Mt CO₂)

emissions_df |>
  filter(year %in% c(2025, 2030, 2035, 2040, 2045, 2050)) |>
  mutate(emissions_Mt = round(emissions_Mt, 2)) |>
  kable(col.names = c("Year", "Emissions (Mt CO₂)"),
        caption   = "Annual grid emissions (selected years)") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Annual grid emissions (selected years)
Year Emissions (Mt CO₂)
2025 136.31
2030 109.04
2035 81.78
2040 54.52
2045 27.26
2050 0.00

8 · Visualisations

# Consistent colour palette across all plots
source_colours <- c(
  coal    = "#4a4a4a",
  gas     = "#f4a261",
  solar   = "#e9c46a",
  wind    = "#57cc99",
  hydro   = "#4cc9f0",
  nuclear = "#9b5de5"
)

# Shared theme
theme_netzero <- theme_minimal(base_size = 12) +
  theme(
    plot.title      = element_text(face = "bold", size = 13),
    plot.subtitle   = element_text(colour = "grey40", size = 10),
    legend.position = "bottom",
    legend.title    = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

8.1 Energy Mix Over Time

demand_line <- data.frame(year = years, demand = D)

ggplot(prod_df, aes(x = year, y = production, fill = source)) +
  geom_area(alpha = 0.9, colour = "white", linewidth = 0.2) +
  geom_line(data = demand_line, aes(x = year, y = demand),
            colour = "black", linewidth = 0.9, linetype = "dashed",
            inherit.aes = FALSE) +
  annotate("text", x = 2043, y = max(D) * 1.04,
           label = "Demand", size = 3.2, colour = "black") +
  scale_fill_manual(values = source_colours, name = "Source") +
  scale_x_continuous(breaks = seq(2025, 2050, 5)) +
  scale_y_continuous(labels = comma) +
  labs(
    title    = "Australia's Optimal Energy Mix, 2025–2050",
    subtitle = "Net-zero constraint enforced in 2050",
    x = "Year", y = "Energy Production (TWh)"
  ) +
  theme_netzero
Stacked area chart showing the optimal energy mix. The dashed line is total demand.

Stacked area chart showing the optimal energy mix. The dashed line is total demand.

8.2 Emissions Trajectory

ggplot(emissions_df, aes(x = year, y = emissions_Mt)) +
  geom_area(fill = "#e76f51", alpha = 0.35) +
  geom_line(colour = "#e76f51", linewidth = 1.2) +
  geom_point(colour = "#e76f51", size = 2.2) +
  geom_hline(yintercept = 0, linetype = "dashed",
             colour = "grey20", linewidth = 0.8) +
  annotate("text", x = 2026, y = max(emissions_df$emissions_Mt) * 0.04,
           label = "Net-zero target", hjust = 0, size = 3.2) +
  scale_x_continuous(breaks = seq(2025, 2050, 5)) +
  scale_y_continuous(labels = comma) +
  labs(
    title    = "CO₂ Emissions Trajectory, 2025–2050",
    subtitle = "Grid emissions under cost-optimal decarbonisation pathway",
    x = "Year", y = "Emissions (Mt CO₂)"
  ) +
  theme_netzero +
  theme(legend.position = "none")
Total annual CO₂ emissions declining to zero in 2050.

Total annual CO₂ emissions declining to zero in 2050.

8.3 Capacity Growth

ggplot(cap_df, aes(x = year, y = capacity, fill = source)) +
  geom_area(alpha = 0.85, colour = "white", linewidth = 0.2) +
  scale_fill_manual(values = source_colours, name = "Source") +
  scale_x_continuous(breaks = seq(2025, 2050, 5)) +
  scale_y_continuous(labels = comma) +
  labs(
    title    = "Installed Capacity Growth by Source, 2025–2050",
    subtitle = "GW of installed generation capacity",
    x = "Year", y = "Installed Capacity (GW)"
  ) +
  theme_netzero
Total installed capacity by source over the planning horizon.

Total installed capacity by source over the planning horizon.

8.4 Annual New Builds

build_df |>
  filter(build > 0.01) |>
  ggplot(aes(x = year, y = build, fill = source)) +
  geom_col(position = "stack", width = 0.8) +
  scale_fill_manual(values = source_colours, name = "Source") +
  scale_x_continuous(breaks = seq(2025, 2050, 5)) +
  labs(
    title    = "Annual New Capacity Built, 2025–2050",
    subtitle = "GW of new generation capacity commissioned each year",
    x = "Year", y = "New Capacity (GW)"
  ) +
  theme_netzero
New capacity commissioned each year by source type.

New capacity commissioned each year by source type.

9 · Scenario Analysis

The parameter block in Section 4 is the single point of modification for alternative scenarios. Adjust any values there and re-knit the document to regenerate all tables and charts automatically.

Scenario Parameter change
A — Accelerated solar B_max["solar"] <- 15
B — Nuclear phase-in K0["nuclear"] <- 3; B_max["nuclear"] <- 2
C — High demand growth D <- 300 + 8 * (yr_idx - 1)
D — Stricter interim target Add interim emissions cap constraint at t = 11 (year 2035)
E — Carbon capture on gas Reduce e_int["gas"] to reflect partial CCS

To compare multiple scenarios programmatically, wrap Sections 4–8 in a named list of parameter sets and bind result data frames prior to plotting.

10 · Session Info

sessionInfo()
#> R version 4.5.1 (2025-06-13)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Tahoe 26.2
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: Australia/Sydney
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] kableExtra_1.4.0      knitr_1.50            scales_1.4.0         
#> [4] tidyr_1.3.1           dplyr_1.1.4           ggplot2_3.5.2        
#> [7] ROI.plugin.glpk_1.0-0 ompr.roi_1.0.2        ompr_1.0.4           
#> 
#> loaded via a namespace (and not attached):
#>  [1] sass_0.4.10         generics_0.1.4      xml2_1.3.8         
#>  [4] slam_0.1-55         stringi_1.8.7       lattice_0.22-7     
#>  [7] digest_0.6.37       magrittr_2.0.3      evaluate_1.0.4     
#> [10] grid_4.5.1          RColorBrewer_1.1-3  fastmap_1.2.0      
#> [13] jsonlite_2.0.0      Matrix_1.7-3        backports_1.5.0    
#> [16] purrr_1.1.0         listcomp_0.4.1      viridisLite_0.4.2  
#> [19] textshaping_1.0.1   numDeriv_2016.8-1.1 lazyeval_0.2.2     
#> [22] jquerylib_0.1.4     cli_3.6.5           registry_0.5-1     
#> [25] rlang_1.1.6         ROI_1.0-2           Rglpk_0.6-5.1      
#> [28] withr_3.0.2         cachem_1.1.0        yaml_2.3.10        
#> [31] tools_4.5.1         checkmate_2.3.3     vctrs_0.6.5        
#> [34] R6_2.6.1            lifecycle_1.0.4     stringr_1.5.1      
#> [37] pkgconfig_2.0.3     bslib_0.9.0         pillar_1.11.0      
#> [40] gtable_0.3.6        data.table_1.17.8   glue_1.8.0         
#> [43] systemfonts_1.3.2   xfun_0.52           tibble_3.3.0       
#> [46] tidyselect_1.2.1    rstudioapi_0.17.1   farver_2.1.2       
#> [49] htmltools_0.5.8.1   labeling_0.4.3      svglite_2.2.2      
#> [52] rmarkdown_2.29      compiler_4.5.1