Calculating retail prices from demand response target schedules to operate domestic electric water heaters

The paper proposes a demand response scheme controlling many domestic electric water heaters (DEWHs) with a price function to consume electric power according to a target schedule. It discusses at length the design of an algorithm to calculate the price function from a target schedule. The price function is used by the control of each DEWH to automatically and optimally minimize its local heating costs. It is demonstrated that the resulting total power consumption approximates the target schedule. The algorithm was successfully validated by simulation with a realistic set of 50 DEWHs assuming perfect knowledge of parameters and water consumption. It is shown that the algorithm is also applicable to clusters of large numbers of DEWHs with statistical knowledge only. However, this leads to a slightly higher deviation from the target schedule.


Introduction
In electricity networks, suppliers have to inject the amount of electricity that is consumed by their customers at the same time. Therefore, the customers' consumption is estimated for each time interval of 15 min of the next day. The supplier then instructs power plants to produce this power, by sending each a so called target schedule containing the amount of power to be produced for each time interval (Konstantin 2013).
This scheme of controlling production according to expected consumption is questioned by increasing amounts of electricity produced by wind and solar farms according to weather conditions and not according to desired schedules. Thus, demand response (DR) approaches were developed to adapt the course of consumption to that of availability. One option is providing customers with time dependent prices for the day ahead. This gives customers an incentive to shift their power consumption to times with lower prices corresponding to times of higher power availability (Vardakas et al. 2015). The approach protects the customer's privacy by design as no central entity has knowledge about the customer's flexibility in power consumption. However, it requires that suppliers can fix prices leading to given consumption schedules.
A domestic electric water heater (DEWH) is good candidate for a load performing DR in a household. Its thermal energy storage allows shifting its electrical consumption in a way that is not noticed and does not require intervention by users. DEWHs consume a major share of domestic electricity. In Germany, it is about 9.3% of the electrical energy of all households, even though only about 15% have one (Stryi-Hipp et al. 2007).
This paper discusses the design of an algorithm that implements DR by calculating a price function for a given target schedule that makes a set of DEWHs to consume power approximately according to the schedule. The algorithm is validated numerically by simulations. In a first step, perfect knowledge about every DEWH is assumed. In a second step this assumption is relaxed by only using statistics about clusters of many similar DEWHs, bringing the approach closer to practical application.

