Capacity Expansion

---

Lecture 16

November 2, 2026

Overview

• Review and Questions
• The Capacity Expansion Problem
• Solving and Reading the Answer
• Adding Renewables
• Carbon Limits and Fuel Prices
• Key Takeaways and Upcoming Schedule
• Appendix

Review and Questions

Economic Dispatch

• Dispatch runs the plants we already have at least cost: the merit order, bent by minimum outputs and ramping.
• The price of electricity is the demand constraint’s shadow price: the cost of the plant that served the last MWh.

Raises the question: what should we build?

Ask whether anything about economic dispatch was unclear, especially the sign of the shadow price (relaxing demand means serving one less MWh, so JuMP reports a negative number).

Questions?

Text: VSRIKRISH to 22333

URL: https://pollev.com/vsrikrish

See Results

The Capacity Expansion Problem

Capacity Expansion

Capacity expansion: which generating resources to build, and how much, to meet future demand (load).

• Typical objective: minimize cost.
• Other goals may enter as constraints: reducing CO₂ emissions, fuel diversity, reliability.

We’ll look at a greenfield problem (no existing plants), not a brownfield one (adding to an existing fleet).

Capacity expansion sits at the long end of the decision time scales: years, not hours. The dispatch decision is nested inside it: whatever we build has to be run every hour, and the running costs are part of the case for building it. Most real planning is brownfield, adding to a fleet that is already there.

Plant Types

• natural gas: combined cycle (CCGT) or combustion turbines (CT);
• coal and nuclear;
• renewables: wind, solar, hydro, geothermal.

They differ in what they cost to build and what they cost to run, and the model has to weigh one against the other.

Preview the result: cheap-to-build plants that are expensive to run (CTs) suit loads that occur only a few hours a year; expensive-to-build plants that are cheap to run suit loads that are there all year. The rest of the lecture makes that precise.

Generator and Demand Data

gens = DataFrame(CSV.File("data/capacity_expansion/generators.csv"))
n_types = 4                                # the thermal types: nuclear, coal, CCGT and CT
gas_price = 7.21                           # delivered to U.S. power plants in 2022 [$/MMBtu]
coal_price = 2.36

# Variable cost of each plant type [$/MWh, to the dollar] at given fuel prices:
# the fuel burned per MWh times its price, plus other running costs
function variable_costs(gas, coal)
    costs = zeros(Int, nrow(gens))
    for g in 1:nrow(gens)
        fuel_price = 0
        if gens[g, :Fuel] == "Gas"
            fuel_price = gas
        elseif gens[g, :Fuel] == "Coal"
            fuel_price = coal
        end
        costs[g] = round(Int, gens[g, :OtherVarCost] + gens[g, :HeatRate] * fuel_price)
    end
    return costs
end
# The data file's variable costs are the 2022 ones
@assert variable_costs(gas_price, coal_price) == gens.VarCost

println("| Generator | Fixed cost (\\\$/MW-yr) | Variable cost (\\\$/MWh) |")
println("|:--|--:|--:|")
for g in 1:n_types
    println("| $(gens[g, :Plant]) | $(commas(gens[g, :FixedCost])) | $(gens[g, :VarCost]) |")
end

Generator	Fixed cost ($/MW-yr)	Variable cost ($/MWh)
Nuclear	790,000	8
Coal	392,000	27
NG CCGT	82,000	48
NG CT	65,000	69

Variable costs at 2022 fuel prices: gas $7.21, coal $2.36 per million Btu (MMBtu).

Demand: NYISO Zone C, every hour of 2020 (peak 2,781 MW).

Fixed costs are annualized: the cost of building one MW of capacity, spread over the plant’s life, per year. Variable costs are per MWh generated, and they are mostly fuel. The fuel prices are what U.S. power plants paid on average in 2022, when natural gas cost about twice what it does now. A CCGT burns about 6.27 million Btu of gas per MWh and a CT about 9.12, so every dollar on the gas price adds about 6 and 9 dollars to their costs per MWh. The last part of the lecture changes that price. Nuclear is very expensive to build and cheap to run; combustion turbines are the reverse. Wind, solar and geothermal return later.

Load Duration Curve

sorted_demand = sort(demand, rev=true)
p_ldc = plot(1:n_hours, sorted_demand, color=:black, xlabel="Hours of the year with at least this load",
    ylabel="Load  [MW]", xlims=(0, n_hours), ylims=(0, 3000))
plot!(p_ldc, size=(1200, 360), left_margin=12mm)

The year’s hourly loads, sorted from highest to lowest, show how many hours each load level is reached.

The chronological demand curve (one point per hour through the year) makes it hard to see how often each level occurs. Sorting turns it into a load duration curve: read any point as “the load is at least this high for this many hours”. The top few hundred MW are needed only a handful of hours a year; the bottom 1,000 MW or so are needed every hour.

Decision Variables

What are our variables?

Variable	Meaning
\(x_g\)	installed capacity (MW) of generator type \(g \in \mathcal{G}\)
\(y_{g,t}\)	production (MWh) of type \(g\) in hour \(t \in \mathcal{T}\)
\(NSE_t\)	non-served energy (MWh) in hour \(t\)

Must we meet the very highest hour? Price non-served energy at $9,000/MWh and let the model decide.

Two kinds of decisions: how much to build (x, one per technology) and how to run it each hour (y, one per technology per hour). Non-served energy is load we choose not to serve; its cost stands in for the damage of a shortfall. 9,000 dollars per MWh is in the range of real market price caps.

The Formulation

Write the objective and constraints.

\[ \begin{aligned} \min_{\mathbf{x}, \mathbf{y}, \mathbf{NSE}} \quad & \sum_{g \in \mathcal{G}} \text{FixedCost}_g \, x_g + \sum_{t \in \mathcal{T}} \sum_{g \in \mathcal{G}} \text{VarCost}_g \, y_{g,t} \\ & \quad + \sum_{t \in \mathcal{T}} \text{NSECost} \, NSE_t \\[0.3em] \text{subject to:} \quad & \sum_{g \in \mathcal{G}} y_{g,t} + NSE_t \geq d_t \qquad \forall t \in \mathcal{T} \\[0.3em] & y_{g,t} \leq x_g \qquad \forall g \in \mathcal{G}, t \in \mathcal{T} \\[0.3em] & x_g, \, y_{g,t}, \, NSE_t \geq 0 \end{aligned} \]

