aggregator.py

#  MIT License
#
#  Copyright (c) 2025 Honda Research Institute Europe GmbH
#
#  Permission is hereby granted, free of charge, to any person obtaining a copy
#  of this software and associated documentation files (the "Software"), to deal
#  in the Software without restriction, including without limitation the rights
#  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#  copies of the Software, and to permit persons to whom the Software is
#  furnished to do so, subject to the following conditions:
#
#  The above copyright notice and this permission notice shall be included in all
#  copies or substantial portions of the Software.
#
#  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
#  SOFTWARE.
#
#    SPDX-License-Identifier: MIT
#

import numpy as np
from dataclasses import dataclass

import cosimos.interfaces.controller
from cosimos.interfaces.model import State as BaseState
from cosimos.interfaces.parameters import Parameters as BaseParameters
from cosimos.data.sources import DataSource
from cosimos.interfaces.controller import MPCControllerOutput
from cosimos.simulator.history import History
from cosimos.simulator.simulation import Simulation
from ..helpers.control import discretize_state_space_model
from ..config import EMSConfig


def state_space_model_thermal(
    config: EMSConfig,
) -> (np.ndarray, np.ndarray, np.ndarray):
    """
    Aggregated thermal state space model

    States x = [theta_b, theta_s]
    Inputs u = [P_chp, Q_rad, Q_cool_b, Q_cool_s]
    Disturbances d = [theta_a, Q_other_b, Q_other_s]

    :param config: EMS parameters
    :return: discrete-time system matrix A, input matrix B, disturbance matrix S
    """

    # coupling building zone <-> server zones
    beta_bs = np.sum(np.sum(config.beta_coupling[1:7, -2:]))

    # heat transfer coefficient of building
    H_air_b = np.sum(config.H_air[:7])
    C_thermal_b = np.sum(config.C_thermal[:7])

    H_air_s = np.sum(config.H_air[7:])
    C_thermal_s = np.sum(config.C_thermal[7:])

    A = np.array(
        [
            [-(H_air_b + beta_bs) / C_thermal_b, beta_bs / C_thermal_b],
            [beta_bs / C_thermal_s, -(H_air_s + beta_bs) / C_thermal_s],
        ]
    )

    B = np.array(
        [
            [1 / (config.c_chp * C_thermal_b), 1 / C_thermal_b, 1 / C_thermal_b, 0],
            [0, 0, 0, 1 / C_thermal_s],
        ]
    )

    S = np.array(
        [
            [H_air_b / C_thermal_b, 1 / C_thermal_b, 0],
            [H_air_s / C_thermal_s, 0, 1 / C_thermal_s],
        ]
    )

    A_dis, B_dis, S_dis = discretize_state_space_model(
        A, B, S, dt=config.sample_time / 60 / 60
    )

    return A_dis, B_dis, S_dis


@dataclass(frozen=True)
class State(BaseState):
    E_bat: float
    """ Energy stored in the stationary battery in kWh """

    theta_b: float
    """ Temperature in the building zone in °C """

    theta_s: float
    """ Temperature in the server zone in °C """

    P_grid_peak_current: float
    """ Grid peak power demand incurred until current time step """


@dataclass(frozen=True)
class Parameters(BaseParameters):
    P_dem: float | np.ndarray | DataSource  # < 0
    """ Building load demand in kW, counted negatively """

    P_ren: float | np.ndarray | DataSource  # > 0
    """ PV power in kW, counted positively """

    theta_air: float | np.ndarray | DataSource
    """ Ambient temperature in °C """

    epsilon: np.ndarray | DataSource = None
    """ Model error / heat disturbances acting on building and server zone """


@dataclass(frozen=True)
class ControllerOutput(MPCControllerOutput):
    P_bat: float | np.ndarray  # >0 -> charge, <0 -> discharge
    """ (Dis)charge power of stationary battery in kW, >0 -> charge, <0 -> discharge """

    P_grid: float | np.ndarray
    """ Power from/to power grid in kW, >0 -> from grid, <0 -> to grid """

    P_chp: float | np.ndarray
    """ Electrical power generated by the CHP in kW """

    Q_chp: float | np.ndarray
    """ Heating power generated by the CHP in kW """

    Q_rad: float | np.ndarray
    """ Heating power generated by the gas boiler in kW """

    Q_cool_b: float | np.ndarray
    """ Cooling power generated by the building cooling machine in kW, <0 """

    P_cool_b: float | np.ndarray
    """ Electrical power used by the building cooling machine in kW, <0 """

    Q_cool_s: float | np.ndarray
    """ Cooling power generated by the server zone cooling machine in kW, <0 """

    P_cool_s: float | np.ndarray  # < 0
    """ Electrical power used by the server zone cooling machine in kW, <0 """


class Model(cosimos.interfaces.model.Model):
    """
    Model of the aggregator-only system.
    I.e. stationary battery + building zone + server zone
    """

    A_dis: np.ndarray
    B_dis: np.ndarray
    S_dis: np.ndarray
    config: EMSConfig

    def __init__(self, config: EMSConfig):
        self.A_dis, self.B_dis, self.S_dis = state_space_model_thermal(config)
        self.config = config

    def step(
        self,
        time: float,
        state: State,
        parameters: Parameters,
        controller_output: ControllerOutput,
        step_size: float,
        history: History,
        simulation: Simulation,
    ) -> (State, ControllerOutput):
        theta = state.to_vector()[1:3]

        # simulate cooling efficiencies
        if self.config.cooling_control_mode is self.config.THERMAL_CONTROL:
            Q_cool_b = controller_output.Q_cool_b
            P_cool_b = 1 / self.config.eps_c * controller_output.Q_cool_b

            Q_cool_s = controller_output.Q_cool_s
            P_cool_s = 1 / self.config.eps_c * controller_output.Q_cool_s
        else:
            P_cool_b = controller_output.P_cool_b
            Q_cool_b = self.config.eps_c * controller_output.P_cool_b

            P_cool_s = controller_output.P_cool_s
            Q_cool_s = self.config.eps_c * controller_output.P_cool_s

        u = np.array(
            [
                controller_output.P_chp,
                controller_output.Q_rad,
                Q_cool_b,
                Q_cool_s,
            ]
        )
        d = np.array(
            [parameters.theta_air, self.config.Q_other_b, self.config.Q_other_s]
        )

        theta_plus = self.A_dis @ theta + self.B_dis @ u + self.S_dis @ d.T
        if parameters.epsilon is not None:
            theta_plus += parameters.epsilon

        # todo: cap P_bat with what's physically possible

        E_bat_plus = state.E_bat + controller_output.P_bat * (step_size / 60 / 60)


        # update P_grid to be balance of all actual electric powers
        P_grid_actual = (
            controller_output.P_bat
            - controller_output.P_chp
            - P_cool_b
            - P_cool_s
            - parameters.P_ren
            - parameters.P_dem
        )

        state_plus = State(
            E_bat=E_bat_plus,
            theta_b=theta_plus[0],
            theta_s=theta_plus[1],
            P_grid_peak_current=max(state.P_grid_peak_current, P_grid_actual),
        )

        updated_output = controller_output.update(
            P_grid=P_grid_actual,
            Q_cool_b=Q_cool_b,
            P_cool_b=P_cool_b,
            Q_cool_s=Q_cool_s,
            P_cool_s=P_cool_s,
        )

        return state_plus, updated_output