MPC-based interval number optimization for electric water heater scheduling in uncertain environments

Abstract

In this paper, interval number optimization and model predictive control are proposed to handle the uncertain-but-bounded parameters in electric water heater load scheduling. First of all, interval numbers are used to describe uncertain parameters including hot water demand, ambient temperature, and real-time price of electricity. Moreover, the traditional thermal dynamic model of electric water heater is transformed into an interval number model, based on which, the day-ahead load scheduling problem with uncertain parameters is formulated, and solved by interval number optimization. Different tolerance degrees for constraint violation and temperature preferences are also discussed for giving consumers more choices. Furthermore, the model predictive control which incorporates both forecasts and newly updated information is utilized to make and execute electric water heater load schedules on a rolling basis throughout the day. Simulation results demonstrate that interval number optimization either in day-ahead optimization or model predictive control format is robust to the uncertain hot water demand, ambient temperature, and real-time price of electricity, enabling customers to flexibly adjust electric water heater control strategy.

Keywords

electric water heater / load scheduling / interval number optimization / model predictive control / uncertainty

Introduction

With the emergence of time-varying retail pricing schemes, more incentives have been provided by utilities to encourage consumers to participate in demand-side management. As an essential enabler for consumers to optimize the use of their smart home appliances, household load scheduling is vital in demand-side management [¹,²]. For consumers, household load scheduling is able to reduce their electricity payment; for utility companies, it can help reduce the peak-to-average ratio of power systems [³].

Many valuable achievements have been made in household load scheduling. In Ref. [⁴], an appliance commitment algorithm has been proposed to schedule thermostatically controlled household loads, and the electric water heater (EWH) studied as a representative. In Ref. [⁵], a traversal-and-pruning (TP) algorithm has been introduced for EWH load scheduling. Compared with the traditional algorithm, this novel TP algorithm is able to significantly reduce the electricity bill. However, in the above mentioned references, all forecast parameters in the scheduling have just been deemed as certain values, i.e., the forecast has been considered totally accurate, which is not true in real world. Generally, forecast errors are inevitable in some parameters such as weather condition, consumers’ behavior, and etc [⁶]. Therefore, it is of great significance to study uncertain optimization for household load scheduling. In Ref. [⁷], a comparison of different linear and nonlinear methods for home energy resource scheduling has been made, taking the presence of data uncertainty into account. In Ref. [⁸], in order to avoid the stochastic optimization from the curse of dimensionality, a two-stage stochastic optimization approach has been introduced. In Ref. [⁹], fuzzy programming has been introduced to tackle the uncertain optimization in home energy management, in which uncertain parameters like electricity prices and outdoor temperature are described by fuzzy parameters.

Although a great amount of achievements have been made in uncertain household load scheduling by using stochastic programming and fuzzy programming, yet further development is restricted by the inherent defects of stochastic programming and fuzzy programming. First, it requires a considerable amount of sampling information on uncertainty in constructing accurate probability distributions or fuzzy membership functions. However, it often seems difficult and expensive to obtain sufficient uncertainty information. Secondly, as both probability distributions and fuzzy membership functions need to be introduced into the optimization problem, the computational complexity will increase more or less and sometimes such an uncertain optimization problem is too hard to be solved [¹⁰].

In comparison to traditional optimization methodologies that deal with uncertainties, the application of interval number optimization, which has attracted increasing attentions recently, in household load scheduling has not yet been extensively studied. As a set-theoretical and non-probabilistic method, interval analysis needs only information about the bounds of the magnitude of uncertainties, not necessarily knowing the specific probabilistic distribution densities. As less uncertainty knowledge is needed, interval number optimization is suitable for addressing the household load scheduling with limited uncertainty information [¹¹]. Thermostatically controlled appliances are the most energy-consuming home appliances [¹²–¹⁴], hence, in order to simplify the statement, the EWH load scheduling is taken as an example to demonstrate household load scheduling using interval number optimization, in which three uncertain parameters (hot water demand, ambient temperature, and real-time price of electricity) are discussed. Note that few studies have been conducted which simultaneously take all the above three uncertain parameters into account in household load scheduling, as it is too complicated to analyze and compute those parameters in traditional methods. However, interval number optimization is able to create suitable load schedules which are robust to the uncertainties and flexible to diversified consumers’ demand.

In practice, because the forecast errors increase with the extension of the forecast time horizon [¹⁵], the uncertain parameters represented by the interval number which are given before the optimization are not accurate enough and the day-ahead scheduling is not good enough. In this paper, model predictive control is introduced to address this problem. In Ref. [¹⁶], a stochastic model predictive control (MPC) has been applied to handle the variability and uncertainties of the power of the PV panels. In Ref. [¹⁷], a closed-loop prediction formulation of robust model predictive control has been used to deal with the problem of robust model predictive control of a building envelope and a HVAC system. In this paper, EWH load scheduling using interval number optimization is formulated in both day-ahead and MPC pattern. It is in nature an integer nonlinear programming problem which can be solved by using the binary particle swarm optimization (BPSO). To be specific, first, the thermal dynamic model of EWH with interval numbers is established. Next, the comparison of interval number by using the possibility degree of interval number is quantitatively described. After that, the binary particle swarm optimization and model predictive control are introduced. Finally, the simulation results get presented and analyzed.

Establishment of the EWH thermal dynamic model with interval numbers

Traditional thermal dynamic model of EWH

The basic working principles of EWH are as follows: when the water temperature decreases to the lower limit of comfort zone, the EWH turns “on” to heat the water; when the water temperature rises to the upper limit, it turns “off.” Figure 1 shows the typical thermal characteristic curve of EWH without considering hot water use.

In essence, the heaters of EWH work because of two reasons: one is cold water inflows due to hot water demand; the other is heat losses from the heat exchange with the environment. In summary, hot water temperature is affected by parameters such as physical characteristics, hot water demand, ambient temperature, and etc. Physical characteristics could be described by equivalent thermal parameters (ETP) which are usually treated as constants. In contrast, the hot water demand and ambient temperature are hard to be predicted accurately in practical situations, as customer behavior is difficult to be foreseen precisely and the weather sometimes changes unexpectedly. If they are only described as constant values, the practical water temperature in the tank will violate the comfort zone frequently under the disturbance of prediction errors. Thus, it is necessary to consider water demand and ambient temperature as uncertain parameters.

The behavior of EWH could be formulated by the thermal dynamic model [⁴]. Without considering hot water use, the thermal dynamic model of EWH is given by

(1){\theta }_{n+\mathrm{1}}={\theta }_{e,n}+X_{n}QR-({\theta }_{e,n}+X_{n}QR-{\theta }_n)\exp [-(t_{n+\mathrm{1}}-t_n)/(RC)].

When the hot water demand is considered, the water temperature should be modified as

(2){\theta }_{n+\mathrm{1}}=[{\theta }_{\mathrm{cur}}(M-d_n)+{\theta }_{e,n}{d}_n]/M,

where θ_cur is the θ_n+1 calculated in Eq. (1).

Equations (1) and (2) can be expressed jointly as a nonlinear function of thermal parameters, ambient temperature, hot water demand, and on/off status in Eq. (3).

(3){\theta }_{n+\mathrm{1}}=f({\theta }_n,t_n,Q,C,R,d_n,X_n,{\theta }_{e,n}).

Interval number model

Interval number, which can be uniquely determined by its upper and lower bounds, represents a bounded set of real number between the bounds. The basic concepts of the interval number are introduced as follows.

For the arbitrary real numbers a^down and a^up which are subject to a^down<a^up, A^I = [a^down, a^up] denotes a normal interval number in which a^up is the upper bound and a^down is the lower bound. Besides, the interval number can also be described by the middle point and radius like A^I = [a^down, a^up] = [a^c, a^w], in which the a^c and a^w respectively represent the middle point and radius of the interval number [¹⁸].

The basic operations of interval numbers have been discussed in Ref. [¹⁹]. Moreover, to quantitatively describe the comparison of two interval numbers, the possibility degree of interval number is proposed [²⁰]. The possibility degree of a real number a and an interval number B^I is given by

(4)

The interval number B^I can be regarded as a random variable with uniform distribution whose upper and low bound is the investigated subject in the interval number model. Thus, the uncertain parameters like hot water demand and ambient temperature can be described as d_n^I=[d_n^{\mathrm{down}},d_n^{\mathrm{up}}] and{\theta }_{e,n}^I=[{\theta }_{e,n}^{\mathrm{down}},{\theta }_{e,n}^{\mathrm{up}}]. Based on this transformation, the thermal dynamic models of EWH can be rewritten with interval numbers [¹¹].

Without considering hot water use, the model is given by

(5){\theta }_{n+\mathrm{1}}^I={\theta }_{e,n}^I+X_{n}QR-({\theta }_{e,n}^I+X_{n}QR-{\theta }_n^I)\exp [-(t_{n+\mathrm{1}}-t_n)/(RC)],

where {\theta }_{n+\mathrm{1}}^I=[{\theta }_{n+\mathrm{1}}^{\mathrm{down}},{\theta }_{n+\mathrm{1}}^{\mathrm{up}}] and {\theta }_n^I=[{\theta }_n^{\mathrm{down}},{\theta }_n^{\mathrm{up}}] denote hot water temperature ranges at times t_{n+\mathrm{1}} and t_n; and the ambient temperature is expressed as {\theta }_{e,n}^I=[{\theta }_{e,n}^{\mathrm{down}},{\theta }_{e,n}^{\mathrm{up}}].

Similarly, the temperature modification equation considering hot water use can be rewritten as

(6){\theta }_{n+\mathrm{1}}^I=[{\theta }_{\mathrm{cur}}^I(M-d_n^I)+{\theta }_{e,n}^I{d}_n^I]/M,

where {\theta }_{\mathrm{cur}}^I=[{\theta }_{\mathrm{cur}}^{\mathrm{down}},{\theta }_{\mathrm{cur}}^{\mathrm{up}}] represents the {\theta }_{n+\mathrm{1}}^I calculated in Eq. (5).

Though Eqs. (1) and (2) look similar to Eqs. (5) and (6), they are quite different in nature. The former are ordinary algebraic operation but the latter belong to interval arithmetic operation. In interval analysis, the interval numbers are calculated through arithmetic operation rules of interval numbers which are essentially different from traditional algebraic operation rules. In the following, the basic operation rules of interval numbers are defined.

For arbitrary a^{\mathrm{down}},a^{\mathrm{up}},b^{\mathrm{down}},b^{\mathrm{up}}\in R is subject to a^down<a^up, b^down<b^up, A^I = [a^down, a^up], B^I = [b^down, b^up], the rules for interval addition, subtraction, multiplication and division are

(7)A^I+B^I=[a^{\mathrm{down}}+b^{\mathrm{down}},a^{\mathrm{up}}+b^{\mathrm{up}}],

(8)A^I-B^I=[a^{\mathrm{down}}-b^{\mathrm{up}},a^{\mathrm{up}}-b^{\mathrm{down}}],

(9)A^I\times B^I=\left[\begin{array}{l}\min (a^{\mathrm{down}}{b}^{\mathrm{down}},a^{\mathrm{down}}{b}^{\mathrm{up}},a^{\mathrm{up}}{b}^{\mathrm{down}},a^{\mathrm{up}}{b}^{\mathrm{up}}),\\ \max (a^{\mathrm{down}}{b}^{\mathrm{down}},a^{\mathrm{down}}{b}^{\mathrm{up}},a^{\mathrm{up}}{b}^{\mathrm{down}},a^{\mathrm{up}}{b}^{\mathrm{up}})\end{array}\right],

(10)A^I/B^I=\left[\begin{array}{l}\min (a^{\mathrm{down}}/b^{\mathrm{down}},a^{\mathrm{down}}/b^{\mathrm{up}},a^{\mathrm{up}}/b^{\mathrm{down}},a^{\mathrm{up}}/b^{\mathrm{up}}),\\ \max (a^{\mathrm{down}}/b^{\mathrm{down}},a^{\mathrm{down}}/b^{\mathrm{up}},a^{\mathrm{up}}/b^{\mathrm{down}},a^{\mathrm{up}}/b^{\mathrm{up}})\end{array}\right].

Formulation of EWH load scheduling

EWH load scheduling considering uncertainties

The primary objective of EWH load scheduling is to minimize the electricity bill band. At the same time, the temperature constraints should be satisfied even under the uncertain hot water demand and ambient temperature. The pricing mechanism of this paper is real-time pricing (RTP). RTP is difficult to be forecasted due to the constantly changing of power supply-demand relationship. Hence RTP is seen as an uncertain parameter which is also described by interval numbers as p_n^I=[p_n^{\mathrm{down}},p_n^{\mathrm{up}}].

Although the electricity payment is important, user comfort should be guaranteed in the first place. Assuming that the comfort zone is {\theta }_{\mathrm{WH}}^I=[{\theta }_{\mathrm{WH}}^{\mathrm{down}},{\theta }_{\mathrm{WH}}^{\mathrm{up}}], in which {\theta }_{\mathrm{WH}}^{\mathrm{down}} and {\theta }_{\mathrm{WH}}^{\mathrm{up}} denote the lower and upper bounds, if the water temperature exceeds the lower or upper bounds, the load schedule is deemed to be unacceptable.

Mathematically, the objective function is described as

(11)\min f_{\mathrm{original}}(X_n,p_n^I)=\min {\displaystyle \underset{n=\mathrm{1}}{\overset{N}{\mathrm{\Sigma }}}{p}_n^I{X}_n{P}_{\mathrm{EWH}}\cdot (t_{n+\mathrm{1}}-t_n)},

subject to the following constraints

(12){\theta }_{n+\mathrm{1}}^I=f({\theta }_n^I,t_n,Q,C,R,d_n^I,X_n,{\theta }_{e,n}^I),

(13){\theta }_{\mathrm{WH}}^{\mathrm{down}}\le {\theta }_n^I\le {\theta }_{\mathrm{WH}}^{\mathrm{up}},

where N is the number of all the time steps over scheduling horizon and P_EWH denotes the rated power of the EWH.

Transformation of EWH load scheduling model

The formulations of objective function and constraints with interval numbers are concise, but cannot be solved directly. Accordingly, the possibility degree of interval numbers is introduced to convert this uncertain optimization problem into a deterministic one.

First, in order to intuitively describe the mean value and uncertain level of interval numbers, the interval numbers with a middle point and radius are used to describe the objective function [¹⁹]

f_{\mathrm{original}}(X_n,p_n^I)=[f^c(X_n),f^w(X_n)],

in which f^c(X_n) and f^w(X_n) are defined as

(14)f^c(X_n)={\displaystyle \frac{\underset{P_n}{\max }f(X_n,P_n)+\underset{P_n}{\min }f(X_n,P_n)}{\mathrm{2}}},

(15)f^w(X_n)={\displaystyle \frac{\underset{P_n}{\max }f(X_n,P_n)-\underset{P_n}{\min }f(X_n,P_n)}{\mathrm{2}}},

where the \underset{P_n}{\max }f(X_n,P_n) and \underset{P_n}{\min }f(X_n,P_n) mean that the uncertain electricity price p_n^I in Eq. (11) is replaced by the maximum and minimum values of the predicted electricity price interval, respectively.

After that, the objective function is transformed into a bi-objective optimization which is min(f^c(X_n), f^w(X_n)) where f^c(X_n) and f^w(X_n) denote the middle point and radius of electricity bill band respectively. Then, the weight coefficient β is introduced to transform this bi-objective optimization into a single-objective optimization which is formulated as

(16)f_{\mathrm{final}}(X_n)=\beta f^c(X_n)/\delta +(\mathrm{1}-\beta )f^w(X_n)/\mathrm{phi;},

where δ and φ represent normalizing factors of the middle point and radius, which are used to unify the magnitudes of two numbers to avoid the elimination of the smaller one. β is determined according to the consumers’ preference for the middle point and radius of electricity bill interval. If β is assigned to be 1, it means that only the middle point of energy cost interval is valued; if β is assigned to be 0, it means that only the radius of energy cost interval is valued.

The possibility degree is utilized to deal with the uncertain constraints by introducing the tolerance degree for constraint violations. σ∈[0,1] quantitatively describes users’ attitude for the violation of water temperature bounds. As the tolerance degree increases, more violation is allowed. But Eq. (13) is a mandatory constraint which means that the water temperature interval {\theta }_{n+\mathrm{1}}^I is required to be strictly controlled in the comfort zone. This is too rigid for the majority of consumers, so Eq. (13) is transformed as

(17)\{\begin{array}{l}\begin{array}{ll}P({\theta }_n^I\ge {\theta }_{\mathrm{WH}}^{\mathrm{down}})\ge \mathrm{1}-{\sigma }_{\mathrm{1}},& {\sigma }_{\mathrm{1}}={\displaystyle \frac{{\theta }_{\mathrm{WH}}^{\mathrm{down}}-{\theta }^{\mathrm{down}}}{{\theta }^{\mathrm{up}}-{\theta }^{\mathrm{down}}}},\end{array}\\ \begin{array}{ll}P({\theta }_n^I\ge {\theta }_{\mathrm{WH}}^{\mathrm{up}})\ge {\sigma }_{\mathrm{2}},& {\sigma }_{\mathrm{2}}={\displaystyle \frac{{\theta }^{\mathrm{up}}-{\theta }_{\mathrm{WH}}^{\mathrm{up}}}{{\theta }^{\mathrm{up}}-{\theta }^{\mathrm{down}}}}.\end{array}\end{array}

The values of σ₁ and σ₂ denote the customers’ preference for the lower and upper limit of water temperature. When σ₁ >0 or σ₂ >0, it means that the customers tolerate a certain extent of comfort violation. If σ₁ equals σ₂, it means that the customers do not have an obvious preference for a higher or lower water temperature. If σ₁ is larger than σ₂, it means the customers prefer a lower water temperature to a higher water temperature. If σ₁ is less than σ₂, it means that the customers have a preference for a higher temperature.

Moreover, the penalty function method is utilized to transform the constraints as a term of the objective function:

(18)\mathrm{phi;}(P(M_{n_m}^I\ge N_{n_m}^I)-{\sigma }_{n_m})={(\min (\mathrm{0},P(M_{n_m}^I\ge N_{n_m}^I)-{\sigma }_{n_m}))}^{\mathrm{2}},

where M_{n_{\mathrm{1}}}^I={\theta }_n^I,N_{n_{\mathrm{1}}}^I={\theta }_{\mathrm{WH}}^{\mathrm{down}},{\sigma }_{n_{\mathrm{1}}}=\mathrm{1}-{\sigma }_{\mathrm{1}};M_{n_{\mathrm{2}}}^I={\theta }_n^I,N_{n_{\mathrm{2}}}^I={\theta }_{\mathrm{WH}}^{\mathrm{up}}, and {\sigma }_{n_{\mathrm{2}}}={\sigma }_{\mathrm{2}}.

The final unconstrained deterministic optimization problem is given by

(19)\min {\displaystyle \underset{n=\mathrm{1}}{\overset{N}{\mathrm{\Sigma }}}{f}_{\mathrm{final}}(X_n)=}\min {\displaystyle \underset{n=\mathrm{1}}{\overset{N}{\mathrm{\Sigma }}}\left[(\mathrm{1}-\beta )f^c(X_n)+\beta f^w(X_n)+\alpha {\displaystyle \underset{m=\mathrm{1}}{\overset{M}{\mathrm{\Sigma }}}\mathrm{phi;}(P(M_{n_m}^I\ge N_{n_m}^I)-{\sigma }_{n_m})}\right]},

where α is a penalty factor which usually has a big numerical value.

BPSO as solution algorithm and implementation of model predictive control

Binary particle swarm optimization (BPSO)

Particle swarm optimization (PSO) is a stochastic optimization methodology which was first proposed by Eberhart and Kennedy [²¹] and Eberhart et al. [²²]. In PSO, each particle indicates a candidate solution in the solution space. The movement and fitness value of particles are determined by the velocity and objective function respectively. The movement of a particle is affected by the best performing particle and the best location discovered so far. The speed and position are given by

(20)v_{i,k}^{t+\mathrm{1}}=w{v}_{i,k}^t+c_{\mathrm{1}}r()(p_{\mathrm{Gb},i}^t-p_{i,k}^t)+c_{\mathrm{2}}r()(p_{\mathrm{Pb},i}^t-p_{i,k}^t),

(21)p_{i,k}^{t+\mathrm{1}}=p_{i,k}^t+v_{i,k}^{t+\mathrm{1}},

where V_k(v_1,k,…,v_n,k) and P_k(p_1,k,…,p_n,k) denote the speed and position of the ith particle, r( ) is a random number in closed set [0,1], w is the inertia weight while c₁ and c₂ denote the weighted coefficients assigned to the global and personal best solutions.

As a modification of basic PSO, the BPSO is introduced to handle the discrete optimization problem. In BPSO, the possibility mapped by the speed of a particle is determined by sigmoid function (22). Based on this possibility, the coordinate of a particle is updated as follows: If rand()>S(v_i^{t+\mathrm{1}}), p_i^{t+\mathrm{1}}=\mathrm{0}; otherwise, p_i^{t+\mathrm{1}}=\mathrm{1} [²³].

(22)S(v_i^{t+\mathrm{1}})={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp (-v_i^{t+\mathrm{1}})}}.

Mathematically, Eq. (16) is a bi-level unconstrained optimization. Its outer layer can be optimized by BPSO. The inner layer optimization is focused on obtaining electricity price interval f_{\mathrm{original}}(X_n,p_n^I) and water temperature interval {\theta }_n^I.

As f_{\mathrm{original}}(X_n,p_n^I) is accumulated linearly and just involved with one interval number p_n^I, f_{\mathrm{original}}(X_n,p_n^I) can be determined easily. With regard to {\theta }_n^I, two extreme schemes, the maximum scheme and minimum scheme, are used to determine the upper and lower bounds of {\theta }_n^I.

(23){\theta }_n^I=\left[\begin{array}{l}(\mathrm{1}-\partial )({\theta }_{e,n}^{\mathrm{down}}+U_{\mathrm{heat},n})+\partial {\theta }_n^{\mathrm{down}}-((\mathrm{1}-\partial )U_{\mathrm{heat},n}+\partial ({\theta }_n^{\mathrm{up}}-{\theta }_{e,n}^{\mathrm{down}}))(d_n^{\mathrm{up}}/M),\\ (\mathrm{1}-\partial )({\theta }_{e,n}^{\mathrm{up}}+U_{\mathrm{heat},n})+\partial {\theta }_n^{\mathrm{up}}-((\mathrm{1}-\partial )U_{\mathrm{heat},n}+\partial ({\theta }_n^{\mathrm{down}}-{\theta }_{e,n}^{\mathrm{up}}))(d_n^{\mathrm{down}}/M)\end{array}\right],

(24){\theta }_n^I=[f({\theta }_n^I,t_n,Q,C,R,d_n^{\mathrm{up}},X_n,{\theta }_{e,n}^{\mathrm{down}}),f({\theta }_n^I,t_n,Q,C,R,d_n^{\mathrm{down}},X_n,{\theta }_{e,n}^{\mathrm{up}})].

The maximum scheme which consists of d_n^{\mathrm{down}}and {\theta }_{e,n}^{\mathrm{up}} could get the upper bound of {\theta }_n^I while the minimum scheme which consists of d_n^{\mathrm{up}} and {\theta }_{e,n}^{\mathrm{down}} could get the lower bound of {\theta }_n^I. The detailed proof process is described as follows:

Proof.

Combining Eqs. (1) and (2), there is

(25)\begin{array}{l}{\theta }_{n+\mathrm{1}}=f({\theta }_n,t_n,Q,C,R,d_n,X_n,{\theta }_{e,n})\\ \iff [({\theta }_{e,n}+X_{n}QR-({\theta }_{e,n}+X_{n}QR-{\theta }_n)\exp [-(t_{n+\mathrm{1}}-t_n)/(RC)])(M-d_n)+{\theta }_{e,n}{d}_n]/M.\end{array}

To simplify this expression, the following equations are set:

(26)U_{\mathrm{heat},n}=X_{n}QR,

(27)\partial =\exp [-(t_{n+\mathrm{1}}-t_n)/(RC)].

After that, there is

(28)\begin{array}{l}{\theta }_{n+\mathrm{1}}=[({\theta }_{e,n}+U_{\mathrm{heat},n}-\partial ({\theta }_{e,n}+U_{\mathrm{heat},n}-{\theta }_n))(M-d_n)+{\theta }_{e,n}{d}_n]/M\\ \iff (\mathrm{1}-\partial )({\theta }_{e,n}+U_{\mathrm{heat},n})+\partial {\theta }_n-((\mathrm{1}-\partial )U_{\mathrm{heat},n}+\partial ({\theta }_n-{\theta }_{e,n}))(d_n/M).\end{array}

Transform the uncertain-but-bounded parameters to intervals, there is

(29)\begin{array}{l}{\theta }_{n+\mathrm{1}}^I=(\mathrm{1}-\partial )({\theta }_{e,n}^I+U_{\mathrm{heat},n})+\partial {\theta }_n^I-((\mathrm{1}-\partial )U_{\mathrm{heat},n}+\partial ({\theta }_n^I-{\theta }_{e,n}^I))(d_n^I/M)\\ \iff [(\mathrm{1}-\partial )({\theta }_{e,n}^{\mathrm{down}}+U_{\mathrm{heat},n})+\partial {\theta }_n^{\mathrm{down}},(\mathrm{1}-\partial )({\theta }_{e,n}^{\mathrm{up}}+U_{\mathrm{heat},n})+\partial {\theta }_n^{\mathrm{up}}]-\\ [(\mathrm{1}-\partial )U_{\mathrm{heat},n}+\partial ({\theta }_n^{\mathrm{down}}-{\theta }_{e,n}^{\mathrm{up}}),(\mathrm{1}-\partial )U_{\mathrm{heat},n}+\partial ({\theta }_n^{\mathrm{up}}-{\theta }_{e,n}^{\mathrm{down}})]\times [d_n^{\mathrm{down}},d_n^{\mathrm{up}}]/M\\ \iff [(\mathrm{1}-\partial )({\theta }_{e,n}^{\mathrm{down}}+U_{\mathrm{heat},n})+\partial {\theta }_n^{\mathrm{down}},(\mathrm{1}-\partial )({\theta }_{e,n}^{\mathrm{up}}+U_{\mathrm{heat,}n})+\partial {\theta }_n^{\mathrm{up}}]-\\ [((\mathrm{1}-\partial )U_{\mathrm{heat},n}+\partial ({\theta }_n^{\mathrm{down}}-{\theta }_{e,n}^{\mathrm{up}}))(d_n^{\mathrm{down}}/M),((\mathrm{1}-\partial )U_{\mathrm{heat},n}+\partial ({\theta }_n^{\mathrm{up}}-{\theta }_{e,n}^{\mathrm{down}})(d_n^{\mathrm{up}}/M))]\end{array}

(30)\begin{array}{l}\iff \left[\begin{array}{l}(\mathrm{1}-\partial )({\theta }_{e,n}^{\mathrm{down}}+U_{\mathrm{heat},n})+\partial {\theta }_n^{\mathrm{down}}-((\mathrm{1}-\partial )U_{\mathrm{heat},n}+\partial ({\theta }_n^{\mathrm{up}}-{\theta }_{e,n}^{\mathrm{down}}))(d_n^{\mathrm{up}}/M),\\ (\mathrm{1}-\partial )({\theta }_{e,n}^{\mathrm{up}}+U_{\mathrm{heat},n})+\partial {\theta }_n^{\mathrm{up}}-((\mathrm{1}-\partial )U_{\mathrm{heat,}n}+\partial ({\theta }_n^{\mathrm{down}}-{\theta }_{e,n}^{\mathrm{up}}))(d_n^{\mathrm{down}}/M)\end{array}\right]\\ \iff [f({\theta }_n^I,t_n,Q,C,R,d_n^{\mathrm{up}},X_n,{\theta }_{e,n}^{\mathrm{down}}),f({\theta }_n^I,t_n,Q,C,R,d_n^{\mathrm{down}},X_n,{\theta }_{e,n}^{\mathrm{up}})].\end{array}

Obviously, when the minimum scheme (d_n^{\mathrm{up}},{\theta }_{e,n}^{\mathrm{down}}) and maximum scheme (d_n^{\mathrm{down}},{\theta }_{e,n}^{\mathrm{up}}) are taken, the bounds of uncertain hot water temperature could be obtained.

Model predictive control

As a robust optimization method, the interval number optimization has many prominent advantages for dealing with the uncertainty problems. For example, the analysis and calculation are quite concise and it does not need any additional auxiliary variables. However, if only the day-ahead optimization is conducted and executed, the interval number optimization cannot handle the unexpected incidents which may occur the next day. For the EWH, the hot water demand is decided by the consumers’ random behavior, so the water demand will have some unexpected extra increase or decrease beyond the prediction interval. Under these circumstances, the model predictive control can be introduced to address the problem.

The receding horizon optimization (RHO) is the core of model predictive control. Earlier in this paper, the day-ahead EWH load scheduling model has been established to make load schedules throughout the whole time horizon. In contrast, the RHO can be seen as a type of local optimization. It just involves a subset of the whole time horizon, which is called as optimization time domain in this paper. Specifically, the local optimization involves the current time step, with accurate information, and some time steps after this step, with forecast information. The RHO could renovate the optimization parameters according to the latest state of the system. Hence the RHO method could obtain better results compared with the single day-ahead optimization.

The time frame of RHO is depicted in Fig. 2 [¹⁵,¹⁷]. The EWH is controlled at each time step in a rolling way. Specifically, at time step n, the EWH is controlled according to the schedule made at time step n−1, and the schedule for future T time steps is made.

Implementation process

Figure 3 describes the specific flowchart of EWH load scheduling using MPC-based interval number optimization, with BPSO as the optimization algorithm. According to time frame of RHO in Fig. 2, the time step n is set and begin to work out the schedule for future T steps beyond this step. First of all, n is equal to 1. The BPSO methodology is used to get the EWH schedule. In Eq. (19), the values of f^c(X_n), f^w(X_n), and

\phi (P(M_{n_m}^I\ge N_{n_m}^I)-{\sigma }_{n_m})

are related to both decision variables X_n and uncertain parameters, which lead to the fact that it cannot be optimized at one shot directly. Therefore, the problem needs to be solved by the iterative interaction between two layers. The outer layer, which takes BPSO as the optimization algorithm, is used to update decision variables X_n. Then X_n is sent to the inner layer, which takes linear interval number programming as the optimization algorithm to compute the corresponding intervals of the electricity bill and hot water temperature. Based on the decision variable X_n and corresponding intervals of hot water temperature, the weighted electricity bill and penalty function can be calculated to get the fitness value of design vector X_n which is utilized to update the decision variables X_n+1 of the next iteration. After the calculation of the outer layer and inner layer, the EWH scheduling results are obtained, which means that when this RHO process ends, a new process will be restarted by resetting the time step n + 1 until the entire time range is completed.

Case studies

Basic parameters

The forecasted RTP interval and actual RTP demonstrated in Fig. 4 is extended from the basic electricity price in Ref. [²⁴]. Similarly, the curves of hot water demand and ambient temperature are shown in Figs. 5 and 6, respectively. The gray lines denote the actual values and the black lines denote the predicted values. The upper and lower bounds are the 110% and 90% of day-ahead predicted value respectively [¹¹].

Moreover, the basic parameters of the EWH can be referred to in 2012 ASHRAE Handbook [²⁵], and are listed in Table 1.

The scheduling horizon begins from 0 a.m. until 12 p.m. and the step size is 15 min, hence N is equal to 96. The power rating of water heater P_EWH is 4.5 kW [²⁶]. The parameters in Eq. (16) are δ=2.2, φ =0.15. The value of weight coefficient β is deemed as a constant value in the analysis. The reason for this is that the main focus of this work is to study the impact of different tolerance degrees of constraint violation and temperature preferences on hot water temperature profile and electricity bill. Taking the value of weight coefficient in Refs. [¹⁸,²⁷] into account, the weight coefficient β is set as 0.5.

Day-ahead optimization

Case 1: zero tolerance degree

Under the disturbances of uncertainty parameters, if the EWH load scheduling is just treated as a deterministic optimization problem, violation may occur frequently in temperature constraints. For example, it is assumed that the low bound of hot water demand and upper bound of ambient temperature are applied in the EWH load scheduling. Due to inevitable forecast errors, the actual values of the hot water demand and ambient temperature will fluctuate in the intervals as shown in Figs. 5 and 6. As a result, the actual interval of water temperature is shown in Fig. 7, in which the black area represents the possible water temperature interval under the fluctuations of uncertain parameters.

In Fig. 7, the minimum of actual water temperature is 55.85°C, but the comfort zone is 59.5°C, 70.5°C. Obviously, such a great error is unacceptable for customers.

To avoid this, the scheme of zero tolerance degree (Scheme 1) is proposed to control the actual water temperature. In this scheme, both σ₁ and σ₂ are assigned to be zero to ensure that there are no violations in the results. At this moment, the penalty function will be very large when the scheduling plan causes violations, and the plan will be forced to adjust. When Scheme 1 is implemented in load scheduling, the water temperature interval is exhibited in Fig. 8.

Compared with Fig. 7, it is obviously noticed that Scheme 1 is robust to variations of the hot water demand and ambient temperature. That is, the hot water temperature always subjects to comfort constraints even under the disturbance of uncertain hot water demand and ambient temperature.

Case 2: different tolerance degrees and temperature preferences