Board item; the fragment is the answer. Build the objective in two pieces: investment (fixed cost times capacity) and operation (variable cost times generation, plus the penalty on unserved energy). Constraints: each hour, generation plus unserved energy covers demand; no plant produces more than its installed capacity; everything is non-negative. Ask: is this an LP? Yes: it’s linear in every variable, capacity is divisible (we model total MW, not individual units), and demand is treated as known. Each of those three is an assumption that later methods relax: discrete units (mixed integer programming) and uncertainty (stochastic programming).

Solving and Reading the Answer

Implementation in JuMP

NSECost = 9000
gencap = Model(HiGHS.Optimizer)
@variables(gencap, begin
    x[g in 1:n_types] >= 0                  # capacity [MW]
    y[g in 1:n_types, t in 1:n_hours] >= 0  # generation [MWh]
    NSE[t in 1:n_hours] >= 0                # non-served energy [MWh]
end)
@objective(gencap, Min, sum(gens[g, :FixedCost] * x[g] for g in 1:n_types) +
    sum(gens[g, :VarCost] * y[g, t] for g in 1:n_types, t in 1:n_hours) + NSECost * sum(NSE))
@constraint(gencap, load[t in 1:n_hours], sum(y[:, t]) + NSE[t] >= demand[t])
@constraint(gencap, availability[g in 1:n_types, t in 1:n_hours], y[g, t] <= x[g])
set_silent(gencap); optimize!(gencap)

8,784 hours times four technologies is about 35,000 generation variables; HiGHS solves it in a few seconds. Name the load and availability constraints so we can ask for their shadow prices later.

Utilization

How much of the year is a plant’s capacity used?

\[\text{utilization} = \frac{\text{energy produced in the year}}{\text{capacity} \times \text{hours in the year}}\]

• A 100 MW plant that produced 219,600 MWh in 2020’s 8,784 hours: 219,600 ÷ (100 × 8,784) = 25%.
• It depends on how the plant is run, which the model decides hour by hour.

A plant that ran at full output all year would be at 100%, and one that never ran at zero. Many sources call this number the plant’s capacity factor. Today we keep the term capacity factor for what the weather allows a wind or solar plant to produce, which comes later, and use utilization for what the model chooses to run.

Results

capacity = value.(x)
generation = [sum(value.(y[g, :])) for g in 1:n_types]
unserved = value.(NSE)
nse_hours = count(unserved .> 1e-6)
# Cost and emissions of this plan, to compare with the plans that follow
thermal_cost = objective_value(gencap)                                            # [$/yr]
thermal_emissions = sum(gens[g, :Emissions] * generation[g] for g in 1:n_types)   # [t CO₂/yr]
# One decimal place, and plain 0 rather than -0.0 for technologies that aren't built
function share(v)
    rounded_v = round(v; digits=1)
    if rounded_v == 0
        return "0"
    end
    return string(rounded_v)
end
# Energy produced as a share of what the capacity could produce in a year [%];
# a dash for a type that isn't built
function utilization(built, produced)
    if built < 0.5
        return "–"
    end
    return string(round(Int, 100 * produced / (built * n_hours)))
end
println("| Resource | Installed (MW) | % of peak | % of energy | Utilization (%) |")
println("|:--|--:|--:|--:|--:|")
for g in 1:n_types
    println("| $(gens[g, :Plant]) | $(commas(capacity[g])) | $(share(100 * capacity[g] / peak_demand)) | ",
        "$(share(100 * generation[g] / sum(demand))) | $(utilization(capacity[g], generation[g])) |")
end

Resource	Installed (MW)	% of peak	% of energy	Utilization (%)
Nuclear	0	0	0	–
Coal	0	0	0	–
NG CCGT	2,154	77.4	99.2	81
NG CT	596	21.4	0.8	2

Up to 31 MW goes unserved, in 7 hours of the year.

No nuclear and no coal: neither is worth building at these costs. Combined-cycle gas carries almost all of the energy, and its utilization is 81%. Combustion turbines are 21% of peak demand but supply only 0.8% of the energy: they’re for the peaks. And there is some unserved energy, despite the 9,000 dollar penalty: 94 MWh over 7 hours. The next slides explain all three facts from the cost data alone. The total cost is $960M a year.

Which Plant Serves Which Load?

# One colour per plant type, in the order of the data file
type_colors = [cb_purple, :gray40, cb_orange, cb_vermillion, cb_blue, cb_sky, cb_yellow]

# The load duration curve, filled by the plant type that serves each level of load.
# Plants are used in order of running cost (`order`, cheapest first), so each type
# fills a band as tall as its capacity.
function band_plot(build, order)
    p = plot(xlabel="Hours of the year with at least this load", ylabel="Load  [MW]",
        xlims=(0, n_hours), ylims=(0, 3000), legend=:outerright)
    lower = zeros(n_hours)
    served = 0
    for g in order
        if build[g] > 0.5
            served += build[g]
            upper = min.(sorted_demand, served)
            plot!(p, 1:n_hours, upper, fillrange=lower, color=type_colors[g], linewidth=0,
                label=gens[g, :Plant])
            lower = upper
        end
    end
    plot!(p, 1:n_hours, sorted_demand, fillrange=lower, color=:black, linewidth=0, label="Unserved")
    plot!(p, 1:n_hours, sorted_demand, color=:black, linewidth=2)
    return p
end

thermal_build = zeros(nrow(gens))
thermal_build[1:n_types] = capacity
p_bands = band_plot(thermal_build, [1, 2, 3, 4])
plot!(p_bands, size=(1200, 330), left_margin=12mm)

The CCGT band is full most of the year (utilization 81%). The CT band is used only at the far left (2%).

The load duration curve again, coloured by the plant that serves each MW of load. The plants run in order of cost, so the CCGTs take every MW of load up to their capacity and the CTs take what is above it. The height of a band is that type’s capacity and its area is its energy, so utilization is how much of its rectangle the band fills. The CCGT rectangle is mostly full. The CT rectangle is nearly empty: 596 MW that are used 2% of the time. The load at the far right, which is there in every hour, is called base load; the load at the far left, there for a few hours, is peak load. Here one type serves both the base and the middle. The unserved load is a sliver above the CTs at the left edge, too thin to see at this scale.