Retail pricing for demand response
In liberalized electricity markets, a supplier is an enterprise that purchases electricity at wholesale (or by operating power plants) and resells it to end customers, using networks operated by other enterprises (Lübkert et al. 2018a). The focus of this paper is the supplier's ability to control the consumption of its customer's DEWHs by using time dependent retail prices. Therefore, the supplier calculates a price function for the day ahead from a target schedule of the desired power consuption. The prices are sent to the customer's DEWHs, which automatically shift heating times to minimize their cost. This makes the resulting total consumption of all DEWHs follow the target schedule. The approach differs from most schemes in literature that maximize the suppliers profit only (Meng and Zeng 2013;Zugno et al. 2013;Wei et al. 2015;Kovacevic et al. 2017;Meng et al. 2017) or do not use retail prices at all (Binding et al. 2013).
In every time interval of 15 min (called a slot), a supplier has to inject the same amount of power into the network that is consumed by its customers. If there is an imbalance, the supplier has to buy balancing energy from the network operator (Konstantin 2013). Regulations force suppliers to minimize imbalances 1 . Therefore, suppliers perform balancing in a planning process on the day before power flow.
Suppliers use two means for planning, these are estimations of uncontrollable power flows and schedules to adjust controllable power flows. Estimations are mainly used for the power consumptions of customers for each slot of the next day. Essentially, measurements from the previous year or general consumption statistics a.k.a. standard load profiles are considered for this purpose (Konstantin 2013). Furthermore, the power produced by wind and solar farms is estimated based on weather forecasts. Schedules are vectors containing the power progression as one value per slot. Suppliers use it to instruct power plants to produce required amounts of power in each slot. In the future, suppliers may also prepare schedules to instruct DR systems to consume a desired consumption.
The use of varying retail prices is a well-known scheme to incentivize customers to participate in DR schemes. The price functions used by suppliers follow specific schemes as documented in Vardakas et al. (2015). With real-time pricing (RTP) the supplier regularly updates the prices and notifies customers some time before these are valid. A subscheme is day-ahead RTP (DA-RTP), where all prices of a day are fixed on the day before power flow.
For suppliers, RTP produces new challenges for estimating consumptions. Traditional estimations based on consumptions of the previous year or general statistics are no more applicable, as consumption progressions differ between years. This is intended for DR, as consumptions shall follow power availability, which also differs between years, e.g. because wind and solar farms lead to a correlation between availability and weather conditions. Hence, suppliers need new means to estimate their customers' consumption for given retail price functions or they need means for controlling DR loads directly, leading to known consumptions. Note that DA-RTP is different from day-ahead auctions, used for trading electricity at wholesale energy exchanges. At exchanges, trades are only concluded for matching bids (Konstantin 2013). As trades are always concluded some time before power flow, the amount of power to be delivered or consumed is always known accurately, not requiring estimations.
If a supplier wants to influence its customer's consumption with RTP, it has to find prices that make its customers shift their consumption as desired. Finding optimal retail prices can be formulated as a bilevel problem (Dempe 2002), in which customers minimize their electricity costs based on received prices and the resulting consumption impacts the supplier's objective and thus the prices. Other works often consider maximizing the supplier's profit resulting from the difference of their income for energy sold and its procurement costs at day-ahead markets (Meng and Zeng 2013;Zugno et al. 2013;Wei et al. 2015;Kovacevic et al. 2017;Meng et al. 2017). To estimate the resulting procurement costs, forecasts of exchange prices can be considered (Zugno et al. 2013). Alternatively, prices can be expressed as a function of time (e.g. for each hour) and the required energy quantity (Meng and Zeng 2013;Kovacevic et al. 2017;Chang et al. 2013). The available power of the supplier may be complemented by additional power sources, e.g. renewable sources (Chang et al. 2013) or battery storage (Wei et al. 2015;Kovacevic et al. 2017). Furthermore, positive or negative amounts of balancing energy have to be purchased during power flow (Zugno et al. 2013;Chang et al. 2013). Chang et al. considers this as part of the supplier's objective, e.g. when trades are already concluded (Chang et al. 2013). If the available power is known in advance, the problem becomes a load shaping problem, where the aggregated load shall follow a target schedule (Lübkert et al. 2018a;Hunziker et al. 2017).
Finding the optimal price function to be sent to all customers is a complex task, unless customers' utility functions are convex and differentiable, which allows to analytically determine strategies (Yu and Hong 2017). A simple solution is to use a price that is inverse to the power availability of the supplier. However, this leads to synchronized load peaks and thus high costs for power balancing during the day (Chang et al. 2013). To reduce the peaks, the loads need to distribute their energy consumption. This can be achieved by exchanging messages between customers (Chang et al. 2013) or with the supplier (Yu and Hong 2017;Safdarian et al. 2014 andSafdarian et al. 2016). Such an approach allows iteratively shaping the aggregated consumption profile directly or indirectly by updating the prices (Kovacevic et al. 2017). Alternatively, households may directly participate in micro markets as prosumers (Horta et al. 2017), buying and selling energy from and to neighbors or operators. When behavioral knowledge of the customers is given, the bilevel problem can be solved by converting it into a single-level problem (Zugno et al. 2013;Wei et al. 2015;Kovacevic et al. 2017) or by heuristic optimization such as genetic algorithms (Meng and Zeng 2013;Meng et al. 2017). None of these approaches constructs prices making the total consumption follow a target schedule.
Load shaping regarding a target schedule becomes simpler when retail prices increase with larger power consumption by using inclining block rates (IBR), as considered in Hunziker et al. (2017). As customers face stepwise increasing costs depending on the power level, they prevent load peaks to reduce costs. Therefore, this approach allows an iterative price adjustment to find a price function inducing the desired load shape. However, it is not necessary to penalize customers' peak load (e.g. cooking) if the peaks of all customers are sufficiently distributed. The approach of this work thus considers the variance of customers' appliances in order to determine prices instead of using IBR.
The effects of changes in retail prices on the power consumption pattern of DEWHs were analyzed in Lübkert et al. (2018a). The analysis used a bilevel model with full knowledge about every DEWH and assumes that each individual DEWH perfectly minimizes its costs with respect to the provided prices. It was shown that resulting price functions cannot operate a single DEWH according to a target schedule. However, well-defined price changes bring DEWHs into one of two states with significantly different power consumption. It was conjectured that this may allow constructing a price function that makes many DEWHs to consume power according to a given target schedule provided that the power consumption in both states can be calculated. This paper discusses the design of an algorithm doing so.

