Introduction

In this project, an energy hub is created to provide electricity, heating and cooling to a residential area. The energy hub consists of different types of energy sources, energy storage and energy conversion devices. It is connected to the grid and the residential area and designed to minimize the total cost while satisfying the energy demand of the residential area. In this project, the energy hub is modeled as a linear programming problem and solved using the COIN-OR Branch and Cut solver (CBC). The restrictions of the energy hub are built on Python.

Modeling

In this section, the energy hub is modeled as a linear programming problem. The order of introduction will align with the order of implementation in Python.

The available devices are defined and divided into two groups: the conversion devices $conv_i$ and the storage devices $stor_i$. For conversion devices, each device $conv_i$ has a base capacity $c_{base,conv_i}$ and a usage quantity $N_{conv_i}$. The maximum capacity of $conv_i$ is defined as,

$$ Cap^{power}{conv_i} = c{base,conv_i} \times N_{conv_i} $$

In the following analysis, the four typical week are merged into a continuous time $t$, where each $t$ represents an hour to be processed. Assume at the $t$-th hour, the conversion device $conv_i$ has a power input $P^{in}_t$ and a power output $P^{out}_t$. The device should satisfy the following constraints,

$$0 \leq P^{in}_t \leq Cap^{power} $$

if the device has multiple outputs,

$$ P^{out,elec}_t = \eta^{elec} P^{in}_t$$

$$ P^{out,heat}_t = \eta^{heat} P^{in}_t$$

if the device only has one output,

$$ P^{out}_t = \eta P^{in}_t$$

For storage devices, each device $stor_i$ has a base power capacity $c_{base,stor_i}$, a stored energy capacity $c_{stor,stor_i}$ and a usage quantity $N_{stor_i}$. The maximum capacity and maximum stored energy of $stor_i$ is defined as,

$$Cap^{power}{stor_i} = c{base,stor_i} \times N_{stor_i}$$

$$ Cap^{stor}{stor_i} = c{stor,stor_i} \times N_{stor_i}$$

The $i$-th storage device $stor_i$ has a stored energy $E_{stor_i,t}$, charging power $P^{in}_{t}$ and discharging power $P^{out}_t$. The charging and discharging efficiency is $\eta$. Besides, $F^{in}_{stor_i,t}=1$ if the device is charging, while $F^{out}_{stor_i,t} = 1$ if the device is discharging. The storage device should satisfy the following constraints,

$$ E_{stor_i,t} = 0, t=0,167,335,503,671$$

$$E_{stor_i,t}=E_{stor_i,t-1}+\eta P^{in}_t-P^{out}_t/\eta$$

$$ 0 \leq E_{stor_i,t} \leq Cap^{stor}_{stor_i}$$

$$0 \leq P^{in}t \leq Cap^{power}{stor_i}, Cap^{power}_{stor_i}$$

$$F^{in}{stor_i,t}+F^{out}{stor_i,t} \leq 1$$

where Eq. \eqref{eq4} is to ensure the stored energy does not diminish after the operation of one day. Eq. \eqref{eq5} is the transition equation of the stored energy. Eq. \eqref{eq6} to Eq. \eqref{eq8} are the constraints of the storage devices. Eq. \eqref{eq9} is to ensure the device do not charge and discharge at the same time.

At each $t$-th hour, the residential area has a demand of electricity $D^{elec}_t$, heating $D^{heat}_t$ and cooling $D^{cool}_t$. The energy hub should satisfy the following constraints, if the supply is not satisfying the requirements, there will be a shedding load $S_t$.

Let the total generation from all electricity-producing devices be $Elec^{gen}_t$, the total discharging electricity from all storage devices be $Elec^{dis}_t$ and total charging electricity from all storage devices be $Elec^{cha}_t$, all the electricity input in conversion devices are $Elec^{in}_t$. And the electricity purchased from the grid is $Elec^{grid}_t$. The electricity balance is,

$$Elec^{grid}_t + Elec^{gen}_t + Elec^{dis}_t = D^{elec}_t + Elec^{cha}_t + Elec^{in}_t - S^{elec}_t$$

$$Heat^{gen}_t + Heat^{dis}_t = D^{heat}_t + Heat^{cha}_t + Heat^{in}_t - S^{heat}_t$$

$$Cool^{gen}_t + Cool^{dis}_t = D^{cool}_t + Cool^{cha}_t + Cool^{in}_t - S^{cool}_t$$

The gas consumption of the energy hub is $Gas^{in}_t$, the gas purchased from the grid is $Gas^{grid}_t$. The gas balance is,