Screening Curves

# Annual cost of one MW of each technology, as a function of how many hours it runs
hours_run = 0:n_hours
annual_cost = zeros(length(hours_run), n_types + 1)
for g in 1:n_types
    annual_cost[:, g] = gens[g, :FixedCost] .+ gens[g, :VarCost] .* hours_run
end
annual_cost[:, end] = NSECost .* hours_run          # not serving: no build cost

ct_vs_ccgt = (gens[3, :FixedCost] - gens[4, :FixedCost]) / (gens[4, :VarCost] - gens[3, :VarCost])
nse_vs_ct = gens[4, :FixedCost] / (NSECost - gens[4, :VarCost])
# The next slide's equations quote these four numbers
@assert gens[3, :FixedCost] == 82000 && gens[3, :VarCost] == 48
@assert gens[4, :FixedCost] == 65000 && gens[4, :VarCost] == 69
p_screen = plot(xlabel="Hours run per year", ylabel="Cost per MW  [\$1000/yr]",
    xlims=(0, n_hours), ylims=(0, 900), legend=:outerright)
for (j, (name, color, style)) in enumerate(zip([gens[1:n_types, :Plant]; "Unserved"],
        [cb_purple, :gray40, cb_orange, cb_vermillion, :black], [:solid, :dash, :solid, :dashdot, :dot]))
    plot!(p_screen, hours_run, annual_cost[:, j] ./ 1000, color=color, linestyle=style, label=name)
end
vline!(p_screen, [ct_vs_ccgt], color=:black, linewidth=1, linestyle=:dash, label=false)
annotate!(p_screen, ct_vs_ccgt + 150, 720, text("$(commas(ct_vs_ccgt)) h", 16, :left, plot_font))
plot!(p_screen, size=(1200, 440), left_margin=12mm)

Each line is what one MW of a technology costs per year if it runs for a given number of hours: fixed cost at zero hours, rising at the variable cost. The cheapest technology for a load that lasts h hours is the lowest line at h. At these prices nuclear and coal are not the lowest at any number of hours in a year. CTs are cheapest below the crossover with CCGT; CCGTs above. Unserved energy (the dotted line, steepest of all) is cheapest only for loads that last a few hours: its crossover with CTs is too close to zero to see at this scale.

Where the Lines Cross

Find the crossovers: below how many hours a year is a CT cheaper than a CCGT? Not serving cheaper than a CT?

\[ \begin{aligned} \text{CT vs. CCGT:} \quad & 65{,}000 + 69h = 82{,}000 + 48h \\ \text{Unserved vs. CT:} \quad & 9{,}000h = 65{,}000 + 69h \end{aligned} \]

A CT is cheaper below 810 hours; not serving is cheaper below 7.28 hours.

Board item; the fragment is the answer. Set the two annual costs equal and solve. CCGT costs $17,000 a year more to build and $21 per MWh less to run, so it pays back after about 810 hours. A CT costs $65,000 a year to build; not serving costs $8,931 more per MWh than running the CT, so a CT pays for itself only if it would run more than about seven hours a year. At these prices coal does not beat CCGT within a year: the crossover is at 14,762 hours. Nuclear’s crossover with CCGT is at about 17,700 hours.

The Optimal Capacity Mix

ccgt_level = sorted_demand[ceil(Int, ct_vs_ccgt)]
ct_top = sorted_demand[ceil(Int, nse_vs_ct)]
p_full = plot(1:n_hours, sorted_demand, color=:black, xlabel="Hours", ylabel="Load  [MW]",
    xlims=(0, n_hours), ylims=(0, 3000), title="Whole year")
hline!(p_full, [ccgt_level], color=cb_orange, linestyle=:dash)
vline!(p_full, [ct_vs_ccgt], color=:gray50, linewidth=1)
annotate!(p_full, 2300, ccgt_level + 200, text("CCGT below $(commas(ccgt_level)) MW", 15, :left, plot_font))
p_top = plot(1:24, sorted_demand[1:24], color=:black, marker=:circle, xlabel="Hours", ylabel="Load  [MW]",
    title="Top 24 hours", ylims=(2600, 2800))
hline!(p_top, [ct_top], color=cb_vermillion, linestyle=:dash)
vline!(p_top, [nse_vs_ct], color=:gray50, linewidth=1)
annotate!(p_top, 11, ct_top + 22, text("CT up to $(commas(ct_top)) MW", 15, :left, plot_font))
plot(p_full, p_top, layout=(1, 2), size=(1200, 350), left_margin=12mm)

CCGT 2,154 MW, CT 596 MW, 31 MW unserved, which is exactly what the LP builds.

Read the crossovers off the load duration curve. Load that lasts more than 810 hours is cheapest served by CCGT: that’s the curve’s height at 810 hours, 2,154 MW. Load lasting between about 7 and 810 hours is cheapest served by CTs: up to the curve’s height at 8 hours. Load lasting fewer hours isn’t worth building for at all. This reproduces the LP’s build to the megawatt, and it answers why there’s unserved energy despite the 9,000 dollar penalty: the last 31 MW are needed for only 7 hours a year, and building capacity for them would cost more than the penalty. The 7.28-hour threshold depends on the penalty: at 10,000 dollars per MWh it’s 6.5 hours.

Adding Renewables

What Changes With Renewables?

Wind and solar can’t produce at full capacity every hour: their capacity factor \(c_{g,t}\) depends on the weather.

The capacity constraint changes from \[y_{g,t} \leq x_g \qquad \text{to} \qquad y_{g,t} \leq x_g \times c_{g,t}\]

For thermal plants, c is 1 in every hour. A MW of solar that produces at midday when load is low is worth less than one that produces at the evening peak. That’s why the screening-curve shortcut breaks: the cost of a renewable MW depends on the hours it produces, not just on how many. The next slides add wind and solar to today’s problem and work through what changes.

Where Capacity Factors Come From

• Resource data: wind, sun, water or heat available each hour.
• Often a typical meteorological year (TMY) or a few actual meteorological years (AMY), sometimes only selected periods.
• Both underrepresent variability, but keep the problem tractable.

