Merge pull request #326 from davidusb-geek/davidusb-geek/dev/thermal_…

…model

Davidusb geek/dev/thermal model

CHANGELOG.md

@@ -1,5 +1,10 @@
 # Changelog
 
+## 0.10.4 - 2024-07-10
+### Improvement
+- Added a new thermal modeling, see the new section in the documentation for help to implement this of model for thermal deferrable loads
+- Improved documentation
+
 ## 0.10.3 - 2024-07-06
 ### Improvement
 - Added improved support for `def_start_penalty` option

docs/conf.py

@@ -22,7 +22,7 @@
 author = 'David HERNANDEZ'
 
 # The full version, including alpha/beta/rc tags
-release = '0.10.3'
+release = '0.10.4'
 
 # -- General configuration ---------------------------------------------------
 

docs/index.md

@@ -23,6 +23,7 @@ lpems.md
 forecasts.md
 mlforecaster.md
 mlregressor.md
+thermal_model.md
 study_case.md
 config.md
 emhass.md

docs/thermal_model.md

@@ -0,0 +1,118 @@
+# Deferrable load thermal model
+
+EMHASS supports defining a deferrable load as a thermal model.
+This is useful to control thermal equipement as: heaters, air conditioners, etc.
+The advantage of using this approach is that you will be able to define your desired room temperature jsut as you will do with your real equipmenet thermostat.
+Then EMHASS will deliver the operating schedule to maintain that desired temperature while minimizing the energy bill and taking into account the forecasted outdoor temperature.
+
+A big thanks to @werdnum for proposing this model and the initial code for implementing this.
+
+## The thermal model
+
+The thermal model implemented in EMHASS is a linear model represented by the following equation:
+
+$$
+    T_{in}^{pred}[k+1] = T_{in}^{pred}[k] + P_{def}[k]\frac{\alpha_h\Delta t}{P_{def}^{nom}}-(\gamma_c(T_{in}^{pred}[k] - T_{out}^{fcst}[k]))
+$$
+
+where $k$ is each time instant, $T_{in}^{pred}$ is the indoor predicted temperature, $T_{out}^{fcst}$ is the outdoor forecasted temperature and $P_{def}$ is the deferrable load power.
+
+In this model we can see two main configuration parameters:
+- The heating rate $\alpha_h$ in degrees per hour.
+- The cooling constant $\gamma_c$ in degrees per hour per degree of cooling.
+
+These parameters are defined according to the thermal characteristics of the building/house.
+It was reported by @werdnum, that values of $\alpha_h=5.0$ and $\gamma_c=0.1$ were reasonable in his case. 
+Of course these parameters should be adapted to each use case. This can be done with with history values of the deferrable load operation and the differents temperatures (indoor/outdoor).
+
+The following diagram tries to represent an example behavior of this model:
+
+![](./images/thermal_load_diagram.svg)
+
+## Implementing the model
+
+To implement this model we need to provide a configuration for the discussed parameters and the input temperatures. You need to pass in the start temperature, the desired room temperature per timestep, and the forecasted outdoor temperature per timestep.
+
+We will control this by using data passed at runtime.
+The first step will be to define a new entry `def_load_config`, this will be used as a dictionary to store any needed special configuration for each deferrable load.
+
+For example if we have just **two** deferrable loads and the **second** load is a **thermal load** then we will define `def_load_config` as for example:
+```
+'def_load_config': {
+    {},
+    {'thermal_config': {
+        'heating_rate': 5.0,
+        'cooling_constant': 0.1,
+        'overshoot_temperature': 24.0,
+        'start_temperature': 20,
+        'desired_temperatures': [...]
+    }}
+}
+```
+
+Here the `desired_temperatures` is a list of float values for each time step.
+
+Now we also need to define the other needed input, the `outdoor_temperature_forecast`, which is a list of float values. The list of floats for `desired_temperatures` and the list in `outdoor_temperature_forecast` should have proper lengths, if using MPC the length should be at least equal to the prediction horizon.
+
+Here is an example modified from a working example provided by @werdnum to pass all the needed data at runtime.
+This example is given for the following configuration: just one deferrable load (a thermal load), no PV, no battery, an MPC application, pre-defined heating intervals times. 
+
+```
+rest_command:
+  emhass_forecast:
+    url: http://localhost:5000/action/naive-mpc-optim
+    method: post
+    timeout: 300
+    payload: >
+      {% macro time_to_timestep(time) -%}
+        {{ (((today_at(time) - now()) / timedelta(minutes=30)) | round(0, 'ceiling')) % 48 }}
+      {%- endmacro %}
+      {%- set horizon = (state_attr('sensor.electricity_price_forecast', 'forecasts')|length) -%}
+      {%- set heated_intervals = [[time_to_timestep("06:30")|int, time_to_timestep("07:30")|int], [time_to_timestep("17:30")|int, time_to_timestep("23:00")|int]] -%}
+      {
+        "prediction_horizon": {{ horizon }},
+        "load_cost_forecast": {{
+          (
+            [states('sensor.general_price')|float(0)]
+            + state_attr('sensor.electricity_price_forecast', 'forecasts')
+            |map(attribute='per_kwh')
+            |list
+          )[:horizon]
+        }},
+        "pv_power_forecast": [
+          {% set comma = joiner(", ") -%}
+          {%- for _ in range(horizon) %}{{ comma() }}0{% endfor %}
+        ],
+        "def_load_config": {
+          "thermal_config": {
+            "heating_rate": 5.0,
+            "cooling_constant": 0.1,
+            "overshoot_temperature": 24.0,
+            "start_temperature": {{ state_attr("climate.living", "current_temperature") }},
+            "heater_desired_temperatures": [
+              {%- set comma = joiner(", ") -%}
+              {%- for i in range(horizon) -%}
+                {%- set timestep = i -%}
+                {{ comma() }}
+                {% for interval in heated_intervals if timestep >= interval[0] and timestep <= interval[1] %}
+                21
+                {%- else -%}
+                0
+                {%- endfor %}
+              {%- endfor %}
+            ]
+          }
+        },
+        "outdoor_temperature_forecast": [
+          {%- set comma = joiner(", ") -%}
+          {%- for fc in state_attr('weather.openweathermap', 'forecast') if (fc.datetime|as_datetime) > now() and (fc.datetime|as_datetime) - now() < timedelta(hours=24) -%}
+            {%- if loop.index0 * 2 < horizon -%}
+              {{ comma() }}{{ fc.temperature }}
+              {%- if loop.index0 * 2 + 1 < horizon -%}
+                {{ comma() }}{{ fc.temperature }}
+              {%- endif -%}
+            {%- endif -%}
+          {%- endfor %}
+        ]
+      }
+```