$$ Gas^{grid}_t = Gas^{in}_t $$

The operation cost of the energy hub is,
$C^{OP} = \sum_{t=0}^{671} (C^{OP,gas}_t + C^{OP,elec}_t + C^{shed}_t)$, where $C^{OP,gas}_t$ is the gas producing cost, $C^{OP,elec}_t$ is the electricity producing cost and $C^{shed}_t=S^{elec}_t+S^{heat}_t+S^{cool}_t$ is the shedding cost.

If one device is selected, the cost of the device is $C^{cap}_i$, and the annual investment cost is $C^{cap} = \sum_{i} C^{cap}_i N_i$.

Finally, the objective function is to minimize the total cost,

$$\min \frac{365}{6*7} C^{cap} + C^{OP}$$

With all the constraints and the objective function mentioned in Eq. \eqref{eq1} to Eq. \eqref{final}, the energy hub planning is modeled as a linear programming problem and solved using the COIN-OR Branch and Cut solver (CBC) on Python.

Experiment Results

The simulation results are shown on appendix. As the results show, the energy hub can provide electricity, heating and cooling to the residential area in all four seasons. The consumption of the residential area is satisfied by the energy hub. But the energy hub always needs to import electricity and gas from the grid, which is the main cost of the energy hub. The lithium battery is the most used storage device in the energy hub, which is used to store electricity in the night and discharge in the day. The internal combustion engine is the most used conversion device in the energy hub, which is used to generate electricity and heat in the peak hours.

It is worth noting that since the internal combustion engine cannot be used solely for generating electricity and the generated heat must also be consumed, the internal combustion engine cannot be used for electricity generation during off-peak hours.

As for the cooling demand, the majority of it is satisfied through the charging and discharging of the ice storage. This is likely because the cooling demand is highly concentrated in certain periods of the day.

Conclusion

In this project, an energy hub is created to provide electricity, heating and cooling to a residential area. The energy hub is modeled as a linear programming problem and solved using the COIN-OR Branch and Cut solver. As the results show, the plan generate by the solver can satisfy the energy demand of the residential area in all four seasons. However, due to the CBC only being able to find a local optimum within a limited time, this solution still can be improved.

Code

import numpy as np
import pandas as pd
from pulp import LpProblem, LpMinimize, LpVariable, lpSum, LpStatus, value, PULP_CBC_CMD
import csv

file_path = 'load_price_data.xlsx'
data = pd.ExcelFile(file_path)
seasons = ['Season1(Winter)', 'Season2(Spring)', 'Season3(Summer)', 'Season4(Autumn)']

hourly_demand = []
hourly_prices = []

for season in seasons:
    season_data = data.parse(season)
    hourly_demand.append(season_data[['elec_load(MW)', 'heating_load(MW)', 'cooling_load(MW)']].values)
    hourly_prices.append(season_data[['elec_price(HKD/MWh)', 'gas_price(HKD/m^3)']].values)

hourly_demand = np.vstack(hourly_demand)
hourly_prices = np.vstack(hourly_prices)

gas_conversion_factor = 100
hourly_gas_prices_mwh = hourly_prices[:, 1] * gas_conversion_factor 

