# A Thermal Load Forecasting Algorithm Based on Trajectory Tracking

Authors: Guo Ling; Du Weiwei; Yang Jiancheng; Wang Liping; Wan Ping; Liu Ling
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The heating load forecast provides a basis for saving heating energy. Considering the nonstationary, nonlinear, and large time-delay characteristics of thermal load, this paper introduces the trajectory tracking stability theory into the field of load forecasting and proposes a heuristic correction that can ensure the convergence of forecast errors and does not depend on the system prediction model algorithm. The Lyapunov method is used to derive an error convergence criterion that has nothing to do with the prediction model, and a heuristic correction algorithm is designed for the predicted value with error divergence trend to ensure the error convergence of the load forecast sequence.


---

## Body

## 1. Introduction

Thermal load forecasting is based on the historical data of the heating process and to analyze and study the inherent changes in the thermal load demand which can be combined with weather, temperature, or holiday types. The purpose of studying it is to estimate and speculate the load demand in advance. It is the basis of thermal scheduling and energy-saving control [1].The heating process is a very complex dynamic system, and the relationship between energy and material transmission is highly nonlinear. It is difficult to establish a mathematical model to predict the heating load through a physical model, so most of the current forecasting methods are based on historical data. Based on statistical analysis [2–9], as early as 1984, Werner conducted tests on multiple district heating systems in Sweden and analyzed the main factors affecting the heat load, outdoor temperature, natural wind, and solar radiation. The results of the study showed that 60% of the total heat load of the heating network can be considered. It is affected by the outdoor temperature. The influence of natural wind will increase the heat load by 1% to 4%. The solar radiation heat can reduce the heat load by 1% to 5%. The heat consumption of domestic hot water will be affected by workdays or weekends. Different, the average domestic hot water load accounts for 30% of the total load, and the pipe network heat loss accounts for 5% to 8% of the total load. Wojdyga analyzed the relevant meteorological factors that affect the heat supply of the central heating system, including outdoor temperature, solar radiation intensity, and wind speed. Caldera et al. analyzed the relationship between the buildingʼs volume, shape factor, window-to-wall ratio, window transmittance, indoor and outdoor temperature, and other geometric parameters and thermal physical parameters and the buildingʼs heating energy consumption. Dotzauer introduced social factors into the heating load prediction and used the outdoor temperature–heating load curve to establish a heating load prediction model. By analyzing outdoor temperature, wind speed, solar radiation intensity, and other meteorological factors that affect heat load, Yetemen et al. established an autoregressive model for heat load forecasting, and Pedersen et al. established a linear regression model for heat load forecasting. Popescu and others further introduced factors such as indoor temperature, circulating water flow, and water supply temperature at the outlet of the heat exchange station into the calculation model of building heat consumption and obtained a multiple regression prediction model of building heat consumption. Nielsen et al. linked heat supply with meteorological data and date information and established a gray box forecast model for heat supply based on statistics-related theories.In addition to the abovementioned prediction models based on classical statistical theory, artificial intelligence methods represented by neural networks have also been widely used in heat load forecasting research [10–14]. Aydinalp et al. established a neural network-based building heating energy consumption and domestic hot water energy consumption model by inputting the characteristic parameters of the building, the performance parameters of the heating system, the average indoor temperature, the number of heating degree days, and the economic situation of the family. To predict the energy consumption of residential buildings in Canada, Kwok et al. used a combination of models and neural network predictors to predict building energy and conducted theoretical analysis. Kusiak et al. used humidity and temperature as inputs and used neural network methods to predict building steam load. Duanmu et al. proposed a novel HCLFM model to predict the cooling load of buildings. In addition, in order to meet the actual needs of small sample modeling, the support vector machine method is also used for building energy consumption prediction [15, 16]. Dong et al. used the support vector machine method to predict building energy consumption in tropical areas. Four commercial buildings in Singapore were randomly selected as the research object. The method selected outdoor monthly average temperature, relative humidity, and sunshine as input variables to predict tropical areas. The relative error of prediction of building energy consumption is less than 4%, which verifies the feasibility and adaptability of the support vector machine method to predict building energy consumption in tropical areas. Li et al., respectively, used support vector machines and neural networks to predict building energy consumption and compared the prediction results. These studies apply intelligent algorithms to heating load forecasting, but the research is conducted offline, staying at the laboratory simulation level, and no application examples of energy-saving control in heating engineering are given.Throughout the research status, the current load forecasting methods are mostly independent or combined use of the following methods: regression analysis [17, 18], time series [19], Kalman filter [20], neural network [21, 22], gray model, wavelet analysis, support vector machine [23], combination method [24], fuzzy processing [25], etc. Among them, regression analysis and time series methods have high requirements for historical data and large forecast errors; neural network convergence speed is slow, generalization ability is poor, and training is easy to overfit; gray system model has better performance when the original data has an exponential growth trend. High prediction accuracy: obviously, if you can combine multiple methods for research and application based on the advantages and disadvantages of various methods, you will get higher prediction accuracy, which also makes the idea of combined prediction that have developed in recent years. In addition, most of the existing forecasting methods are point forecasting methods, which cannot determine the range of possible fluctuations in the forecast results, and does not reflect the uncertain factors in the system. In these studies, the traditional heating load forecasting method is difficult to adapt to the heating load nonlinearity, resulting in low forecast accuracy, poor accuracy, and even failure, which cannot meet the needs of heating energy saving. The existing intelligent heating load forecasting method is researched for offline, staying at the laboratory simulation level, and no application examples of energy saving in heating engineering are given, and heating load forecasting software is not developed, making the theoretical results unable to be applied to actual projects. . However, with the continuous improvement of load forecasting requirements for heat dispatching and heating energy-saving control, it is necessary to further study the accuracy and accuracy of heating load forecasting, while ensuring the universality and stability of the algorithm.According to the actual demand of the heating system for the prediction accuracy of the heat load prediction, this paper uses the traditional ARMA time series method [26] for linear load prediction and BP neural network method for nonlinear load prediction according to the characteristics of the heating load. In view of the shortcomings of existing prediction methods, the idea of trajectory tracking control is introduced into the field of load forecasting [27]. A heuristic load forecasting algorithm based on trajectory tracking is proposed, which can ensure the convergence of forecast errors and does not depend on the system prediction model. The final simulation results analyze the performance comparison of several algorithms.


## 2. Classic Thermal Load Forecasting Method

### 2.1. Time Series Forecasting Method

The historical data of the heating load is an ordered collection sampled and recorded at a certain time interval, so it is a time series. The time series method is a more mature algorithm developed in the current heating system load forecast. According to the historical data of load, establish a mathematical model describing the change of the heating load with time, establish the expression of load prediction on the basis of this model, and predict the future load.As a typical representative of the time series model, the ARMA model is expressed as(1)Yt=ω1Yt−1+ω2Yt−2+⋯+ωpYt−p+et−φ1et−1−φ2et−2−⋯−φqet−q,where Yt is the value of the stationary time series at time t, et is the random interference acting on the time series at time t, the actual parameters ω1, ω2, …, ωp are the autoregressive coefficients, and the actual parameters φ1, φ2, …, φq are the movements average coefficients.By performing statistical processing on the above model, the autocorrelation coefficient can be obtained as(2)ρ1=1−ϕ1θ1ϕ1−θ11+θ12−2ϕ1θ1,ρ2=ϕ1ρ1,…ρk=ϕ1ρk−1.Therefore, the autocorrelation coefficientρ1 is a function of θ1 and ϕ1, and the autocorrelation function starts from ρ1 and decays exponentially. If ϕ1>0, the exponential decay of the autocorrelation function is smooth; if ϕ1<0, the decay of the autocorrelation function is alternating and oscillates between positive and negative values. The sign of ρ1 is determined by ϕ1−θ1, which determines the direction in which the exponential decay tends to zero. The initial value of its partial autocorrelation function is ϕ11=ρ1, which will decay exponentially thereafter. If ϕ1>0, ϕkk is smooth exponential decay; if ϕ1<0, ϕkk exponential decay is oscillating. The sign of ϕ11 is the same as ρ1, which is also determined by ϕ1−θ1.The prediction of the ARMA model requires that the sequence must be stable and that the influencing factors on the heating load must be basically the same within the time frame studied. If the given sequence is not a stationary sequence, you must preprocess the given sequence to make it stable, and then use the ARMA model to model.The heat load time series modeling process is shown in Figure1. The main steps are as follows:Step 1. Preprocess the sequence to determine whether the time series of heat load is a stationary impure random sequence. If it is a nonstationary sequence, process the sequence to meet the conditions of ARMA model modeling; that is, the processed sequence is a stationary non–white noise sequence.Step 2. Calculate the values of the sample autocorrelation coefficient (ACF) and sample partial autocorrelation coefficient (PACF) of the observation sequence.Step 3. According to the sample autocorrelation coefficient and partial autocorrelation coefficient, select the appropriate ARMA model to fit.Step 4. Estimate the value of unknown parameters in the model.Step 5. Test the validity of the model. If the fitting model fails the test, go to Step 3; reselect the model and then fit.Step 6. Model optimization, if the fitting model passes the test, still turn to Step 2, fully consider various possibilities, establish multiple fitting models, and select the optimal model from all the models that pass the test.Step 7. Use the fitted ARMA model to predict the future trend of the sequence, that is, the load value at the future moment.Figure 1
Flowchart of modeling based on time series method.The advantage of the time series method is that it requires less data, less workload, and faster calculation speed, which can reflect the continuity of the recent changes in load. The disadvantage of the time series method is that the modeling process is more complicated and requires higher theoretical knowledge, and the model has higher requirements for the stability of the original time series. It is only suitable for the prediction of relatively uniform load changes. These factors often lead to the model prediction error being large.


### 2.2. Neural Network Prediction Method

The artificial neural network method selects the load of the past period as a training sample, constructs a suitable network structure, and uses a certain training algorithm to train the network to meet the accuracy requirements. This neural network is used as a load prediction model in neural network prediction. In general, an optimized BP neural network is selected to establish a prediction model. The main work of modeling is as follows.


#### 2.2.1. Selection of Model Input and Output

The choice of predictive model input will have a great influence on the learning of the network. Regarding the load forecasting model of BP neural network, the difference mainly lies in the selection of input elements and the determination of the number of hidden layers. Since the current neural networkʼs ability to promote prediction is mainly determined through experiments, the selection of input variables is also based on empirical judgment. Generally, the real load data of 24 or 48 points per day in the days before the prediction day is selected as the model input, and the number of nodes in the input layer of the neural network is determined accordingly. In addition, in terms of the output of the model, the output of the model should be the load of the cycle corresponding to the predicted day.


#### 2.2.2. Determination of the Number of Network Layers

When each hidden node adopts the S-type compression function shown in the following formula, a hidden layer is enough to realize any decision classification problem, so a three-layer BP neural network is often used in heating load prediction:(3)y=11+e−x.


#### 2.2.3. Determination of the Number of Hidden Layers

To train a three-layer BP network with a single hidden layer, the number of hidden layer nodes can be selected according to the following empirical formula:(4)n1=n+m+a.In the formula,n is the number of input nodes, m is the number of output nodes, and a is a constant between 1 and 10. Because the selection of the number of nodes in each layer in the neural network based on the BP algorithm has a great influence on the performance of the network, the number of nodes in the hidden layer needs to be properly selected. After many times of training the data, and according to the training results, when the hidden layer node is about 15, the results obtained are more reasonable.


#### 2.2.4. Normalization Processing of Data

Before the neural network is trained and tested, in order to meet the network input and output data requirements and improve the convergence speed, in order to avoid the saturation of neurons, the network input data and actual output data need to be calculated according to equation (5) to [−1, 1]. In the interval, after using ANN prediction, the output layer is converted back to the load value according to equation (6):(5)y=2×x−xminxmax−xmin−1,(6)x=0.5y+1xmax−xmin+xmin.In the formula,x and y are the values before and after normalization of the input sample, and xmin and xmax are the minimum and maximum values of the load in the sample collection, respectively.


#### 2.2.5. Weight Initial Value Setting

The selection of the initial weight is one of many factors that affect the training of the network tending to have the smallest acceptable error. The initialization of the weight strongly affects the final solution. Generally, the initial weight of the network is selected aswmax=1.0.The advantages of artificial neural networks are the ability to adapt to a large number of unstructured and nonprecision laws, the characteristics of information memory, autonomous learning, knowledge reasoning, and optimized calculation, as well as strong calculation ability, complex mapping ability, and fault tolerance. And various intelligent processing capabilities, especially its learning and adaptive functions, are not available in other algorithms. The shortcomings are that the number of neural network layers and the number of neurons are mostly determined based on subjective experience, it is difficult to determine the network structure scientifically, the learning speed is slow, and there are local minimum points.


## 2.1. Time Series Forecasting Method

The historical data of the heating load is an ordered collection sampled and recorded at a certain time interval, so it is a time series. The time series method is a more mature algorithm developed in the current heating system load forecast. According to the historical data of load, establish a mathematical model describing the change of the heating load with time, establish the expression of load prediction on the basis of this model, and predict the future load.As a typical representative of the time series model, the ARMA model is expressed as(1)Yt=ω1Yt−1+ω2Yt−2+⋯+ωpYt−p+et−φ1et−1−φ2et−2−⋯−φqet−q,where Yt is the value of the stationary time series at time t, et is the random interference acting on the time series at time t, the actual parameters ω1, ω2, …, ωp are the autoregressive coefficients, and the actual parameters φ1, φ2, …, φq are the movements average coefficients.By performing statistical processing on the above model, the autocorrelation coefficient can be obtained as(2)ρ1=1−ϕ1θ1ϕ1−θ11+θ12−2ϕ1θ1,ρ2=ϕ1ρ1,…ρk=ϕ1ρk−1.Therefore, the autocorrelation coefficientρ1 is a function of θ1 and ϕ1, and the autocorrelation function starts from ρ1 and decays exponentially. If ϕ1>0, the exponential decay of the autocorrelation function is smooth; if ϕ1<0, the decay of the autocorrelation function is alternating and oscillates between positive and negative values. The sign of ρ1 is determined by ϕ1−θ1, which determines the direction in which the exponential decay tends to zero. The initial value of its partial autocorrelation function is ϕ11=ρ1, which will decay exponentially thereafter. If ϕ1>0, ϕkk is smooth exponential decay; if ϕ1<0, ϕkk exponential decay is oscillating. The sign of ϕ11 is the same as ρ1, which is also determined by ϕ1−θ1.The prediction of the ARMA model requires that the sequence must be stable and that the influencing factors on the heating load must be basically the same within the time frame studied. If the given sequence is not a stationary sequence, you must preprocess the given sequence to make it stable, and then use the ARMA model to model.The heat load time series modeling process is shown in Figure1. The main steps are as follows:Step 1. Preprocess the sequence to determine whether the time series of heat load is a stationary impure random sequence. If it is a nonstationary sequence, process the sequence to meet the conditions of ARMA model modeling; that is, the processed sequence is a stationary non–white noise sequence.Step 2. Calculate the values of the sample autocorrelation coefficient (ACF) and sample partial autocorrelation coefficient (PACF) of the observation sequence.Step 3. According to the sample autocorrelation coefficient and partial autocorrelation coefficient, select the appropriate ARMA model to fit.Step 4. Estimate the value of unknown parameters in the model.Step 5. Test the validity of the model. If the fitting model fails the test, go to Step 3; reselect the model and then fit.Step 6. Model optimization, if the fitting model passes the test, still turn to Step 2, fully consider various possibilities, establish multiple fitting models, and select the optimal model from all the models that pass the test.Step 7. Use the fitted ARMA model to predict the future trend of the sequence, that is, the load value at the future moment.Figure 1
Flowchart of modeling based on time series method.The advantage of the time series method is that it requires less data, less workload, and faster calculation speed, which can reflect the continuity of the recent changes in load. The disadvantage of the time series method is that the modeling process is more complicated and requires higher theoretical knowledge, and the model has higher requirements for the stability of the original time series. It is only suitable for the prediction of relatively uniform load changes. These factors often lead to the model prediction error being large.


## 2.2. Neural Network Prediction Method

The artificial neural network method selects the load of the past period as a training sample, constructs a suitable network structure, and uses a certain training algorithm to train the network to meet the accuracy requirements. This neural network is used as a load prediction model in neural network prediction. In general, an optimized BP neural network is selected to establish a prediction model. The main work of modeling is as follows.


### 2.2.1. Selection of Model Input and Output

The choice of predictive model input will have a great influence on the learning of the network. Regarding the load forecasting model of BP neural network, the difference mainly lies in the selection of input elements and the determination of the number of hidden layers. Since the current neural networkʼs ability to promote prediction is mainly determined through experiments, the selection of input variables is also based on empirical judgment. Generally, the real load data of 24 or 48 points per day in the days before the prediction day is selected as the model input, and the number of nodes in the input layer of the neural network is determined accordingly. In addition, in terms of the output of the model, the output of the model should be the load of the cycle corresponding to the predicted day.


### 2.2.2. Determination of the Number of Network Layers

When each hidden node adopts the S-type compression function shown in the following formula, a hidden layer is enough to realize any decision classification problem, so a three-layer BP neural network is often used in heating load prediction:(3)y=11+e−x.


### 2.2.3. Determination of the Number of Hidden Layers

To train a three-layer BP network with a single hidden layer, the number of hidden layer nodes can be selected according to the following empirical formula:(4)n1=n+m+a.In the formula,n is the number of input nodes, m is the number of output nodes, and a is a constant between 1 and 10. Because the selection of the number of nodes in each layer in the neural network based on the BP algorithm has a great influence on the performance of the network, the number of nodes in the hidden layer needs to be properly selected. After many times of training the data, and according to the training results, when the hidden layer node is about 15, the results obtained are more reasonable.


### 2.2.4. Normalization Processing of Data

Before the neural network is trained and tested, in order to meet the network input and output data requirements and improve the convergence speed, in order to avoid the saturation of neurons, the network input data and actual output data need to be calculated according to equation (5) to [−1, 1]. In the interval, after using ANN prediction, the output layer is converted back to the load value according to equation (6):(5)y=2×x−xminxmax−xmin−1,(6)x=0.5y+1xmax−xmin+xmin.In the formula,x and y are the values before and after normalization of the input sample, and xmin and xmax are the minimum and maximum values of the load in the sample collection, respectively.


### 2.2.5. Weight Initial Value Setting

The selection of the initial weight is one of many factors that affect the training of the network tending to have the smallest acceptable error. The initialization of the weight strongly affects the final solution. Generally, the initial weight of the network is selected aswmax=1.0.The advantages of artificial neural networks are the ability to adapt to a large number of unstructured and nonprecision laws, the characteristics of information memory, autonomous learning, knowledge reasoning, and optimized calculation, as well as strong calculation ability, complex mapping ability, and fault tolerance. And various intelligent processing capabilities, especially its learning and adaptive functions, are not available in other algorithms. The shortcomings are that the number of neural network layers and the number of neurons are mostly determined based on subjective experience, it is difficult to determine the network structure scientifically, the learning speed is slow, and there are local minimum points.


## 2.2.1. Selection of Model Input and Output

The choice of predictive model input will have a great influence on the learning of the network. Regarding the load forecasting model of BP neural network, the difference mainly lies in the selection of input elements and the determination of the number of hidden layers. Since the current neural networkʼs ability to promote prediction is mainly determined through experiments, the selection of input variables is also based on empirical judgment. Generally, the real load data of 24 or 48 points per day in the days before the prediction day is selected as the model input, and the number of nodes in the input layer of the neural network is determined accordingly. In addition, in terms of the output of the model, the output of the model should be the load of the cycle corresponding to the predicted day.


## 2.2.2. Determination of the Number of Network Layers

When each hidden node adopts the S-type compression function shown in the following formula, a hidden layer is enough to realize any decision classification problem, so a three-layer BP neural network is often used in heating load prediction:(3)y=11+e−x.


## 2.2.3. Determination of the Number of Hidden Layers

To train a three-layer BP network with a single hidden layer, the number of hidden layer nodes can be selected according to the following empirical formula:(4)n1=n+m+a.In the formula,n is the number of input nodes, m is the number of output nodes, and a is a constant between 1 and 10. Because the selection of the number of nodes in each layer in the neural network based on the BP algorithm has a great influence on the performance of the network, the number of nodes in the hidden layer needs to be properly selected. After many times of training the data, and according to the training results, when the hidden layer node is about 15, the results obtained are more reasonable.


## 2.2.4. Normalization Processing of Data

Before the neural network is trained and tested, in order to meet the network input and output data requirements and improve the convergence speed, in order to avoid the saturation of neurons, the network input data and actual output data need to be calculated according to equation (5) to [−1, 1]. In the interval, after using ANN prediction, the output layer is converted back to the load value according to equation (6):(5)y=2×x−xminxmax−xmin−1,(6)x=0.5y+1xmax−xmin+xmin.In the formula,x and y are the values before and after normalization of the input sample, and xmin and xmax are the minimum and maximum values of the load in the sample collection, respectively.


## 2.2.5. Weight Initial Value Setting

The selection of the initial weight is one of many factors that affect the training of the network tending to have the smallest acceptable error. The initialization of the weight strongly affects the final solution. Generally, the initial weight of the network is selected aswmax=1.0.The advantages of artificial neural networks are the ability to adapt to a large number of unstructured and nonprecision laws, the characteristics of information memory, autonomous learning, knowledge reasoning, and optimized calculation, as well as strong calculation ability, complex mapping ability, and fault tolerance. And various intelligent processing capabilities, especially its learning and adaptive functions, are not available in other algorithms. The shortcomings are that the number of neural network layers and the number of neurons are mostly determined based on subjective experience, it is difficult to determine the network structure scientifically, the learning speed is slow, and there are local minimum points.


## 3. Introduction of Trajectory Tracking Stability Theory

### 3.1. Trajectory Tracking Error Model

Considering the existence of a moving target in two-dimensional space, its ideal trajectory is shown in curveC in Figure 2 in the inertial coordinate system.Figure 2
Schematic diagram of target trajectory tracking.Suppose there is a moving pointPx,yT to track the trajectory C, and the angle between the velocity v of point P and the positive direction of the y-axis is θ; that is,(7)θ=arctanx˙y˙.Then, the position motion model of pointP can be expressed as(8)P˙=x˙y˙=vsinθvcosθ=Rθ0v.Among them,(9)Rθ=cosθsinθ−sinθcosθ.At any time, the trajectory coordinate system is defined as follows: taking the desired position of pointP (that is, the true position of the target) Pa=xa,yaT as the origin, the ya-axis is along the tangent of Pa, and the xa-axis is perpendicular to ya and pointing to the right, as shown in Figure 3.Figure 3
Schematic diagram of coordinated system and tracking error.Then, the rotation angleθa of the trajectory coordinate system relative to the inertial coordinate system is(10)θa=arctanx˙ay˙a.Referring to equation (8), the motion model of the desired position Pa can be expressed as(11)P˙a=x˙ay˙a=vasinθavacosθa=Rθa0va.Among them,va is the desired speed of point P (that is, the true speed of the target), and the direction is consistent with the ya-axis, and Rθa and Rθ have the same form.Therefore, the tracking errorε of P and Pa in the trajectory coordinate system is(12)ε=es=ΔXcosθa−ΔsinθaΔXsinθa+ΔYcosθa=cosθasinθa−sinθacosθaTΔXΔY=RTθaP−Pa.Among them,e and s are the lateral and forward tracking errors of P and Pa in the trajectory coordinate system, respectively, and ΔX and ΔY are the lateral and longitudinal tracking errors of P and Pa in the inertial coordinate system, respectively.From equation (12), we can know that εTε=P−PaTRRTP−Pa=P−Pa2. Obviously, there is only if ε ⟶ 0, so the fundamental goal of P to track the ideal trajectory is to make ε ⟶ 0.


### 3.2. Trajectory Tracking Stability Criterion

In control theory, the Lyapunov criterion [28] described in Lemma 1 is commonly used to judge the stability of the tracking system and to check the error convergence trend.Lemma 1.
For continuous nonlinear time-varying systemsχ˙=fχ,t≥0, if a scalar function Vχ with continuous first-order partial derivatives for χ can be constructed, V0=0 satisfies the condition for all nonzero state points χ in the state space ℜ: ① Vχ is positive definite; ② V˙χ=dVχ/dt is negative definite; and ③when χ⟶∞, then Vχ⟶∞. Then, the system equilibrium state χ=0 is asymptotically stable in a large range; that is, χ asymptotically converges to zero.Based on the above lemma to analyze the tracking error convergence, the key is to construct the Lyapunov candidate functionVχ, and different Lyapunov candidate functions may obtain different stability criteria. Considering that the main purpose of this article is to verify the applicability of the trajectory tracking stability theory in load forecasting, no special research has been conducted on the selection of Lyapunov candidate functions, so the most commonly used Lyapunov function type is directly used:(13)V=12εTε.Then,V derivates t:(14)V˙=εTε˙=εTR˙TθaP−Pa+RTθaP˙−P˙a=εT0−θ˙aθ˙a0RTθaP−Pa+RTθaP˙−P˙a=εT0−θ˙aθ˙a0ε+εTRTθaP˙−P˙a=εTRTθaP˙−P˙a.Therefore, if equation (14) is negatively set at ε≠0, the tracking error of the moving point P to the target trajectory will converge to zero asymptotically, which can ensure the stability of trajectory tracking.If the actual load sequence is regarded as a trajectory on the time axis, high-precision load prediction is to realize the high-precision tracking of the trajectory composed of predicted points and the actual load trajectory. Therefore, the convergence analysis of the load forecast error is essentially equivalent to the convergence analysis of the tracking error of the predicted point.


## 3.1. Trajectory Tracking Error Model

Considering the existence of a moving target in two-dimensional space, its ideal trajectory is shown in curveC in Figure 2 in the inertial coordinate system.Figure 2
Schematic diagram of target trajectory tracking.Suppose there is a moving pointPx,yT to track the trajectory C, and the angle between the velocity v of point P and the positive direction of the y-axis is θ; that is,(7)θ=arctanx˙y˙.Then, the position motion model of pointP can be expressed as(8)P˙=x˙y˙=vsinθvcosθ=Rθ0v.Among them,(9)Rθ=cosθsinθ−sinθcosθ.At any time, the trajectory coordinate system is defined as follows: taking the desired position of pointP (that is, the true position of the target) Pa=xa,yaT as the origin, the ya-axis is along the tangent of Pa, and the xa-axis is perpendicular to ya and pointing to the right, as shown in Figure 3.Figure 3
Schematic diagram of coordinated system and tracking error.Then, the rotation angleθa of the trajectory coordinate system relative to the inertial coordinate system is(10)θa=arctanx˙ay˙a.Referring to equation (8), the motion model of the desired position Pa can be expressed as(11)P˙a=x˙ay˙a=vasinθavacosθa=Rθa0va.Among them,va is the desired speed of point P (that is, the true speed of the target), and the direction is consistent with the ya-axis, and Rθa and Rθ have the same form.Therefore, the tracking errorε of P and Pa in the trajectory coordinate system is(12)ε=es=ΔXcosθa−ΔsinθaΔXsinθa+ΔYcosθa=cosθasinθa−sinθacosθaTΔXΔY=RTθaP−Pa.Among them,e and s are the lateral and forward tracking errors of P and Pa in the trajectory coordinate system, respectively, and ΔX and ΔY are the lateral and longitudinal tracking errors of P and Pa in the inertial coordinate system, respectively.From equation (12), we can know that εTε=P−PaTRRTP−Pa=P−Pa2. Obviously, there is only if ε ⟶ 0, so the fundamental goal of P to track the ideal trajectory is to make ε ⟶ 0.


## 3.2. Trajectory Tracking Stability Criterion

In control theory, the Lyapunov criterion [28] described in Lemma 1 is commonly used to judge the stability of the tracking system and to check the error convergence trend.Lemma 1.
For continuous nonlinear time-varying systemsχ˙=fχ,t≥0, if a scalar function Vχ with continuous first-order partial derivatives for χ can be constructed, V0=0 satisfies the condition for all nonzero state points χ in the state space ℜ: ① Vχ is positive definite; ② V˙χ=dVχ/dt is negative definite; and ③when χ⟶∞, then Vχ⟶∞. Then, the system equilibrium state χ=0 is asymptotically stable in a large range; that is, χ asymptotically converges to zero.Based on the above lemma to analyze the tracking error convergence, the key is to construct the Lyapunov candidate functionVχ, and different Lyapunov candidate functions may obtain different stability criteria. Considering that the main purpose of this article is to verify the applicability of the trajectory tracking stability theory in load forecasting, no special research has been conducted on the selection of Lyapunov candidate functions, so the most commonly used Lyapunov function type is directly used:(13)V=12εTε.Then,V derivates t:(14)V˙=εTε˙=εTR˙TθaP−Pa+RTθaP˙−P˙a=εT0−θ˙aθ˙a0RTθaP−Pa+RTθaP˙−P˙a=εT0−θ˙aθ˙a0ε+εTRTθaP˙−P˙a=εTRTθaP˙−P˙a.Therefore, if equation (14) is negatively set at ε≠0, the tracking error of the moving point P to the target trajectory will converge to zero asymptotically, which can ensure the stability of trajectory tracking.If the actual load sequence is regarded as a trajectory on the time axis, high-precision load prediction is to realize the high-precision tracking of the trajectory composed of predicted points and the actual load trajectory. Therefore, the convergence analysis of the load forecast error is essentially equivalent to the convergence analysis of the tracking error of the predicted point.


## 4. Process of Heat Load Prediction Algorithm Based on Trajectory Tracking

### 4.1. Condition of Error Convergence for Heating Load Forecast

Theorem 1.
At any time, ifθ ＞ θa holds when ΔY > 0, or θ ＜ θa holds when ΔY < 0, the prediction error will converge to zero asymptotically.In a two-dimensional coordinate system with thex-axis as the time and the y-axis as the load value, the load time series can be regarded as a trajectory. If moving point P=x,yT represents the predicted load value y at time x and Pa=xa,yaT represents the true load value ya at time x, since both have the same time axis, the lateral error of point P and its expected position Pa in the inertial coordinate system must meet ΔX = 0, as shown in Figure 4.Therefore, we can further expand equation (14) to obtain(15)V˙=escosθa−sinθasinθacosθavsinθ−vasinθavcosθ−vacosθa=svcosθ−θa−va+evsinθ−θa=svcosθr+evsinθr−sva.When ΔY = 0, it is known from ΔX = 0 and equation (12) that ε = 0; that is, the prediction error converges.When ΔY ≠ 0, considering thatP and Pa have the same time axis, there are only two possible geometric relationships; that is, point P is directly above or below point Pa, so we can solve the prediction error according to the following two cases. The sufficient condition for the convergence of the prediction error can be expressed as follows:Figure 4
Convergence analysis of forecasting error. (a) ΔY > 0, θa < π/2; (b) ΔY > 0, θa > π/2; (c) ΔY < 0, θa < π/2; and (d) ΔY < 0, θa > π/2.
(a)(b)(c)(d)Case 1.
WhenP is directly above Pa, that is, ΔY > 0.
From the trajectory coordinate system in Figure4, the following geometric features can be obtained intuitively:① Ifθa<π/2, then e < 0, s > 0, as shown in Figure 4(a).② Ifθa>π/2, then e < 0, s < 0, as shown in Figure 4(b).③ For anyθa, there are(16)cosθa=ss2+e2,sinθa=−es2+e2.
Then, by equation (15),(17)V˙=svcosθr+evsinθr−sva=vs2+e2scosθrs2+e2−−esinθrs2+e2−svavs2+e2=vs2+e2cosθ−svavs2+e2.
In the actual load curve, the real value sequence adds a data point every time ΔT, and the predicted value sequence also adds a data point every time ΔT, so v ＝ van and v>0; that is, in a two-dimensional coordinate system with the x-axis as the time axis, the actual load value and the predicted value extend toward the positive direction of the x-axis at the same speed. From equation (17) and Lemma 1, if and only if cosθ<sva/vs2+e2=cosθa, that is, when θ>θa, V˙<0, the prediction error converges to zero.Case 2.
WhenP is directly below Pa, that is, ΔY < 0.
From the track coordinate system of Figure4, the following geometric features can also be obtained intuitively:① Ifθa<π/2, then e > 0, s < 0, as shown in Figure 4(c).② Ifθa>π/2, then e > 0, s ＞ 0, as shown in Figure 4(d).③ For anyθa, there are(18)cosθa=−ss2+e2,sinθa=es2+e2.
Then, by equation (15),(19)V˙=svcosθr+evsinθr−sva=vs2+e2scosθrs2+e2+esinθrs2+e2−svavs2+e2=vs2+e2cosθ+svavs2+e2.
An analysis similar to Case 1 can be made. Therefore, from equation (19) and Lemma 1, it can also be seen that, if and only if cosθ>−sva/vs2+e2=cosθa, that is, when θ<θa, V˙<0, the prediction error converges to zero.
Noticed that the above analysis process is completely derived from trajectory tracking and its stability theory and has nothing to do with the specific prediction model used. From this, we have obtained the conclusions about the error convergence conditions that have nothing to do with the prediction model in this paper.


### 4.2. Error Correction Strategy for Heating Load Forecast

At any sampling timex, the historical predicted value yx and the real values yax and yax−1 are known, and the predicted value yx+1 at the next time can be calculated by any existing prediction model (defined as the preliminary prediction model); then the value Δy, the angle θ between the ray from point x,yx to point x+1,yx+1 and the positive direction of the y-axis, and the angle θa between the ray from point x−1,yax−1 to point x,yax and the positive direction of the y-axis can be calculated. Then, according to Theorem 1, it can be judged whether the prediction error convergence condition is satisfied.Obviously, if the error convergence condition is satisfied, it indicates that the predicted valueyx+1 solved by the preliminary prediction model has a tendency to converge to zero, and no correction is needed at this time; if the error convergence condition is not satisfied, it indicates that the predicted value yx+1 solved by the preliminary prediction model is implicit. The trend of error divergence should be corrected at this time.For the predicted value with error divergence trend, this paper further proposes an adaptive correction algorithm for the predicted value. The basic idea of the algorithm is to calculateΔY and tanθa from known conditions and then use the monotonic feature of the tangent function to adjust the value of tanθ based on Theorem 1, so that the predicted value yx+1 at any time meets the error convergence condition. The method of recalculating the predicted value yx+1 is as follows:① IfΔY>0 and yax−yax−1=0, then(20)yx+1=k1yx,0<k1<1.② IfΔY>0 and yax−yax−1>0, then(21)yx+1=yx+k2yax−yax−1,0<k2<1.③ IfΔY>0 and yax−yax−1<0, then(22)yx+1=yx+k3yax−yax−1,k3>1.④ IfΔY<0 and yax−yax−1=0, then(23)yx+1=k4yx,k4>1.⑤ IfΔY<0 and yax−yax−1>0, then(24)yx+1=yx+k5yax−yax−1,k5>1.⑥ IfΔY<0 and yax−yax−1<0, then(25)yx+1=yx+k6yax−yax−1,0<k6<1.The correction coefficientsk1–k6 are selected by the user. The following theoretically proves the effectiveness of this revised algorithm.Proof.
TakeX > 0 as an example.
Whenyax−yax−1, then θa=π/2. From equations (7) and (20), we can get(26)tanθ=1yx+1−yx=1k1−1yx.
Since0<k1<1 and yx>0, then tanθ<0 and θ>π/2, so θ>θa holds. When yax−yax−1>0, from equations (7) and (10), we can find that equation (7) is equivalent to(27)tanθa=1yax−yax−1=k2yx+1−yx=k2tanθ.
Since0<k2<1 and tanθa>0, according to the monotonic characteristics of tangent function, θ>θa was established.
Whenyax−yax−1<0, from equations (7) and (10), we can find that equation (22) is equivalent to(28)tanθa=1yax−yax−1=k3yx+1−yx=k3tanθ.
Sincek3>1 and tanθa<0, so tanθa<tanθ<0, according to the monotonic characteristics of tangent function, θ>θa was established.
In summary,ΔY>0 holds when θ>θa, and the revised predicted value yx+1 meets the error convergence condition of Theorem 1. Similarly, it can be proved that ΔY<0, so the correction algorithm is theoretically effective.


### 4.3. Full Description Based on Trajectory Tracking

According to the previous analysis and derivation process, the load forecasting algorithm proposed in this paper is described as follows:Step 1. Choose an existing prediction model (such as ARMA or BPNN model) to establish a preliminary prediction model of load time series.Step 2. At time x, use the preliminary prediction model to calculate the preliminary prediction value yx+1 at time.Step 3. From the predicted values yx and yx−1 at times x and x+1, use equation (7) to calculate the angle θ between the ray from point x,yx to point x+1,yx+1 and the positive direction of the y-axis.Step 4. From the predicted values yax and yax−1 at times x and x−1, use equation (9) to calculate the angle θa between the ray from point x−1,yax−1 to point x,yax and the positive direction of the y-axis.Step 5. From the predicted value yx and the real value yax at time x, we can calculate ΔY=yx−yax.Step 6. According to Theorem 1 of this paper, judge whether ΔY, θ, and θa meet the convergence condition of the prediction error. If it is satisfied, the predicted value yx+1 solved by the preliminary prediction model is considered valid; otherwise, go to Step 7.Step 7. According to the prediction value correction algorithm proposed in this paper, solve the load prediction value yx+1 at time x + 1 again.From the previous analysis and algorithm description, this method is actually a heuristic correction algorithm for the existing prediction model. The heuristic rules of the algorithm are derived based on the theory of trajectory tracking stability and have nothing to do with the specific prediction model selected, so it is a heuristic correction algorithm that can ensure the convergence of prediction errors and does not depend on the system model. In addition, during the iterative process along the sampling time, the algorithm feeds back the inherent laws contained in the historical prediction sequence to the future prediction, which improves the prediction accuracy.


## 4.1. Condition of Error Convergence for Heating Load Forecast

Theorem 1.
At any time, ifθ ＞ θa holds when ΔY > 0, or θ ＜ θa holds when ΔY < 0, the prediction error will converge to zero asymptotically.In a two-dimensional coordinate system with thex-axis as the time and the y-axis as the load value, the load time series can be regarded as a trajectory. If moving point P=x,yT represents the predicted load value y at time x and Pa=xa,yaT represents the true load value ya at time x, since both have the same time axis, the lateral error of point P and its expected position Pa in the inertial coordinate system must meet ΔX = 0, as shown in Figure 4.Therefore, we can further expand equation (14) to obtain(15)V˙=escosθa−sinθasinθacosθavsinθ−vasinθavcosθ−vacosθa=svcosθ−θa−va+evsinθ−θa=svcosθr+evsinθr−sva.When ΔY = 0, it is known from ΔX = 0 and equation (12) that ε = 0; that is, the prediction error converges.When ΔY ≠ 0, considering thatP and Pa have the same time axis, there are only two possible geometric relationships; that is, point P is directly above or below point Pa, so we can solve the prediction error according to the following two cases. The sufficient condition for the convergence of the prediction error can be expressed as follows:Figure 4
Convergence analysis of forecasting error. (a) ΔY > 0, θa < π/2; (b) ΔY > 0, θa > π/2; (c) ΔY < 0, θa < π/2; and (d) ΔY < 0, θa > π/2.
(a)(b)(c)(d)Case 1.
WhenP is directly above Pa, that is, ΔY > 0.
From the trajectory coordinate system in Figure4, the following geometric features can be obtained intuitively:① Ifθa<π/2, then e < 0, s > 0, as shown in Figure 4(a).② Ifθa>π/2, then e < 0, s < 0, as shown in Figure 4(b).③ For anyθa, there are(16)cosθa=ss2+e2,sinθa=−es2+e2.
Then, by equation (15),(17)V˙=svcosθr+evsinθr−sva=vs2+e2scosθrs2+e2−−esinθrs2+e2−svavs2+e2=vs2+e2cosθ−svavs2+e2.
In the actual load curve, the real value sequence adds a data point every time ΔT, and the predicted value sequence also adds a data point every time ΔT, so v ＝ van and v>0; that is, in a two-dimensional coordinate system with the x-axis as the time axis, the actual load value and the predicted value extend toward the positive direction of the x-axis at the same speed. From equation (17) and Lemma 1, if and only if cosθ<sva/vs2+e2=cosθa, that is, when θ>θa, V˙<0, the prediction error converges to zero.Case 2.
WhenP is directly below Pa, that is, ΔY < 0.
From the track coordinate system of Figure4, the following geometric features can also be obtained intuitively:① Ifθa<π/2, then e > 0, s < 0, as shown in Figure 4(c).② Ifθa>π/2, then e > 0, s ＞ 0, as shown in Figure 4(d).③ For anyθa, there are(18)cosθa=−ss2+e2,sinθa=es2+e2.
Then, by equation (15),(19)V˙=svcosθr+evsinθr−sva=vs2+e2scosθrs2+e2+esinθrs2+e2−svavs2+e2=vs2+e2cosθ+svavs2+e2.
An analysis similar to Case 1 can be made. Therefore, from equation (19) and Lemma 1, it can also be seen that, if and only if cosθ>−sva/vs2+e2=cosθa, that is, when θ<θa, V˙<0, the prediction error converges to zero.
Noticed that the above analysis process is completely derived from trajectory tracking and its stability theory and has nothing to do with the specific prediction model used. From this, we have obtained the conclusions about the error convergence conditions that have nothing to do with the prediction model in this paper.


## 4.2. Error Correction Strategy for Heating Load Forecast

At any sampling timex, the historical predicted value yx and the real values yax and yax−1 are known, and the predicted value yx+1 at the next time can be calculated by any existing prediction model (defined as the preliminary prediction model); then the value Δy, the angle θ between the ray from point x,yx to point x+1,yx+1 and the positive direction of the y-axis, and the angle θa between the ray from point x−1,yax−1 to point x,yax and the positive direction of the y-axis can be calculated. Then, according to Theorem 1, it can be judged whether the prediction error convergence condition is satisfied.Obviously, if the error convergence condition is satisfied, it indicates that the predicted valueyx+1 solved by the preliminary prediction model has a tendency to converge to zero, and no correction is needed at this time; if the error convergence condition is not satisfied, it indicates that the predicted value yx+1 solved by the preliminary prediction model is implicit. The trend of error divergence should be corrected at this time.For the predicted value with error divergence trend, this paper further proposes an adaptive correction algorithm for the predicted value. The basic idea of the algorithm is to calculateΔY and tanθa from known conditions and then use the monotonic feature of the tangent function to adjust the value of tanθ based on Theorem 1, so that the predicted value yx+1 at any time meets the error convergence condition. The method of recalculating the predicted value yx+1 is as follows:① IfΔY>0 and yax−yax−1=0, then(20)yx+1=k1yx,0<k1<1.② IfΔY>0 and yax−yax−1>0, then(21)yx+1=yx+k2yax−yax−1,0<k2<1.③ IfΔY>0 and yax−yax−1<0, then(22)yx+1=yx+k3yax−yax−1,k3>1.④ IfΔY<0 and yax−yax−1=0, then(23)yx+1=k4yx,k4>1.⑤ IfΔY<0 and yax−yax−1>0, then(24)yx+1=yx+k5yax−yax−1,k5>1.⑥ IfΔY<0 and yax−yax−1<0, then(25)yx+1=yx+k6yax−yax−1,0<k6<1.The correction coefficientsk1–k6 are selected by the user. The following theoretically proves the effectiveness of this revised algorithm.Proof.
TakeX > 0 as an example.
Whenyax−yax−1, then θa=π/2. From equations (7) and (20), we can get(26)tanθ=1yx+1−yx=1k1−1yx.
Since0<k1<1 and yx>0, then tanθ<0 and θ>π/2, so θ>θa holds. When yax−yax−1>0, from equations (7) and (10), we can find that equation (7) is equivalent to(27)tanθa=1yax−yax−1=k2yx+1−yx=k2tanθ.
Since0<k2<1 and tanθa>0, according to the monotonic characteristics of tangent function, θ>θa was established.
Whenyax−yax−1<0, from equations (7) and (10), we can find that equation (22) is equivalent to(28)tanθa=1yax−yax−1=k3yx+1−yx=k3tanθ.
Sincek3>1 and tanθa<0, so tanθa<tanθ<0, according to the monotonic characteristics of tangent function, θ>θa was established.
In summary,ΔY>0 holds when θ>θa, and the revised predicted value yx+1 meets the error convergence condition of Theorem 1. Similarly, it can be proved that ΔY<0, so the correction algorithm is theoretically effective.


## 4.3. Full Description Based on Trajectory Tracking

According to the previous analysis and derivation process, the load forecasting algorithm proposed in this paper is described as follows:Step 1. Choose an existing prediction model (such as ARMA or BPNN model) to establish a preliminary prediction model of load time series.Step 2. At time x, use the preliminary prediction model to calculate the preliminary prediction value yx+1 at time.Step 3. From the predicted values yx and yx−1 at times x and x+1, use equation (7) to calculate the angle θ between the ray from point x,yx to point x+1,yx+1 and the positive direction of the y-axis.Step 4. From the predicted values yax and yax−1 at times x and x−1, use equation (9) to calculate the angle θa between the ray from point x−1,yax−1 to point x,yax and the positive direction of the y-axis.Step 5. From the predicted value yx and the real value yax at time x, we can calculate ΔY=yx−yax.Step 6. According to Theorem 1 of this paper, judge whether ΔY, θ, and θa meet the convergence condition of the prediction error. If it is satisfied, the predicted value yx+1 solved by the preliminary prediction model is considered valid; otherwise, go to Step 7.Step 7. According to the prediction value correction algorithm proposed in this paper, solve the load prediction value yx+1 at time x + 1 again.From the previous analysis and algorithm description, this method is actually a heuristic correction algorithm for the existing prediction model. The heuristic rules of the algorithm are derived based on the theory of trajectory tracking stability and have nothing to do with the specific prediction model selected, so it is a heuristic correction algorithm that can ensure the convergence of prediction errors and does not depend on the system model. In addition, during the iterative process along the sampling time, the algorithm feeds back the inherent laws contained in the historical prediction sequence to the future prediction, which improves the prediction accuracy.


## 5. Simulation and Result Analysis

It is of no practical significance to compare this method with other autonomous prediction algorithms alone. The effectiveness of this method is reflected in the revision and improvement of existing prediction models. In the experimental part, the classic autoregressive moving average (ARMA) model [29] and Back Propagation Neural Network (BPNN) model [30] are taken as examples to compare and analyze the correction effect of the method under different prediction models and prediction field of view conditions. Use the load sample data GEFCom2012 published by the Global Energy Forecasting Competition to simulate on the MATLAB platform. In all experiments, the parameters of the modified algorithm in this paper are selected as k1=k2=k6=0.9, k3=k4=k5=1.1.


### 5.1. Comparison Experiment with the ARMA Prediction Model