Always check what weather data a study used. Ours is for a good wind and solar site, not New York.

A typical meteorological year stitches together typical months from many years, which smooths out the extremes: the windless week, the cloudy month. Planning for a typical year can under-build for the bad ones. That’s a structural choice in the model, in the sense of October 7. Today’s profiles are for a much better site than upstate New York, so the plan they give is for that kind of resource. It is not a forecast for New York.

Wind and Solar Capacity Factors

profiles = DataFrame(CSV.File("data/capacity_expansion/wind_solar_capacity_factors.csv"))
# The profiles cover a 365-day year and 2020 has 366 days, so February 28 is used again for February 29
before_leap_day = 1:(59 * 24)
feb28 = (58 * 24 + 1):(59 * 24)
after_leap_day = (59 * 24 + 1):nrow(profiles)
wind_cf = [profiles.Wind[before_leap_day]; profiles.Wind[feb28]; profiles.Wind[after_leap_day]]
solar_cf = [profiles.Solar[before_leap_day]; profiles.Solar[feb28]; profiles.Solar[after_leap_day]]
wind_average = sum(wind_cf) / n_hours
solar_average = sum(solar_cf) / n_hours

# Average capacity factor in each hour of the day
n_days = div(n_hours, 24)
wind_by_hour = zeros(24)
solar_by_hour = zeros(24)
for h in 1:24
    wind_by_hour[h] = sum(wind_cf[h:24:end]) / n_days
    solar_by_hour[h] = sum(solar_cf[h:24:end]) / n_days
end
p_profiles = plot(1:24, wind_by_hour, color=cb_sky, label="wind", xlabel="Hour of the day",
    ylabel="Average capacity factor", ylims=(0, 1), xticks=0:6:24, legend=:topleft)
plot!(p_profiles, 1:24, solar_by_hour, color=cb_orange, linestyle=:dash, label="solar")
plot!(p_profiles, size=(1200, 350), left_margin=12mm)

Over the year, wind averages 34% and solar 35%.

The figure averages each hour of the day over the whole year. Solar produces only in daylight and peaks at midday. Wind is steadier across the day and a little stronger at night. The two annual averages are close, but the hours are very different, and that is what will matter.

Geothermal, Wind and Solar

geothermal_row = n_types + 1
wind_row = n_types + 2
solar_row = n_types + 3
n_all = nrow(gens)
# Average share of capacity that is available over the year, by type
average_available = ones(n_all)
average_available[wind_row] = wind_average
average_available[solar_row] = solar_average
println("::: {.small}\n")
println("| Generator | Fixed cost (\\\$/MW-yr) | Variable cost (\\\$/MWh) | Average capacity factor (%) |")
println("|:--|--:|--:|--:|")
for g in (n_types + 1):n_all
    println("| $(gens[g, :Plant]) | $(commas(gens[g, :FixedCost])) | $(gens[g, :VarCost]) | ",
        "$(round(Int, 100 * average_available[g])) |")
end
println("\n:::")

Generator	Fixed cost ($/MW-yr)	Variable cost ($/MWh)	Average capacity factor (%)
Geothermal	470,000	0	100
Wind	91,000	0	34
Solar	70,000	0	35

These three join the four thermal types. Geothermal is available in every hour, but only where the resource exists.

None of the three has a fuel cost. Wind and solar are cheap to build and available only when the weather allows. Geothermal costs more to build than coal and runs whenever it is wanted. Its cost is the Energy Information Administration’s estimate for a 50 MW plant. Real sites are small and need the right geology, so letting the model build as much as it wants is generous to geothermal. Because geothermal is always available, the screening rule still applies to it: it beats a CCGT for loads that last more than 8,083 hours a year.

Do Screening Curves Still Work?

# Wind and solar on the screening chart: a fixed cost each year, and nothing per MWh
solar_vs_ct = (gens[solar_row, :FixedCost] - gens[4, :FixedCost]) / gens[4, :VarCost]
wind_vs_ccgt = (gens[wind_row, :FixedCost] - gens[3, :FixedCost]) / gens[3, :VarCost]
# What a MWh costs once the capacity factor is counted [$/MWh]
wind_energy_cost = gens[wind_row, :FixedCost] / (wind_average * n_hours)
solar_energy_cost = gens[solar_row, :FixedCost] / (solar_average * n_hours)
geothermal_energy_cost = gens[geothermal_row, :FixedCost] / n_hours
# A CCGT's whole cost per MWh when it runs all year [$/MWh]
ccgt_all_year_cost = gens[3, :FixedCost] / n_hours + gens[3, :VarCost]

p_renewable_screen = plot(xlabel="Hours run per year", ylabel="Cost per MW  [\$1000/yr]",
    xlims=(0, n_hours), ylims=(0, 520), legend=:outerright)
plot!(p_renewable_screen, hours_run, annual_cost[:, 3] ./ 1000, color=cb_orange, label="NG CCGT")
plot!(p_renewable_screen, hours_run, annual_cost[:, 4] ./ 1000, color=cb_vermillion, linestyle=:dashdot, label="NG CT")
hline!(p_renewable_screen, [gens[wind_row, :FixedCost] / 1000], color=cb_sky, linestyle=:dash, label="Wind")
hline!(p_renewable_screen, [gens[solar_row, :FixedCost] / 1000], color=cb_green, linestyle=:dot, label="Solar")
plot!(p_renewable_screen, size=(1200, 370), left_margin=12mm)

By the screening rule, solar serves every load over 72 hours a year and no CCGT is built. What is missing?

Board item. Wind and solar cost the same each year however much they run, so their lines are flat. Solar’s line is below CCGT’s at every number of hours, and below CT’s beyond 72 hours. Wind’s is below CCGT’s beyond 188 hours. Let the room say what is wrong before moving on. The chart assumes a plant can run in whichever hours it is needed. A MW of solar runs only when the sun is up: 35% of its capacity on average, mostly around midday, and never on request.

The Cost of a Solar MWh

A MW of solar costs $70,000 a year and averages 35% of its capacity. What does a MWh cost?

• Solar: 70,000 ÷ (0.346 × 8,784 h) = $23/MWh
• Wind: 91,000 ÷ (0.342 × 8,784 h) = $30/MWh
• CCGT fuel: $48/MWh

