Endogenous Learning Curves in Multi-Horizon Dynamic Investment Optimisation

Consider a long-term multi-year investment problem where CSP (Concentrated Solar Power) has a learning curve such that

$$LCOE = c_0 \left(\frac{x_t}{x_0}\right)^{-\gamma} + c_1$$

where $c_0$ is cost at start, $c_1$ is material cost and $x_t$ is cumulative capacity in the investment interval $t$. Thus, $x_0$ is the initial cumulative CSP capacity.

Additionally, there are nuclear and coal generators for which there is no potential for reducing their LCOE.

We build an optimisation to minimise the cost of supplying a flat demand $d=100$ GW with the given technologies between 2020 and 2050, where a CO$_2$ budget cap is applied.

Hint: Problem formulation is to be found further along this notebook.

Task: Explore different discount rates, learning rates, CO$_2$ budgets. For instance

  • No learning for CSP and no CO$_2$ budget would result in a coal-reliant system.
  • A CO$_2$ budget and no learning prefers a system built on nuclear.
  • A CO$_2$ budget and learning results in a system with CSP.

NB The learning curve coupling makes the problem non-linear, so you need to install the non-linear interior-point solver ipopt:

conda install -c conda-forge ipopt

Licence

Copyright 2019 Tom Brown (KIT)

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

Imports

In [1]:
from pyomo.environ import ConcreteModel, Var, Objective, NonNegativeReals, Constraint, Suffix, exp
from pyomo.opt import SolverFactory
import pandas as pd
import matplotlib.pyplot as plt
plt.style.use('bmh')
%matplotlib inline

Parameters

In [2]:
techs = ["coal","nuclear","CSP"]
colors = ["#707070","#ff9000","#f9d002"]
parameters = pd.DataFrame(columns=techs)
parameters.loc["current annuity"] = [15.*8760,65.*8760,150.*8760] # EUR/MW/a
parameters.loc["potential annuity"] = [15.*8760,65.*8760,35.*8760] # EUR/MW/a
parameters.loc["learning parameter"] = [0.,0.,0.333]
parameters.loc["marginal cost"] = [35.,10.,0.] #EUR/MWhel
parameters.loc["specific emissions"] = [1.,0.,0.] #tcO2/MWhel
parameters.loc["lifetime"] = [40,40,30]  #years
parameters.loc["existing age"] = [20,0,0] #years
parameters.loc["existing capacity"] = [100,0,0] #GW

parameters.loc["current LCOE"] = parameters.loc["current annuity"]/8760 + parameters.loc["marginal cost"]
parameters.loc["potential LCOE"] = parameters.loc["potential annuity"]/8760 + parameters.loc["marginal cost"]

parameters
Out[2]:
coal nuclear CSP
current annuity 131400.0 569400.0 1314000.000
potential annuity 131400.0 569400.0 306600.000
learning parameter 0.0 0.0 0.333
marginal cost 35.0 10.0 0.000
specific emissions 1.0 0.0 0.000
lifetime 40.0 40.0 30.000
existing age 20.0 0.0 0.000
existing capacity 100.0 0.0 0.000
current LCOE 50.0 75.0 150.000
potential LCOE 50.0 75.0 35.000
In [70]:
#discount rate
rate = 0.05

#demand in GW
demand = 100.

# considered years
years = list(range(2020,2070))


scenario = "no_co2-no_learning"
scenario = "co2-0p2-no_learning"
scenario = "co2-0p2-learning"


if "no_learning" in scenario:
    parameters.loc["learning parameter"] = 0
else:
    parameters.at["learning parameter","CSP"] = 0.333

    
# carbon budget in average tCO2/MWh_el    
if "no_co2" in scenario:
    co2_budget = 2.
else:
    co2_budget = 0.2

Build Model

Note: We use pyomo for building optimisation problems in python. This is also what pypsa uses under the hood.

In [71]:
model = ConcreteModel("discounted total costs")

Generator capacity available for tech $s$ in year $a$ $$G_{s,a}$$

In [72]:
model.generators = Var(techs, years, within=NonNegativeReals)

Generator dispatch for tech $s$ in year $a$ $$g_{s,a}$$

In [73]:
model.generators_dispatch = Var(techs, years, within=NonNegativeReals)

New capacity built for tech $s$ in year $a$ $$Q_{s,a}$$

In [74]:
model.generators_built = Var(techs, years, within=NonNegativeReals)
$$c_{s,a}$$
In [75]:
model.fixed_costs = Var(techs, years, within=NonNegativeReals)

The objective is to minimise the system costs:

$$\min \quad \sum_{s\in S, a\in A} \frac{1}{10^6\cdot (1+r)^{a}} \left( o_{s,a} \cdot g_{s,a} \cdot 8760 + \sum_{b} c_{s,b} Q_{s,b} \mathbb{I}(a \geq b) \mathbb{I}(a < b+L_s) \right) $$
In [76]:
# in billion EUR

# annuities from existing generators
# in billion (MW to GW *1e3, then devide by 1e9)
constant =sum(parameters.at["existing capacity",tech]*parameters.at["current annuity",tech]/1e6/(1+rate)**(year-years[0]) for tech in techs for year in years if year < years[0] + parameters.at["lifetime",tech] - parameters.at["existing age",tech])
print(constant)