devices = [
    # 以下为能源转换设备
    {"name": "CHP A", "input": "Gas", "output": ["Electricity", "Heat"], "efficiency": [0.22, 0.5], "cost": 8.5, "lifetime": 30, "base_capacity": 2},
    {"name": "CHP B", "input": "Gas", "output": ["Electricity", "Heat"], "efficiency": [0.25, 0.55], "cost": 7.9, "lifetime": 30, "base_capacity": 8},
    {"name": "Internal Combustion Engine", "input": "Gas", "output": ["Electricity", "Heat"], "efficiency": [0.35, 0.5], "cost": 4.4, "lifetime": 30, "base_capacity": 8},
    {"name": "Electric Boiler", "input": "Electricity", "output": ["Heat"], "efficiency": [0.9], "cost": 0.6, "lifetime": 20, "base_capacity": 2},
    {"name": "Heat Pump A", "input": "Electricity", "output": ["Heat"], "efficiency": [5], "cost": 3.8, "lifetime": 20, "base_capacity": 2},
    {"name": "Heat Pump B", "input": "Electricity", "output": ["Heat"], "efficiency": [6], "cost": 3.5, "lifetime": 20, "base_capacity": 10},
    {"name": "Compression Refrigerator A", "input": "Electricity", "output": ["Cooling"], "efficiency": [4], "cost": 0.4, "lifetime": 10, "base_capacity": 0.5},
    {"name": "Compression Refrigerator B", "input": "Electricity", "output": ["Cooling"], "efficiency": [5], "cost": 0.6, "lifetime": 20, "base_capacity": 12},
    {"name": "Water Absorption Refrigerator", "input": "Heat", "output": ["Cooling"], "efficiency": [0.8], "cost": 0.6, "lifetime": 20, "base_capacity": 12},
    {"name": "Gas Boiler", "input": "Gas", "output": ["Heat"], "efficiency": [0.95], "cost": 0.4, "lifetime": 30, "base_capacity": 2},
    # 以下为储能设备
    {"name": "Lithium Battery", "input": "Electricity", "output": ["Electricity"], "efficiency": [0.92], "cost": 4.0, "lifetime": 10, "base_capacity": 1, "storage_capacity": 2},
    {"name": "Hot Water Storage", "input": "Heat", "output": ["Heat"], "efficiency": [0.95], "cost": 0.09, "lifetime": 20, "base_capacity": 8, "storage_capacity": 40},
    {"name": "Ice Storage", "input": "Cooling", "output": ["Cooling"], "efficiency": [0.95], "cost": 0.09, "lifetime": 20, "base_capacity": 8, "storage_capacity": 40}
]

devices_dict = {device['name']: device for device in devices}

prob = LpProblem("EnergyHubOptimization", LpMinimize)

conversion_devices = [device for device in devices if 'storage_capacity' not in device]
storage_devices = [device for device in devices if 'storage_capacity' in device]

device_capacity_vars = {}
device_num_vars = {}
max_units = 50 

for device in conversion_devices:
    num_var = LpVariable(f"{device['name']}_num", lowBound=0, upBound=max_units, cat='Integer')
    capacity_var = LpVariable(f"{device['name']}_capacity", lowBound=0, upBound=device['base_capacity'] * max_units)
    device_capacity_vars[device['name']] = capacity_var
    device_num_vars[device['name']] = num_var
    prob += capacity_var == device['base_capacity'] * num_var, f"{device['name']}_Capacity_Definition"

# 储能设备变量设置
storage_charge_vars, storage_discharge_vars, storage_soc_vars = {}, {}, {}
storage_power_capacity_vars, storage_energy_capacity_vars = {}, {}
storage_eta = {}
storage_charge_status = {}
storage_discharge_status = {}

for device in storage_devices:
    name = device['name']
    storage_charge_vars[name] = [LpVariable(f"{name}_charge_{t}", lowBound=0) for t in range(672)]
    storage_discharge_vars[name] = [LpVariable(f"{name}_discharge_{t}", lowBound=0) for t in range(672)]
    storage_soc_vars[name] = [LpVariable(f"{name}_soc_{t}", lowBound=0) for t in range(672)]

    # 分别定义充放电功率和储能容量变量，并设定上限
    num_var = LpVariable(f"{name}_num", lowBound=0, upBound=max_units, cat='Integer')
    storage_power_capacity_vars[name] = LpVariable(f"{name}_power_capacity", lowBound=0, upBound=device['base_capacity'] * max_units)
    storage_energy_capacity_vars[name] = LpVariable(f"{name}_energy_capacity", lowBound=0, upBound=device['storage_capacity'] * max_units)

    storage_eta[name] = device['efficiency'][0]
    device_num_vars[name] = num_var

    prob += storage_power_capacity_vars[name] == device['base_capacity'] * num_var, f"{name}_Power_Capacity_Definition"
    prob += storage_energy_capacity_vars[name] == device['storage_capacity'] * num_var, f"{name}_Energy_Capacity_Definition"
    device_capacity_vars[name] = storage_energy_capacity_vars[name]  # 储能设备的容量变量

    storage_charge_status[name] = [LpVariable(f"{name}_charge_status_{t}", cat='Binary') for t in range(672)]
    storage_discharge_status[name] = [LpVariable(f"{name}_discharge_status_{t}", cat='Binary') for t in range(672)]

    device['P_max'] = device['base_capacity'] * max_units
    device['E_max'] = device['storage_capacity'] * max_units