scripts/script_debug_optim.py

@@ -64,24 +64,31 @@
 
     df_input_data = fcst.get_load_cost_forecast(df_input_data)
     df_input_data = fcst.get_prod_price_forecast(df_input_data)
-    input_data_dict = {
-        'retrieve_hass_conf': retrieve_hass_conf,
-    }
+    input_data_dict = {'retrieve_hass_conf': retrieve_hass_conf}
 
     # Set special debug cases
+
+    # Solver configurations
     optim_conf.update({'lp_solver': 'PULP_CBC_CMD'}) # set the name of the linear programming solver that will be used. Options are 'PULP_CBC_CMD', 'GLPK_CMD' and 'COIN_CMD'. 
     optim_conf.update({'lp_solver_path': 'empty'})  # set the path to the LP solver, COIN_CMD default is /usr/bin/cbc
+
+    # Semi continuous and constant values
     optim_conf.update({'treat_def_as_semi_cont': [True, False]})
     optim_conf.update({'set_def_constant': [True, False]})
+
+    # A sequence of values
     # optim_conf.update({'P_deferrable_nom': [[500.0, 100.0, 100.0, 500.0], 750.0]})
 
+    # Using a battery
     optim_conf.update({'set_use_battery': False})
     optim_conf.update({'set_nocharge_from_grid': False})
     optim_conf.update({'set_battery_dynamic': True})
     optim_conf.update({'set_nodischarge_to_grid': True})
 
+    # A hybrid inverter case
     plant_conf.update({'inverter_is_hybrid': False})
 