model.objective = Objective(expr=constant +
                                 sum(model.generators_built[tech,year]*model.fixed_costs[tech,year]/1e6*sum(1/(1+rate)**(yearb-years[0]) for yearb in years if ((yearb>= year) and (yearb < year + parameters.at["lifetime",tech])))
                                     for year in years
                                     for tech in techs) + 
                                 sum(model.generators_dispatch[tech,year]*parameters.at["marginal cost",tech]*8760/1e6/(1+rate)**(year-years[0])
                                     for year in years
                                     for tech in techs))

171.94111609602413

Add a constraint such that demand is met by generator dispatch:

$$\forall a\in A: \quad d = \sum_{s \in S} g_{s,a}$$
In [77]:
def balance_constraint(model, year):
    return demand == sum(model.generators_dispatch[tech,year] for tech in techs)
model.balance_constraint = Constraint(years, rule=balance_constraint)
$$ g_{s,a} \leq G_{s,a} $$
In [78]:
def generator_constraint(model, tech, year):
    return model.generators_dispatch[tech,year] <= model.generators[tech,year]
model.generator_constraint = Constraint(techs, years, rule=generator_constraint)

Add a constraint on carbon dioxide emissions:

$$\sum_{s\in S, a\in A} G_{s,a} \cdot e_{t} \leq \hat{e} \cdot |A| \cdot d$$
In [79]:
def co2_constraint(model):
    return co2_budget*len(years)*demand >= sum(model.generators_dispatch[tech,year]*parameters.at["specific emissions",tech] for tech in techs for year in years)
model.co2_constraint = Constraint(rule=co2_constraint)
In [80]:
def fixed_cost_constraint(model,tech,year):
    if parameters.at["learning parameter",tech] == 0:
        return model.fixed_costs[tech,year] == parameters.at["current annuity",tech]
    else:
        return model.fixed_costs[tech,year] == parameters.at["potential annuity",tech] + (parameters.at["current annuity",tech]-parameters.at["potential annuity",tech])*(1+sum(model.generators[tech,yeart] for yeart in years if yeart < year))**(-parameters.at["learning parameter",tech])
model.fixed_cost_constraint = Constraint(techs, years, rule=fixed_cost_constraint)
In [81]:
def build_years(model,tech,year):
    if year < years[0] + parameters.at["lifetime",tech] - parameters.at["existing age",tech]:
        constant = parameters.at["existing capacity",tech]
    else:
        constant = 0.
    
    return model.generators[tech,year] == constant + sum(model.generators_built[tech,yearb] for yearb in years if ((year>= yearb) and (year < yearb + parameters.at["lifetime",tech])))
model.build_years = Constraint(techs, years, rule=build_years)

Hint: You can print the model formulation with model.pprint()

Solve Model

In [82]:
opt = SolverFactory("ipopt")
In [83]:
results = opt.solve(model,suffixes=["dual"],keepfiles=False)

Results

Optimised cost (in billion euros NPV):

In [84]:
print(model.objective())

1019.7276625042617

The unoptimized cost (where everything is supplied by coal) is:

In [85]:
print(8760*demand*parameters.at["current LCOE","coal"]*len(years)/1e6)

2190.0
In [86]:
dispatch = pd.DataFrame(0.,index=years,columns=techs)
for year in years:
    for tech in techs:
        dispatch.at[year,tech] = model.generators_dispatch[tech,year].value
In [87]:
fig, ax = plt.subplots()

fig.set_size_inches((10,6))

dispatch.plot(kind="area",stacked=True,color=colors,ax=ax,linewidth=0)
ax.set_xlabel("year")
ax.set_ylabel("dispatch [GW]")

fig.tight_layout()

fig.savefig("{}-dispatch.pdf".format(scenario),transparent=True)

Plotting the development of the technology mix of the optimal solution over time:

In [88]:
capacities = pd.DataFrame(0.,index=years,columns=techs)
for year in years:
    for tech in techs:
        capacities.at[year,tech] = model.generators[tech,year].value
In [89]:
fig, ax = plt.subplots()

fig.set_size_inches((10,6))

capacities.plot(kind="area",stacked=True,color=colors,ax=ax,linewidth=0)
ax.set_xlabel("year")
ax.set_ylabel("capacity [GW]")

fig.tight_layout()

fig.savefig("{}-capacity.pdf".format(scenario),transparent=True)
In [90]:
build_years = pd.DataFrame(0.,index=years,columns=techs)
for year in years:
    for tech in techs:
        build_years.at[year,tech] = model.generators_built[tech,year].value
In [91]:
fig, ax = plt.subplots()

fig.set_size_inches((10,6))

build_years.plot(kind="area",stacked=True,color=colors,ax=ax,linewidth=0)
ax.set_xlabel("year")
ax.set_ylabel("new capacity built [GW]")

fig.tight_layout()

fig.savefig("{}-new_capacity.pdf".format(scenario),transparent=True)

Plotting the development of the costs of the technology over time:

In [92]:
costs = pd.DataFrame(0.,index=years,columns=techs)
for year in years:
    for tech in techs:
        costs.at[year,tech] = model.fixed_costs[tech,year].value/8760. + parameters.at["marginal cost",tech]
In [93]:
fig, ax = plt.subplots()

fig.set_size_inches((10,6))

costs.plot(color=colors,ax=ax,linewidth=3)
ax.set_xlabel("year")
ax.set_ylabel("LCOE [EUR/MWh]")
ax.set_ylim([0,160])


fig.tight_layout()

fig.savefig("{}-lcoe.pdf".format(scenario),transparent=True)
In [ ]: