Electricity Consumption Prediction of Solid Electric Thermal Storage with a Cyber–Physical Approach

: This paper proposes a cyber–physical approach to enhance the prediction accuracy of electricity consumption of solid electric thermal storage (SETS) system, which integrates a physical model and a data-based cyber model. In the cyber–physical model, the prediction error of the physical model is used as an input of the cyber model to further calibrate the prediction error. Firstly, customers’ behavior characteristics are extracted by the integration of K-means and one-versus-one support vector machine. Secondly, based on the behavior characteristics and ambient temperature, the physical model is developed to predict daily electricity consumption. Finally, the error levels of physical model are classiﬁed, together with the temperature and prediction values of the physical model, are selected as the inputs of the cyber model using the back propagation (BP) neural network to calibrate the results of the physical model. The effectiveness of the proposed cyber–physical model (CPM) is veriﬁed by a 1 MW SETS system. The simulation results show that, compared with the physical model (PM) and cyber model (CM), the maximum relative errors (MRE) with the CPM are reduced to 25.4% and 4.8%, respectively.


Introduction
Thermal energy storage is considered as one of the advanced energy technologies [1]. Electric energy can be stored in the form of heat during off-peak demand periods and used for heating of rooms during peak demand periods. The improvement of thermal storage is useful to reasonably arrange the electricity consumption of thermal storage loads and promote the thermal storage peak shaving incentive mechanism realization. Therefore, solid electric thermal storage (SETS) has become one of the most promising solutions as a flexible demand response (DR) in demand-side management (DSM) [2]. Consequently, SETS prediction is an important precondition for peak shaving in DSM. Due to the electric-thermal time shift characteristics of SETS, the electricity prediction can provide multiple options for peak shaving and dispatching of the power system. However, there are few studies on the SETS prediction, and the prediction accuracy needs to be improved. Therefore, if the SETS prediction is applied to the actual operation of the power system, it is necessary to enhance the prediction accuracy of SETS based on advanced algorithms.
There are few studies on the electricity prediction of those devices. A physical model (PM) of SETS is integrated into an energy management system for isolated microgrids [3]. SETS is used to accommodate wind power to supply heat load in isolated microgrids [4]. A PM for residential

•
To the best of the authors' knowledge, this is the first work to use the cyber-physical approach to predict the SETS' load change. The physical and cyber components of SETS are integrated. • Using the existing knowledge of thermodynamics, a SETS PM is developed by considering the customers' behavior characteristics.

•
The load data of 1MW SETS is used to validate the CPM, and the results show that, compared with the PM and the CM, the maximum relative errors (MRE) with the CPM are reduced to 25.4% and 4.8%, respectively.
The rest of this paper is organized as follows: The PM of SETS, including its structure, principle, formula derivation, the customers' behavior characteristics of SETS, and influencing factors are presented in Section 2. The cyber-physical approach combining the PM and CM of SETS is proposed in Section 3. Simulations are conducted to validate the effectiveness of the cyber-physical approach in Section 4. The conclusion is drawn in Section 5.

PM of SETS
In the PM of SETS, the structure of the temperature-based energy flow of the SETS is proposed in Section 2.2. (Principle). The heat conduction of SETS is described in Section 2.3. (Thermal Energy Storage). The heat transfer formulas of the heat exchanger are given in Section 2.4.1 (Heat Transfer). The Formula (14) is developed by the authors in Section 2.4.2 (Customers Heating). The behavior extraction of SETS is firstly considered in the model. Figure 1 shows a schematic diagram of SETS. It contains electric heating wires, magnesia bricks, the temperature sensor of bricks, two layers of thermal isolation (including perlite and ceramic fiber) wrapped in steel plates, heat exchanger, and fan. When the SETS operates, the resistive heating wires generate thermal energy, which will be stored in the bricks for discharging. The cold air is blown into the interior by the fan, and transferred into hot air by the heating wires for changing the entrance temperature of the heat exchanger. The cold water is circulated into the heat exchanger through the return water pipe and heated by the hot air to improve the customers' room interior temperature. SETS becomes a convenient DSM tool for utility companies.

