2.1. Mathematical Model
A grey-box energy model of a residential house was created. The model contains several energy sources and appliances—see
Figure 1. Firstly, the model contains appliances that have a fixed power consumption when turned on. Although these appliances might not necessarily have a fixed consumption on start-up, such dynamics were not included in the simulation, as the timescale of the simulations was several days. These appliances were built into the model by specifying a timeframe for when they should operate on any given day of the week. Examples of such appliances are kitchen appliances, a washing machine, a TV, a PC, and lighting fixtures.
Secondly, the model contains electric heating and cooling systems that most affect the overall power consumption. Predominantly, the outside temperature and the heating effect of solar irradiation forecasts have the biggest effect on indoor temperature. The heating system respects the demands of the user by keeping the desired temperature.
Thirdly, the model contains a photovoltaic system. Its expected energy generation forecast profile is taken from the average values of solar irradiation in the last several years [
36]. The weather forecast also contains information about expected cloudiness. This expected cloudiness was then used to bring variability to the forecasted solar irradiation. Lastly, a model of a battery of an electric vehicle was built into the model but it was only used as an appliance in the current version of the model.
The model was primarily designed for a mid-latitude residential house in Europe, where the most important appliance is heating in winter and cooling in summer. However, it can easily be parameterized for other regions.
2.1.1. Appliances with Constant Power Consumption
Such appliances in the model are kitchen appliances, a washing machine, a TV, a PC, and lighting. All of these are considered to have constant power consumption, so they are defined by their timing during the week.
2.1.2. Sunshine Irradiation Forecast
For the irradiation forecast, a system from the European Union Science Hub called PVGIS (Photovoltaic Geographical Information System) was used [
36]. It contains historical data on irradiation in different areas around the world measured for several years. The monthly irradiation data show a typical day’s irradiation of a given month for a given place. In this case, Palermo (Italy) was chosen as the city of interest. Such an irradiation profile was then recalculated based on the expected cloudiness in the upcoming days by:
where
Icloud is the irradiation with cloudiness taken into account, W/m2;
Iforecasted is the typical irradiation on a clear day, W/m2;
C is the level of cloudiness, %.
Forecasted irradiation was calculated for each timestep (5 min) using linear interpolation for both the cloudiness forecast and the typical day’s irradiation course.
The attenuation of the irradiation by clouds is considered to be 90%. That means that even when the expected cloudiness level is 100% there is still some part of the irradiation that gets to the surface of the Earth.
The main idea behind the use of (1) was to obtain an approximation of variation in the irradiation on the forecast horizon and to see how such variation would affect the model’s output.
2.1.3. Heating and Cooling Subsystem
Both heating and cooling are connected to the weather forecast. The model takes into account the desired temperature indoors and the forecasted data. The inside temperature is considered to change dynamically. Such dynamics are defined by a constant mass M, which describes heat accumulation. The energy balance is described by the discrete-time differential equation of the simplified zone thermal model [
37,
38]:
where
k is the current time step;
Pheating is the heating power, W;
PSun is the heating power of the Sun, W;
Pcooling is the cooling power, W;
K is the thermal loss coefficient, W/°C;
Toutside is the outside temperature, °C;
Tinside is the inside temperature, °C;
M is a constant that consists of mass and specific heat;
Ts is the sampling period, s.
Constants M and K were chosen to reflect experience in a real house; Moreover, the thermal loss coefficient K would change as a function of solar radiation, and convection due to wind speed and conduction.
The simplified zone thermal model provides a simple way for modeling the house (especially in the case of simulation of a large aggregate of buildings) for which more detailed and complicated models cannot be used for real-time optimization. However, more complex methods such as ASHRAE [
39] would be more suitable, for example, for comfort evaluations and for the sizing of the HVAC system in the presence of more thermal zones. There are also software programs, such as TRNSYS, that provide the possibility to have detailed simulations. However, the main idea behind this paper was to consider the parts of the model that might affect the resulting energy demands the most.
In order to achieve the desired temperature inside the home, a PI controller was implemented into the model. It was considered that a controller would be present in the building so that the effects of changing the desired temperature could be seen. The PI controller used in the model is described by:
where
u is the control output for heating or cooling, W;
k is the current timestep;
kp is the proportional gain of the controller;
Ts is the sampling period, s;
Ti is the integration time constant, s;
e is the control error, °C.
The controller’s output u is calculated in each timestep, defined by the sampling period, from the control error e, which is the difference between the desired indoor temperature and the actual indoor temperature. The control error used in the PI controller is formulated as:
where
wT is the desired indoor temperature, °C;
Tinside is the measured indoor temperature, °C.
The calculated control output is then used as the heating or cooling power needed for the current time step in order to follow the desired indoor temperature in the next time step. Both heating and cooling have a defined maximal power consumption and efficiency of transforming electrical energy into heating or cooling power, so the control output cannot be infinite but is bounded.
In the first step of the simulation, the actual indoor temperature is considered to be equal to the desired temperature unless the heating or cooling maximal power is exceeded. The amount of heating or cooling power needed to achieve the desired temperature is calculated from the static equation considering the outside temperature and the sunshine power:
In case the calculated power exceeds the limits of heating and cooling, then the actual indoor temperature is calculated from the static equation by the following formula, which is taken from Equation (5):
Then, the desired temperature in the first step is considered to be the same as the calculated indoor temperature in order to start the simulation with a control deviation equal to zero. In the following time steps, the indoor temperature is calculated by:
The actual heating and cooling power used for indoor temperature calculation is divided by the efficiency value. Each system has its own efficiency, and its range is between 0 and 1, where 0 means 0% efficiency and 1 means 100% efficiency.
The resulting forecast, shown in
Figure 2, was obtained for the building’s thermal loss coefficient of 0.25 kW/°C, M representing mass and specific heat of 3000, heating maximal power consumption of 3 kW with heating efficiency 90%, maximal power consumption of cooling of 4 kW and the efficiency of cooling 60%. These values were chosen based on experience in terms of the energy consumption values of a real house.
The sampling period and the simulation timestep were set to 300 s. The controller’s proportional gain kP was set to 4 and the parameter of the integration time constant Ti was set to 12,000 s.
The gain of the controller corresponds to the inverse of the system gain, and the integration time constant of the controller is consistent with the dominant time constant of the system. This setting leads approximately to aperiodic control response.
The heating power of sunshine is defined by the area affected by the sun during the day and the solar irradiation forecast by the formula:
where
PSun is the heating power from the Sun that is transferred indoors, W;
Aexposed is the area that is exposed to the Sun, m2;
ksolar_heat is the coefficient of transformation of the irradiation to heating power.
2.1.4. Photovoltaic Subsystem
The photovoltaic subsystem is built into the model in a similar way to the heating power of the sun. It is described by the area of the photovoltaic panels, the irradiation forecast with cloudiness levels, and the efficiency of turning the solar power into electrical power. Electrical power obtained from the photovoltaic panels is described by:
where
Ppv is the electrical power generated by the photovoltaics, W;
Apv is the area of the photovoltaic panels, m2;
kpv_efficiency is the efficiency factor of photovoltaic panels.
The first part of
Figure 3 and
Figure 4 shows the overall predicted power consumption of the home and the predicted power generation from the photovoltaic panels. By combining both terms, the overall power consumption of the home is lowered by the power generated by the photovoltaic panels, as shown by the energy balance on the second graph. The surplus of generated energy results in the possibility of selling it back to the grid. These results were obtained for two different areas covered by photovoltaic panels-namely 10 m
2 and 30 m
2.
Figure 5 and
Figure 6 show the amount of money that is forecasted to be paid in each step as well as a cumulated sum of the money that will be paid on the forecast horizon.
The efficiency of the panels was considered to be constantly 20%. The pricing was considered to be 0.185 EUR per kWh consumed from the grid and 0.018 EUR per kWh sold back to the grid. The tariffs were considered to be fix-priced; however, in future studies, it will be important to consider flexible price tariffs.
2.1.5. Electric Vehicle (EV) Battery Subsystem
The EV battery subsystem in this model is considered only as an appliance. The state of charge (SOC) of the battery in charging mode is given by:
where
SOC(k) is the state of charge in the current step, %;
SOC(k − 1) is the state of charge in the previous step, %;
Pcharging is the charging power, kW;
εch is the charging efficiency;
Cmax is the maximum possible battery capacity, kWh.
The discharging mode is described by:
where
Pdischarged is the discharged power on the output of the battery, kW;
εd is the discharging efficiency.
The energy from the battery in the model is only used for driving the electric vehicle and not as a supporting energy system for the home, which could be the next feature implemented in the model.
The first part of
Figure 7 shows the charging cycles that were fixed in the week. The second part shows the discharging, which reflects the car being driven. These timings were fixed for testing purposes. The third part shows the State of charge of the battery on the forecast horizon.
These results were obtained for the battery defined with a maximum capacity of 24 kWh, maximal charging power of 3 kW, charging efficiency of 80%, maximal discharging power of 6 kW and discharging efficiency of 60%.
2.1.6. Weather Forecast
Weather forecast uses the OpenWeatherMap API, which has different types of forecasts available via a personal key that is obtained after a free registration [
40]. The 5-day forecast with a 3 h timestep was used in this work. The forecast contains information about outside temperature, pressure, humidity, cloudiness, wind speed and direction, rain, and snow. The first part of
Figure 8 shows the expected outside temperature for the upcoming 5 days. The second part shows the solar irradiation forecast that takes into account the expected levels of cloudiness. However, at the end of July 2020, there was no cloudiness expected, so the irradiation was not affected.
Figure 9 shows the overall power consumption prediction of each appliance as a stacked area plot.
The prediction of the model’s behavior depends on the precision of the forecast. Forty weather forecasts were saved during the summer in order to make an analysis of the forecast error. The forecast with a minimal future prediction (the first forecast) was used as the real value for all other forecasts, as it was not possible to make real measurements to make the comparisons. Such error is important to understand the possible precision of the predictions and whether a 5-day weather forecast might provide any substantial valuable information. The forecasts were saved between the 22nd of July and the 30th of July 2020. The forecasts were downloaded every 3 h during these days. The total number of forecasts for comparison was 39. Because any real measurement was not available, the first forecasted value was considered the real value.
The prediction error tends to be quite large for the five-day forecasts but then the error stabilizes at a value around 0.6 °C. The changes in the error can be seen in
Figure 10. which is a box plot where each box represents an error distribution for a specific time sample in all 40 forecasts.
Figure 11 shows the distribution of errors. Finally, what different forecasts looked like in comparison to the first forecast (as described above) temperature can be seen in
Figure 12. It shows that the forecasts have a tendency to overestimate the maximums during the day as well as underestimate the minimums during the night.
2.3. Energy Demand Optimization
The biggest part of the overall energy consumption of the building is caused by heating and cooling, which correspond with the desired indoor temperature that the controller tries to achieve. An optimization strategy based on the compromise between comfort and price was implemented into the model.
The cost function is composed of two terms, namely, the deviation from the comfort temperature and the price paid for the electricity. It is expected that the owner would define the desired indoor temperature. One constant desired temperature was set for the whole horizon in the model. This defined desired temperature was then considered as the best achievable comfort. This comfort has some value for the owner, and deviating from it has a ‘price’ weighted by parameter r.
When the main aim of the optimization is to lower the price paid for the electricity consumed by the HVAC, the best result is obtained by turning off the HVAC. When comfort is of the utmost importance, then the price paid for the HVAC is defined only by the circumstances such as outdoor temperature, solar irradiance and the quality of control by the controller. If the owner is willing to make a compromise with the comfort, then the optimization finds the balance between saving money and keeping the comfort within some bounds. These bounds are defined by the user as an acceptable deviation from the desired temperature.
Such an optimization approach falls into the category of multi-objective optimization [
41]. The problem of multi-objective optimization can be described as finding an optimum for a set of objective functions that are to be optimized simultaneously. This optimization problem can be generally described as
where k is the number of objective functions and X is the search space or solution space defined by constraint functions such as inequalities or equalities.
These objective functions represent conflicting objectives because improving one objective leads to degrading the other. This means that there is no single solution that would be optimal for all objective functions. A solution then consists of some trade-offs. Such a solution is called the Pareto optimal, when moving away from this solution might lead to improvement in one objective function but degrade the other [
41].
It is thus clear that there is no single solution to multi-objective optimization, but rather a set of equally optimal solutions. Such a set is called the Pareto front [
42]. Finding a solution to a multi-objective optimization problem can be achieved in different ways as overviewed in [
41]. One possibility is to have no information about the preferences among the objectives. In such cases, the method is about finding a set of neutral compromises between the objectives. Another possibility is to subject the solutions to preferences, which can be set by the user or the decision maker before the optimization runs. These methods are called a priori. The optimization runs first and finds a set of Pareto optimal solutions from which the user or the decision maker then chooses one with regard to the preferences. This method is called a posteriori.
This work uses a weighted sum method that transforms the formulation (8) into a composite objective function [
41]:
where
wi is the weights;
fi is the objective functions;
k is the number of objective functions;
x is the vector of variables.
In this paper, the composite objective function is thus defined by a weighted sum formulation:
where
r is the weight of comfort;
N is the number of samples from the weather forecast;
Tdesired_by_user(k) is the desired temperature defined by the user, °C, in the k-th time step;
Tdesired_delta(k) is the optimized desired temperature, °C in the k-th time step;
P(k) is the price paid for heating and cooling, in EUR, in the k-th time step.
In this case, the variable is the real-valued vector containing the differences of the desired temperature, T
desired_delta, in the time steps of the forecasting horizon:
The price for heating and cooling is calculated by:
where
P(k) is the price paid for heating and cooling, in EUR, in the k-th time step;
pconsumed is the power consumed by heating and cooling, kW;
pkWs is the price per kWs consumed, EUR/kWs;
Ts is the sampling period, s.
The price paid for heating and cooling does not have a corresponding weight, or rather, the weight is considered to be equal to 1. This approach would fall into the category of ranking methods, a subset of rating methods, where one of the objectives is given a fixed weight of 1 [
41].
As also mentioned in [
41], the absolute values of the weights might not represent the actual preferences of the user because the magnitudes of each objective function might vary greatly, thus affecting the resulting composite objective function differently. In such cases, the weights might compensate for these different magnitudes of each objective function more than reflecting the preferences of the user.
The optimization calculates Tdesired_delta so that the objective function in (13) is minimized. Two scenarios have been considered: one in which the value of Tdesired_delta is constant for the whole forecast horizon (one variable for the optimization problem), and another in which there are three values of Tdesired_delta used at different times each day (three variables for the optimization problem). These times reflect the presence of the residents inside the home. The value of f is calculated for the whole forecast horizon. If r is equal to 0 then the resulting change in the desired temperature corresponds to minimizing the price on the forecast horizon, so the optimization becomes a parametric optimization instead of a multi-objective optimization. Such a case shows the maximal potential for a decrease in price on the forecast horizon while the optimization is constrained to finding one or three values of Tdesired_delta_k.
Each T
desired_delta_k is considered to be bounded. In this model scenario, they are bounded differently.
These bounds are used at different times during the days of the forecast. One is used when the user is present, another when he is not, and the third is used at night. The bounds, as well as the timings, can be set by the user and reflect the maximal distance from the desired temperature that the user is willing to accept. These bounds were considered to be fixed, but it is possible that the user’s preferences or willingness to accept discomfort might change during the day. Such changes in the preferences might be reflected in changes to the bounds as well as the weights of the weighted sum (13).
The weighted sum (13) acts as the utility function, which is supposed to approximate the user’s preference function. Such a function would represent the user’s satisfaction with changes to the value of each objective function. As mentioned in [
41], the utility function used in the form of a weighted sum is a linear approximation of the preference function. However, the preference function might be non-linear, meaning the preferences might change in the criterion space, for example, if one of the objective functions reaches a certain value.
As the objective function is calculated using an iterative algorithm with non-linearities (i.e.,: constraints on the value of u), an iterative method is used for the solution of the optimization problem.
The function fminbnd in MATLAB has been used for the optimization [
43]. The algorithm is based on a golden section search and parabolic interpolation [
43]. This function is a one-dimensional minimizer that finds a minimum for a problem specified by
where
x1 is the lower bound;
x2 is the higher bound.
The main idea behind the optimization in this paper was to show how the predicted energy price changes with respect to different values of the desired indoor temperature. This reflects a compromise between comfort and price. The pricing range can best be shown by calculating the optimization for different values of the weight of comfort.
The resulting pricing range can be seen in
Figure 17. It is divided into three parts. The first graph shows how the price changes with regard to the value of the comfort weight. The second graph represents a distance from the ideal comfort for different values of comfort weight. It is calculated as the sum of the difference between the desired temperature originally set by the user and the optimized desired temperature. Moreover, the third graph shows how each of the three optimized desired temperature deltas changes for different values of the comfort weight.
It can be seen that there is a compromise between the price paid and comfort. By increasing the value of the weight, comfort gains more importance, and the price rises. These graphs also show that there are 4 or 5 main areas the optimization leads to.
The desired temperature deltas in the last graph were defined as three deltas used at different times of the day. The different timings represent the times when the user is at home, and when the user is not present, and the last timing represents night-time.
The temperature delta, when the user is not present, was used during the day from 11.00 to 17.00. Its range was between −4 and +4 °C. The delta did not change very much, which was caused by the high temperatures and sunshine during the day. So, the optimization found that changing this temperature delta did not lead to a significant change in the power consumed by cooling since when increasing the desired temperature, the cooling still ran at the maximal level.
The temperature delta when the user was present was used at hours from 9.00 to 10.00 and then from 18.00 to 23.00. Its range was set to boundaries at −2 and +2 °C. When the weight was low, indicating that comfort was not important, this delta was set to −2 °C. This means that when the desired temperature set by the user was lowered to −2 °C, that would be 20 °C in the model scenario. As the weight of comfort increased, this delta lowered until it was 0 °C. In this case, the desired temperature was 22 °C, which was the same as defined by the user.
The third temperature delta was in the night-time frame and was used between 00.00 and 8.00. Its range was set to boundaries at −3 and +3 °C. This delta stayed at the value of −3 °C, resulting in the desired temperature of 19 °C for the first half of the values of the weight.
As the weight increased, comfort gained importance, which led to all desired temperature deltas being 0 °C. This is exactly the temperature that the user set as the ideal one.
The final price range obtained was from 33.44 EUR to 42.36 EUR for heating and cooling over 5 days for the modeled house.
So, such space of possible compromises could then be used for the a posteriori evaluation of a given weather forecast. The possibilities would change with the weather forecast.
A clearer example from the results would be the following. The user can take a look at the graph of the prices and the desired temperature changes for that given price. In this scenario, the user could choose to pay 35.40 EUR with the compromise of changing the desired temperature when he is present by −0.32 °C when he is not present by +0.15 °C and by −3 °C at night.
The model used in the optimization process had 10 m2 of photovoltaic panels. Increasing their area or also considering flexible energy tariffs might lead to different results.