citylearn.py

import gym
from gym.utils import seeding
import numpy as np
import pandas as pd
import json
from gym import spaces
from energy_models import Battery, HeatPump, ElectricHeater, EnergyStorage, Building
from reward_function import reward_function_sa, reward_function_ma
from pathlib import Path

gym.logger.set_level(40)


# Reference Rule-based controller. Used as a baseline to calculate the costs in CityLearn
# It requires, at least, the hour of the day as input state
class RBC_Agent:
    def __init__(self, actions_spaces):
        self.actions_spaces = actions_spaces
        self.action_tracker = []

    def select_action(self, states):
        hour_day = states[0]

        # multiplier = 0.4
        # # Daytime: release stored energy
        # a = [[0.0 for _ in range(len(self.actions_spaces[i].sample()))] for i in range(len(self.actions_spaces))]
        # if 7 <= hour_day <= 11:
        #     a = [[-0.05 * multiplier for _ in range(len(self.actions_spaces[i].sample()))] for i in
        #          range(len(self.actions_spaces))]
        # elif 12 <= hour_day <= 15:
        #     a = [[-0.05 * multiplier for _ in range(len(self.actions_spaces[i].sample()))] for i in
        #          range(len(self.actions_spaces))]
        # elif 16 <= hour_day <= 18:
        #     a = [[-0.11 * multiplier for _ in range(len(self.actions_spaces[i].sample()))] for i in
        #          range(len(self.actions_spaces))]
        # elif 19 <= hour_day <= 22:
        #     a = [[-0.06 * multiplier for _ in range(len(self.actions_spaces[i].sample()))] for i in
        #          range(len(self.actions_spaces))]
        #
        # # Early nightime: store DHW and/or cooling energy
        # if 23 <= hour_day <= 24:
        #     a = [[0.085 * multiplier for _ in range(len(self.actions_spaces[i].sample()))] for i in
        #          range(len(self.actions_spaces))]
        # elif 1 <= hour_day <= 6:
        #     a = [[0.1383 * multiplier for _ in range(len(self.actions_spaces[i].sample()))] for i in
        #          range(len(self.actions_spaces))]

        # Daytime: release stored energy
        a = [[0.0 for _ in range(len(self.actions_spaces[i].sample()))] for i in range(len(self.actions_spaces))]
        if 9 <= hour_day <= 21:
            a = [[-0.08 for _ in range(len(self.actions_spaces[i].sample()))] for i in range(len(self.actions_spaces))]

        # Early nighttime: store DHW and/or cooling energy
        if (1 <= hour_day <= 8) or (22 <= hour_day <= 24):
            a = [[0.091 for _ in range(len(self.actions_spaces[i].sample()))] for i in range(len(self.actions_spaces))]

        return np.array(a, dtype='object')


def auto_size(buildings):
    for building in buildings.values():

        # Autosize guarantees that the DHW device is large enough to always satisfy the maximum DHW demand
        if building.dhw_heating_device.nominal_power == 'autosize':

            # If the DHW device is a HeatPump
            if isinstance(building.dhw_heating_device, HeatPump):

                # We assume that the heat pump is always large enough to meet the highest heating or cooling demand
                # of the building
                building.dhw_heating_device.nominal_power = np.array(
                    building.sim_results['dhw_demand'] / building.dhw_heating_device.cop_heating).max()

            # If the device is an electric heater
            elif isinstance(building.dhw_heating_device, ElectricHeater):
                building.dhw_heating_device.nominal_power = (
                        np.array(building.sim_results['dhw_demand']) / building.dhw_heating_device.efficiency).max()

        # Autosize guarantees that the cooling device is large enough to always satisfy the maximum DHW demand
        if building.cooling_device.nominal_power == 'autosize':
            building.cooling_device.nominal_power = (
                    np.array(building.sim_results['cooling_demand']) / building.cooling_device.cop_cooling).max()

        # Defining the capacity of the storage devices as a number of times the maximum demand
        building.dhw_storage.capacity = max(building.sim_results['dhw_demand']) * building.dhw_storage.capacity
        building.cooling_storage.capacity = max(
            building.sim_results['cooling_demand']) * building.cooling_storage.capacity

        # Done in order to avoid dividing by 0 if the capacity is 0
        if building.dhw_storage.capacity <= 0.00001:
            building.dhw_storage.capacity = 0.00001
        if building.cooling_storage.capacity <= 0.00001:
            building.cooling_storage.capacity = 0.00001


