distributor.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
import control as ct
from dataclasses import dataclass

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


def state_space_model(config: EMSConfig) -> (np.ndarray, np.ndarray, np.ndarray):
    """
    State space model of the aggregator system, i.e. the detailed thermal behavior of the 9 zones

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

    C_thermal_matrix = np.tile(config.C_thermal, (len(config.C_thermal), 1)).T

    A = (-np.diag(config.H_air) + config.beta_coupling) / C_thermal_matrix

    # heating only in zones 1-7
    B_cooling = np.diag(1 / config.C_thermal)
    B_heating = B_cooling[:, :7]

    B = np.hstack((B_heating, B_cooling))

    S_ambient = config.H_air / config.C_thermal
    S_ambient = np.expand_dims(S_ambient, axis=1)   # turn into column vector

    S_Q_other = np.diag(1 / config.C_thermal)
    S = np.hstack((S_ambient, S_Q_other))

    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):
    theta_1: float
    """ Temperature in zone 1 in °C """

    theta_2: float
    """ Temperature in zone 2 in °C """

    theta_3: float
    """ Temperature in zone 3 in °C """

    theta_4: float
    """ Temperature in zone 4 in °C """

    theta_5: float
    """ Temperature in zone 5 in °C """

    theta_6: float
    """ Temperature in zone 6 in °C """

    theta_7: float
    """ Temperature in zone 7 in °C """

    theta_8: float
    """ Temperature in zone 8 in °C """

    theta_9: float
    """ Temperature in zone 9 in °C """


@dataclass(frozen=True)
class Parameters(BaseParameters):
    theta_ref: float | np.ndarray | DataSource
    """ Reference temperature for building zones 1-7 in °C """

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

    Q_heat: float | np.ndarray
    """ Building zone heating budget to be distributed in kW """

    Q_cool_b: float | np.ndarray
    """ Building zone cooling budget to be distributed in kW """

    Q_cool_s: float | np.ndarray
    """ Server zone cooling budget to be distributed in kW """

    epsilon: np.ndarray | DataSource = None
    """ Model error / heat disturbances acting on zones 1-9 """


@dataclass(frozen=True)
class ControllerOutput(MPCControllerOutput, Vectorizable):
    Q_heat_1: float | np.ndarray
    """ Heating power allocated to zone 1 in kW """

    Q_heat_2: float | np.ndarray
    """ Heating power allocated to zone 2 in kW """

    Q_heat_3: float | np.ndarray
    """ Heating power allocated to zone 3 in kW """

    Q_heat_4: float | np.ndarray
    """ Heating power allocated to zone 4 in kW """

    Q_heat_5: float | np.ndarray
    """ Heating power allocated to zone 5 in kW """

    Q_heat_6: float | np.ndarray
    """ Heating power allocated to zone 6 in kW """

    Q_heat_7: float | np.ndarray
    """ Heating power allocated to zone 7 in kW """

    # ---------------------------------------------

    Q_cool_1: float | np.ndarray
    """ Cooling power allocated to zone 1 in kW """

    Q_cool_2: float | np.ndarray
    """ Cooling power allocated to zone 2 in kW """

    Q_cool_3: float | np.ndarray
    """ Cooling power allocated to zone 3 in kW """

    Q_cool_4: float | np.ndarray
    """ Cooling power allocated to zone 4 in kW """

    Q_cool_5: float | np.ndarray
    """ Cooling power allocated to zone 5 in kW """

    Q_cool_6: float | np.ndarray
    """ Cooling power allocated to zone 6 in kW """

    Q_cool_7: float | np.ndarray
    """ Cooling power allocated to zone 7 in kW """

    Q_cool_8: float | np.ndarray
    """ Cooling power allocated to zone 8 in kW """

    Q_cool_9: float | np.ndarray
    """ Cooling power allocated to zone 9 in kW """

    # ---------------------------------------------

    P_cool_1: float | np.ndarray
    """ Cooling power allocated to zone 1 in kW """

    P_cool_2: float | np.ndarray
    """ Cooling power allocated to zone 2 in kW """

    P_cool_3: float | np.ndarray
    """ Cooling power allocated to zone 3 in kW """

    P_cool_4: float | np.ndarray
    """ Cooling power allocated to zone 4 in kW """

    P_cool_5: float | np.ndarray
    """ Cooling power allocated to zone 5 in kW """

    P_cool_6: float | np.ndarray
    """ Cooling power allocated to zone 6 in kW """

    P_cool_7: float | np.ndarray
    """ Cooling power allocated to zone 7 in kW """

    P_cool_8: float | np.ndarray
    """ Cooling power allocated to zone 8 in kW """

    P_cool_9: float | np.ndarray
    """ Cooling power allocated to zone 9 in kW """

    @property
    def Q_heat(self) -> np.ndarray:
        """ Vector containing Q_heat 1 - Q_heat 7"""
        return np.array([getattr(self, f"Q_heat_{i}") for i in range(1, 8)])

    @property
    def Q_cool(self) -> np.ndarray:
        """ Vector containing Q_cool_1 - Q_cool_9 """
        return np.array([getattr(self, f"P_cool_{i}") for i in range(1, 10)])

    @property
    def P_cool(self) -> np.ndarray:
        """ Vector containing P_cool_1 - P_cool_9 """
        return np.array([getattr(self, f"P_cool_{i}") for i in range(1, 10)])


class Model(cosimos.interfaces.model.Model):
    """
    Model of the distributor-only system.
    I.e. 9 temperature zones with heating/cooling in each zone,
    without coupling to the electrical system.
    """

    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(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()
        u = controller_output.to_vector()[:-9]
        d = np.concatenate(([parameters.theta_air], self.config.Q_other))

        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

        state_plus = State.from_vector(theta_plus)

        return state_plus, controller_output