Based on the actual measurement of a certain residential area in Jinan on December 9, 2018, the load value after a single prediction field of view of 0.5 h was predicted and experimentally measured, using a separate ARMA model (ARMA method) and based on the ARMA model. The error-corrected method (ARMA + EC method) obtained the forecast and measured results at 48 oʼclock on the day, and a comparative analysis was conducted.Based on the historical load value of 120 points totaling 60 hours before the current day as the data source, an ARMA prediction model was established. The classic (p, q)-order ARMA prediction model can be expressed by equation (1).The main work of modeling with ARMA lies in model order determination and coefficient estimation. Firstly, use Autocorrelation Function (ACF) to analyze and remove the fixed trend items in the original data source. Then, the Bayesian Information Criterion (BIC), Yule–Walker, and Newton–Raphson algorithms are used to determine the model order, autoregressive coefficient, and moving average coefficient, respectively. In the experiment, the BIC values corresponding to each order ARMA candidate model are shown in Table1.Table 1
BIC values for ARMA model candidate.
Orderq=1q=2q=3q=4q=5q=6p=15.955.985.975.995.935.95p=25.915.955.966.055.955.94p=35.975.955.996.025.876.00p=45.965.995.865.895.775.81p=55.996.046.015.965.755.77p=66.026.036.005.955.795.80It can be seen that the ARMA (5, 5) model will get the minimum BIC value of 5.7496, so ARMA (5, 5) is used as the prediction model in the experiment. At the same time, the order coefficients of the ARMA (5, 5) prediction model areω1 ＝ 0.9189, ω2 ＝ −0.2617, ω3 ＝ 0.3394, ω4 ＝ 0.7065, ω5 ＝ 0.2710, φ1 ＝ −0.1748, φ2 ＝ −0.3256, φ3 ＝ 0.0260, φ4 ＝ −0.5668, and φ5 ＝ −0.7124.Two methods are used to obtain the forecast result of 48 oʼclock on the day, as shown in Figure5. Analysis of the data in the figure shows the following:① At 8:00, the predicted value is greater than the true value, that is,ΔY>0; the ray connected by the true value from 7:30 to 8:00 forms an angle θa with the positive direction of the y-axis; from 8:00 to 8:30, the ray connected by the predicted value forms an angle θ with the positive direction of the y-axis. Obviously, at this time, ΔY>0 and θ>θa. Theorem 1 shows that the load value predicted by the classic ARMA model at 8:30 already has the error convergence trend, and the method in this paper does not need to be corrected, so the prediction results of the two methods are completely equal.② At 19:30, the predicted value is less than the true value; that is,ΔY<0. The rays formed by the true value from 19:00 to 19:30 form an angle θa with the positive direction of the y-axis, from 19:30 to 20:00 The rays formed by the predicted values form an angle θ with the positive direction of the y-axis. Obviously, ΔY<0 and θ>θa at this time. Theorem 1 shows that the load value predicted by the classic ARMA model at 20:00 does not have the error convergence trend, and it needs to be recalculated using the revised algorithm in this paper, so as to obtain the load forecast value at 20:00 with less error.Figure 5
Forecasting result comparison between ARMA method and ARMA + EC method.Similar predictions can be made at other times. The results show that, in the stage where the error of the classic ARMA prediction model is gradually increasing, the method in this paper can effectively evaluate the convergence of the prediction error and make adaptive corrections to ensure the convergence of the prediction error and improve the prediction accuracy.


### 5.2. Comparison Experiment with the BPNN Prediction Model

Taking December 21, 2018 as the forecasting comparison date, the load value after a single forecast horizon of 0.5 h is predicted, and the separate BPNN model (BPNN method) and the error correction method in this paper (BPNN + EC method) based on the BPNN model are, respectively, adopted. We can get the forecast result at 48 oʼclock on the forecast day and make a comparative analysis.In the past 60 days, a total of 2880 historical load values were used as data sources to establish a BPNN prediction model. As a very widely used neural network prediction model, the main job of BPNN modeling is to determine the number of nodes in each layer. In the experiment of this paper, the input layer is designed as 8 nodes, which correspond to the 4 load data closest to the current time and the load data at the same time in the last 4 days; and according to Kolmogorovʼs theorem and empirical knowledge, the hidden layer is selected as 17 intermediate nodes; Bayesian regularization (BR) and Levenberg–Marquardt (LM) rules are used to learn and train the network model, and all input data during training are normalized to [0, 1].Two methods are used to obtain the comparison result between the forecast and actual measurement at 48 oʼclock on the day, as shown in Figure6. The data in the figure is compared and analyzed using the same method as the comparison experiment of the ARMA prediction model. The results show that for the BPNN prediction model, the correction method can also obtain a similar improvement effect; that is, when the prediction value of the BPNN model shows an error divergence trend, this method can effectively evaluate the convergence of the prediction error and perform adaptive correction to ensure the error convergence and improve prediction accuracy.Figure 6
Forecasting result comparison between BPNN method and BPNN + EC method.


### 5.3. Comparative Experiments in Different Prediction Fields

Starting from November 2nd, 2018, 60 consecutive days as the comparison period for forecasting and measuring, the load values after different forecast nodes 0.5 h, 1 h, 2 h, 4 h, and 8 h are predicted, using separate ARMA method and ARMA + EC in this paper The EC method, as well as the separate BPNN method and the BPNN + EC method in this article, obtained the comparison results of 60 days of multipoint prediction and actual measurement for comparative analysis.In order to verify the improvement effect of the algorithm in this paper, the most commonly used statistical indicators Mean Absolute Error (MAE) and Root Mean Square Error (RMSE) are selected to evaluate the prediction accuracy, and two new indicators the Reduction Rate of MAE (RRMAE) and the Reduction Rate of RMSE (RRRMSE) are defined to evaluate the degree of correction of the prediction error of this method. The mathematical expressions of the indicators are as follows:(29)MAE=1N∑t=1Nyt−yat,RMSE=∑t=1Nyt−yat2N,RRMAE=MAEtypical−MAEproposedMAEtypical×100%,RRRMSE=RMSEtypical−RMSEproposedRMSEtypical×100%.N is the total prediction points, yt and yat are the predicted value and the true value at the time of t,MAEtypical, MAEtypical, and RMSEtypical are the MAE and RMSE values of the existing prediction model, and MAEproposed and RMSEproposed are the MAE and RMSE values of the model modified by the method in this paper.The above two experimental methods were used to establish ARMA and BPNN prediction models for different prediction nodes (0.5 h, 1 h, 2 h, 4 h, and 8 h). For different prediction nodes, the ARMA method and the ARMA + EC method are used to obtain the comparison value between the multipoint prediction and the actual measurement for 60 consecutive days from November 2, 2018 to February 1, 2019; BPNN method and BPNN + EC are used methods to obtain the comparison value between the multipoint prediction and the actual measurement for 60 consecutive days in the same period. Calculating the statistical indicators, the results are shown in Tables2 and 3.Table 2
Comparison between ARMA and ARMA + EC methods with multiple forecasting nodes.
Prediction nodesMAE (MW)RRMAE (%)RMSE (MW)RRRMSE (%)ARMAARMA + ECARMAARMA + EC0.522.4618.7216.6825.3321.4815.19123.3919.5616.3526.5522.5914.89225.4021.4115.7328.8124.7114.23429.2224.9814.5033.4829.1213.00836.9932.4012.3943.1638.5810.60Table 3
Comparison between BPNN and BPNN + EC method with multiple forecasting nodes.
Prediction nodes (h)MAE (MW)RRMAE (%)RMSE (MW)RRRMSE (%)BPNNBPNN + ECBPNNARMA + EC0.519.3915.8418.3123.2819.4916.30120.1916.5318.1224.2720.3816.05221.7717.9117.7226.3122.2215.55424.9520.8216.5730.3025.9014.52831.3126.6115.0038.3533.4812.68Experimental results show the following:① Compared with the single ARMA and BPNN prediction models, the prediction error of the method in this paper is significantly reduced. Compared with the ARMA model, the maximum reduction rates of the two indexes of MAE and RMSE reached 16.6759% and 15.1931%. Compared with the BPNN model, the maximum reduction rates of the two indexes of MAE and RMSE reached 18.31149% and 16.2956%.② For different prediction models and prediction nodes, the method is independent of the selected system model, and the improvement of prediction performance is both effective and robust. In all cases, MAE and RMSE both have a performance improvement of at least 10%.③ As the prediction node increases, the load value fluctuation between two adjacent samples increases, resulting in the corresponding increase in the MAE and RMSE indicators of the two algorithms; that is, the prediction accuracy decreases with the increase of the prediction field of view, so this method is more suitable for short-term load forecasting


## 5.1. Comparison Experiment with the ARMA Prediction Model

