how to add emission related global constraint in pypsa?

Question

I have created a very small program to understand the workings of Pypsa. It has an hourly time interval load profile with random values between 2000MW-3000MW with 1-month total values. To serve this load profile I have added 11 generators:

3 Coal
2 Gas
3 Nuclear
2 Solar
1 Backup

There are 5 basic properties of generators that I have added to the model:

Capacity (MW)
Variable cost (Rupee/MWh)
Ramp-up rate (%)
Ramp-down rate (%)
CO2 Emission rate (Ton CO2/MWh)

I was able to add the manual constraint into the model where I can restrict the overall CO2 emission from the coal and gas plants. This is the code that I have written:

import pypsa
import numpy as np
import pandas as pd
from pyomo.environ import Constraint
from pyomo.environ import value

start_mt = 1
start_yr = 2022
end_mt = 1
end_yr = 2022
end_day = 31
frequency = 60

snapshots = pd.date_range("{}-{}-01".format(start_yr, start_mt), "{}-{}-{} 23:59".format(end_yr, end_mt, end_day),
                          freq=str(frequency) + "min")
np.random.seed(len(snapshots))

# Create a PyPSA network
network = pypsa.Network()

# Add a load bus
network.add("Bus", "Bus")
network.set_snapshots(snapshots)

load_profile = np.random.randint(2000, 3000, len(snapshots))

# Add the load to the network
network.add("Load", "Load profile", bus="Bus", p_set=load_profile)

# Define the generator data dictionary
generator_data = {
    'coal1': {'capacity': 800, 'ramp up': 0.6, 'ramp down': 0.6, 'variable cost': 10, 'co2_emission_factor': 0.95},
    'coal2': {'capacity': 600, 'ramp up': 0.6, 'ramp down': 0.6, 'variable cost': 11, 'co2_emission_factor': 0.95},
    'coal3': {'capacity': 500, 'ramp up': 0.6, 'ramp down': 0.6, 'variable cost': 11, 'co2_emission_factor': 0.95},
    'gas1': {'capacity': 600, 'ramp up': 0.45, 'ramp down': 0.45, 'variable cost': 12, 'co2_emission_factor': 0.45},
    'gas2': {'capacity': 600, 'ramp up': 0.45, 'ramp down': 0.45, 'variable cost': 13, 'co2_emission_factor': 0.45},
    'nuclear1': {'capacity': 300, 'ramp up': 0.01, 'ramp down': 0.01, 'variable cost': 4, 'co2_emission_factor': 0.03},
    'nuclear2': {'capacity': 400, 'ramp up': 0.01, 'ramp down': 0.01, 'variable cost': 3, 'co2_emission_factor': 0.03},
    'nuclear3': {'capacity': 250, 'ramp up': 0.01, 'ramp down': 0.01, 'variable cost': 3, 'co2_emission_factor': 0.03},
    'solar1': {'capacity': 150, 'ramp up': 1, 'ramp down': 1, 'variable cost': 1, 'co2_emission_factor': 0.0},
    'solar2': {'capacity': 200, 'ramp up': 1, 'ramp down': 1, 'variable cost': 2, 'co2_emission_factor': 0.0},
    'backup': {'capacity': 2000, 'ramp up': 1, 'ramp down': 1, 'variable cost': 20, 'co2_emission_factor': 0.0},
}


# Add generators to the network
for name, data in generator_data.items():
    network.add("Generator", name,
                bus="Bus",
                p_nom=data['capacity'],
                marginal_cost=data['variable cost'],
                ramp_limit_up=data['ramp up'],
                ramp_limit_down=data['ramp down'],
                )


def extra_functionality(network, snaphots):
    def co2_limiter(model):
        coalGen = sum([model.generator_p[name, i] for i in list(network.snapshots) for name in ['coal1', 'coal2', 'coal3']])
        gasGen = sum([model.generator_p[name, i] for i in list(network.snapshots) for name in ['gas1', 'gas2']])
        nuclearGen = sum([model.generator_p[name, i] for i in list(network.snapshots) for name in ['nuclear1', 'nuclear2', 'nuclear3']])

        coal_co2_emissions = coalGen * 0.95
        gas_co2_emissions = gasGen * 0.45
        nuclear_co2_emissions = nuclearGen * 0.03
        total_co2 = (coal_co2_emissions + gas_co2_emissions + nuclear_co2_emissions) / 1000000
        expr = total_co2 <= 0.04
        return expr

    network.model.co2_limiter = Constraint(rule=co2_limiter)



solver_name = "gurobi"
network.lopf(network.snapshots, solver_name=solver_name, extra_functionality=extra_functionality)


dispatch = network.generators_t.p
total_gen = dispatch.sum()

co2 = sum([total_gen[gen] * data['co2_emission_factor'] for gen, data in generator_data.items()])
print('CO2 Emission = ',co2)
dispatch['load profile'] = load_profile

dispatch.to_excel('fuel wise dispatch.xlsx')

When you run this program with and without the co2_limiter constraint you ll notice that without this constraint the total CO2 Emission ll be around 860287 Ton CO2/ MWh and when you activate the co2_limiter constraint then the total emission ll be around 399999 Ton CO2/MWh ( because we have set the limit to <= 0.04 Million Ton CO2/MWh) and in dispatch you call also see although backup generator is the costliest one but since it has no co2 emission value it ll use that generator. My question is in Pypsa we do have some internal Global constraints to control this CO2 emission. I have gone through the documentation but still, I was not able to get how to use this in my case. Can anyone please help?

Answer 1

You can indeed use a global constraint for that without writing any extra functionality.

https://pypsa.readthedocs.io/en/latest/components.html#global-constraints

It's an interaction between carriers, which store data on the specific emissions in t/MWh_th, the generators which have an associated carrier, and a global constraint that limits the total amount of emissions, in this example to 100t CO2.

n.add("Carrier", "gas", co2_emissions=0.2)

n.add("Generator", "gas power plant", bus="bus0", carrier="gas", ...)

n.add(
    "GlobalConstraint",
    "CO2Limit",
    carrier_attribute="co2_emissions",
    sense="<=",
    constant=100,
)