+    # Setting some negative values on production prices
     df_input_data.loc[df_input_data.index[25:30],'unit_prod_price'] = -0.07
     df_input_data['P_PV_forecast'] = df_input_data['P_PV_forecast']*2
     P_PV_forecast = P_PV_forecast*2

scripts/script_thermal_model_optim.py

@@ -0,0 +1,137 @@
+# -*- coding: utf-8 -*-
+import pickle
+import random
+import numpy as np
+import pandas as pd
+import pathlib
+import plotly.express as px
+import plotly.subplots as sp
+import plotly.io as pio
+pio.renderers.default = 'browser'
+pd.options.plotting.backend = "plotly"
+
+from emhass.retrieve_hass import RetrieveHass
+from emhass.optimization import Optimization
+from emhass.forecast import Forecast
+from emhass.utils import get_root, get_yaml_parse, get_days_list, get_logger
+
+# the root folder
+root = str(get_root(__file__, num_parent=2))
+emhass_conf = {}
+emhass_conf['config_path'] = pathlib.Path(root) / 'config_emhass.yaml'
+emhass_conf['data_path'] = pathlib.Path(root) / 'data/'
+emhass_conf['root_path'] = pathlib.Path(root)
+
+# create logger
+logger, ch = get_logger(__name__, emhass_conf, save_to_file=False)
+
+if __name__ == '__main__':
+    get_data_from_file = True
+    params = None
+    show_figures = True
+    template = 'presentation'
+
+    retrieve_hass_conf, optim_conf, plant_conf = get_yaml_parse(emhass_conf, use_secrets=False)
+    retrieve_hass_conf, optim_conf, plant_conf = \
+        retrieve_hass_conf, optim_conf, plant_conf
+    rh = RetrieveHass(retrieve_hass_conf['hass_url'], retrieve_hass_conf['long_lived_token'], 
+                            retrieve_hass_conf['freq'], retrieve_hass_conf['time_zone'],
+                            params, emhass_conf, logger)
+    if get_data_from_file:
+        with open(emhass_conf['data_path'] / 'test_df_final.pkl', 'rb') as inp:
+            rh.df_final, days_list, var_list = pickle.load(inp)
+        retrieve_hass_conf['var_load'] = str(var_list[0])
+        retrieve_hass_conf['var_PV'] = str(var_list[1])
+        retrieve_hass_conf['var_interp'] = [retrieve_hass_conf['var_PV'], retrieve_hass_conf['var_load']]
+        retrieve_hass_conf['var_replace_zero'] = [retrieve_hass_conf['var_PV']]
+    else:
+        days_list = get_days_list(retrieve_hass_conf['days_to_retrieve'])
+        var_list = [retrieve_hass_conf['var_load'], retrieve_hass_conf['var_PV']]
+        rh.get_data(days_list, var_list,
+                        minimal_response=False, significant_changes_only=False)
+    rh.prepare_data(retrieve_hass_conf['var_load'], load_negative = retrieve_hass_conf['load_negative'],
+                            set_zero_min = retrieve_hass_conf['set_zero_min'], 
+                            var_replace_zero = retrieve_hass_conf['var_replace_zero'], 
+                            var_interp = retrieve_hass_conf['var_interp'])
+    df_input_data = rh.df_final.copy()
+
+    fcst = Forecast(retrieve_hass_conf, optim_conf, plant_conf,
+                    params, emhass_conf, logger, get_data_from_file=get_data_from_file)
+    df_weather = fcst.get_weather_forecast(method='csv')
+    P_PV_forecast = fcst.get_power_from_weather(df_weather)
+    P_load_forecast = fcst.get_load_forecast(method=optim_conf['load_forecast_method'])
+    df_input_data = pd.concat([P_PV_forecast, P_load_forecast], axis=1)
+    df_input_data.columns = ['P_PV_forecast', 'P_load_forecast']
+
+    df_input_data = fcst.get_load_cost_forecast(df_input_data)
+    df_input_data = fcst.get_prod_price_forecast(df_input_data)
+    input_data_dict = {'retrieve_hass_conf': retrieve_hass_conf}
+
+    # Set special debug cases
+
+    # Solver configurations
+    optim_conf.update({'lp_solver': 'PULP_CBC_CMD'}) # set the name of the linear programming solver that will be used. Options are 'PULP_CBC_CMD', 'GLPK_CMD' and 'COIN_CMD'. 
+    optim_conf.update({'lp_solver_path': 'empty'})  # set the path to the LP solver, COIN_CMD default is /usr/bin/cbc
+
+    # Semi continuous and constant values
+    optim_conf.update({'treat_def_as_semi_cont': [True, False]})
+    optim_conf.update({'set_def_constant': [True, False]})
+
+    # Thermal modeling
+    df_input_data['outdoor_temperature_forecast'] = [random.normalvariate(10.0, 3.0) for _ in range(48)]
+
+    runtimeparams = {
+        'def_load_config': [
+            {},
+            {'thermal_config': {
+                'heating_rate': 5.0,
+                'cooling_constant': 0.1,
+                'overshoot_temperature': 24.0,
+                'start_temperature': 20,
+                'desired_temperatures': [21]*48,
+                }
+            }
+        ]
+    }
+    if 'def_load_config' in runtimeparams:
+        optim_conf["def_load_config"] = runtimeparams['def_load_config']
+
+    costfun = 'profit'
+    opt = Optimization(retrieve_hass_conf, optim_conf, plant_conf, 
+                       fcst.var_load_cost, fcst.var_prod_price,  
+                       costfun, emhass_conf, logger)
+    # opt_res_dayahead = opt.perform_dayahead_forecast_optim(
+    #     df_input_data, P_PV_forecast, P_load_forecast)
+    unit_load_cost = df_input_data[opt.var_load_cost].values # €/kWh
+    unit_prod_price = df_input_data[opt.var_prod_price].values # €/kWh
+    opt_res_dayahead = opt.perform_optimization(df_input_data, P_PV_forecast.values.ravel(), 
+                                                P_load_forecast.values.ravel(), 
+                                                unit_load_cost, unit_prod_price, debug=True)
+
+    # Let's plot the input data
+    fig_inputs_dah = df_input_data.plot()
+    fig_inputs_dah.layout.template = template
+    fig_inputs_dah.update_yaxes(title_text = "Powers (W) and Costs(EUR)")
+    fig_inputs_dah.update_xaxes(title_text = "Time")
+    if show_figures:
+        fig_inputs_dah.show()
+
+    vars_to_plot = ['P_deferrable0', 'P_deferrable1','P_grid', 'P_PV', 'predicted_temp_heater1', 'target_temp_heater1', 'P_def_start_1', 'P_def_bin2_1']
+    if plant_conf['inverter_is_hybrid']:
+        vars_to_plot = vars_to_plot + ['P_hybrid_inverter']
+    if plant_conf['compute_curtailment']:
+        vars_to_plot = vars_to_plot + ['P_PV_curtailment']
+    if optim_conf['set_use_battery']:
+        vars_to_plot = vars_to_plot + ['P_batt'] + ['SOC_opt']
+    fig_res_dah = opt_res_dayahead[vars_to_plot].plot() # 'P_def_start_0', 'P_def_start_1', 'P_def_bin2_0', 'P_def_bin2_1'
+    fig_res_dah.layout.template = template
+    fig_res_dah.update_yaxes(title_text = "Powers (W)")
+    fig_res_dah.update_xaxes(title_text = "Time")
+    if show_figures:
+        fig_res_dah.show()
+
+    print("System with: PV, two deferrable loads, dayahead optimization, profit >> total cost function sum: "+\
+        str(opt_res_dayahead['cost_profit'].sum())+", Status: "+opt_res_dayahead['optim_status'].unique().item())
+
+    print(opt_res_dayahead[vars_to_plot])
+

setup.py

@@ -19,7 +19,7 @@
 
 setup(
     name='emhass',  # Required
-    version='0.10.3',  # Required
+    version='0.10.4',  # Required
     description='An Energy Management System for Home Assistant',  # Optional
     long_description=long_description,  # Optional
     long_description_content_type='text/markdown',  # Optional (see note above)