Modular Expansion¶
In this example, we demonstrate how to handle the modular expansion feature in PyPSA. Modular expansion allows you to specify discrete steps for the expansion of components. This is particularly useful for technologies that can only be built in fixed block sizes, such as nuclear generators.
We start by loading the 3-hourly resolved model.energy example and adding a nuclear generator as expansion option (with ramp limits of 1%/h and a minimum part load of 70%). Initially, we allow the nuclear generator to be built in continuous sizes.
import pypsa
from pypsa.costs import annuity
n = pypsa.examples.model_energy()
n.add(
"Generator",
"nuclear",
bus="electricity",
p_nom_extendable=True,
marginal_cost=15,
capital_cost=annuity(0.07, 50) * 8_000_000,
p_min_pu=0.7,
ramp_limit_up=0.03,
ramp_limit_down=0.03,
)
n.optimize(solver_name="gurobi")
n.generators.p_nom_opt
INFO:pypsa.network.io:Imported network 'Model-Energy' has buses, carriers, generators, links, loads, storage_units, stores
/tmp/ipykernel_5193/1646313121.py:18: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable. n.optimize(solver_name="gurobi")
INFO:linopy.model: Solve problem using Gurobi solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io:Writing objective.
Writing constraints.: 0%| | 0/23 [00:00<?, ?it/s]
Writing constraints.: 87%|████████▋ | 20/23 [00:00<00:00, 192.54it/s]
Writing constraints.: 100%|██████████| 23/23 [00:00<00:00, 170.47it/s]
Writing continuous variables.: 0%| | 0/11 [00:00<?, ?it/s]
Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 378.64it/s]
INFO:linopy.io: Writing time: 0.18s
Set parameter WLSAccessID
Set parameter WLSSecret
Set parameter LicenseID to value 2537914
Academic license 2537914 - for non-commercial use only - registered to l.___@tu-berlin.de
Read LP format model from file /tmp/linopy-problem-rui80sqa.lp
Reading time = 0.08 seconds
obj: 75925 rows, 32127 columns, 156259 nonzeros
Set parameter LogToConsole to value 0
Warning: environment still referenced so free is deferred (Continue to use WLS)
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 32127 primals, 75925 duals Objective: 6.54e+09 Solver model: available Solver message: 2
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Generator-p-ramp_limit_up, Generator-p-ramp_limit_down, Link-ext-p-lower, Link-ext-p-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance, Store-energy_balance were not assigned to the network.
name load shedding 10901.160000 wind 4208.294998 solar 5280.709789 nuclear 7539.778273 Name: p_nom_opt, dtype: float64
As optimal capacity for the nuclear generator, we obtain 7,639 MW.
Now, let's assume that the nuclear generator can only be built in discrete steps of 1,000 MW. We can set the p_nom_mod attribute to 1,000 MW, which introduces an integer variable that constraints the optimised capacity to be a multiple of 1,000 MW. To constrain the option space for the integer variable, and help the solver a bit, we can also set the upper limit p_nom_max to 10,000 MW.
This problem is solved as a mixed-integer linear program (MILP), which will take longer to solve than the previous continuous linear program (LP).
n.generators.loc["nuclear", "p_nom_mod"] = 1000
n.generators.loc["nuclear", "p_nom_max"] = 10000
n.optimize(solver_name="gurobi")
n.generators.p_nom_opt
/tmp/ipykernel_5193/3454292476.py:3: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable. n.optimize(solver_name="gurobi")
INFO:linopy.model: Solve problem using Gurobi solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io:Writing objective.
Writing constraints.: 0%| | 0/25 [00:00<?, ?it/s]
Writing constraints.: 88%|████████▊ | 22/25 [00:00<00:00, 209.59it/s]
Writing constraints.: 100%|██████████| 25/25 [00:00<00:00, 187.46it/s]
Writing continuous variables.: 0%| | 0/12 [00:00<?, ?it/s]
Writing continuous variables.: 100%|██████████| 12/12 [00:00<00:00, 410.32it/s]
Writing integer variables.: 0%| | 0/1 [00:00<?, ?it/s]
Writing integer variables.: 100%|██████████| 1/1 [00:00<00:00, 362.11it/s]
INFO:linopy.io: Writing time: 0.18s
Set parameter WLSAccessID
Set parameter WLSSecret
Set parameter LicenseID to value 2537914
Academic license 2537914 - for non-commercial use only - registered to l.___@tu-berlin.de
Read LP format model from file /tmp/linopy-problem-w5aa8zja.lp
Reading time = 0.08 seconds
obj: 75927 rows, 32128 columns, 156262 nonzeros
Set parameter LogToConsole to value 0
WARNING:linopy.solvers:Dual values of MILP couldn't be parsed
Warning: environment still referenced so free is deferred (Continue to use WLS)
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 32128 primals, 0 duals Objective: 6.55e+09 Solver model: available Solver message: 2
INFO:pypsa.optimization.optimize:No shadow prices were assigned to the network.
name load shedding 10901.160000 wind 2978.500684 solar 4365.968538 nuclear 8000.000000 Name: p_nom_opt, dtype: float64
From the re-optimised results, we can see that the optimal capacity for the nuclear generator is now 8,000 MW, a multiple of 1,000 MW.
As a short concluding excursion into the results, we can also illustrate how the optimal dispatch of this nuclear generator is constrained by ramp limits and minimum part loads in January 2019.
n.generators_t.p.loc["2019-01"].plot(figsize=(6, 3))
<Axes: xlabel='snapshot'>