Both cost less than the gas they save, so some of each is worth building.

Board item; the fragments carry the answer. Divide the annual cost of a MW by the MWh it actually produces. A CCGT that is already built saves $48 of fuel for every MWh something else supplies. A MWh from wind or from solar costs less than that, so replacing gas energy with either saves money. Neither can be counted on in a given hour, so neither replaces the gas capacity itself. Comparing a new plant’s cost per MWh with the fuel it would replace is the comparison to make whenever the capacity it competes with is needed anyway.

Geothermal is a different case, because it is always available. Its energy costs $54 per MWh, more than a CCGT’s fuel. But it can take the place of a CCGT altogether, and a CCGT that runs all year costs $57 per MWh once its fixed cost is counted.

Adding Renewables in JuMP

# Share of each type's capacity that is available each hour: all of it for the
# thermal types and geothermal, the capacity factor for wind and solar
available = ones(n_all, n_hours)
available[wind_row, :] = wind_cf
available[solar_row, :] = solar_cf

# The least-cost plan for the listed plant types, with an optional cap on annual
# emissions [t CO₂]. Costs default to the data file's; pass others to change them.
function solve_plan(types; cap=Inf, fixed_cost=gens.FixedCost, var_cost=gens.VarCost)
    model = Model(HiGHS.Optimizer)
    @variables(model, begin
        x[g in types] >= 0                  # capacity [MW]
        y[g in types, t in 1:n_hours] >= 0  # generation [MWh]
        nse[t in 1:n_hours] >= 0            # non-served energy [MWh]
    end)
    @expression(model, emitted, sum(gens[g, :Emissions] * y[g, t] for g in types, t in 1:n_hours))
    @objective(model, Min, sum(fixed_cost[g] * x[g] for g in types) +
        sum(var_cost[g] * y[g, t] for g in types, t in 1:n_hours) + NSECost * sum(nse))
    @constraint(model, [t in 1:n_hours], sum(y[g, t] for g in types) + nse[t] >= demand[t])
    @constraint(model, [g in types, t in 1:n_hours], y[g, t] <= x[g] * available[g, t])
    if isfinite(cap)
        @constraint(model, emissions_cap, emitted <= cap)
    end
    set_silent(model)
    optimize!(model)

    build = zeros(n_all)              # capacity by type [MW]; zero for types left out
    hourly = zeros(n_all, n_hours)    # generation by type and hour [MWh]
    energy = zeros(n_all)             # generation by type over the year [MWh]
    for g in types
        build[g] = value(x[g])
        for t in 1:n_hours
            hourly[g, t] = value(y[g, t])
        end
        energy[g] = sum(hourly[g, :])
    end
    carbon_price = 0
    if isfinite(cap)
        carbon_price = shadow_price(emissions_cap)
    end
    return (types=types, build=build, hourly=hourly, energy=energy, unserved=value.(nse),
        cost=objective_value(model), emissions=value(emitted), carbon_price=carbon_price)
end

@constraint(model, [g in types, t in 1:n_hours],
    y[g, t] <= x[g] * available[g, t])

Only the availability constraint changes: available[g, t] is 1 for the thermal types and geothermal, and the capacity factor for wind and solar.

The folded cell holds the whole model as a function, so the same code can be solved for any set of plant types, with or without an emissions cap, and later at other fuel prices. Everything in it is the model from earlier except this one line. Writing the constraint with indexed sets keeps it one line however many technologies there are.

Results With Renewables

with_renewables = solve_plan(1:n_all)
# Totals of the whole-MW values the table shows, so the text and the table agree
renewables_added = round(with_renewables.build[wind_row]) + round(with_renewables.build[solar_row])
gas_before = round(capacity[3]) + round(capacity[4])
gas_after = round(with_renewables.build[3]) + round(with_renewables.build[4])
# The notes say geothermal is not built
@assert with_renewables.build[geothermal_row] < 0.5

# For the notes: the plan if geothermal were the only addition, and the plan at
# overnight costs with no tax credits (EIA: wind $153,000, solar $141,000 per MW-yr)
with_geothermal = solve_plan(1:geothermal_row)
no_credit_cost = copy(gens.FixedCost)
no_credit_cost[wind_row] = 153000
no_credit_cost[solar_row] = 141000
no_credit = solve_plan(1:n_all; fixed_cost=no_credit_cost)

println("::: {.small}\n")
println("| | Thermal types only | With renewables |")
println("|:--|--:|--:|")
for g in 1:n_all
    if g <= n_types
        before = commas(capacity[g])
        built_before = capacity[g]
    else
        before = "–"
        built_before = 0
    end
    if built_before > 0.5 || with_renewables.build[g] > 0.5 || g == geothermal_row
        println("| $(gens[g, :Plant]) (MW) | $(before) | $(commas(with_renewables.build[g])) |")
    end
end
println("| Cost (\\\$M/yr) | $(commas(thermal_cost / 1e6)) | $(commas(with_renewables.cost / 1e6)) |")
println("| CO~2~ (kt/yr) | $(commas(thermal_emissions / 1000)) | $(commas(with_renewables.emissions / 1000)) |")
println("\n:::")

	Thermal types only	With renewables
NG CCGT (MW)	2,154	1,579
NG CT (MW)	596	650
Geothermal (MW)	–	0
Wind (MW)	–	1,842
Solar (MW)	–	1,817
Cost ($M/yr)	960	729
CO2 (kt/yr)	6,659	2,325

Wind and solar add 3,659 MW; gas falls by only 521 MW.

Three things change. Wind and solar are built: 3,659 MW between them. The gas plants run less: the CCGTs’ utilization falls from 81% to 37%, which is where the saving of $231M a year and the 65% cut in emissions come from. And the gas capacity falls much less, from 2,750 to 2,229 MW, because demand after dark still has to be met.

Nuclear, coal and geothermal are not built. Geothermal’s energy costs more than wind’s and solar’s, and gas supplies the capacity for the evenings at a lower cost. If wind and solar could not be built, the plan would include 1,386 MW of geothermal, running all year.

This result depends on the costs. These wind and solar costs give energy at about $30 and $23 per MWh at this site, close to recent U.S. estimates that include tax credits. At overnight costs without credits ($153,000 and $141,000 per MW-year, from the Energy Information Administration’s capital cost estimates), the same model builds 681 MW of wind and solar and 1,135 MW of geothermal.

Load Net of Wind and Solar

# What is left for the other plants each hour, once wind and solar have produced
# all that the weather allows
net_load = demand .- with_renewables.build[wind_row] .* wind_cf .-
    with_renewables.build[solar_row] .* solar_cf
sorted_net = sort(net_load, rev=true)
surplus_hours = count(net_load .< 0)
# Share of what wind and solar could have produced that goes unused [%]
could_produce = with_renewables.build[wind_row] * sum(wind_cf) + with_renewables.build[solar_row] * sum(solar_cf)
unused_share = 100 * (1 - (with_renewables.energy[wind_row] + with_renewables.energy[solar_row]) / could_produce)

p_net = plot(1:n_hours, sorted_demand, color=:black, label="Load",
    xlabel="Hours of the year with at least this load", ylabel="Load  [MW]",
    xlims=(0, n_hours), legend=:outerright)
plot!(p_net, 1:n_hours, sorted_net, color=cb_orange, linestyle=:dash, label="Net of wind and solar")
hline!(p_net, [0], color=:gray50, linewidth=1)
plot!(p_net, size=(1200, 320), left_margin=12mm)

The other plants’ peak falls only from 2,781 to 2,316 MW. In 2,466 hours, wind and solar alone could exceed the load.

The load duration curve again, with a second curve for what is left after wind and solar. Each curve is sorted separately, so a position along the axis is not the same hour on both. Read three things. The left end, the highest net load, falls by only about a sixth, because the highest net loads come on evenings with little wind. That is why the gas capacity falls so little. The middle of the curve drops by more than half, so the gas plants that remain are used far less. And the right end goes below zero: in those hours wind and solar could produce more than the whole load, and 9% of what they could produce over the year goes unused.

Carbon Limits and Fuel Prices

Capping Emissions

half_cap = 0.5 * thermal_emissions          # [t CO₂/yr]
capped_thermal = solve_plan(1:n_types; cap=half_cap)
capped_all = solve_plan(1:n_all; cap=half_cap)
capped_quarter = solve_plan(1:n_all; cap=0.25 * thermal_emissions)
capped_tenth = solve_plan(1:n_all; cap=0.1 * thermal_emissions)
# For the notes: the tightest cap where there is no geothermal resource
no_geothermal = [1, 2, 3, 4, wind_row, solar_row]
capped_tenth_no_geothermal = solve_plan(no_geothermal; cap=0.1 * thermal_emissions)

# The next slide's title and formula quote these numbers
@assert round(Int, -capped_thermal.carbon_price) == 94
@assert gens[1, :FixedCost] == 790000 && gens[1, :VarCost] == 8 && gens[3, :Emissions] == 0.43
# The table leaves coal out, which is right only if no plan builds it
@assert max(with_renewables.build[2], capped_thermal.build[2], capped_all.build[2],
    capped_quarter.build[2], capped_tenth.build[2]) < 0.5

# One table row for a plan: MW by type (a dash where the plan could not build that
# type), then its cost and the cap's shadow price (0 when the cap does not bind)
shown_types = [1, 3, 4, geothermal_row, wind_row, solar_row]
function plan_row(label, plan)
    row = "| $(label) |"
    for g in shown_types
        if g in plan.types
            row *= " $(commas(plan.build[g])) |"
        else
            row *= " – |"
        end
    end
    price_text = "0"
    if -plan.carbon_price > 0.005
        price_text = "−" * string(round(-plan.carbon_price; digits=2))
    end
    return row * " $(commas(plan.cost / 1e6)) | $(price_text) |"
end

println("::: {.small}\n")
println("| Types,\u00a0cap | Nuclear | CCGT | CT | Geothermal | Wind | Solar | Cost | Shadow\u00a0price |")
println("|:--|--:|--:|--:|--:|--:|--:|--:|--:|")
println(plan_row("Thermal,\u00a050%", capped_thermal))
println(plan_row("All,\u00a050%", capped_all))
println(plan_row("All,\u00a025%", capped_quarter))
println(plan_row("All,\u00a010%", capped_tenth))
println("\n:::")

Types, cap	Nuclear	CCGT	CT	Geothermal	Wind	Solar	Cost	Shadow price
Thermal, 50%	880	1,345	525	–	–	–	1,275	−94.42
All, 50%	0	1,579	650	0	1,842	1,817	729	0
All, 25%	0	1,486	581	156	2,210	1,941	759	−78.46
All, 10%	0	708	451	1,102	1,000	1,027	845	−109.7

Capacities in MW, cost in $M a year, shadow price in $/t. Emission factors \(e_g\): coal 1, CCGT 0.43, CT 0.55 t/MWh; the other types 0.

Each cap is a share of the thermal-only plan’s emissions, 6,659 kt of CO₂ a year: \(\sum_{g,t} e_g \, y_{g,t} \leq E\).

Without a cap the thermal-only plan emits 6,659 thousand tonnes a year. Halving that with thermal types only means building 880 MW of nuclear in place of CCGT, and the plan’s cost goes from $960M to $1,275M a year.

With renewables available the same cap changes nothing. The uncapped plan with wind and solar already emits only 35% of the thermal-only total, so the cap does not bind and its shadow price is zero. Tighter caps do bind. At 25% the plan adds wind, solar and a little geothermal. At 10% geothermal serves the base of the load, and the plan builds less wind and solar than at 25%, because a plant that is always available leaves fewer hours for them to fill.

Nuclear is not built in any plan that can build geothermal, which costs less to build and nothing to run. Where there is no geothermal resource, nuclear is the alternative, and the 10% cap then costs $570 a tonne against $110.

Why $94 a Tonne?

Swap a MW of CCGT for nuclear, both running all year:

\[ \frac{\overbrace{(790{,}000 / 8{,}784 + 8)}^{\text{nuclear, \$/MWh}} - \overbrace{(82{,}000 / 8{,}784 + 48)}^{\text{CCGT, \$/MWh}}}{0.43 \text{ t/MWh}} \]

That is $94.4 per tonne, the thermal-only shadow price.

This carbon price is the cost of abating the last tonne. With renewables it is $0 at this cap, $78 at 25% and $110 at 10%.

Board item; the fragments carry the answer. With thermal types only, the LP reported a shadow price of −$94.42 a tonne on the last slide. The shadow price tells a utility what one more tonne of allowed emissions is worth, or what the last tonne of abatement costs. With thermal types only there is one way to abate: replace a CCGT with a nuclear plant, both running all year, and the arithmetic gives the shadow price. With renewables there are cheaper ways to avoid a tonne, and at a 50% cap they are in the plan already, so the price is zero. The price rises as the cap tightens because the cheapest abatement is used first. A cap’s shadow price depends on what the system is able to build.

Natural Gas Prices Since 2019

gas_history = DataFrame(CSV.File("data/capacity_expansion/gas_prices_monthly.csv";
    types=Dict(:Month => String)))
n_months = nrow(gas_history)
# Position of each month on the time axis: January 2019 is 2019.0
month_position = zeros(n_months)
for i in 1:n_months
    month_position[i] = 2019 + (i - 1) / 12
end
invasion = 2022 + 1 / 12          # February 2022: Russia invades Ukraine
hormuz = 2026 + 2 / 12            # March 2026: no tankers pass the Strait of Hormuz

# Averages the next slides' prices are based on [$/MMBtu]
months_2020 = 13:24
months_2022 = 37:48
months_2025 = 73:84
months_2026 = 85:n_months          # January to July 2026
@assert gas_history.Month[73] == "2025-01" && gas_history.Month[85] == "2026-01"
us_2020 = sum(gas_history.UnitedStates[months_2020]) / length(months_2020)
us_2022 = sum(gas_history.UnitedStates[months_2022]) / length(months_2022)
us_2025 = sum(gas_history.UnitedStates[months_2025]) / length(months_2025)
europe_2026 = sum(gas_history.Europe[months_2026]) / length(months_2026)
europe_peak, europe_peak_index = findmax(gas_history.Europe)
# The month of that peak in words: the file starts in January 2019
month_names = ["January", "February", "March", "April", "May", "June", "July", "August",
    "September", "October", "November", "December"]
europe_peak_month = month_names[mod(europe_peak_index - 1, 12) + 1] * " " *
    string(2019 + div(europe_peak_index - 1, 12))

p_gas = plot(month_position, gas_history.Europe, color=cb_vermillion, linestyle=:dash, label="Europe",
    ylabel="Gas price  [\$/MMBtu]", xticks=2019:2026, xlims=(2019, 2026.75), ylims=(0, 75),
    legend=:outerright)
plot!(p_gas, month_position, gas_history.UnitedStates, color=cb_blue, label="United States")
vline!(p_gas, [invasion, hormuz], color=:gray40, linewidth=1, linestyle=:dot)
annotate!(p_gas, invasion - 0.08, 66, text("Russia invades\nUkraine", 15, :right, plot_font))
annotate!(p_gas, hormuz - 0.08, 66, text("Strait of Hormuz\ncloses to tankers", 15, :right, plot_font))
plot!(p_gas, size=(1200, 300), left_margin=12mm)

U.S. gas averaged $3.53 in 2025 and Europe’s $15.06 so far in 2026. Our base is the $7.21 U.S. power plants paid in 2022.

Monthly prices in dollars per million Btu: the Henry Hub spot price for the United States, and the International Monetary Fund’s European price. Two shocks are marked. Before 2022 Europe bought much of its gas from Russia by pipeline. After the invasion of Ukraine it had to replace that gas with liquefied gas shipped from elsewhere, and its price reached $70 in August 2022. The U.S. price rose too, from $2.03 in 2020 to $6.42 in 2022: U.S. gas is exported as LNG, which ties the two markets together loosely. In 2026 tankers stopped passing the Strait of Hormuz at the end of February, which cut off Qatar’s LNG, and Europe’s price rose from $11 in February to $18 in March. The spike in the U.S. line in January 2026 is a winter storm.

Our base case uses what U.S. power plants paid for delivered gas in 2022, which is a little above that year’s Henry Hub average. The next slides solve the plan again at two other prices: $3.50, about the U.S. price in 2025, and $15, about Europe’s price in 2026.

Results at Other Gas Prices

today_gas = 3.5        # about the U.S. average for 2025 [$/MMBtu]
europe_gas = 15        # about Europe's average for January to July 2026
today_costs = variable_costs(today_gas, coal_price)
europe_costs = variable_costs(europe_gas, coal_price)

thermal_today = solve_plan(1:n_types; var_cost=today_costs)
all_today = solve_plan(1:n_all; var_cost=today_costs)
thermal_2022 = solve_plan(1:n_types)
thermal_europe = solve_plan(1:n_types; var_cost=europe_costs)
all_europe = solve_plan(1:n_all; var_cost=europe_costs)
# The figure leaves nuclear out, which is right only if none of these plans builds it
@assert max(thermal_today.build[1], all_today.build[1], thermal_2022.build[1],
    thermal_europe.build[1], all_europe.build[1]) < 0.5

# Capacity by type [MW] for the six plans: one row per plan, one column per type.
# The types are stacked in this order, from the bottom of each bar up.
plans = [thermal_today, all_today, thermal_2022, with_renewables, thermal_europe, all_europe]
stack_types = [2, geothermal_row, 3, 4, wind_row, solar_row]
n_plans = length(plans)
plan_capacity = zeros(n_plans, length(stack_types))
for i in 1:n_plans
    for j in 1:length(stack_types)
        plan_capacity[i, j] = plans[i].build[stack_types[j]]
    end
end

# Two bars for each gas price, with a gap between prices
bar_positions = [1, 2, 3.5, 4.5, 6, 7]
p_prices = plot(ylabel="Capacity  [MW]", ylims=(0, 8000), xlims=(0.3, 7.7), legend=:outerright,
    xticks=(bar_positions, ["thermal", "all", "thermal", "all", "thermal", "all"]))
bar_bottom = zeros(n_plans)
for j in 1:length(stack_types)
    g = stack_types[j]
    bar_top = bar_bottom .+ plan_capacity[:, j]
    bar!(p_prices, bar_positions, bar_top, fillrange=bar_bottom, bar_width=0.7, linewidth=0,
        color=type_colors[g], label=gens[g, :Plant])
    bar_bottom = bar_top
end
# Each plan's emissions above its bar, and the gas price above each pair
for i in 1:n_plans
    annotate!(p_prices, bar_positions[i], bar_bottom[i] + 330,
        text("$(commas(plans[i].emissions / 1000)) kt", 14, plot_font))
end
annotate!(p_prices, 1.5, 7250, text("\$3.50 gas", 16, plot_font))
annotate!(p_prices, 4, 7250, text("\$7.21 gas", 16, plot_font))
annotate!(p_prices, 6.5, 7250, text("\$15 gas", 16, plot_font))
plot!(p_prices, size=(1200, 290), left_margin=12mm)

A CCGT and a CT cost $25 and $35 per MWh to run at $3.50 gas, $97 and $140 at $15.

With thermal types only, $15 gas makes coal the largest plant and CO₂ (above each bar) 2.2 times the base case.

Each bar is a new solve with only the gas price changed. The six plans cost $603M, $570M, $960M, $729M, $1,253M and $922M a year, left to right. At $3.50 the plans look like the base case with less wind and solar, because the fuel they save is worth less. At $15 the thermal-only plan changes character: coal, which no plan has built so far, becomes the largest plant at 1,720 MW, and emissions more than double. A high gas price with no limit on carbon raises emissions here. With renewables the high price does the opposite. Geothermal, which wind and solar kept out at $7.21, is now built (1,155 MW), because the gas plants that covered the evenings have become expensive to run, and emissions fall to 613 thousand tonnes.

Coal is worth building for load that lasts all year once gas costs more than about $9.5 per million Btu. That is the price at which coal’s line on the screening chart meets the CCGT’s at the end of the year.

Base, Intermediate and Peak Load

# The thermal-only plan at $15 gas, and the two crossovers of its screening curves [hours]
coal_vs_ccgt_europe = (gens[2, :FixedCost] - gens[3, :FixedCost]) / (europe_costs[3] - europe_costs[2])
ct_vs_ccgt_europe = (gens[3, :FixedCost] - gens[4, :FixedCost]) / (europe_costs[4] - europe_costs[3])
p_three = band_plot(thermal_europe.build, [1, 2, 3, 4])
plot!(p_three, size=(1200, 280), left_margin=12mm)

• Coal serves the base load, lasting over 4,429 hours.
• CCGT serves the intermediate load, 395 to 4,429 hours.
• CT serves the peak load, under 395 hours.

The thermal-only plan at $15 gas, drawn the same way as the first plan. Now there are three bands, and each crossover of the screening curves marks where one ends: coal costs $310,000 a year more than a CCGT to build and $70 less per MWh to run, so it wins for load lasting more than 4,429 hours. Utilization falls from band to band: 95% for coal, 23% for CCGT and 1% for CT. A plant that is expensive to build and cheap to run is worth it only if it runs most of the year; that is a baseload plant. A peaker is the reverse. At the 2022 U.S. price and at today’s, a CCGT is both the baseload and the intermediate plant, and coal is not built.

Key Takeaways and Upcoming Schedule

Key Takeaways

• Capacity expansion chooses what to build and how to run it together, at least total cost.
• Screening curves read a thermal build off the load duration curve: each type serves the hours in which its total cost per MW is lowest.
• The value of wind and solar depends on when they produce.
• A CO₂ cap’s shadow price is a carbon price: the cost of abating the last tonne.
• The plan depends on fuel prices: expensive gas moves the base load to other plants.

Next Classes

Wednesday, November 4: Quiz 4, then mixed integer programming.

Monday, November 9: Solid waste and network models.

Wednesday, November 11: Lab 3, fixed costs and facility siting in JuMP.

Assessments

Mini-Project 2: Due Thursday, November 19. Its final part (build decisions and a CO₂ cap) can now be done.

Homework 7: Assigned Wednesday, due Thursday, November 12.

Quiz 4: Wednesday, November 4.

Appendix

A Week With Wind and Solar

# The week around the year's peak hour, under the plan with renewables
peak_hour = argmax(demand)
week = max(1, peak_hour - 84):min(n_hours, peak_hour + 83)
built = findall(with_renewables.build .> 0.5)
week_mix = zeros(length(week), length(built) + 1)
for (j, g) in enumerate(built)
    week_mix[:, j] = with_renewables.hourly[g, week]
end
week_mix[:, end] = with_renewables.unserved[week]
areaplot(1:length(week), week_mix, label=permutedims([gens[built, :Plant]; "Unserved"]),
    color=permutedims([type_colors[built]; :black]),
    xlabel="Hour of the peak week", ylabel="Generation  [MW]", legend=:outerright,
    size=(1200, 350), left_margin=12mm)

Gas still covers the evenings, so its capacity falls much less than its energy.

One week, centred on the year’s peak hour. Solar fills the middle of each day and wind comes and goes. Every evening the gas plants climb back to nearly the output they would have without wind and solar. This is the duck curve, seen from the planner’s side: the capacity is still needed, but it runs fewer hours.

Hourly Prices

hourly_price = round.(Int, [-shadow_price(load[t]) for t in 1:n_hours])
# Hours whose price is not a plant's running cost or the penalty
between_hours = count((hourly_price .!= gens[3, :VarCost]) .& (hourly_price .!= gens[4, :VarCost]) .&
    (hourly_price .!= NSECost))
println("| Price (\\\$/MWh) | Hours |")
println("|--:|--:|")
for p in sort(unique(hourly_price))
    println("| $(commas(p)) | $(commas(count(==(p), hourly_price))) |")
end

Price ($/MWh)	Hours
48	7,974
59	1
69	801
2,552	1
9,000	7

Scarcity hours ($9,000) pay back the CTs’ fixed costs.

The price each hour is the load constraint’s shadow price, as in dispatch. Most hours the CCGT or CT is marginal. In the unserved hours the price is the penalty itself. In 2 hours the price sits in between: those hours are degenerate (one type’s capacity exactly meets demand), so the shadow price there isn’t unique. A CT earns its fixed cost back only in the scarcity hours: 7.28 hours times $8,931 of margin is $65,000, its annual fixed cost.