Introduction to optihood

optihood is based on oemof (open energy modelling framework), in particular oemof.solph and oemof.thermal packages. In oemof, an energy system is defined as a combination of linear component models of different types namely, sources, transformers, sinks, and buses. The components models are linked together using flow objects which may have associated costs (optional) to formulate the target function to optimize. Then, oemof makes use of pyomo, an open-source modelling language for optimization solutions, to define a MILP optimization problem. The framework optihood offers a flexible and easy-to-use environment with several useful features for a neighborhood beyond what oemof can provide. The basic construct of oemof is modified in a sense that the components are grouped by buildings and the results (time series of energy production/consumption, costs, emissions, etc.) can be obtained per building.

Grid electricity, natural gas, or any other form of energy consumed by the system can be considered as energy sources. An energy source is modelled simply as a “Source” from oemof solph. An energy demand can be related to electricity, space heat and domestic hot water, and is modelled as a “Sink” from oemof solph. In terms of energy conversion and storage technologies, the following are presently implemented: air-source heat pump (ASHP), ground-source heat pump (GSHP), combined heating and power (CHP), solar thermal collector, PV, electric heating rod, gas boiler, electrical battery and thermal storage.

optihood_architecture

To summarize, the technologies are classified in three categories:

  • Energy inputs: all energy vectors that are energy sources for the systems, e.g. fuels, grid electricity, energy from the environment like solar radiation, etc.

  • Energy converters: equipment that converts inputs into usable energy (heat or electricity) for the end consumers

  • Energy storages: device to store electricity or heat locally in order to be consumed later.

energy_types

Energy system component models

The energy system components can be classified into energy converters and storages. We use constant efficiency models for CHP, gas boiler and electric heating rods, where a fixed efficiency is pre-defined. These fixed efficiencies are defined by the user in the input scenario file.

Heat pumps

Heat pumps (ASHP and GSHP) are modelled based on a bi-quadratic polynomial fit of the condenser heating power (\(\dot{ q }_c\)) and the electrical consumption power of the compressor (\(\dot{w}_{cp}\)):

\[\begin{split}&\dot{q}_c = bq_1 + bq_2 \cdot \bar{T}_{e,in} + bq_3 \cdot \bar{T}_{c,out} + bq_4 \cdot \bar{T}_{e,in} \cdot{\bar{T}_c,out} + bq_5 \cdot \bar{T}^2_{e,in} + bq_6 \\ &\dot{w}_{cp} = bp_1 + bp_2 \cdot \bar{T}_{e,in} + bp_3 \cdot \bar{T}_{c,out} + bp_4 \cdot \bar{T}_{e,in} \cdot \bar{T}_{c,out} + bp_5 \cdot \bar{T}^2_{e,in} + bp_6 \cdot \bar{T}^2_{c,out}\end{split}\]

where, \(T_{e,in}\) and \(T_{c,out}\) are fluid temperatures at the inlet of the evaporator and the outlet of the condenser, respectively. \(\bar{T}\) denotes the normalized temperature and is defined as \(\bar{T} = \frac{T[^{\circ} \text{C}]}{273.15}\). For the solution of the system of equations the Brent solver is used [1]. The polynomial coefficients \(b_{qi}\) and \(b_{pi}\) are calculated from the catalog heat pump data using the multidimensional least square fitting algorithm of Scipy [2] in Python.

Table 1: Parameters, inputs and outputs of heat pump model.

HP_model_param

Polynomial fit analysis for heat pump model

R410A-predict-Cop

Figure 1: Typical coefficient of performance map (COP) for a R410A heat pump obtained using the two equations above.

R410A-Qcond

Figure 2: Differences between experimental and fitted data using the full polynomial formulation from the two equations above for condenser heat.

R410A-COP

Figure 3: Differences between experimental and fitted data using the full polynomial formulation from the two equations above for coefficient of performance (COP).

However, this model is non-linear. A way to overcome the non-linearity would be to fix the \(\bar{T}_{c,out}\) to 35 °C and 65 °C, respectively, for space heating (SH) and domestic hot water (DHW). Thus we would use for example:

\[\dot{q}_c = bq_1 + bq_2 \cdot \bar{T}_{e,in} + bq_3 \cdot \frac{35}{273.15} + bq_4 \cdot \bar{T}_{e,in} \cdot \frac{35}{273.15} + bq_5 \cdot \bar{T}_{e,in}^2 + bq_6 \cdot \frac{35}{273.15}^2\]
\[\dot{w}_{cp} = bp_1 + bp_2 \cdot \bar{T}_{e,in} + bp_3 \cdot \frac{35}{273.15} + bp_4 \cdot \bar{T}_{e,in} \cdot \frac{35}{273.15} + bp_5 \cdot \bar{T}_{e,in}^2 + bp_6 \cdot \frac{35}{273.15}^2\]

The fitted data for the HP08L-M-BC air/water heat pump using the proposed approach described by the two equations above are provided in Fig. 4-5 and Table 2, while the fitted heat pump coefficients are given in Table 3. While, the fitted data for the ProDomo13-R410A brine/water heat pump using the proposed approach described by the two equations above are provided in Fig. 6-7 and Table 4, while the fitted heat pump coefficients are given in Table 5.

Table 2: Differences between experiments and fitted data for the HP08L-M-BC air/water heat pump using the two equations above. \(error=100 \cdot |\frac{Q_{exp}-Q_{num}}{Q_{exp}}|\) and \(RMS = \sqrt { \sum{\frac{(Q_{exp}-Q_{num})^2}{n_p}} }\) where \(n_p\) is the number of data points.

HP_table2

Table 3: Fitted coefficients for the HP08L-M-BC air/water heat pump using the two equations above.

HP_table3 HP08L-M-BC-COP-1

Figure 4: Differences between experimental and fitted data of HP08L-M-BC air/water heat pump using the proposed approach from the two equations above for coefficient of performance (COP).

HP08L-M-BC-Qcond-1

Figure 5: Differences between experimental and fitted data of HP08L-M-BC air/water heat pump using the proposed approach from the two equations above for condenser heat.

ProDomo13-R410A-COP-1

Figure 6: Differences between experimental and fitted data of ProDomo13-R410A brine/water heat pump using the proposed approach from the two equations above for coefficient of performance (COP).

ProDomo13-R410A-Qcond-1

Figure 7: Differences between experimental and fitted data of ProDomo13-R410A brine/water heat pump using the proposed approach from the two equations above for condenser heat.

Table 4: Differences between experiments and fitted data for the ProDomo13-R410A brine/water heat pump using the two equations above. \(error=100 \cdot |\frac{Q_{exp}-Q_{num}}{Q_{exp}}|\) and \(RMS = \sqrt { \sum{\frac{(Q_{exp}-Q_{num})^2}{n_p}} }\) where \(n_p\) is the number of data points.

HP_table4

Table 5: Fitted coefficients for the ProDomo13-R410A brine/water heat pump using the two equations above.

HP_table5

Solar thermal collector

A module to calculate the usable heat of a flat plate collector is described in details in Solar thermal collector. The model for solar thermal collector is taken from the oemof thermal package.

PV

The installed PV provides electricity to the building during the irradiation hours. Along with the battery, the usual strategy is to store the PV surplus power in the battery to be consumed at later hours of the planning horizon. The maximum available power \(pv_t^{avail}\) of the PV is a built function that depends on the PV cell temperature, the ambient temperature and the total solar horizontal irradiation. These formulas, as well as the decision variables and the characteristics of the PV are stated in the next Table. PV modules production profiles are pre-calculated before the optimization.

Two-zone thermal energy storage

A simplified 2-zone-model of a stratified thermal energy storage is implemented and described indetails in Stratified thermal storage. The model for stratified thermal storage is taken from the oemof thermal package.

Combined production transformer

A new transformer called combined production transformer which extends the features of oemof “Transformer” was defined. Since some transformers like HP can have different efficiencies for SH and DHW production (DHW needs a higher temperature than SH), this transformer offers the possibility to consider those different efficiencies. It allows to produce both space heating (SH) and domestic hot water (DHW) during the same timestep while respecting the input/output balance constraint.

\[P_{input}(t) = \frac{P_{DHW}(t)}{\eta_{DHW}} + \frac{P_{SH}(t)}{\eta_{SH}}, \forall t\]