Principle
The basic assumption of the proposed PM is that the temperature of the bricks is spatially uniform in any transient process. This assumption implies that the temperature gradients within the bricks are negligible. The proposed model is based on the following assumptions: all the electric energy consumed by the electric heating wires is stored in the bricks in the form of thermal energy. The thermal energy stored in the SETS is mainly released by thermal radiation and convection. The isolation layer can prevent the thermal energy of hot air in the SETS from flowing outside. In heat transfer of the heat exchanger, the heat loss between the heat exchanger and the ambient environment, the pipe thermal resistance and fouling effect are ignored. In Figure 2, the PM of SETS in this paper is based on the basic principle of heat transfer instead of entransy dissipation-based thermal resistance theory presented in [22]. The initial and final temperatures of the bricks are assumed to be T ini and T fin , respectively. The internal and external temperatures of the isolation layer are T int and T ext , respectively. T ini is equal to T int before the electric heating wires are energized. The T ini is increased as the heating wires are energized. In the heating process, thermal energy always transfers heat to the surroundings through the isolation layer, so T int can be calculated by the average temperature of the bricks. According to a large number of survey and analysis of SETS operation, the T ext approximation is 50 • C. The boundary and initial conditions of the PM are described in Equation (1).
The inlet and outlet temperatures of the hot air through the heat exchanger are T a,in and T a,o , respectively. The inlet and outlet temperatures of the circulating water through the heat exchanger are T w,in and T w,o , respectively. T cri a,room is the criterion minimum temperature of customers' room. The room temperature T a,room is influenced by the environmental conditions (e.g., ambient temperature T a,amb . In Section 2.7. (Influencing Factors) explains why ambient temperature is selected). In order to keep the room temperature close to T cri a,room , the supply and return water temperatures T w,in and T w,o can be increased correspondingly. In the inlet of the heat exchanger, a high-temperature air is required to be input, so the stored thermal energy final temperature T fin needs to be higher (about 700 • C). T fin can be deduced by the change of ambient temperature T a,amb , and the average power consumption of the SETS can be obtained by the T fin . According to engineering experience, the variation ranges of T ini and T a,o are very little that they are assumed to be constant. T w,o is affected by the ambient temperature T a,amb , so the input and control variables of the heat exchanger T w,o and T w,in are deduced based on Equation (12).
Without considering the temperature gradients within the bricks, the initial temperature T ini and final temperature T fin changes in the SETS are used to predict the thermal energy consumption. The proposed SETS PM prediction of the consumed thermal energy is equal to the electricity by the ambient temperature change.

Thermal Energy Storage
According to the above principle, the electric energy consumption E pro,t is equal to the sum of thermal energy E sto,t stored in the bricks and the thermal energy loss of the SETS E los,t as shown in Equation (2).
The E los,t is not considered in the modeling. So the Equation (2) can be a simplify as E pro,t = E sto,t . The thermal energy E sto,t is transferred to the hot air by thermal convection φ conv , thermal radiation φ rad , and thermal conduction φ cond . The three heat transfer formulas are introduced as follows:

Thermal Convection
Based on the initial temperature T ini and final temperature T fin of the thermal storage, the thermal convection from bricks into the isolation layer is given by Equation (3). The thermal convection equation is conveniently expressed by Newton's law of cooling [23].
where h conv is the surface heat transfer coefficient of the bricks, A cha is the surface area of the round holes in the bricks as shown in Figure 3. Many bricks are superimposed to form the core frame of the thermal storage. Two heat bricks are stacked together to form a circular hole in the center of the bricks as shown in the zoom figure. Electric heating wires are placed in the middle of the hole. L, W, H, and D are the length, width, height, and diameter of the brick, respectively. A cha is the area of hot airflow channel which is stacked by two bricks together, A cha = L × D × N, where N is the number of bricks in a row. Based on the Nusselt-Number N u , the heat transfer coefficient h conv can be calculated by the Zukauskas formula [24]. λ is the thermal conductivity. In heat convection, multiple round holes can be approximately regarded as the tube bundle. The average heat transfer performance of the tube bundle is related to the flow Reynolds-number (R e ) [24].
where u a is the hot air flow rate, D is bricks round hole diameter, and ν T ave is the aerodynamic viscosity at the temperature of T ave = T fin /2. According to the range of R e , the Nusselt-number (N u ) can be calculated by Zukauskas formula in Equation (5).
where P r a is the Prandtl-number (PN) of the fluid average temperature, PN reflects the contrast between momentum diffusion and thermal diffusion in fluid. P r c is a constant. P r c = P r a /P r bri 0.25 , P r bri is the PN of the bricks average surface temperature. According to the temperature T int and T fin , the P r a and P r bri can be queried. When 10 2 < R e ≤ 20 6 , the corresponding N u can be selected. In most cases, the calculation results satisfy the condition that 10 3 < R e ≤ 20 5 . This paper uses N u = 0.27R 0.63 e P 0.36 r g P r c as an example. Substituting Equation (5) into Equation (3), the thermal energy φ conv is obtained in Equation (6).