System model preserving privacy
This paper uses a bilevel model to control DEWHs with DA-RTP prices, which is based on the model from Lübkert et al. (2018a). The supplier is the leader. In its planning process, the supplier determines a target schedule of the power to be consumed by all DEWHs on the next day. From the schedule, it fixes one retail price for each slot. The DEWHs are the followers automatically minimizing their costs while considering their users' constraints. Therefore, the local control of each DEWH receives the prices and solves a linear program (LP) to find the heating times leading to its lowest costs for its known hot water consumptions. The supplier's goal is to minimize the mean square error between target schedule and total consumption of all DEWHs.
Adjusting the power consumption of DR loads on the basis of target schedules fits well with the common planning process used by suppliers. As discussed in the previous section, schedules are commonly used to adjust controllable power flows such as the power produced by power plants. Suppliers could use DR loads for controllability in a similar way, just with power flows having a reversed sign. However, this requires that suppliers can make the loads follow the target schedule accurately.
The following assumes DEWHs to be controlled with prices, because this protects the privacy of its users. A common alternative is the control of the DEWHs' heating directly by the supplier or some operator of a virtual power plant (You et al. 2009). As the supplier can then choose all heating times, it can ensure that the schedule is adhered to. However, to comply with the requirements of each user, the supplier needs to know each user's flexibilities. For DEWHs it needs to know the expected hot water consumptions of the user that is related to its personal behavioral pattern. This is private information that should be protected as much as possible. On the other hand, if decisions are made locally based on changing retails prices, knowledge about the user's flexibility is needed in the DEWH only. The supplier is still informed about the heating decisions of each DEWH, as the power consumption per slot is needed for accounting. However, details of the times of hot water consumption can be concealed.

Follower's problem of DEWHs
Today, DEWHs permanently keep the stored water at a preset temperature by using a thermostat and a heating element at the bottom of the tank. Whenever the temperature falls below the thermostat's dead band due to losses to the environment or drawing water, the heater is turned on until the upper temperature bound of the dead band is reached.
In order to shift the DEWH's energy consumption for DR, the water can be preheated to higher temperatures by using a more sophisticated controller. However, higher temperatures lead to increased energy losses to the environment due to the imperfect insulation. Preheating thus requires a sufficient price difference to financially compensate additional losses.
The simplified thermal model for a DEWH described by differential Eq.
(1) allows formulating the LP in (5)-(7) to minimize the operational costs with time-varying electricity prices p i as developed in Lübkert et al. (2018a) and Kepplinger et al. (2015). The LP requires considering discrete time, by dividing the optimization horizon into time slots of equal length t and solving (1) for each time slot with constant heating power P, heat conductivity G, water draw power equivalent HW and environmental temperature T env . The heat capacity C results from the product of the volume and the specific heat capacity of water. The discrete solution is given by (2)-(4), which describes the temperature change within a time slot. The water temperature T i converges exponentially to T ∞ , with the temperature change rate τ defining the duration of changing about 63.2% towards T ∞ (see (2)). Thus, the temperature changes faster for smaller τ -values and vice versa.
The decision variables h i determine the portion of each slot used for heating. It allows high precision of the solution while keeping the complexity low by using a LP with less time slots compared to an integer linear program (ILP) considering binary heating decisions. However, the heating element can only be turned on or off. The resulting average power h i · P heater can thus only be achieved with many DEWHs equally distributing their heating times in the slot. Furthermore, the solution must consider the constraints of the user regarding minimum T min and maximum temperature T max . A detailed description of the LP is given in Lübkert et al. (2018a).

Leader's problem of supplier
To be able to produce day-ahead real-time prices, the supplier has forecasts of the available power and the demand of price insensitive loads for the next day. Additionally, it has a lower estimation of the total energy demand required by all DEWHs. The energy demand can be distributed over the day within the limits of the user's constraints. The supplier can choose the target schedule such that the remaining available power is utilized in a beneficial way, i.e. by preventing down-regulation of renewable energy sources or buying additional energy at high costs. The procedure of selecting an optimal target schedule is out of scope of this work. The supplier induces the DEWHs to follow the target schedule by sending a price function to all customers. The price function consists of discrete time-dependent prices and is implemented as price vector p. To protect the customers' privacy, the supplier is assumed to have statistical knowledge about hot water usage and device parameters only. However, the subsequent discussion initially requires the supplier to have full knowledge of hot water usage and device parameters. This is relaxed to statistical knowledge later. The supplier applies the pricing algorithm proposed in this work on the selected target schedule to calculate the closest possible aggregated consumption profile and the related price vector p. The supplier uses unidirectional communication to the customers sending p as only shared information. The customers do not reply with resulting schedules.
The objective of the supplier is to find the price vector p minimizing the difference of the target scheduleĥ to the aggregated heating profile of all DEWHs (8). The objective extends (Lübkert et al. 2018a), where a single DEWH is considered. Prices are constant during each slot. The individual heating profile vector h k of DEWH with number k results from minimizing its operational cost for the given price vector (9)-(11).
Simple pricing schemes The model was validated for simple pricing schemes, which are constant prices and prices that are inverse proportional to the target schedule. Both schemes lead to huge deviations from the target schedule as shown in Fig. 1. This motivates the need for a more suitable pricing algorithm as proposed in this work. Using a constant price makes the heating profile independent of the desired target schedule h. The cost-optimal heating profile h const maintains the minimum temperature and heats as late as possible. As shown in Fig. 1, the energy consumption is lower than desired at night and higher during the evening. The realistic target schedule h as well as the DEWH setup are the same as used for validating the proposed pricing algorithm later in this paper. With inverse proportional prices, loads are shifted to times of high availability. However, only few hours are used by all DEWHs to preheat to maximum temperature when prices reach a local minimum. Thus, the heating profile h inv is far lower than desired at most times and exceeds the desired demand by a factor of about four in three hours. Due to the increased temperature, the total energy demand also increases resulting in an average daily heating duration of 3.44 h instead of 3.3 h at minimum temperature of 40°C. The increase is comparable to the energy demand which is required when maintaining the average temperature of 55°C.