Scheme 1 is robust to all possible situations, however it is also expensive. For the majority of consumers, they prefer to make a trade-off between the electricity bill and comfort. Therefore, several schemes with different tolerance degrees and temperature preferences are taken into account in this part. There are two assessment factors included in a scheme. One is the tolerance degree for comfort violation; the other one is the temperature preference which describes the preference for a higher or lower water temperature.

Conservative schemes discussed here includes a low tolerance degree with no temperature preference (Scheme 2: σ₁ = σ₂ = 0.2), a low tolerance degree with a higher temperature preference (Scheme 3:σ₁ = 0.1, σ₂ = 0.3), a low tolerance degree with a lower temperature preference (Scheme 4:σ₁ = 0.3, σ₂ = 0.1), a high tolerance degree with no temperature preference (Scheme 5:σ₁ = σ₂ = 0.4), a high tolerance degree with a higher temperature preference (Scheme 6:σ₁ = 0.3, σ₂ = 0.5), and a high tolerance degree with a lower temperature preference (Scheme 7:σ₁ = 0.5, σ₂ = 0.3).

As discussed in Section 3.2, the uncertain electricity bill can be expressed as p=[p^c, p^w], in which p^c is the middle point of electricity bill band and p^w is the radius of electricity bill band. Moreover, the weighted mean value of electricity bill band is described as \bar{p}=\beta \cdot p^c/\delta +(\mathrm{1}–\beta )\cdot p^w/\phi. The middle point, radius and weighted mean value are expressed in Table 2.

When it comes to electricity bills, Scheme 1 is the most expensive strategy as it shows an absolute demand for comfort. And with the increase of σ, the costs will reduce (Scheme 3 is cheaper than Scheme 1 but more expensive than Scheme 5). At the same time, the preference for a higher water temperature is more expensive than the preference for a lower water temperature (Scheme 2 is cheaper than Scheme 3 but more expensive than Scheme 4). The interval curves of hot water temperature from Scheme 2 to Scheme 7 are displayed in Fig. 9.

To quantize constraint violations in each scheme, the violation degree is described in nth step with \Delta {\theta }_{\mathrm{degree},n} as

(31)\Delta {\theta }_{\mathrm{degree},n}={\displaystyle \frac{{\theta }_{\mathrm{violation},n}-{\theta }_{\mathrm{set}}}{{\theta }_n^{\mathrm{up}}-{\theta }_n^{\mathrm{down}}}}\times \mathrm{100}\%,

(32)

(33)

Correspondingly, the violation degrees of different schemes are illustrated in Fig. 10.

As shown in Figs. 9 and 10, the violation degrees of Scheme 2–4, are less than those of Scheme 5–7, as they possess lower tolerance degrees for constraint violation; in Schemes 2 and 5, constraint violations on the upper limit and the lower limit share similar numbers, as no temperature preferences are required in the above two schemes; in Schemes 4 and 6, more constraint violations occur on the upper limit, as they both show a preferences for a higher temperature; in Schemes 5 and 7, the situation is contrary. More constraint violations occur on the lower limit, as they both show a preference for a lower temperature.

The above analysis indicates that the proposed method is able to flexibly control the water temperature according to different consumers’ demands.

Model prediction optimization

Case 3: model prediction optimization

In actual operation, if an unexpected situation occurs, the day-ahead plan cannot be corrected in time. Taking the scheduling plan in Section 5.2.1 as an example, it is assumed that the water consumption is greatly reduced between period 24 and period 40, and the actual temperature curve is as shown in Fig. 11. There are many violations of comfort constraints, which means that the interval number optimization could not solve the sudden incident.

The simulation results only considering RHO is shown in Fig. 12, whose optimization time domain T is 16. The single receding horizon optimization strategy cannot ensure that there is no comfort violation under the environment setting of this paper.

The simulation results after MPC-based interval number optimization of which the optimization time domain T is 16 is shown in Fig. 13. The hot water temperature basically satisfies the comfort constraints for all scheduling periods, which means the schedules have a strong robustness. In Fig. 14, the simulation process diagram in the receding horizon optimization is revealed which explains the principle of RHO clearly. The RHO could response to the sudden incident of water demand efficiently and avoid the actual water temperature of EWH being out of control.

Case 4: optimization time domain

RHO is the local optimization, which means the selection of the optimization time domain is vitally important. To study the influence of the optimization time domain T reasonably and accurately, the optimization time domain T can be set from 2 to 48, corresponding to 0.5–12 h in reality. The simulation results of different Ts in Scheme 1 are presented in Fig. 15.

It can be seen from Fig. 15 that when T is small, the electricity cost is high due to not including much future information. As the time horizon T increases, the electricity cost decreases. When T is set as 16, the electricity cost is minimum. Then, when the time horizon T is set from 16 to 48, the electricity cost increases as T increases. The reason for this is that the uncertainties of electricity price and ambient temperature forecasts increases. The uncertainties of electricity price and ambient temperature influence the choice for the time horizon T. A bigger T could get more decision information but could also bring more uncertainties. Therefore, when the uncertainty of the predicted data increases, the optimum value of T decreases correspondingly, and when the uncertainty decreases, the change is reversed. The best choice of T in different situations can be determined by optimization to obtain the minimum electricity cost.

Conclusions

This paper introduces the interval number optimization and MPC into EWH load scheduling. Interval number optimization is robust to the uncertainties and flexible to different consumers’ demands in electricity bill and comfort. Several schemes with different tolerance degrees and temperature preferences are utilized to quantifiably describe diversified requirements for the comfort zone. MPC can handle the sudden incident and compensate for deficiencies of interval number optimization. The simulation results demonstrate that as the tolerance degree increases, the electricity bill decreases. Besides, the preference for low water temperature is cheaper than that for high water temperature. Moreover, both the midpoint and the radius of electricity bill interval would reduce with the increase of tolerance degree. Furthermore, MPC could handle the uncertainty factors and possesses a stronger robustness to customer behavior.