where, \(P\) denotes the operating power for inputs (for example, electricity used by HP) and outputs (SH and DHW), \(\eta\) denotes efficiency of the transformer and \(t\) denotes the time step. Physically the converters cannot supply both SH and DHW at the same time. However, if we consider a timestep of 1 hour it can be considered to be sub-divided into smaller intervals to produce SH and DHW both within 1 hour. The combined production transformer was used for the implementation of heat pumps (ASHP, GSHP), CHP, gas boiler and electric heating rod.

PVT collector

PVT class was implemented within the converters module, which defines the energy conversion technologies supported by optihood. The collector output is modelled based on the characteristic curve model reported in the SwissEnergy sponsored project PVT Wrap-Up (Zenhäusern et al. (2017)). The thermal output of a PVT collector, \(\dot{Q}\), highly depends on the surrounding environment and the operating conditions. The most significant influencing factors are the solar irradiation per collector surface area (\(G\)), ambient air temperature (\(T_{amb}\)) and the mean temperature of the collector fluid (\(T_m\)). The characteristic equation of thermal output of the PVT collector is given by:

\[\frac{\dot Q}{A} =(G - \frac{P_{el}^{DC}}{(\alpha \tau) \cdot A}) \cdot \eta_0 - a_1(T_m - T_{amb}) - a_2 (T_m - T_{amb})^2\]

where A stands for the gross area of the collector surface, \(P_{el}^{DC}\) stands for the DC electrical output of the collector, (alpha tau) is the transmission absorption product of the collector, \(\eta_0\) is the maximum thermal efficiency, \(a_1\) is the linear heat loss coefficient and \(a_2\) is the quadratic heat loss coefficient of the collector. A corresponding label \(PVT\) was added to the energy conversion technology processing function, to allow the definition of a PVT collector in the input excel/config file while preparing the optimization problem.

Layered thermal energy storage and discrete temperature levels

A discretized thermal energy storage with several predefined discrete temperature levels was implemented. Moreover, the heat production technologies such as heat pumps, CHP, solar thermal collectors, etc. were extended to allow multiple output flows (at different temperature levels). It should be noted that the temperature levels are predefined and each heat production technology, therefore, has a predefined hourly efficiency related to a specific temperature level. The number of discrete temperature levels is parameterized and can be defined in the input scenario excel file. In order to use discrete temperature levels, the temperatureLevels parameters has to be True when the EnergyNetwork class is instantiated:

import EnergyNetworkGroup as EnergyNetwork