Impact of prices on cost-optimal heating profiles
Creating price functions to induce an energy consumption following a target schedule requires understanding the impact of prices on the cost-optimal solution and predicting the resulting heating profile of a DEWH. This section describes how prices can be calculated such that a specific DEWH changes its target temperature and thus the heating state at a certain time. Furthermore, it is explained how the heating rate is chosen in each individual hour, allowing to predict the total heating profile.

Changes of heating state
In order to understand how price functions impact the cost-optimal heating profile, the bilevel problem (8) of finding a price to shape the heating profile of a single DEWH was analyzed in Lübkert et al. (2018a). It was understood that the resulting optimal price function allows the solver to select many heating profiles resulting at the same costs. The reason is that the price compensates the costs for earlier heating at any time slot. For this purpose, the price can be formulated as an exponentially increasing function (12) using the temperature change rate τ of the considered DEWH as a parameter called τ p = (τ p 1 , . . . , τ p n ). Such a price determined with τ p i = τ for all i using (12) leads to a constant effective price determined with (13) of p eff,i = p ref e n t/τ . The effective price represents a modified price information incorporating additional costs for higher standby losses due to preheating to higher temperatures (Lübkert et al. 2017). A constant effective price allows arbitrary heating profiles fulfilling the constraints and leading to minimum temperature at the end of the last slot.
In contrast, if parameter τ p i is smaller or larger than the τ -value of the DEWH, the heating profile will deterministically either heat to the maximum or it will always maintain minimum temperature respectively. Preheating is always beneficial when prices are generated with τ p i < τ, as this leads to a monotonic increasing effective price, whereas τ p i > τ leads to a monotonic decreasing effective price and thus prevents any preheating. Thus, the heating state H i in slot i can be determined directly as shown in (14).
A change of the heating state from maintaining minimum (H i = MIN) to heating to maximum temperature (H i+1 = MAX) leads to a local minimum in the effective price (13) and thus to a peak in the resulting heating profile as a rise from minimum to maximum temperature is required to fulfill the upcoming demand. Changing the heating state from MAX to MIN generates a local maximum in the effective price, which stops the heating until the temperature reaches the minimum again.

Prediction of cost-optimal heating profiles
Price functions generated with (12) cause changes in the heating state H i of a costoptimizing DEWH depending on the temperature change rate τ , as discussed in the previous section. In order to induce an aggregated heating profile of many DEWHs, the impact of such state changes on the heating profile needs to be known in detail. This section extends the findings of Lübkert et al. (2018a), which only considered monotonic decreasing τ p i values to induce peak loads at certain times. In principle, whenever the heating state changes, the DEWH will either start or stop preheating. Before the start of a preheating phase (H i = MIN), the DEWH will heat as late as possible and thus heats the demand of the slot to maintain minimum temperature and only preheats when heating the full slot is insufficient. At the start of a preheating phase (H i = MAX) the heating rate will be at a maximum level, e.g. heating the full slot (h i = 1). However, the actual heating rate depends on the maximum required temperature, which depends on the demand of later slots and also on later heating states as well as on the effective price (13).
The required temperature is selected such that the minimum temperature will be reached at the next local minimum effective price which is the start of next heating phase. When the maximum temperature is reached during the preheating phase, the temperature will be maintained until the effective price is below the next local minimum. Thus, any drawn energy will be restored by heating, similar to maintaining minimum temperature, but with higher standby losses. If the required temperature cannot be reached within a single slot, the heating rate will continue to be higher than the demand to reach the temperature. A different behavior occurs, if the effective price of the slot before the heating phase starts is below the price of the slot after the start (p eff,i−1 < p eff,i+1 ). Then, the heating rate is selected such that the resulting temperature allows to reach the required temperature within the next slot. Figure 2 shows an example LP solution for a price, which was generated to induce two changes from heating state MIN to MAX. The resulting heating and temperature profile is shown at the top and the selected τ p as well as the resulting effective price at the bottom of the figure. The local minima of the effective prices are at slot 5, 16, and 24. At both state changes (first two minima), the heating rate is set to maximum value (h i = 1). A preheating before the state change occurs in both cases, as the preceding effective price is slightly lower than the following price. The maximum temperature is maintained in slot 6 and 7 until it is sufficient to fulfill the energy demand until the next state change. The second preheating phase heats below the maximum allowed level and thus may turn off heating until the very last slot, where the effective price is slightly smaller than in slot 15 and thus the demanded energy will be restored by heating as it occurs to maintain minimum temperature.