# 储能设备的约束
for device in storage_devices:
    name = device['name']
    eta = storage_eta[name]
    P_max = device['P_max']
    E_max = device['E_max']

    for t in range(672):
        if t == 0:
            prob += (
                storage_soc_vars[name][t] == storage_charge_vars[name][t] * eta - storage_discharge_vars[name][t] * (1 / eta),
                f"{name}_SOC_{t}"
            )
        else:
            prob += (
                storage_soc_vars[name][t] == storage_soc_vars[name][t-1]
                + storage_charge_vars[name][t] * eta
                - storage_discharge_vars[name][t] * (1 / eta),
                f"{name}_SOC_{t}"
            )
        prob += (
            storage_soc_vars[name][t] <= storage_energy_capacity_vars[name],
            f"{name}_SOC_Capacity_{t}"
        )
        prob += (
            storage_soc_vars[name][t] >= 0,
            f"{name}_SOC_NonNegative_{t}"
        )
        prob += (
            storage_charge_vars[name][t] <= storage_power_capacity_vars[name],
            f"{name}_Charge_Power_Capacity_Limit_{t}"
        )
        prob += (
            storage_charge_vars[name][t] <= P_max * storage_charge_status[name][t],
            f"{name}_Charge_Power_Status_Limit_{t}"
        )
        prob += (
            storage_discharge_vars[name][t] <= storage_power_capacity_vars[name],
            f"{name}_Discharge_Power_Capacity_Limit_{t}"
        )
        prob += (
            storage_discharge_vars[name][t] <= P_max * storage_discharge_status[name][t],
            f"{name}_Discharge_Power_Status_Limit_{t}"
        )
        # 防止同时充放电
        prob += (
            storage_charge_status[name][t] + storage_discharge_status[name][t] <= 1,
            f"{name}_No_Simultaneous_Charge_Discharge_{t}"
        )

    # 每周结束时的SOC清零约束
    for t in [167, 335, 503, 671]:
        prob += (
            storage_soc_vars[name][t] == 0,
            f"{name}_End_of_Week_SOC_Reset_{t}"
        )

device_input_vars = {}
device_output_vars = {}
device_operation_vars = {}
max_capacity = {} 

for device in conversion_devices:
    device_name = device['name']
    device_input_vars[device_name] = [LpVariable(f"{device_name}_input_{t}", lowBound=0) for t in range(672)]
    device_output_vars[device_name] = {output: [LpVariable(f"{device_name}_{output}_output_{t}", lowBound=0) for t in range(672)] for output in device['output']}
    device_operation_vars[device_name] = [LpVariable(f"{device_name}_operating_{t}", cat='Binary') for t in range(672)]
    max_capacity[device_name] = device['base_capacity'] * max_units

electricity_purchase = [LpVariable(f"Electricity_Purchase_{t}", lowBound=0) for t in range(672)]
gas_purchase = [LpVariable(f"Gas_Purchase_{t}", lowBound=0) for t in range(672)]

# 添加设备的输入、输出和运行状态约束
for device in conversion_devices:
    device_name = device['name']
    M = max_capacity[device_name] 
    for t in range(672):
        prob += (
            device_input_vars[device_name][t] <= device_capacity_vars[device_name],
            f"{device_name}_Capacity_Limit_{t}"
        )
        prob += (
            device_input_vars[device_name][t] <= M * device_operation_vars[device_name][t],
            f"{device_name}_Operation_Limit_{t}"
        )

        if 'Electricity' in device['output'] and 'Heat' in device['output']:
            prob += (
                device_output_vars[device_name]['Electricity'][t] == device_input_vars[device_name][t] * device['efficiency'][0],
                f"{device_name}_Electricity_Output_{t}"
            )
            prob += (
                device_output_vars[device_name]['Heat'][t] == device_input_vars[device_name][t] * device['efficiency'][1],
                f"{device_name}_Heat_Output_{t}"
            )
        else:
            for j, output in enumerate(device['output']):
                prob += (
                    device_output_vars[device_name][output][t] == device_input_vars[device_name][t] * device['efficiency'][j],
                    f"{device_name}_{output}_Output_{t}"
                )

shed_load_vars = {energy_type: [LpVariable(f"shed_{energy_type}_{t}", lowBound=0) for t in range(672)] for energy_type in ['Electricity', 'Heat', 'Cooling']}

# 供需平衡约束
for t in range(672):
    elec_demand, heat_demand, cool_demand = hourly_demand[t]

    total_electricity_generation = lpSum(
        device_output_vars[device['name']]['Electricity'][t]
        for device in conversion_devices if 'Electricity' in device['output']
    )
    total_electricity_storage_discharge = lpSum(
        storage_discharge_vars[storage['name']][t]
        for storage in storage_devices if storage['output'][0] == 'Electricity'
    )
    total_electricity_storage_charge = lpSum(
        storage_charge_vars[storage['name']][t]
        for storage in storage_devices if storage['input'] == 'Electricity'
    )
    total_electricity_input_in_devices = lpSum(
        device_input_vars[device_name][t]
        for device_name in device_input_vars
        if devices_dict[device_name]['input'] == 'Electricity'
    )

    prob += (
        electricity_purchase[t] + total_electricity_generation + total_electricity_storage_discharge
        == elec_demand - shed_load_vars['Electricity'][t] + total_electricity_input_in_devices + total_electricity_storage_charge,
        f"Electricity_Balance_{t}"
    )

    total_heat_generation = lpSum(
        device_output_vars[device['name']]['Heat'][t]
        for device in conversion_devices if 'Heat' in device['output']
    )
    total_heat_storage_discharge = lpSum(
        storage_discharge_vars[storage['name']][t]
        for storage in storage_devices if storage['output'][0] == 'Heat'
    )
    total_heat_storage_charge = lpSum(
        storage_charge_vars[storage['name']][t]
        for storage in storage_devices if storage['input'] == 'Heat'
    )
    total_heat_input_in_devices = lpSum(
        device_input_vars[device_name][t]
        for device_name in device_input_vars
        if devices_dict[device_name]['input'] == 'Heat'
    )

    prob += (
        total_heat_generation + total_heat_storage_discharge
        == heat_demand - shed_load_vars['Heat'][t] + total_heat_input_in_devices + total_heat_storage_charge,
        f"Heat_Balance_{t}"
    )

    total_cooling_generation = lpSum(
        device_output_vars[device['name']]['Cooling'][t]
        for device in conversion_devices if 'Cooling' in device['output']
    )
    total_cooling_storage_discharge = lpSum(
        storage_discharge_vars[storage['name']][t]
        for storage in storage_devices if storage['output'][0] == 'Cooling'
    )
    total_cooling_storage_charge = lpSum(
        storage_charge_vars[storage['name']][t]
        for storage in storage_devices if storage['input'] == 'Cooling'
    )
    total_cooling_input_in_devices = lpSum(
        device_input_vars[device_name][t]
        for device_name in device_input_vars
        if devices_dict[device_name]['input'] == 'Cooling'
    )

    prob += (
        total_cooling_generation + total_cooling_storage_discharge
        == cool_demand - shed_load_vars['Cooling'][t] + total_cooling_input_in_devices + total_cooling_storage_charge,
        f"Cooling_Balance_{t}"
    )

    total_gas_consumption = lpSum(
        device_input_vars[device_name][t]
        for device_name in device_input_vars
        if devices_dict[device_name]['input'] == 'Gas'
    )
    prob += (
        gas_purchase[t] == total_gas_consumption,
        f"Gas_Balance_{t}"
    )

interest_rate = 0.04
total_investment_cost = lpSum(
    (interest_rate / (1 - (1 + interest_rate) ** -device['lifetime'])) * device['cost'] * device_num_vars[device['name']] 
    for device in conversion_devices
) * 1e6 + lpSum(
    (interest_rate / (1 - (1 + interest_rate) ** -device['lifetime'])) * device['cost'] * device_num_vars[device['name']]
    for device in storage_devices
) * 1e6  

shed_penalty_cost = 1e6
total_shed_cost = lpSum(
    shed_penalty_cost * (shed_load_vars['Electricity'][t] + shed_load_vars['Heat'][t] + shed_load_vars['Cooling'][t])
    for t in range(672)
)

variable_operation_costs = lpSum(
    electricity_purchase[t] * hourly_prices[t, 0] + gas_purchase[t] * hourly_gas_prices_mwh[t]
    for t in range(672)
)

prob += total_investment_cost + (variable_operation_costs + total_shed_cost) * 365 / (4 * 7), "Total_Cost"
solver = PULP_CBC_CMD(timeLimit=1200)

prob.solve(solver)


print("优化状态:", LpStatus[prob.status])
if prob.status == 1:
    for device in conversion_devices + storage_devices:
        num_devices = int(round(value(device_num_vars[device['name']])))
        print(f"{device['name']}: {num_devices}")

    total_cost = value(total_investment_cost + variable_operation_costs * 365 / (4 * 7) + total_shed_cost)
    print(f"\n总年化投资成本: {value(total_investment_cost):.2f} HKD")
    print(f"运行成本: {value(variable_operation_costs * 365 / (4 * 7)):.2f} HKD")
    print(f"\n总成本 (百万 HKD): {total_cost / 1e6:.2f}")