def building_loader(data_path, building_attributes, weather_file, solar_profile, carbon_intensity, electricity_price,
                    building_ids, buildings_states_actions, save_memory=True):
    with open(building_attributes) as json_file:
        data = json.load(json_file)

    buildings, observation_spaces, action_spaces = {}, [], []
    s_low_central_agent, s_high_central_agent, appended_states = [], [], []
    a_low_central_agent, a_high_central_agent, appended_actions = [], [], []
    for uid, attributes in zip(data, data.values()):
        if uid in building_ids:

            battery = Battery(capacity=attributes['Battery']['capacity'],
                              capacity_loss_coef=attributes['Battery']['capacity_loss_coefficient'],
                              loss_coef=attributes['Battery']['loss_coefficient'],
                              efficiency=attributes['Battery']['efficiency'],
                              nominal_power=attributes['Battery']['nominal_power'],
                              power_efficiency_curve=attributes['Battery']['power_efficiency_curve'],
                              capacity_power_curve=attributes['Battery']['capacity_power_curve'],
                              save_memory=save_memory)

            heat_pump = HeatPump(nominal_power=attributes['Heat_Pump']['nominal_power'],
                                 eta_tech=attributes['Heat_Pump']['technical_efficiency'],
                                 t_target_heating=attributes['Heat_Pump']['t_target_heating'],
                                 t_target_cooling=attributes['Heat_Pump']['t_target_cooling'], save_memory=save_memory)

            electric_heater = ElectricHeater(nominal_power=attributes['Electric_Water_Heater']['nominal_power'],
                                             efficiency=attributes['Electric_Water_Heater']['efficiency'],
                                             save_memory=save_memory)

            chilled_water_tank = EnergyStorage(capacity=attributes['Chilled_Water_Tank']['capacity'],
                                               loss_coef=attributes['Chilled_Water_Tank']['loss_coefficient'],
                                               save_memory=save_memory)

            dhw_tank = EnergyStorage(capacity=attributes['DHW_Tank']['capacity'],
                                     loss_coef=attributes['DHW_Tank']['loss_coefficient'], save_memory=save_memory)

            building = Building(buildingId=uid, dhw_storage=dhw_tank, cooling_storage=chilled_water_tank,
                                electrical_storage=battery, dhw_heating_device=electric_heater,
                                cooling_device=heat_pump, save_memory=save_memory)

            data_file = str(uid) + '.csv'
            simulation_data = data_path / data_file
            with open(simulation_data) as csv_file:
                data = pd.read_csv(csv_file)

            building.sim_results['cooling_demand'] = list(data['Cooling Load [kWh]'])
            building.sim_results['dhw_demand'] = list(data['DHW Heating [kWh]'])
            building.sim_results['non_shiftable_load'] = list(data['Equipment Electric Power [kWh]'])
            building.sim_results['month'] = list(data['Month'])
            building.sim_results['day'] = list(data['Day Type'])
            building.sim_results['hour'] = list(data['Hour'])
            building.sim_results['daylight_savings_status'] = list(data['Daylight Savings Status'])
            building.sim_results['t_in'] = list(data['Indoor Temperature [C]'])
            building.sim_results['avg_unmet_setpoint'] = list(data['Average Unmet Cooling Setpoint Difference [C]'])
            building.sim_results['rh_in'] = list(data['Indoor Relative Humidity [%]'])

            with open(weather_file) as csv_file:
                weather_data = pd.read_csv(csv_file)

            building.sim_results['t_out'] = list(weather_data['Outdoor Drybulb Temperature [C]'])
            building.sim_results['rh_out'] = list(weather_data['Outdoor Relative Humidity [%]'])
            building.sim_results['diffuse_solar_rad'] = list(weather_data['Diffuse Solar Radiation [W/m2]'])
            building.sim_results['direct_solar_rad'] = list(weather_data['Direct Solar Radiation [W/m2]'])

            # Reading weather forecasts
            building.sim_results['t_out_pred_6h'] = list(weather_data['6h Prediction Outdoor Drybulb Temperature [C]'])
            building.sim_results['t_out_pred_12h'] = list(
                weather_data['12h Prediction Outdoor Drybulb Temperature [C]'])
            building.sim_results['t_out_pred_24h'] = list(
                weather_data['24h Prediction Outdoor Drybulb Temperature [C]'])

            building.sim_results['rh_out_pred_6h'] = list(weather_data['6h Prediction Outdoor Relative Humidity [%]'])
            building.sim_results['rh_out_pred_12h'] = list(weather_data['12h Prediction Outdoor Relative Humidity [%]'])
            building.sim_results['rh_out_pred_24h'] = list(weather_data['24h Prediction Outdoor Relative Humidity [%]'])

            building.sim_results['diffuse_solar_rad_pred_6h'] = list(
                weather_data['6h Prediction Diffuse Solar Radiation [W/m2]'])
            building.sim_results['diffuse_solar_rad_pred_12h'] = list(
                weather_data['12h Prediction Diffuse Solar Radiation [W/m2]'])
            building.sim_results['diffuse_solar_rad_pred_24h'] = list(
                weather_data['24h Prediction Diffuse Solar Radiation [W/m2]'])

            building.sim_results['direct_solar_rad_pred_6h'] = list(
                weather_data['6h Prediction Direct Solar Radiation [W/m2]'])
            building.sim_results['direct_solar_rad_pred_12h'] = list(
                weather_data['12h Prediction Direct Solar Radiation [W/m2]'])
            building.sim_results['direct_solar_rad_pred_24h'] = list(
                weather_data['24h Prediction Direct Solar Radiation [W/m2]'])

            # Reading the building attributes
            building.building_type = attributes['Building_Type']
            building.climate_zone = attributes['Climate_Zone']
            building.solar_power_capacity = attributes['Solar_Power_Installed(kW)']

            with open(solar_profile) as csv_file:
                data = pd.read_csv(csv_file)

            building.sim_results['solar_gen'] = list(
                attributes['Solar_Power_Installed(kW)'] * data['Hourly Data: AC inverter power (W)'] / 1000)

            with open(carbon_intensity) as csv_file:
                data = pd.read_csv(csv_file)

            building.sim_results['carbon_intensity'] = list(data['kg_CO2/kWh'])

            with open(electricity_price) as csv_file:
                data = pd.read_csv(csv_file)

            building.sim_results['electricity_price'] = list(data['Price'])

            # Finding the max and min possible values of all the states, which can then be used by the RL agent to scale the states and train any function approximators more effectively
            s_low, s_high = [], []
            for state_name, value in zip(buildings_states_actions[uid]['states'],
                                         buildings_states_actions[uid]['states'].values()):
                if value == True:
                    if state_name == "net_electricity_consumption":
                        # lower and upper bounds of net electricity consumption are rough estimates and may not be completely accurate. Scaling this state-variable using these bounds may result in normalized values above 1 or below 0.
                        _net_elec_cons_upper_bound = max(
                            np.array(building.sim_results['non_shiftable_load']) - np.array(
                                building.sim_results['solar_gen']) + np.array(
                                building.sim_results['dhw_demand']) / .8 + np.array(building.sim_results[
                                                                                        'cooling_demand']) + building.dhw_storage.capacity / .8 + building.cooling_storage.capacity / 2)
                        s_low.append(0.)
                        s_high.append(_net_elec_cons_upper_bound)
                        s_low_central_agent.append(0.)
                        s_high_central_agent.append(_net_elec_cons_upper_bound)

                    elif (state_name != 'cooling_storage_soc') and (state_name != 'dhw_storage_soc') and (
                            state_name != 'electrical_storage_soc'):
                        s_low.append(min(building.sim_results[state_name]))
                        s_high.append(max(building.sim_results[state_name]))

                        # Create boundaries of the observation space of a centralized agent (if a central agent is
                        # being used instead of decentralized ones). We include all the weather variables used as
                        # states, and use the list appended_states to make sure we don't include any repeated states
                        # (i.e. weather variables measured by different buildings)
                        if state_name in ['t_in', 'avg_unmet_setpoint', 'rh_in', 'non_shiftable_load', 'solar_gen']:
                            s_low_central_agent.append(min(building.sim_results[state_name]))
                            s_high_central_agent.append(max(building.sim_results[state_name]))

                        elif state_name not in appended_states:
                            s_low_central_agent.append(min(building.sim_results[state_name]))
                            s_high_central_agent.append(max(building.sim_results[state_name]))
                            appended_states.append(state_name)
                    else:
                        s_low.append(0.0)
                        s_high.append(1.0)
                        s_low_central_agent.append(0.0)
                        s_high_central_agent.append(1.0)

            '''The energy storage (tank) capacity indicates how many times bigger the tank is compared to the maximum hourly energy demand of the building (cooling or DHW respectively), which sets a lower bound for the action of 1/tank_capacity, as the energy storage device can't provide the building with more energy than it will ever need for a given hour. The heat pump is sized using approximately the maximum hourly energy demand of the building (after accounting for the COP, see function autosize). Therefore, we make the fair assumption that the action also has an upper bound equal to 1/tank_capacity. This boundaries should speed up the learning process of the agents and make them more stable rather than if we just set them to -1 and 1. I.e. if Chilled_Water_Tank.Capacity is 3 (3 times the max. hourly demand of the building in the entire year), its actions will be bounded between -1/3 and 1/3'''
            a_low, a_high = [], []
            for action_name, value in zip(buildings_states_actions[uid]['actions'],
                                          buildings_states_actions[uid]['actions'].values()):
                if value == True:
                    if action_name == 'cooling_storage':

                        # Avoid division by 0
                        if attributes['Chilled_Water_Tank']['capacity'] > 0.000001:
                            a_low.append(max(-1.0 / attributes['Chilled_Water_Tank']['capacity'], -1.0))
                            a_high.append(min(1.0 / attributes['Chilled_Water_Tank']['capacity'], 1.0))
                            a_low_central_agent.append(max(-1.0 / attributes['Chilled_Water_Tank']['capacity'], -1.0))
                            a_high_central_agent.append(min(1.0 / attributes['Chilled_Water_Tank']['capacity'], 1.0))
                        else:
                            a_low.append(-1.0)
                            a_high.append(1.0)
                            a_low_central_agent.append(-1.0)
                            a_high_central_agent.append(1.0)
                    elif action_name == 'dhw_storage':
                        if attributes['DHW_Tank']['capacity'] > 0.000001:
                            a_low.append(max(-1.0 / attributes['DHW_Tank']['capacity'], -1.0))
                            a_high.append(min(1.0 / attributes['DHW_Tank']['capacity'], 1.0))
                            a_low_central_agent.append(max(-1.0 / attributes['DHW_Tank']['capacity'], -1.0))
                            a_high_central_agent.append(min(1.0 / attributes['DHW_Tank']['capacity'], 1.0))
                        else:
                            a_low.append(-1.0)
                            a_high.append(1.0)
                            a_low_central_agent.append(-1.0)
                            a_high_central_agent.append(1.0)

                    elif action_name == 'electrical_storage':
                        a_low.append(-1.0)
                        a_high.append(1.0)
                        a_low_central_agent.append(-1.0)
                        a_high_central_agent.append(1.0)

            building.set_state_space(np.array(s_high), np.array(s_low))
            building.set_action_space(np.array(a_high), np.array(a_low))

            observation_spaces.append(building.observation_space)
            action_spaces.append(building.action_space)

            buildings[uid] = building

    observation_space_central_agent = spaces.Box(low=np.float32(np.array(s_low_central_agent)),
                                                 high=np.float32(np.array(s_high_central_agent)), dtype=np.float32)
    action_space_central_agent = spaces.Box(low=np.float32(np.array(a_low_central_agent)),
                                            high=np.float32(np.array(a_high_central_agent)), dtype=np.float32)

    for building in buildings.values():

        # If the DHW device is a HeatPump
        if isinstance(building.dhw_heating_device, HeatPump):
            # Calculating COPs of the heat pumps for every hour
            building.dhw_heating_device.cop_heating = building.dhw_heating_device.eta_tech * (
                    building.dhw_heating_device.t_target_heating + 273.15) / (
                                                              building.dhw_heating_device.t_target_heating -
                                                              weather_data['Outdoor Drybulb Temperature [C]'])
            building.dhw_heating_device.cop_heating[building.dhw_heating_device.cop_heating < 0] = 20.0
            building.dhw_heating_device.cop_heating[building.dhw_heating_device.cop_heating > 20] = 20.0
            building.dhw_heating_device.cop_heating = building.dhw_heating_device.cop_heating.to_numpy()

        building.cooling_device.cop_cooling = building.cooling_device.eta_tech * (
                building.cooling_device.t_target_cooling + 273.15) / (weather_data[
                                                                          'Outdoor Drybulb Temperature [C]'] - building.cooling_device.t_target_cooling)
        building.cooling_device.cop_cooling[building.cooling_device.cop_cooling < 0] = 20.0
        building.cooling_device.cop_cooling[building.cooling_device.cop_cooling > 20] = 20.0
        building.cooling_device.cop_cooling = building.cooling_device.cop_cooling.to_numpy()

        building.reset()

    auto_size(buildings)

    return buildings, observation_spaces, action_spaces, observation_space_central_agent, action_space_central_agent