τ -price algorithm
Results from the previous section allow constructing an algorithm, which generates a price function that induces an aggregated heating profile h of a known set of DEWHs following a target scheduleĥ. It exploits the effect that price functions resulting from (12) determine the heating state of a DEWH depending on the relation of the τ -value of the DEWH and of the prices. Thus, the solution space is reduced to prices corresponding to (12) and the decision variable p of the leader's problem is replacedwith τ p .
In principle the proposed Algorithm 1 utilizes a hill climbing approach by increasing or decreasing the number of DEWHs which are in the heating state MAX or MIN at each time slot i. The direction of change depends on whether the heating profile resulting from the current price is below or above the desired profile. As a change of the heating state at one time slot may impact the heating profile in a later or earlier slot, the algorithm iteratively tries to improve the aggregated profile and only considers changes that reduce the total error within the regarded horizon. The algorithm terminates when the error converges with respect to some precision (see lines 5 and 14). In Algorithm 1, planned state changes are represented by the peak plan vector , which is initially zero for each slot and thus all DEWHs remain in heating state MIN. A positive value i will induce additional i DEWHs to change from MIN to MAX heating state in slot i, whereas a negative value will do the opposite. The sum of all -values until time slot k thus determines the number of DEWHs which are in heating state MAX at that moment and must be positive at all time slots.
For each time slot in the horizon, the algorithm determines the total error resulting from the current peak plan and adjusts the value at time slot i in the direction of the error at slot i (line 8-10). Because of the iterative approach, not only current bounds (line 11) but also the subsequent bounds have to be ensured (line 12). For this purpose, the cumulative sum of the peak plan is determined. The minimum and maximum number of devices is applied in each time slot and then converted back into the adjusted peak plan by taking the difference (see line 35-44). If the new resulting total error is below the previous achieved error, the adjusted peak plan˜ is accepted or otherwise rejected (line 14-18). Incrementing the peak plan only by one step per time slot avoids an early depletion of flexibility, as once all DEWHs are in MAX heating state, the aggregated heating profile cannot be further increased in the following slots.
In order to determine the aggregated heating profile and the corresponding error (line 30-34) regarding the desired profile, the peak plan is converted into a price vector p (line 23-29) which is applied to a method for approximating the resulting heating profile (line 32), e.g. by solving the LP for each DEWH or using some heuristic. The vector p results from applying vector τ p to (12). The values τ p i are determined from the vector τ containing the available values in descending order. To set m DEWHs into heating state MAX at slot i, τ p i is selected between τ m and τ m+1 (line 25). The border values τ 0 and τ d+1 are defined by some value significantly outside of the considered range of τ -values.
The available values within the vector τ need to be unique, because otherwise an increment of i would have no effect. In practice, adjacent values must have an adequate minimum difference, due to a limited precision of LP solvers. Otherwise a τ p value between two very close τ -values may induce an undefined behavior of two DEWHs instead of changing the heating state to MIN or MAX. Thus, considering DEWHs with non-unique or similar τ -values leads to the effect, that incrementing i actually changes the heating states of multiple DEWHs.
Algorithm 2 creates a vector τ regarding this limitation, given any vectorτ of sorted τ -values of considered DEWHs and minimum difference of τ . This is done by ignoring everyτ k which has a smaller difference to its preceding value. The drawback of this method is that the reduction may be unevenly distributed. This leads to a few larger gaps, were a change in the peak plan changes the heating state of many DEWHs at once. However, an alternative approach, where instead of the difference to the previous value inτ , the difference to the last value in τ is considered, may again lead to the case that the selected τ p is too close to the τ -value of a DEWH.

Evaluation for a single day
The performance of the proposed algorithm is assessed with a realistic setup of DEWHs for a fixed horizon of one day in this section as well as with a receding horizon over 6 days in the next section. The results are compared to alternative price scenarios of poor performance, using constant or inverse proportional prices shown in Fig. 1. For this purpose, the algorithm was implemented in Octave/Matlab using the Yalmip toolbox (Löfberg 2004;Löfberg 2017).