#set a time period for the optimization problem
TimePeriod = np.date_range("2021-01-01 00:00:00". "2021-12-31 23:00:00", freq = "60min"

#create an energy network and set the network parameters from an excel file
network = EnergyNetwork(timeperiod, temperatureLevels=True)

The discrete temperature levels defined in the input scenario file, set the temperatures of the output flows of the heat conversion technologies. Depending on the time resolution of the optimization problem, it may not be acceptable for a heat conversion technology to produce heat at more than one temperature levels in a single time step. Therefore, limit_active_flow_count constraint of oemof solph package (Hilpert et al. (2018)) was used to permit only one of the heat output flows to remain active at a given time step. A class ThermalStorageTemperatureLevels was developed to represent a discretized thermal energy storage. The model of a layered thermal energy storage is a combination of dual temperature zone storages from oemof thermal python package (Hilpert et al. (2018)). The dual temperature zone storages include predefined calculations for top/bottom and lateral surface losses. While the lateral surface losses are preserved for the storage layers at each temperature level, the top and bottom surface losses should only be considered for the topmost (i.e. at the highest temperature level) and the lowest (i.e. at the lowest temperature) layers. The fixed one-time investment cost of the discretized thermal energy storage should be added to the objective function only once (instead of being added for each layer separately). These functionalities are implemented within the ThermalStorageTemperatureLevels class. Moreover, the total storage volume \(V_{stor}\) is calculated as the sum of individual layer volumes (\(v_i\)), as follows:

\[\sum_{i=1}^n v_i = V_{stor}\]

where \(n\) denotes the number of discrete temperature levels.

A constraint called multiTemperatureStorageCapacityConstaint was developed to implement the following rule on the storage volume capacity:

\[V_{stor,min} \leq V_{stor} \leq V_{stor,max}\]

where \(V_{stor,max}\) and \(V_{stor,min}\) represent the minimum and the maximum limits for the storage volume. The Figure below shows a graphical representation of a layered thermal energy storage with three discrete temperature levels. The DHW demand is met using the topmost temperature level at 65 °C i.e. highest temperature, while the lowest temperature level at 35 °C is used to cover the SH demand. A rule for charging the thermal energy storage was implemented, such that the energy inflow at a given storage layer (except the lowest layer), equals the energy outflow from the preceding storage layer. Therefore, in order to supply thermal energy at 50 °C to the storage, the same volume added at the 50 °C layer should be displaced from layer below, i.e. from the 35 °C storage level (as shown in Figure 11). This means that the energy conversion technologies can heat water from 35 °C to 50 °C and from 50 °C to 65 °C, in that order.

layered_storage layered_storage_1

Building model

A linear RC building model is presently under-development to replace the static space heating demand profiles. A building model is a grey-box model which is often used to depict the thermal behaviour of a building in a simplified manner. It is implemented as a custom sink component along with a set of new constraints.

building_model_oemof

The specific building model implemented in optihood was proposed and validated in [3] and is characterized by three thermal spaces:

  • wall and building mass

  • indoor air

  • distribution system

Each thermal space is at a certain temperature at a particular timestep. Moreover, each thermal space has two key parameters which represent the thermal resistance and thermal capacity. The temperature of each thermal space is influenced by the temperature of adjascent thermal spaces, heat flow, internal heat gains and ambient weather conditions.

building_model

The parameters and variables of the RC model are described below:

Parameters

\(R_{ind}\)

Thermal resistance between indoor and wall states [K/kW]

\(R_{wall}\)

Thermal resistance between wall state and outside [K/kW]

\(R_{dis}\)

Thermal resistance between indoor and distribution system states [K/kW]

\(C_{ind}\)

Thermal capacity of the indoor air state [kWh/K]

\(C_{wall}\)

Thermal capacity of the wall state [kWh/K]

\(C_{dis}\)

Thermal capacity of the distribution system state [kWh/K]

\(gA\)

Aperture area of the windows [\(m^2\)]

\(Q^{dis}_{min}\)

Minimum operating power from the tank to the distribution system [kW]

\(Q^{dis}_{max}\)

Maximum operating power from the tank to the distribution system [kW]

\(T^{ind}_{min}\)

Indoor minimum comfort temperature [°C]

\(T^{ind}_{max}\)

Indoor maximum comfort temperature [°C]

Exogenous input parameters

\(T^{amb}_{t}\)

Ambient outside air temperature at \(t^{th}\) timestep [°C]

\(I^{H}_{t}\)

Total horizontal irradiation at \(t^{th}\) timestep [kW/\(m^2\)]

\(Q^{occ}_{t}\)

Internal heat gains from occupants at \(t^{th}\) timestep [kW]

Boundary parameters

\(T^{ind}_{init}\)

Indoor initial temperature [°C]

\(T^{wall}_{init}\)

Wall initial temperature [°C]

\(T^{dis}_{init}\)

Distribution system initial temperature [°C]

State variables

\(T^{ind}_t\)

Indoor temperature at \(t^{th}\) timestep [°C]

\(T^{wall}_t\)

Wall temperature at \(t^{th}\) timestep [°C]

\(T^{dis}_t\)

Distribution system temperature at \(t^{th}\) timestep [°C]

\(\epsilon^{ind}_t\)

Violation of indoor comfort temperature range at \(t^{th}\) timestep [°C]

\(\delta^{ind}_t\)

Violation of indoor final temperature requirement [°C]

\(P^{dis}_t\)

Electric consumption of the distribution system

Decision variable

\(Q^{dis}_t\)

Heating power from the tank to the distribution system at \(t^{th}\) timestep [kW]

The state space equations of the building model are:

state_space_eq
The final constraints of the building model are:
constraint1 constraint2 constraint3 constraint4 constraint5
[1] M Galassi, J Davies, B Gough, G Jungman, M Booth, and F Rossi. GNU Scientific Library Reference Manual. John Wiley & Sons, 2nd edition, 2003.
[2] Eric Jones, Travis Oliphant, Pearu Peterson, et al. SciPy: Open source scientific tools for Python, 2001.
[3] T. Péan, R. Costa Castelló y J. Salom, Price and carbon-based energy flexibility of residential heating and cooling loads using model predictive control, Sustainable Cities and Society, vol. 50, 2019.