Based on the actual measurement of a certain residential area in Jinan on December 9, 2018, the load value after a single prediction field of view of 0.5 h was predicted and experimentally measured, using a separate ARMA model (ARMA method) and based on the ARMA model. The error-corrected method (ARMA + EC method) obtained the forecast and measured results at 48 oʼclock on the day, and a comparative analysis was conducted.Based on the historical load value of 120 points totaling 60 hours before the current day as the data source, an ARMA prediction model was established. The classic (p, q)-order ARMA prediction model can be expressed by equation (1).The main work of modeling with ARMA lies in model order determination and coefficient estimation. Firstly, use Autocorrelation Function (ACF) to analyze and remove the fixed trend items in the original data source. Then, the Bayesian Information Criterion (BIC), Yule–Walker, and Newton–Raphson algorithms are used to determine the model order, autoregressive coefficient, and moving average coefficient, respectively. In the experiment, the BIC values corresponding to each order ARMA candidate model are shown in Table1.Table 1
BIC values for ARMA model candidate.
Orderq=1q=2q=3q=4q=5q=6p=15.955.985.975.995.935.95p=25.915.955.966.055.955.94p=35.975.955.996.025.876.00p=45.965.995.865.895.775.81p=55.996.046.015.965.755.77p=66.026.036.005.955.795.80It can be seen that the ARMA (5, 5) model will get the minimum BIC value of 5.7496, so ARMA (5, 5) is used as the prediction model in the experiment. At the same time, the order coefficients of the ARMA (5, 5) prediction model areω1 ＝ 0.9189, ω2 ＝ −0.2617, ω3 ＝ 0.3394, ω4 ＝ 0.7065, ω5 ＝ 0.2710, φ1 ＝ −0.1748, φ2 ＝ −0.3256, φ3 ＝ 0.0260, φ4 ＝ −0.5668, and φ5 ＝ −0.7124.Two methods are used to obtain the forecast result of 48 oʼclock on the day, as shown in Figure5. Analysis of the data in the figure shows the following:① At 8:00, the predicted value is greater than the true value, that is,ΔY>0; the ray connected by the true value from 7:30 to 8:00 forms an angle θa with the positive direction of the y-axis; from 8:00 to 8:30, the ray connected by the predicted value forms an angle θ with the positive direction of the y-axis. Obviously, at this time, ΔY>0 and θ>θa. Theorem 1 shows that the load value predicted by the classic ARMA model at 8:30 already has the error convergence trend, and the method in this paper does not need to be corrected, so the prediction results of the two methods are completely equal.② At 19:30, the predicted value is less than the true value; that is,ΔY<0. The rays formed by the true value from 19:00 to 19:30 form an angle θa with the positive direction of the y-axis, from 19:30 to 20:00 The rays formed by the predicted values form an angle θ with the positive direction of the y-axis. Obviously, ΔY<0 and θ>θa at this time. Theorem 1 shows that the load value predicted by the classic ARMA model at 20:00 does not have the error convergence trend, and it needs to be recalculated using the revised algorithm in this paper, so as to obtain the load forecast value at 20:00 with less error.Figure 5
Forecasting result comparison between ARMA method and ARMA + EC method.Similar predictions can be made at other times. The results show that, in the stage where the error of the classic ARMA prediction model is gradually increasing, the method in this paper can effectively evaluate the convergence of the prediction error and make adaptive corrections to ensure the convergence of the prediction error and improve the prediction accuracy.


## 5.2. Comparison Experiment with the BPNN Prediction Model

Taking December 21, 2018 as the forecasting comparison date, the load value after a single forecast horizon of 0.5 h is predicted, and the separate BPNN model (BPNN method) and the error correction method in this paper (BPNN + EC method) based on the BPNN model are, respectively, adopted. We can get the forecast result at 48 oʼclock on the forecast day and make a comparative analysis.In the past 60 days, a total of 2880 historical load values were used as data sources to establish a BPNN prediction model. As a very widely used neural network prediction model, the main job of BPNN modeling is to determine the number of nodes in each layer. In the experiment of this paper, the input layer is designed as 8 nodes, which correspond to the 4 load data closest to the current time and the load data at the same time in the last 4 days; and according to Kolmogorovʼs theorem and empirical knowledge, the hidden layer is selected as 17 intermediate nodes; Bayesian regularization (BR) and Levenberg–Marquardt (LM) rules are used to learn and train the network model, and all input data during training are normalized to [0, 1].Two methods are used to obtain the comparison result between the forecast and actual measurement at 48 oʼclock on the day, as shown in Figure6. The data in the figure is compared and analyzed using the same method as the comparison experiment of the ARMA prediction model. The results show that for the BPNN prediction model, the correction method can also obtain a similar improvement effect; that is, when the prediction value of the BPNN model shows an error divergence trend, this method can effectively evaluate the convergence of the prediction error and perform adaptive correction to ensure the error convergence and improve prediction accuracy.Figure 6
Forecasting result comparison between BPNN method and BPNN + EC method.


## 5.3. Comparative Experiments in Different Prediction Fields

Starting from November 2nd, 2018, 60 consecutive days as the comparison period for forecasting and measuring, the load values after different forecast nodes 0.5 h, 1 h, 2 h, 4 h, and 8 h are predicted, using separate ARMA method and ARMA + EC in this paper The EC method, as well as the separate BPNN method and the BPNN + EC method in this article, obtained the comparison results of 60 days of multipoint prediction and actual measurement for comparative analysis.In order to verify the improvement effect of the algorithm in this paper, the most commonly used statistical indicators Mean Absolute Error (MAE) and Root Mean Square Error (RMSE) are selected to evaluate the prediction accuracy, and two new indicators the Reduction Rate of MAE (RRMAE) and the Reduction Rate of RMSE (RRRMSE) are defined to evaluate the degree of correction of the prediction error of this method. The mathematical expressions of the indicators are as follows:(29)MAE=1N∑t=1Nyt−yat,RMSE=∑t=1Nyt−yat2N,RRMAE=MAEtypical−MAEproposedMAEtypical×100%,RRRMSE=RMSEtypical−RMSEproposedRMSEtypical×100%.N is the total prediction points, yt and yat are the predicted value and the true value at the time of t,MAEtypical, MAEtypical, and RMSEtypical are the MAE and RMSE values of the existing prediction model, and MAEproposed and RMSEproposed are the MAE and RMSE values of the model modified by the method in this paper.The above two experimental methods were used to establish ARMA and BPNN prediction models for different prediction nodes (0.5 h, 1 h, 2 h, 4 h, and 8 h). For different prediction nodes, the ARMA method and the ARMA + EC method are used to obtain the comparison value between the multipoint prediction and the actual measurement for 60 consecutive days from November 2, 2018 to February 1, 2019; BPNN method and BPNN + EC are used methods to obtain the comparison value between the multipoint prediction and the actual measurement for 60 consecutive days in the same period. Calculating the statistical indicators, the results are shown in Tables2 and 3.Table 2
Comparison between ARMA and ARMA + EC methods with multiple forecasting nodes.
Prediction nodesMAE (MW)RRMAE (%)RMSE (MW)RRRMSE (%)ARMAARMA + ECARMAARMA + EC0.522.4618.7216.6825.3321.4815.19123.3919.5616.3526.5522.5914.89225.4021.4115.7328.8124.7114.23429.2224.9814.5033.4829.1213.00836.9932.4012.3943.1638.5810.60Table 3
Comparison between BPNN and BPNN + EC method with multiple forecasting nodes.
Prediction nodes (h)MAE (MW)RRMAE (%)RMSE (MW)RRRMSE (%)BPNNBPNN + ECBPNNARMA + EC0.519.3915.8418.3123.2819.4916.30120.1916.5318.1224.2720.3816.05221.7717.9117.7226.3122.2215.55424.9520.8216.5730.3025.9014.52831.3126.6115.0038.3533.4812.68Experimental results show the following:① Compared with the single ARMA and BPNN prediction models, the prediction error of the method in this paper is significantly reduced. Compared with the ARMA model, the maximum reduction rates of the two indexes of MAE and RMSE reached 16.6759% and 15.1931%. Compared with the BPNN model, the maximum reduction rates of the two indexes of MAE and RMSE reached 18.31149% and 16.2956%.② For different prediction models and prediction nodes, the method is independent of the selected system model, and the improvement of prediction performance is both effective and robust. In all cases, MAE and RMSE both have a performance improvement of at least 10%.③ As the prediction node increases, the load value fluctuation between two adjacent samples increases, resulting in the corresponding increase in the MAE and RMSE indicators of the two algorithms; that is, the prediction accuracy decreases with the increase of the prediction field of view, so this method is more suitable for short-term load forecasting


## 6. Summary

Starting from the historical load sequence to study how to reduce the prediction error is an important means to improve the accuracy of heating load prediction. This paper proposes a correction method for load prediction error based on trajectory tracking. This method is completely based on the theory of trajectory tracking stability. There are no special requirements for the prediction model itself; that is, it has nothing to do with the system model, which can theoretically ensure error convergence, and simulation experiments also show its effectiveness.



---

*Source: 5919238-2020-11-17.xml*