Scenario configuration
The proposed algorithm requires a set of DEWHs with a sufficient diversity of τ -values in order to produce reasonable results. The hypothesis is that the existing DEWHs in German households fulfill this requirement. For this purpose a sampling of devices available on the market from selected manufacturers (Löfberg 2014 shows that DEWHs of significantly different volume may have similar τ -values due to the diversity of the heat conductivity. However, although only a limited range of volumes is available, the considered device base shows that it is likely to have an adequate distribution of clustered τ -values for DEWHs on the market. This allows concluding that the existing DEWHs in Germany fulfill the requirement of a sufficient diversity in their τ -values. In this study a sample of 50 devices is constructed assuming a cropped discrete normal distribution of the volume as shown in Fig. 3b, with mean of 70 l and standard deviation of 70 l. The corresponding heat conductivity is selected consecutively from the available values. Additionally, an error between −2 to 2% is added to the volume with uniform distribution. The resulting setup is available at Lübkert et al. (2018b).
Each DEWH's daily energy demand is defined as twice its volume at minimum temperature T min = 40°C, which leads to an average heating duration of 3.3 h/day considering heating elements with rated power of P = 2kW and an inlet water temperature of T cold = 10°C. The demand is distributed throughout the day according to a representative average water consumption profile (Defra 2008). The allowed temperature range is from T min = 40°C to T max = 70°C and the environment temperature is considered as T env = 20°C. To simplify the model, the temperature parameters are assumed to be constant. The algorithm terminates when no further reduction of the error is possible ( = 0). A day is divided into n = 24 time slots of length t = 1h. The initial price p ref is set to 1$/MWh to simplify interpretation of resulting price values, as only the relation of prices is important, not their absolute value.
The desired aggregated heating profileĥ is constructed with (17) incorporating the minimum total energy demand (15) of all DEWHs required for maintaining the minimum temperature. The shape ofĥ follows a virtual power availabilityP with (16), by inverting and normalizing electricity exchange prices, as in Lübkert et al. (2018a). Exemplary prices p ex are taken from Phelix at EPEX SPOT (Energy Monitoring Company and Energy Saving Trust) for January 1, 2016.

Numerical results
With the prices generated by the proposed algorithm, the aggregated energy consumption follows the desired shape of the availability profile. The energy consumption is shown in Fig. 4a, with the aggregated heating profiles (target schedule and LP solution) at the top, the resulting p with corresponding τ p in the center, and the individual heating profiles h k of all considered DEWHs in a heat map at the bottom. The considered 50 τ -values of the DEWHs are reduced to 34 clusters by Algorithm 2 using a minimum difference of τ = 1%, as shown in Fig. 4b.
The algorithm terminates after 6 iterations. The most significant error reduction is achieved within the first two iterations. As shown in Fig. 4c, the root-mean-square  (18), which is a significantly improvement compared to an MAE r of 63.6% for constant and 113% for inverse proportional prices p ex .
During the first seven hours of the day, the algorithm increases the energy consumption following the increasing availability by changing τ p such that in each hour a few DEWHs change to heating state MAX. As the DEWHs have different volumes and water demand, the resulting consumption peaks differ in length and height. As the availability decreases from slot seven, the algorithm increases τ p again. As a result, some DEWHs rely on preheated water without need for further heating for few hours. Others continue heating to maintain either minimum or maximum temperature.
The example availability increases again in the afternoon, where the first devices already stop heating as the maximum temperature now is sufficient to serve the remaining demand of the day. Thus, τ p is decreased again for a short time. Afterwards, the availability quickly drops to a minimum value in slot 18 and the algorithm increases τ p accordingly. However, the aggregated consumption cannot be reduced sufficiently, because many devices did not preheat and still require to heat to satisfy the demand.
In the very last hour, the algorithms fails to increase the energy consumption, as there is no remaining demand to preheat and it is always optimal to finish at minimum temperature. In order to increase the flexibility and include the demand of the next day, a receding horizon can be used. This is done in the next section by sending new prices for the next day during the current day.
In this scenario about half of the DEWHs remain in heating state MIN all the time. These have the smallest considered τ -values, due to smaller volumes or a less efficient insulation. With an average heating duration of 3.38 h, the total energy demand is lower than in the scenario using exchange prices or constant prices when maintaining medium temperature. Although DEWHs with better insulation or larger volume preheat to higher temperatures, they still have lower operational costs for the same energy demand. The average price paid for consumed hot water energy (19) increases monotonically with τ . The DEWH having the largest τ -value has to pay about 1.13 $/MWh, whereas the smallest τ -value leads to costs of 1.26 $/MWh in relation to the arbitrary reference price p ref .
At most times of the day, the resulting heating profile follows the availability favorable, except at slot six, where a sudden decrease occurs. Although an increase of the peak plan in that slot could reduce the local error, it would not reduce but increase the total error. However, it is possible to find a combination of increase and later decrease which would improve the total error. Thus, results could be further improved by a more complex postoptimization method, e.g. using a recursive approach which first reduces the local error and then tries to find an additional change of the peak plan to reduce the total error again. Alternatively, only the local error could be considered to modify the peak plan, which may lead to better results for the first part, but much higher errors at the end of the day. In combination with a receding horizon, this still can lead to better results, but requires some mechanism to ensure convergence. Further alternatives are using a more sophisticated logic, e.g. modifying around minimum and maximum error in one step, or a general heuristic optimization method like genetic algorithm or particle swarm optimization.

Receding horizon and price reset
The results from the previous section suggest that a receding horizon can improve the ability of the proposed algorithm to reduce the error at the end of a day, as the water demand of the next day can be considered for further preheating. Otherwise, all DEWHs will always end at minimum temperature at the same time, independent of the desired heating profile. However, by introducing a receding horizon, the prices would continue increasing exponentially, which is impractical as customers would not accept such prices. Thus, the price needs to be reset to a lower value after some time.
The price reset hinders the benefit of a receding horizon, as again all DEWHs will end at minimum temperature at that time. However, the moment of the reset can be selected such that its impact on the resulting heating profile is small. Because the price reset causes the same behavior as the end of the scheduling horizon, the price finding problem can be split into two independent horizons before and after the reset. Algorithm 3 utilizes this effect to calculate a price including a reset to a lower value.

Algorithm 3 τ -Price Heuristic with Price Reset
r ← index of last minimum value in cumulative sum of p 1..r+1 ← FIND_PRICE(ĥ 1..r+1 ) p r+2..n ← FIND_PRICE(ĥ r+2..n ) After determining the price for the whole horizon, the time slot r is selected to reset the price afterwards. A reasonable time to reset the price is where most DEWHs are at heating state MIN and thus the peak plan is at a minimum level. The last occasion of this level is selected to ensure sufficient flexibility in the time before. However, the reset must occur before the start of the next horizon, as otherwise it could be removed again. Alternatively, the reset needs to be enforced in the next schedule. The example shown in Table 1 leads to a price reset index r = 4, as this is the last minimum value in the cumulative sum of the peak plan . In the next step the price is determined again as a combination of two computed price functions for consecutive sub-horizons. The first sub-horizon ends at slot r + 1 and the second starts at r + 2.
The results for a receding horizon with price reset are shown in Fig. 5 considering power availability from January 1 to 6 respectively. Prices are updated every day at noon (12 p.m.) covering a horizon of 36 hours until midnight. The aggregated heating profile still follows the shape of the availability. However, there are some noticeable errors especially around the price reset where availability is low. Due to this effect, the error reaches an MAE r -value of about 20.8%, which is slightly higher compared to the performance for the first day, but significantly better than using constant prices (MAE r = 75.7%) or exchange prices (MAE r = 102%).
The price resets always occur close to midnight, due to the nature of the virtual power availability constructed from exchange prices, leading to minimum availability in the late afternoon followed by an increase at night. A high availability in a larger period or a very high availability in a short period leads the algorithm to select more devices to create peaks. This increases the overall energy demand and also the maximum price, as the price gradient needs to be steeper for DEWHs with smaller τ -values. Due to the price reset, the maximum price difference is about 25% for the considered DEWH setup.

Average optimization model vs. statistical behavior
In the previous sections the behavior of a set of known DEWHs was analyzed. However, in practice, each DEWH considered for price construction would represent the average behavior of a group of similar DEWHs for which only a stochastic water demand profile is known. In order to assess the applicability of this approach, a Monte Carlo simulation is conducted to analyze the additional error due to a variance of the water consumption patterns. For this purpose 1000 optimal heating profiles are calculated for a single DEWH cluster with a volume of V = 65 l considering different water demand profiles. The total water demand is randomly chosen with 2V ·f N [l at T min ], where f N is normally distributed with mean μ = 1 and standard deviation σ = 0.25. The total demand is distributed into 10 randomly chosen slots regarding the average profile. Figure 6a shows three example profiles on the left and a histogram of the total demand on the right.
The results of the previous section have shown that the most common individual heating profiles include a single pre-heating at beginning, two pre-heating phases, or maintaining minimum temperature during all slots of the horizon. The resulting heating profiles of the Monte Carlo simulation for these three scenarios are shown in Fig. 6. In Table 1 Example peak plan and cumulative sum leading to a price reset index r = 4 the case of maintaining minimum temperature (Fig. 6b), the average of all heating profiles is equal to the result of the average optimization model. Only small deviations occur, due to the randomization. Thus, in the considered scenario of the previous section, about 10-50% of the profiles can be considered as correct.
In the other two cases slightly larger deviations occur. In the case of continuous heating to maximum temperature (Fig. 6c) both results begin with a similar profile, but differ at the end. This is because, some of the consumption patterns have a focus in the morning, which allows to turn off heating earlier. Others have a major consumption in the evening, which still requires some heating at later times. In total the resulting profiles differ from the average model by stretching the demand in the second part of the day, which causes a shortfall first and then a surplus. In the case of a second pre-heating peak in the afternoon (Fig. 6d), this leads to significant reduction of the peak height and a low additional load before and after the peak. The first pre-heating peak remains equal to the average model, but also distributes the demand afterwards into the time between both peaks.
Overall, the resulting heating profiles from the average optimization model are a good estimate of many different optimized profiles. Especially the first part of the horizon fits a c b d Fig. 6 Comparison of averaged optimization and Monte Carlo simulation. a MC HW profiles, b maintain minimum temperature, c heat to maximum, d heat to maximum twice well, but at later slots additional errors are introduced as the energy consumption is higher than expected at times where it is intended to be reduced.
Deviations in volume or heat conductivity would lead to a different τ -value and thus are crucial, as this would change the number of devices considered for a cluster in the averaged model. Different temperature ranges have an impact on the peak height as DEWHs may either heat more or less in a time slot, which reduces the main peak, but increases the consumption before the peak. The rated power has a similar impact, as some devices require more time to heat the same amount of energy. However, in practice rated power is restricted to a few common discrete values.

Discussion
With the proposed algorithm, prices increase exponentially with a single reset to a minimum value per day. This successfully shapes the energy consumption of many DEWHs. However, it is unclear if the prices generated by the algorithm are applicable in practice.
Customers might accept such prices. They lower the operating costs of the DEWHs. The ratio between maximum and minimum price is only about 25%, which is far less compared to other approaches, e.g. up to roughly 400% in Hunziker et al. (2017). Furthermore, prices are similar every day and change very little between the hours. However, the same price function needs to be applied for many different appliances, unless each DEWH has its own electricity meter. In most households with a single meter measuring the entire household consumption, sending different prices to individual appliances is unrealistic. The integration of price insensitive loads (e.g. cooking, lights, multimedia) is still not a problem, as these behave the same for any price function. Its consumption can be forecasted as common today and subtracted from the target schedule. A bigger issue is the integration of cost optimizing shiftable loads (e.g. washing machines, dryer, charging electric vehicles). These devices will run as early as possible after the price reset to exploit the lowest prices.
The effort for integrating other thermostatically controllable loads (TCL) may be low due to their similar thermal model. In fact, the proposed algorithm only requires that TCLs know their τ -value and optimally minimize their operational costs according to it. However, if the loads use a different control mechanism, not fully optimizing for lowest costs, the resulting total consumption is different and deviates from the target schedule. Thus, loads may react with undesirable behavior, e.g. causing synchronized load peaks at lowest prices.
The major challenge with non-DEWH-loads is likely to occur around the price reset. Even though this point in time is selected to have low effect on the consumption of DEWHs, other loads will react on this large price change. Their consumption will be low before and large after the reset. In many cases this prevents following a given target schedule, even though DEWHs account for a large part of the domestic energy demand.

Conclusion
The paper showed that it is possible to construct an algorithm that calculates a price function from a target schedule operating many DEWHs to approximately consume power according to a schedule. With the setup used for validation, prices vary by about 25%. Prices increase exponentially but can be reset to some smaller value in the late evening after a period of low power availability. The prices can also be used for the more realistic case of clusters consisting of many similar DEWHs only described by statistics. However, this leads to a larger deviation from the target schedule.
Future work should quantify the deviation from the target schedule for multiple sets of DEWH clusters. Results need to be compared to the use of mathematical optimization including heuristic approaches. Moreover, numerical requirements should be examined in more detail, such as the minimum distance of the τ -values between clusters, the required precision for prices, and the smallest possible length of slots. Finally, the approach should be extended to further types of loads -thermostatically controllable loads and others -to study the feasibility to control all of these with the same price function.
Endnote 1 In Germany, this is regulated in the contract every supplier has to conclude with its transmission system operator. Its Section 5.2 obligates the supplier to keep deviations from the balance as low as possible by taking reasonable measures (EPEX SPOT).