Thermal Radiation
The thermal energy φ rad released from the bricks to the flow channel hot air, is expressed by Equation (7), which is based on the Stefan-Boltzmann law [23].
where ε is the emissivity of the bricks' surface, σ is the blackbody radiation constant (Stefan-Boltzmann constant is 5.67 × 10 −8 W/m 2 K 4 ), and A bri is the heat transfer surface area of the bricks.

Thermal Conduction
The isolation layer consists of fiberglass felt and perlite coated. The thermal conduction φ cond between the internal hot air and the external ambient environment is given by the Fourier's law of heat conduction in (8) [23].
where δ f δ p and λ f λ p are the respective thickness, and thermal conductivity of fiberglass and perlite, and A iso is the isolation layer area. The prediction of average power consumptionQ sto can be obtained by Equations (6)-(8) as shown in Equation (9).
where T fin is unknown which can be obtained in Section 2.4 (Thermal Bricks Energy Release).

Thermal Bricks Energy Release
Thermal bricks energy release includes heat transfer and customers' heating.

Heat Transfer
When the SETS releases thermal energy, the inlet hot air temperature of the heat exchanger is regulated by the circulating wind speed of the fan. In the heat exchanger, according to the conservation of energy, the thermal energy released Q a by the inlet hot air is equal to the thermal energy absorbed Q w by the outlet circulating water. Therefore, the relationship between the inlet hot air temperature and the outlet circulating water temperature is given in (10) [25].
where,ṁ a andṁ w are the flow rate of hot air and water through the heat exchanger, respectively. c a and c w are the specific heat of hot air and water through the heat exchanger, respectively. Since Q a = Q w , T a,i can be deduced as shown in Equation (11). Ignoring the heat loss from the hot air at the outlet of the SETS to the heat exchanger inlet, and assuming that T a,i ≈ T fin , then, taking this into Equation (9), the average power consumptionQ sto can be obtained.
where, T w,o − T w,in is determined by the Section 2.4.2 (Customers Heating).

Customers Heating
The thermal energy storage is intended to exchange heat with the circulating water to maintain the customers' room temperature greater than or equal to the designed heating temperature T cri a,room . The temperatures of supply and return water are obtained by the room and the ambient temperature in Equation (12).
T cri a,room −T a,amb T cri a,room −T cri a,amb T cri a,room −T a,amb T cri a,room −T cri a,amb 1 1+η where T cri w,o and T cri w,in are the respective temperatures of water supply and return designed for SETS. η is the heat transfer coefficient of the radiator. According to the standard of heating in the north of China, the minimum value of T cri a,room is 18 • C. During the winter heating period, the average ambient temperature of Shenyang is −16.9 • C.
The overshoot of water temperature and the frequent start of the fan are easily caused by the use of T w,o and T w,in owing to the inertia of temperature change. Therefore, in this paper, the highest supply water temperature T w,o,max and lowest return water temperature T w,i,min during the heating period are designed. Taking Shenyang city as an example, the minimum and maximum temperatures of T a,amb are taken into Equation (12) to obtain the temperature curves of supply and return water as shown in Figure 4. The minimum temperature of return water T w,i,min is selected as the benchmark so that the change of temperature T w,o − T w,in = γ can be deduced as shown in Equation (13). Substituting γ into Equation (11), T a,in is obtained as shown in Equation (14).
Based on the above analysis, the thermal energy store in bricks at time t can be obtained as shown in Equation (15).
where,Q sto,t is the average power consumption of SETS at time t. E pro,t is the electricity consumption by the heating wires, and E sto,t is the stored heat energy in the bricks. The working time t is related to the behavior characteristics, and the superscript i is the time index.