In summary, the methods proposed have some advantages over traditional optimization methods. However, there are still some deficiencies and limitations in the above discussion. More appliances should be included in the analysis, such as interruptible loads, non-interruptible loads, and even distributed energy resource. Moreover, the receding horizon optimization is also an influencing factor which remains to be studied.

References

[1] Ruiz N, Claessens B, Jimeno J, López J A, Six D. Residential load forecasting under a demand response program based on economic incentives. International Transactions on Electrical Energy Systems, 2015, 25(8): 1436–1451

[2] Sisodiya S, Shejul K, Kumbhar G B. Scheduling of demand-side resources for a building energy management system. International Transactions on Electrical Energy Systems, 2017, 27(9): e2369

[3] Rastegar M, Fotuhi-Firuzabad M, Aminifar F. Load commitment in a smart home. Applied Energy, 2012, 96: 45–54

[4] Du P, Lu N. Appliance commitment for household load scheduling. IEEE Transactions on Smart Grid, 2011, 2(2): 411–419

[5] Wang C, Zhou Y, Wang J, Peng P. A novel traversal-and-pruning algorithm for household load scheduling. Applied Energy, 2013, 102: 1430–1438

[6] Beaudin M, Zareipour H. Home energy management systems: a review of modelling and complexity. Renewable & Sustainable Energy Reviews, 2015, 45: 318–335

[7] Vivekananthan C, Mishra Y, Li F. Real-time price based home energy management scheduler. IEEE Transactions on Power Systems, 2015, 30(4): 2149–2159

[8] Sucar I, Enrique L. Decision Theory Models for Applications in Artificial Intelligence: Concepts and Solutions. Hershey: IGI Global Press, 2012: 97–143

[9] Wu Z, Zhang X P, Brandt J, Zhou S Y, Li J N. Three control approaches for optimized energy flow with home energy management system. IEEE Power and Energy Technology Systems Journal, 2015, 2(1): 21–31

[10] Shahgoshtasbi D, Jamshidi M M. A new intelligent neuro-fuzzy paradigm for energy-efficient homes. IEEE Systems Journal, 2014, 8(2): 664–673

[11] Wang J, Li Y, Zhou Y. Interval number optimization for household load scheduling with uncertainty. Energy and Building, 2016, 130: 613–624

[12] Wang D, Jia H, Wang C, Lu N, Fan M, Zhou Y, Qi Y. Voltage stability enhancement using thermostatically controlled appliances as a comfort-constrained virtual generator. International Transactions on Electrical Energy Systems, 2015, 25(12): 3509–3522

[13] Kara E C, Bergés M, Hug G. Modeling thermostatically controlled loads to engage households in the smart grid: lessons learned from residential refrigeration units. In: International Conference on Computing in Civil and Building Engineering, Orlando, USA, 2014: 2032–2039

[14] Zhang J, Dominguez-Garcia A. Evaluation of demand response resource aggregation system capacity under uncertainty. IEEE Transactions on Smart Grid, 2017, 9(5): 4577–4586

[15] Sharma I, Dong J, Malikopoulos A A, Street M, Ostrowski J, Kuruganti T, Jackson R. A modeling framework for optimal energy management of a residential building. Energy and Building, 2016, 130: 55–63

[16] Rahmani-Andebili M. Scheduling deferrable appliances and energy resources of a smart home applying multi-time scale stochastic model predictive control. Sustainable Cities & Society, 2017, 32: 338–347

[17] Maasoumy M, Sangiovanni-Vincentelli A. Optimal control of building HVAC systems in the presence of imperfect predictions. In: Dynamics and Systems and Control Conference, Fort Lauderdale, USA, 2012: 257–266

[18] Sengupta A, Pal T K. On comparing interval numbers. European Journal of Operational Research, 2000, 127(1): 28–43

[19] Andrew A M.Applied interval analysis: with examples in parameter and state estimation, robust control and robotics. Kybernetes, 2001, 31(5): 117–123

[20] Jiang C, Han X, Liu G R, Liu G P. A nonlinear interval number programming method for uncertain optimization problems. European Journal of Operational Research, 2008, 188(1): 1–13

[21] Eberhart R, Kennedy J. A new optimizer using particle swarm theory. In: IEEE Proceedings of the 6th International Symposium on Micro Machine and Human Science, New York, USA, 2002: 39–43

[22] Eberhart R, Simpson P, Dobbins R. Computational Intelligence PC Tools. Salt Lake City: Academic Press Professional, Inc. 1996

[23] Kennedy J, Eberhart R C. A discrete binary version of the particle swarm algorithm. In: 1997 IEEE International Conference on Systems, Man, and Cybernetics. Computational Cybernetics and Simulation, Orlando, USA, 1997, 5: 4104–4108

[24] Alrumayh O, Bhattacharya K. Model predictive control-based home energy management system in smart grid. In: IEEE Electrical Power and Energy Conference, Piscataway, USA, 2016: 152–157

[25] American Society of Heating, Refrigerating and Air-Conditioning Engineers, Inc. 2012 ASHRAE Handbook–Heating, Ventilating, and Air-Conditioning Systems and Equipment. 2013, available at app.knovel.com website

[26] Kondoh J, Lu N, Hammerstrom D J. An evaluation of the water heater load potential for providing regulation service. IEEE Transactions on Power Systems, 2011, 26(3): 1309–1316

[27] Jiang C, Han X, Liu G P. A sequential nonlinear interval number programming method for uncertain structures. Computer Methods in Applied Mechanics and Engineering, 2008, 197(49-50): 4250–4265

RIGHTS & PERMISSIONS

Higher Education Press