class CityLearn(gym.Env):
    def __init__(self, data_path, building_attributes, weather_file, solar_profile, building_ids, carbon_intensity=None, electricity_price=None,
                 buildings_states_actions=None, simulation_period=(0, 8759),
                 cost_function=['ramping', '1-load_factor', 'average_daily_peak', 'peak_demand',
                                'net_electricity_consumption', 'electricity_costs'], central_agent=False, save_memory=True, verbose=0):
        with open(buildings_states_actions) as json_file:
            self.buildings_states_actions = json.load(json_file)

        self.data_path = data_path
        self.buildings_states_actions_filename = buildings_states_actions
        self.buildings_net_electricity_demand = []
        self.building_attributes = building_attributes
        self.solar_profile = solar_profile
        self.carbon_intensity = carbon_intensity
        self.electricity_price = electricity_price
        self.building_ids = building_ids
        self.cost_function = cost_function
        self.cost_rbc = None
        self.weather_file = weather_file
        self.central_agent = central_agent
        self.loss = []
        self.verbose = verbose

        params_loader = {'data_path': data_path,
                         'building_attributes': self.data_path / self.building_attributes,
                         'weather_file': self.data_path / self.weather_file,
                         'solar_profile': self.data_path / self.solar_profile,
                         'carbon_intensity': self.data_path / self.carbon_intensity,
                         'electricity_price': self.data_path / self.electricity_price,
                         'building_ids': building_ids,
                         'buildings_states_actions': self.buildings_states_actions,
                         'save_memory': save_memory}

        self.buildings, self.observation_spaces, self.action_spaces, self.observation_space, self.action_space = building_loader(
            **params_loader)

        self.simulation_period = simulation_period
        self.uid = None
        self.n_buildings = len([i for i in self.buildings])
        self.reset()

    def get_state_action_spaces(self):
        return self.observation_spaces, self.action_spaces

    def next_hour(self):
        self.time_step = next(self.hour)
        for building in self.buildings.values():
            building.time_step = self.time_step

    def get_building_information(self):

        np.seterr(divide='ignore', invalid='ignore')
        # Annual DHW demand, Annual Cooling Demand, Annual Electricity Demand
        building_info = {}
        n_years = (self.simulation_period[1] - self.simulation_period[0] + 1) / 8760
        for uid, building in self.buildings.items():
            building_info[uid] = {}
            building_info[uid]['building_type'] = building.building_type
            building_info[uid]['climate_zone'] = building.climate_zone
            building_info[uid]['solar_power_capacity (kW)'] = round(building.solar_power_capacity, 3)
            building_info[uid]['Annual_DHW_demand (kWh)'] = round(sum(building.sim_results['dhw_demand']) / n_years, 3)
            building_info[uid]['Annual_cooling_demand (kWh)'] = round(
                sum(building.sim_results['cooling_demand']) / n_years, 3)
            building_info[uid]['Annual_nonshiftable_electrical_demand (kWh)'] = round(
                sum(building.sim_results['non_shiftable_load']) / n_years, 3)

            building_info[uid]['Correlations_DHW'] = {}
            building_info[uid]['Correlations_cooling_demand'] = {}
            building_info[uid]['Correlations_non_shiftable_load'] = {}

            for uid_corr, building_corr in self.buildings.items():
                if uid_corr != uid:
                    building_info[uid]['Correlations_DHW'][uid_corr] = round((np.corrcoef(
                        np.array(building.sim_results['dhw_demand']),
                        np.array(building_corr.sim_results['dhw_demand'])))[0][1], 3)
                    building_info[uid]['Correlations_cooling_demand'][uid_corr] = round((np.corrcoef(
                        np.array(building.sim_results['cooling_demand']),
                        np.array(building_corr.sim_results['cooling_demand'])))[0][1], 3)
                    building_info[uid]['Correlations_non_shiftable_load'][uid_corr] = round((np.corrcoef(
                        np.array(building.sim_results['non_shiftable_load']),
                        np.array(building_corr.sim_results['non_shiftable_load'])))[0][1], 3)

        return building_info

    def step(self, actions):

        self.buildings_net_electricity_demand = []
        self.current_carbon_intensity = list(self.buildings.values())[0].sim_results['carbon_intensity'][self.time_step]
        self.current_electricity_price = list(self.buildings.values())[0].sim_results['electricity_price'][self.time_step]
        electric_demand = 0
        elec_consumption_electrical_storage = 0
        elec_consumption_dhw_storage = 0
        elec_consumption_cooling_storage = 0
        elec_consumption_dhw_total = 0
        elec_consumption_cooling_total = 0
        elec_consumption_appliances = 0
        elec_generation = 0

        if self.central_agent:
            # If the agent is centralized, all the actions for all the buildings are provided as an ordered list of
            # numbers. The order corresponds to the order of the buildings as they appear on the file
            # building_attributes.json, and only considering the buildings selected for the simulation by the user (
            # building_ids).
            for uid, building in self.buildings.items():

                if self.buildings_states_actions[uid]['actions']['cooling_storage']:
                    # Cooling
                    _electric_demand_cooling = building.set_storage_cooling(actions[0])
                    actions = actions[1:]
                    elec_consumption_cooling_storage += building._electric_consumption_cooling_storage
                else:
                    _electric_demand_cooling = 0

                if self.buildings_states_actions[uid]['actions']['dhw_storage']:
                    # DHW
                    _electric_demand_dhw = building.set_storage_heating(actions[0])
                    actions = actions[1:]
                    elec_consumption_dhw_storage += building._electric_consumption_dhw_storage
                else:
                    _electric_demand_dhw = 0

                # Total heating and cooling electrical loads
                elec_consumption_cooling_total += _electric_demand_cooling
                elec_consumption_dhw_total += _electric_demand_dhw

                # Electrical appliances
                _non_shiftable_load = building.get_non_shiftable_load()
                elec_consumption_appliances += _non_shiftable_load

                # Solar generation
                _solar_generation = building.get_solar_power()
                elec_generation += _solar_generation

                # Adding loads from appliances and subtracting solar generation to the net electrical load of each
                # building
                building_electric_demand = round(
                    _electric_demand_cooling + _electric_demand_dhw + _non_shiftable_load - _solar_generation, 4)

                # Electricity consumed by every building
                building.current_net_electricity_demand = building_electric_demand
                self.buildings_net_electricity_demand.append(
                    -building_electric_demand)  # >0 if solar generation > electricity consumption

                # Total electricity consumption
                electric_demand += building_electric_demand

            assert len(actions) == 0, 'Some of the actions provided were not used'

        else:

            assert len(
                actions) == self.n_buildings, "The length of the list of actions should match the length of the list " \
                                              "of buildings. "

            for a, (uid, building) in zip(actions, self.buildings.items()):

                assert sum(self.buildings_states_actions[uid]['actions'].values()) == len(
                    a), "The number of input actions for building " + str(
                    uid) + " must match the number of actions defined in the list of building attributes."

                # Getting input actions and establishing associations between the components of the action array and
                # their corresponding actions.
                if self.buildings_states_actions[uid]['actions']['electrical_storage']:

                    if self.buildings_states_actions[uid]['actions']['cooling_storage']:
                        # Cooling
                        _electric_demand_cooling = building.set_storage_cooling(a[0])
                        elec_consumption_cooling_storage += building._electric_consumption_cooling_storage

                        # 'Electrical Storage' & 'Cooling Storage' & 'DHW Storage'
                        if self.buildings_states_actions[uid]['actions']['dhw_storage']:
                            # DHW
                            _electric_demand_dhw = building.set_storage_heating(a[1])
                            elec_consumption_dhw_storage += building._electric_consumption_dhw_storage

                            # Electrical
                            _electric_demand_electrical_storage = building.set_storage_electrical(a[2])
                            elec_consumption_electrical_storage += _electric_demand_electrical_storage

                        # 'Electrical Storage' & 'Cooling Storage'
                        else:
                            _electric_demand_dhw = building.set_storage_heating(0.0)
                            # Electrical
                            _electric_demand_electrical_storage = building.set_storage_electrical(a[1])
                            elec_consumption_electrical_storage += _electric_demand_electrical_storage
                    else:
                        _electric_demand_cooling = building.set_storage_cooling(0.0)
                        # 'Electrical Storage' & 'DHW Storage'
                        if self.buildings_states_actions[uid]['actions']['dhw_storage']:
                            # DHW
                            _electric_demand_dhw = building.set_storage_heating(a[0])
                            elec_consumption_dhw_storage += building._electric_consumption_dhw_storage

                            # Electrical
                            _electric_demand_electrical_storage = building.set_storage_electrical(a[1])
                            elec_consumption_electrical_storage += _electric_demand_electrical_storage

                        # 'Electrical Storage'
                        else:
                            _electric_demand_dhw = building.set_storage_heating(0.0)
                            # Electrical
                            _electric_demand_electrical_storage = building.set_storage_electrical(a[0])
                            elec_consumption_electrical_storage += _electric_demand_electrical_storage

                else:

                    _electric_demand_electrical_storage = 0.0

                    if self.buildings_states_actions[uid]['actions']['cooling_storage']:
                        # Cooling
                        _electric_demand_cooling = building.set_storage_cooling(a[0])
                        elec_consumption_cooling_storage += building._electric_consumption_cooling_storage

                        if self.buildings_states_actions[uid]['actions']['dhw_storage']:
                            # DHW
                            _electric_demand_dhw = building.set_storage_heating(a[1])
                            elec_consumption_dhw_storage += building._electric_consumption_dhw_storage

                        else:
                            _electric_demand_dhw = building.set_storage_heating(0.0)

                    else:
                        _electric_demand_cooling = building.set_storage_cooling(0.0)
                        # DHW
                        _electric_demand_dhw = building.set_storage_heating(a[0])
                        elec_consumption_dhw_storage += building._electric_consumption_dhw_storage

                # Total heating and cooling electrical loads
                elec_consumption_cooling_total += _electric_demand_cooling
                elec_consumption_dhw_total += _electric_demand_dhw

                # Electrical appliances
                _non_shiftable_load = building.get_non_shiftable_load()
                elec_consumption_appliances += _non_shiftable_load

                # Solar generation
                _solar_generation = building.get_solar_power()
                elec_generation += _solar_generation

                # Adding loads from appliances and subtracting solar generation to the net electrical load of each
                # building
                building_electric_demand = round(_electric_demand_electrical_storage + _electric_demand_cooling +
                                                 _electric_demand_dhw + _non_shiftable_load - _solar_generation, 4)

                # Electricity consumed by every building
                building.current_net_electricity_demand = building_electric_demand
                self.buildings_net_electricity_demand.append(-building_electric_demand)

                # Total electricity consumption
                electric_demand += building_electric_demand

        self.next_hour()

        if self.central_agent:
            s, s_appended = [], []
            for uid, building in self.buildings.items():

                # If the agent is centralized, we append the states avoiding repetition. I.e. if multiple buildings
                # share the outdoor temperature as a state, we only append it once to the states of the central
                # agent. The variable s_appended is used for this purpose.
                for state_name, value in self.buildings_states_actions[uid]['states'].items():
                    if value == True:
                        if state_name not in s_appended:
                            if state_name in ['t_in', 'avg_unmet_setpoint', 'rh_in', 'non_shiftable_load', 'solar_gen']:
                                s.append(building.sim_results[state_name][self.time_step])
                            elif state_name == 'net_electricity_consumption':
                                s.append(building.current_net_electricity_demand)
                            elif state_name != 'cooling_storage_soc' and state_name != 'dhw_storage_soc':
                                s.append(building.sim_results[state_name][self.time_step])
                                s_appended.append(state_name)
                            elif state_name == 'cooling_storage_soc':
                                s.append(building.cooling_storage._soc / building.cooling_storage.capacity)
                            elif state_name == 'dhw_storage_soc':
                                s.append(building.dhw_storage._soc / building.dhw_storage.capacity)
            self.state = np.array(s)
            rewards = reward_function_sa(self.buildings_net_electricity_demand)
            self.cumulated_reward_episode += rewards

        else:
            # If the controllers are decentralized, we append all the states to each associated agent's list of states.
            self.state = []
            for uid, building in self.buildings.items():
                s = []
                for state_name, value in self.buildings_states_actions[uid]['states'].items():
                    if value == True:
                        if state_name == 'net_electricity_consumption':
                            s.append(building.current_net_electricity_demand)
                        elif (state_name != 'cooling_storage_soc') and (state_name != 'dhw_storage_soc') and (
                                state_name != 'electrical_storage_soc'):
                            s.append(building.sim_results[state_name][self.time_step])
                        elif state_name == 'cooling_storage_soc':
                            s.append(building.cooling_storage._soc / building.cooling_storage.capacity)
                        elif state_name == 'dhw_storage_soc':
                            s.append(building.dhw_storage._soc / building.dhw_storage.capacity)
                        elif state_name == 'electrical_storage_soc':
                            s.append(building.electrical_storage._soc / building.electrical_storage.capacity)

                self.state.append(np.array(s))
            self.state = np.array(self.state, dtype='object')

            rewards = self.reward_function.get_rewards(self.buildings_net_electricity_demand,
                                                       self.current_carbon_intensity, self.current_electricity_price)
            self.cumulated_reward_episode += sum(rewards)

        # Control variables which are used to display the results and the behavior of the buildings at district level
        self.carbon_emissions.append(np.float32(max(0, electric_demand) * self.current_carbon_intensity))
        self.electricity_costs.append(np.float32(max(0, electric_demand) * self.current_electricity_price))
        self.net_electric_consumption.append(np.float32(electric_demand))
        self.electric_consumption_electric_storage.append(np.float32(elec_consumption_electrical_storage))
        self.electric_consumption_dhw_storage.append(np.float32(elec_consumption_dhw_storage))
        self.electric_consumption_cooling_storage.append(np.float32(elec_consumption_cooling_storage))
        self.electric_consumption_dhw.append(np.float32(elec_consumption_dhw_total))
        self.electric_consumption_cooling.append(np.float32(elec_consumption_cooling_total))
        self.electric_consumption_appliances.append(np.float32(elec_consumption_appliances))
        self.electric_generation.append(np.float32(elec_generation))
        self.net_electric_consumption_no_storage.append(np.float32(
            electric_demand - elec_consumption_cooling_storage - elec_consumption_dhw_storage - elec_consumption_electrical_storage))
        self.net_electric_consumption_no_pv_no_storage.append(np.float32(
            electric_demand + elec_generation - elec_consumption_cooling_storage - elec_consumption_dhw_storage - elec_consumption_electrical_storage))

        terminal = self._terminal()
        return self._get_ob(), rewards, terminal, {}

    def reset_baseline_cost(self):
        self.cost_rbc = None

    def reset(self):

        # Initialization of variables
        self.hour = iter(np.array(range(self.simulation_period[0], self.simulation_period[1] + 1)))
        self.next_hour()

        self.carbon_emissions = []
        self.electricity_costs= []
        self.net_electric_consumption = []
        self.net_electric_consumption_no_storage = []
        self.net_electric_consumption_no_pv_no_storage = []
        self.electric_consumption_electric_storage = []
        self.electric_consumption_dhw_storage = []
        self.electric_consumption_cooling_storage = []
        self.electric_consumption_electrical_storage = []
        self.electric_consumption_dhw = []
        self.electric_consumption_cooling = []
        self.electric_consumption_appliances = []
        self.electric_generation = []

        self.cumulated_reward_episode = 0
        self.current_carbon_intensity = 0

        if self.central_agent:
            s, s_appended = [], []
            for uid, building in self.buildings.items():
                building.reset()
                for state_name, value in self.buildings_states_actions[uid]['states'].items():
                    if state_name not in s_appended:
                        if value == True:
                            if state_name in ['t_in', 'avg_unmet_setpoint', 'rh_in', 'non_shiftable_load', 'solar_gen']:
                                s.append(building.sim_results[state_name][self.time_step])
                            elif state_name == 'net_electricity_consumption':
                                s.append(building.current_net_electricity_demand)
                            elif state_name != 'cooling_storage_soc' and state_name != 'dhw_storage_soc':
                                s.append(building.sim_results[state_name][self.time_step])
                                s_appended.append(state_name)
                            elif state_name == 'cooling_storage_soc':
                                s.append(0.0)
                            elif state_name == 'dhw_storage_soc':
                                s.append(0.0)
            self.state = np.array(s)
        else:
            self.reward_function = reward_function_ma(len(self.building_ids), self.get_building_information())

            self.state = []
            for uid, building in self.buildings.items():
                building.reset()
                s = []
                for state_name, value in zip(self.buildings_states_actions[uid]['states'],
                                             self.buildings_states_actions[uid]['states'].values()):
                    if value == True:
                        if state_name == 'net_electricity_consumption':
                            s.append(building.current_net_electricity_demand)
                        elif (state_name != 'cooling_storage_soc') and (state_name != 'dhw_storage_soc') and (
                                state_name != 'electrical_storage_soc'):
                            s.append(building.sim_results[state_name][self.time_step])
                        elif state_name == 'cooling_storage_soc':
                            s.append(0.0)
                        elif state_name == 'dhw_storage_soc':
                            s.append(0.0)
                        elif state_name == 'electrical_storage_soc':
                            s.append(0.0)

                self.state.append(np.array(s, dtype=np.float32))

            self.state = np.array(self.state, dtype='object')

        return self._get_ob()

    def _get_ob(self):
        return self.state

    def _terminal(self):
        is_terminal = bool(self.time_step >= self.simulation_period[1])
        if is_terminal:
            for building in self.buildings.values():
                building.terminate()

            # When the simulation is over, convert all the control variables to numpy arrays so they are easier to plot.
            self.carbon_emissions = np.array(self.carbon_emissions)
            self.electricity_costs = np.array(self.electricity_costs)
            self.net_electric_consumption = np.array(self.net_electric_consumption)
            self.net_electric_consumption_no_storage = np.array(self.net_electric_consumption_no_storage)
            self.net_electric_consumption_no_pv_no_storage = np.array(self.net_electric_consumption_no_pv_no_storage)
            self.electric_consumption_electric_storage = np.array(self.electric_consumption_electric_storage)
            self.electric_consumption_dhw_storage = np.array(self.electric_consumption_dhw_storage)
            self.electric_consumption_cooling_storage = np.array(self.electric_consumption_cooling_storage)
            self.electric_consumption_dhw = np.array(self.electric_consumption_dhw)
            self.electric_consumption_cooling = np.array(self.electric_consumption_cooling)
            self.electric_consumption_appliances = np.array(self.electric_consumption_appliances)
            self.electric_generation = np.array(self.electric_generation)
            #             self.loss.append([i for i in self.get_baseline_cost().values()])

            if self.verbose == 1:
                print('Cumulated reward: ' + str(self.cumulated_reward_episode))

        return is_terminal

    def seed(self, seed=None):
        self.np_random, seed = seeding.np_random(seed)
        return [seed]

    def get_buildings_net_electric_demand(self):
        return self.buildings_net_electricity_demand

    def cost(self):

        # Running the reference rule-based controller to find the baseline cost
        if self.cost_rbc is None:
            env_rbc = CityLearn(self.data_path, self.building_attributes, self.weather_file, self.solar_profile,
                                self.building_ids, carbon_intensity=self.carbon_intensity, electricity_price=self.electricity_price,
                                buildings_states_actions=self.buildings_states_actions_filename,
                                simulation_period=self.simulation_period, cost_function=self.cost_function,
                                central_agent=False)
            _, actions_spaces = env_rbc.get_state_action_spaces()

            # Instantiatiing the control agent(s)
            agent_rbc = RBC_Agent(actions_spaces)

            state = env_rbc.reset()
            done = False
            while not done:
                action = agent_rbc.select_action(
                    [list(env_rbc.buildings.values())[0].sim_results['hour'][env_rbc.time_step]])
                next_state, rewards, done, _ = env_rbc.step(action)
                state = next_state

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                self.cost_rbc, self.cost_rbc_last_yr = env_rbc.get_baseline_cost()
            else:
                self.cost_rbc = env_rbc.get_baseline_cost()

        # Compute the costs normalized by the baseline costs
        cost, cost_last_yr, c_score, c_score_last_yr = {}, {}, [], []
        if 'ramping' in self.cost_function:
            cost['ramping'] = np.abs(
                (self.net_electric_consumption - np.roll(self.net_electric_consumption, 1))[1:]).sum() / self.cost_rbc[
                                  'ramping']
            c_score.append(cost['ramping'])

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['ramping_last_yr'] = np.abs(
                    (self.net_electric_consumption[-8760:] - np.roll(self.net_electric_consumption[-8760:], 1))[
                    1:]).sum() / self.cost_rbc_last_yr['ramping_last_yr']
                c_score_last_yr.append(cost_last_yr['ramping_last_yr'])

        # Finds the load factor for every month (average monthly demand divided by its maximum peak), and averages
        # all the load factors across the 12 months. The metric is one minus the load factor.
        if '1-load_factor' in self.cost_function:
            cost['1-load_factor'] = np.mean([1 - np.mean(self.net_electric_consumption[i:i + int(8760 / 12)]) / np.max(
                self.net_electric_consumption[i:i + int(8760 / 12)]) for i in
                                             range(0, len(self.net_electric_consumption), int(8760 / 12))]) / \
                                    self.cost_rbc['1-load_factor']
            c_score.append(cost['1-load_factor'])

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['1-load_factor_last_yr'] = np.mean([1 - np.mean(
                    self.net_electric_consumption[-8760:][i:i + int(8760 / 12)]) / np.max(
                    self.net_electric_consumption[-8760:][i:i + int(8760 / 12)]) for i in
                                                                 range(0, len(self.net_electric_consumption[-8760:]),
                                                                       int(8760 / 12))]) / self.cost_rbc_last_yr[
                                                            '1-load_factor_last_yr']
                c_score_last_yr.append(cost_last_yr['1-load_factor_last_yr'])

        # Average of all the daily peaks of the 365 day of the year. The peaks are calculated using the net energy
        # demand of the whole district of buildings.
        if 'average_daily_peak' in self.cost_function:
            cost['average_daily_peak'] = np.mean([np.array(self.net_electric_consumption[i:i + 24]).max() for i in
                                                  range(0, len(self.net_electric_consumption), 24)]) / self.cost_rbc[
                                             'average_daily_peak']
            c_score.append(cost['average_daily_peak'])

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['average_daily_peak_last_yr'] = np.mean(
                    [self.net_electric_consumption[-8760:][i:i + 24].max() for i in
                     range(0, len(self.net_electric_consumption[-8760:]), 24)]) / self.cost_rbc_last_yr[
                                                                 'average_daily_peak_last_yr']
                c_score_last_yr.append(cost_last_yr['average_daily_peak_last_yr'])

        # Peak demand of the district for the whole year period.
        if 'peak_demand' in self.cost_function:
            cost['peak_demand'] = np.array(self.net_electric_consumption).max() / self.cost_rbc['peak_demand']
            c_score.append(cost['peak_demand'])

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['peak_demand_last_yr'] = np.array(self.net_electric_consumption[-8760:]).max() / \
                                                    self.cost_rbc_last_yr['peak_demand_last_yr']
                c_score_last_yr.append(cost_last_yr['peak_demand_last_yr'])

        # Positive net electricity consumption for the whole district. It is clipped at a min. value of 0 because the
        # objective is to minimize the energy consumed in the district, not to profit from the excess generation. (
        # Island operation is therefore incentivized)
        if 'net_electricity_consumption' in self.cost_function:
            cost['net_electricity_consumption'] = np.array(self.net_electric_consumption).clip(min=0).sum() / \
                                                self.cost_rbc['net_electricity_consumption']

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['net_electricity_consumption_last_yr'] = np.array(
                    self.net_electric_consumption[-8760:]).clip(
                    min=0).sum() / self.cost_rbc_last_yr['net_electricity_consumption_last_yr']

        if 'carbon_emissions' in self.cost_function:
            cost['carbon_emissions'] = np.array(self.carbon_emissions).sum() / self.cost_rbc['carbon_emissions']

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['carbon_emissions_last_yr'] = np.array(self.carbon_emissions[-8760:]).sum() / \
                                                            self.cost_rbc_last_yr['carbon_emissions_last_yr']
        if 'electricity_costs' in self.cost_function:
            cost['electricity_costs'] = np.array(self.electricity_costs).sum() / self.cost_rbc['electricity_costs']

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['electricity_costs_last_yr'] = np.array(self.electricity_costs[-8760:]).sum() / \
                                                            self.cost_rbc_last_yr['electricity_costs_last_yr']

        # Not used for the challenge
        if 'quadratic' in self.cost_function:
            cost['quadratic'] = (np.array(self.net_electric_consumption).clip(min=0) ** 2).sum() / self.cost_rbc[
                'quadratic']
            c_score.append(cost['quadratic'])

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['quadratic_last_yr'] = (self.net_electric_consumption[-8760:].clip(min=0) ** 2).sum() / \
                                                    self.cost_rbc_last_yr['quadratic_last_yr']
                c_score_last_yr.append(cost_last_yr['quadratic_last_yr'])

        cost['total'] = np.mean([c for c in cost.values()])

        if c_score != []:
            cost['coordination_score'] = np.mean(c_score)

        if c_score_last_yr != []:
            cost_last_yr['coordination_score_last_yr'] = np.mean(c_score_last_yr)

        if self.simulation_period[1] - self.simulation_period[0] > 8760:
            cost_last_yr['total_last_yr'] = np.mean([c for c in cost_last_yr.values()])
            return cost, cost_last_yr

        return cost

    def get_baseline_cost(self):

        # Computes the costs for the Rule-based controller, which are used to normalized the actual costs.
        cost, cost_last_yr = {}, {}
        if 'ramping' in self.cost_function:
            cost['ramping'] = np.abs(
                (self.net_electric_consumption - np.roll(self.net_electric_consumption, 1))[1:]).sum()

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['ramping_last_yr'] = np.abs(
                    (self.net_electric_consumption[-8760:] - np.roll(self.net_electric_consumption[-8760:], 1))[
                    1:]).sum()

        if '1-load_factor' in self.cost_function:
            cost['1-load_factor'] = np.mean([1 - np.mean(self.net_electric_consumption[i:i + int(8760 / 12)]) / np.max(
                self.net_electric_consumption[i:i + int(8760 / 12)]) for i in
                                             range(0, len(self.net_electric_consumption), int(8760 / 12))])

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['1-load_factor_last_yr'] = np.mean([1 - np.mean(
                    self.net_electric_consumption[-8760:][i:i + int(8760 / 12)]) / np.max(
                    self.net_electric_consumption[-8760:][i:i + int(8760 / 12)]) for i in
                                                                 range(0, len(self.net_electric_consumption[-8760:]),
                                                                       int(8760 / 12))])

        if 'average_daily_peak' in self.cost_function:
            cost['average_daily_peak'] = np.mean([self.net_electric_consumption[i:i + 24].max() for i in
                                                  range(0, len(self.net_electric_consumption), 24)])

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['average_daily_peak_last_yr'] = np.mean(
                    [self.net_electric_consumption[-8760:][i:i + 24].max() for i in
                     range(0, len(self.net_electric_consumption[-8760:]), 24)])

        if 'peak_demand' in self.cost_function:
            cost['peak_demand'] = self.net_electric_consumption.max()

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['peak_demand_last_yr'] = self.net_electric_consumption[-8760:].max()

        if 'net_electricity_consumption' in self.cost_function:
            cost['net_electricity_consumption'] = self.net_electric_consumption.clip(min=0).sum()

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['net_electricity_consumption_last_yr'] = self.net_electric_consumption[-8760:].clip(
                    min=0).sum()

        if 'carbon_emissions' in self.cost_function:
            cost['carbon_emissions'] = self.carbon_emissions.sum()

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['carbon_emissions_last_yr'] = self.carbon_emissions[-8760:].sum()

        if 'electricity_costs' in self.cost_function:
            cost['electricity_costs'] = self.electricity_costs.sum()

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['electricity_costs_last_yr'] = self.electricity_costs[-8760:].sum()

        if 'quadratic' in self.cost_function:
            cost['quadratic'] = (self.net_electric_consumption.clip(min=0) ** 2).sum()

            if self.simulation_period[1] - self.simulation_period[0] > 8760:
                cost_last_yr['quadratic_last_yr'] = (self.net_electric_consumption[-8760:].clip(min=0) ** 2).sum()

        if self.simulation_period[1] - self.simulation_period[0] > 8760:
            return cost, cost_last_yr

        return cost