Customers' Behavior Characteristics Extraction
The working time t of the SETS is determined by the different electricity consumption behavior characteristics of customers. For example, two types of behavior characteristics of SETS customer within one month is shown in Figure 5. One type starts from the previous night 22:00 and stops heating at 1:00. Then, SETS starts again at 3:00 and ends the reservoir at 5:00. The other type starts from the previous night 22:00 and then continues to work until the end at 5:00. The consumption time t of one type is t = 3.5h on the left curve, and the other type is t = 7h on the right curve.
Since the load curves of the same electrical behavior have regular working time, it is necessary to cluster the same electricity consumption behavior characteristics with the heat load data. Then, the customers' behavior types are used to determine the working time t of SETS. Also, the E pro,t daily electricity consumption of SETS is predicted by Equation (15). The method of K-means is used to identify customers' behavior types. After all SETS data clustering, the customers' behavior types are obtained. The load data with the same type has a similar working time, the same type is identified using the support vector machine that has been widely used to learn target class [26]. When the target types are larger than binary, the one-versus-one is one of the effective methods to find the discriminative binary classifiers to solve a multiclass problem [27]. Thus, the one-versus-one method is adopted in behavior extraction. The overall process of the behavior characteristics extraction is shown in Algorithm 1.

Algorithm 1
The behavior characteristics extraction algorithm process.
1: Input SETS' load data. 2: Label the data of the same customer from 1, 2, 3 to n; 3: Using unsupervised clustering method (K-Means) to identify target data types c 1 , c 2 ; 4: Statistics time labels y 1 , y 2 , · · · , y m of the same type c j ; 5: i f (m > 2) then 6: Using the one-versus-one method, divide the time labels into binary pairs as the support vector machine output time labels; 7: else 8: Using the support vector machine directly; 9: end i f 10: Output the working time t.

Summary of the PM Prediction
Based on the above analysis, the proposed PM of SETS can be used to predict the daily electricity consumption. The overall PM flowchart is shown in Figure 6. According to the initial values, the range of Reynolds-Number (R e ) is estimated. If R e > 20 6 or R e ≤ 10 2 , return the initialization parameters, otherwise, input the ambient temperature T a,amb to obtain temperatures of T w,o and T w,in . Select the maximum value of supply water T w,o,max and the minimum value of return water T w,in,min . Then, obtain the inlet temperature T a,in . N u is selected based on the range of R e , so input the N u and T a,in ≈ T fin to predict the average power consumptionQ sto . TheQ sto,t is integrated over working hours t. According to the stored thermal energyQ sto,t in the bricks, the predicted value of SETS heating wires electricity consumption E pro,t can be obtained.

Influencing Factors
Since influencing factors (including humidity, wind speed, and temperature) influence the prediction of the CPM, it is necessary to select the strong correlation factors as the input of the cyber and PMs. Thus, the Pearson correlation analysis in Equation (16) is adopted to analyze which one is strongly correlated. The influencing factors and the SETS load are represented as a and b, respectively.
The statistics on the correlation between the influencing factors and SETS load within 24 h is as shown in Figure 7. The histogram explains which influencing factors are strongly correlated with the SETS load. The influencing factors (including humidity, wind speed, and temperature) are negatively correlated with the load, and humidity is positively correlated with the SETS load only in the range of 8:00 to 10:00. The ranges of correlation coefficient of humidity (blue), wind speed (green) and temperature (yellow) are [−0.24, 0.21], [−0.47, 0.13], and [−0.82, −0.66], respectively. It can be seen that the correlation between the temperature and SETS load is strongly negative, while humidity and wind speed are weakly correlated with SETS load. Thus, the temperature is selected as the input of the cyber and PMs.

Cyber-Physical Approach
The parameters of SETS and ambient temperature are the inputs of the SETS PM. The PM prediction result is shown in Figure 8. It can be observed that January has the best fitting with the least error. The test result with the PM still has an irregular error as shown in Table 1. The root mean square error (RMSE), mean absolute error (MAE), and mean absolute percentage error (MAPE) are used to measure the prediction performance of the PM with a deviation between the predicted value and the real value.
where m is the number of data. y i is the real value of SETS load,ŷ i is the prediction value of SETS PM.  The PM predicts error according to the result of RMSE, MAE, and MAPE, and obtains a consistent evaluation. Selecting the month of January which has the smallest RMSE in five months, the predicted result RMSE error is 167.36 kWh. The electricity consumption of real value is 5616 kWh, the prediction accuracy is approximately 29.8%. Therefore, the prediction accuracy of the PM still needs to be further improved.
To improve the prediction accuracy, other prediction methods are considered. The CM is widely used in load prediction because it can accurately model the past. Hence, the application of the CM is considered to predict the electricity consumption of SETS. However, the electric behavior of SETS is different from the other conventional electric loads, it is difficult to predict a sudden increase and decrease in heat load, that is, the predictability of the CM is not good enough. On the contrary, the PM is sensitive to a sudden increase and decrease in the heat load. Therefore, the CM and PM can be combined to improve the prediction accuracy with mutual verification. CPS is widely used in many fields, for example, measurement recovery in false data injection attacks [28], smart home cyberattack detection [29], energy theft detection in multi-tenant data centers [30], and inter-well connectivity estimation in petroleum [31], etc. Hence, the cyber-physical approach is developed to predict electricity consumption of SETS. The CM established by the load data of SETS is used to calibrate the output error of the PM.
The cyber components (including influencing factors) and the physical components (including behavior characteristics and average power consumption) of SETS are integrated using the cyber-physical approach. The cyber-physical approach is used to predict the electricity consumption of SETS. The Pearson correlation analysis formula is used to solve the correlation coefficients ρ a,b between power consumption of SETS and influencing factors. The ambient temperature has the greatest impact on the power consumption of SETS, which is considered as an input of the PM flow.
Firstly, the PM of SETS is utilized to predict the daily average power consumption based on ambient temperature. Secondly, the integration of K-means and one-versus-one support vector machine is applied to extract the behavior characteristics of SETS to obtain the work time t. Then the electricity consumption of SETS is obtained by integrating average power consumption. The performance of the prediction is tested by the error calculation Equation (17). The error levels are divided according to the size of the errors.
Since RMSE is a good reflection of the precision of the measurement, it is used to design the error levels of the PM. According to the result of RMSE, which is represented by s , the subscript indicates s = 1, 2, 3, 4, 5 which represents the months from November 2017 to March 2018, respectively. The number 1 to 5 are 372.05, 310.29, 167.36, 227.1, 389.45, respectively. Therefore, according to the order of s ( 5 > 1 > 2 > 4 > 3 ), the error is divided into 5 levels, and then the data is input in the CM which is used in BP neural network.
Subsequently, the output of PM is calibrated by the error of CM output to obtain the electricity consumption of the SETS. The specific method of the CM is introduced as follows.
The error levels, ambient temperature, and PM prediction values are taken as input vectors x n of BP neural network. Output vectors o n is the error as shown in Figure 9. The number of hidden neurons is 13, the sigmoid function is selected as the activation function for the hidden layer and the output layer, and the gradient descent method is used to update the parameters. In the CM, the influencing factors and the power consumption of SETS are the input vectors x n of BP neural network, and the electricity is the output vectors o n . The number of hidden neurons is set as 25. During training, the CM of the BP neural network can automatically extract the "reasonable rules" between the input layer and output layer data and memorize the learning content in the weight of the network adaptively. The CM can predict the load deviation after training even if there is unknown noise-pollution. BP neural network is used to train the CM. The CM prediction result plus the PM predicted value is the electricity consumption of SETS.

Validation
In order to verify the effectiveness of the proposed CPM method, its prediction result is compared with PM, CM and real value. MRE is used to evaluate the performance of the prediction model.
The proposed method to predict the electricity consumption of SETS is verified using real operation data, the data is obtained from the power grid valley thermal storage system monitoring platform of State Grid Company. This paper uses the data of one customer from the monitoring platform. The SETS is located in the third printing plant of Heping District, Shenyang City, Liaoning Province. The scene device of SETS is shown in Figure 10. The nominal power of the SETS is 1 MW, and the maximum thermal storage capacity is 10 MWh.  Table 2. The data of 120 days are used as the training set, and the data of the later 10 days are used for the testing. The simulation results are shown in Figure 11. The four curves are the individual SETS real values, CPM, PM, and the CM predicted value. The trained model is used to test the load condition in the following 10 days. From the training results, the CPM predicted value coincides with the real value. In the testing results, the comparison data in the orange dotted frame is enlarged for 121-130 days as seen in Figures 12-14. The result after the application of the CPM is significantly better than that individually predicted by the PM and CM alone.

Comparison of CPM with Real Value
The simulation results are shown in Figure 12. The prediction of the CPM coincides with the real value v. The prediction data are shown in Table 3. Comparing the real value v with the CPM, the predicted absolute error is e 1 . Selecting three days (123, 124, and 125) for calculation gives the CPM prediction with the MRE (∆ 1 = e 1 v × 100%) of 2.4%, 3.3%, and 5%, respectively. Therefore, the MRE (∆ 1 ) of the prediction by the CPM is not more than 5%.

Comparison of CPM with PM
The simulation results are shown in Figure 13. Compared with the real value v, the prediction error of the PM is larger than in the CPM. The prediction data are shown in Table 3. Comparing the real value v with the PM, the predicted absolute error is e 2 . Selecting three days (121, 125, and 130) for calculation gives the PM prediction with the MRE (∆ 2 = e 2 v × 100%) of 16.7%, 30.4%, and 20.6%, respectively, while the CPM prediction gives the MRE (∆ 1 ) of 1%, 5%, and 2.3%, respectively. Therefore, the prediction result of the CPM when compared with the PM has the MRE reduced by 25.4%. Figure 13. Comparison of cyber-physical model (CPM) with PM testing curves.

Comparison of CPM with CM
The simulation results are shown in Figure 14. Compared with the real value v, the prediction error of the CM is larger than in the CPM. The prediction data are shown in Table 3. Comparing the real value of v with the CM, the predicted absolute error is e 3 . Selecting three days (122, 123, and 125) for calculation gives the CM prediction with the MRE (∆ 3 = e 3 v × 100%) of 5.4%, 5.6%, and 7.2%, respectively, while the CPM prediction gives the MRE (∆ 1 ) of 0.6%, 2.4%, and 5%, respectively. Therefore, the prediction result of the CPM when compared with the CM has the MRE reduced by 4.8%.

Conclusions
Dispatching electricity consumption of SETS is an effective method for peak shaving in the power system. The prediction of thermal storage is a prerequisite for dynamic optimal dispatching, and the thermal storage can be used to shift times of resources. The CPS approach is proposed to predict power load, which promotes the application of the CPS in smart grids and ubiquitous power internet of things. This paper proposes a cyber-physical approach to predict the electricity consumption of SETS with consideration of the ambient temperature and electric behavior of the SETS into the PM. The CM adopts the BP neural network to calibrate the errors obtained in the PM. The 1MW SETS is established to validate the proposed cyber-physical approach. The simulation results show that when the CPM is compared with the PM, the MRE is reduced by 25.4%, and when compared with the CM is reduced by 4.8%. Using the CPM to calibrate the PM effectively improves the prediction accuracy. Conclusively, the prediction of SETS using the CPM is better than the individual PM and the CM alone.
The recommendations for future research are as follows.
• Thermal storage is an effective method for peak shaving and dispatching in power system. The electricity consumption prediction of SETS is worth to explore the combination with the heat and power generation unit.

•
The application of the proposed cyber-physical model in other resources of the power system is recommended to be studied.