1. Introduction

2. Materials and Methods

2.2. Analysis of Multi-Energy Load Correlation Characteristics Based on the Copula Method

To improve the final prediction accuracy of the short-term multi-energy load forecasting, several aspects must be considered, including the characteristics of the internal load represented by the historical development laws of the multi-energy loads, coupling conversion relationship of the load, characteristics of the external load, and meteorological aspects. Compared to the number of studies on the historical development law of the multi-energy loads and correlation between the multi-energy loads and meteorological factors, only few studies address the characteristics of the complex and flexible coupling conversion of the multi-energy loads. Hence, more accurate characterization methods are needed to effectively optimize the overall performance of the proposed short-term multi-energy load forecasting model. In this study, copula theory is used to model and analyze the correlation characteristics of the non-linear coupling conversion of the time series of the multi-energy loads.

2.2.1. Definition of the Copula Function

The copula theory was developed to solve a joint distribution problem of random variables when their marginal distributions are known. The formulation of Sklar’s theorem introduced the concept of the copula function. Further, models based on the copula function have been widely used in finance and economics, new energy output characteristic analysis, and other fields owing to their ability to explain complex nonlinearities among variables [21,22].

Sklar’s theorem proves that the joint distribution function of multiple n-dimensional variables can be constructed by combining the marginal distribution function of these variables and associated copula function, which represents the complex correlation among variables. Considering an n-dimensional variable $x_1,x_2,\cdots ,x_n$with a marginal distribution function in each dimension $F_1\left(x_1\right),F_2\left(x_2\right),\cdots ,F_n\left(x_n\right)$, the joint distribution and density functions of the n-dimensional variable can be expressed as follows:

$$\mathit{F}(x_1,x_2,\cdots ,x_n)=\mathit{C}[F_1(x_1),F_2(x_2),\cdots ,F_n(x_n)],$$

(10)

$$\mathit{f}(x_1,x_2,\cdots ,x_n)=c(f(x_1),f(x_2),\cdots ,f(x_n))\cdot {\displaystyle \prod _{i=1}^n{f}_i(x_i)}$$

(11)

where $F$ is the joint distribution function of the n-dimensional variable, $C$ the copula distribution function representing the complex correlation between n-dimensional variables, $F_i\left(x_i\right)$ the density function of the variable in the $i^{th}$ dimension, $f$ the joint density function of the n-dimensional variable, and $c$ the copula density function.

Based on equations (10) and (11), the copula density function can be obtained by taking the derivative of the copula distribution function, as expressed in equation (12).

$$c(f(x_1),f(x_2),\cdots ,f(x_n))=\frac{{\partial }^n\mathit{C}[F_1(x_1),F_2(x_2),\cdots ,F_n(x_n)]}{\partial F_1(x_1)\partial F_2(x_2)\cdots \partial F_n(x_n)}$$

(12)

By calculating the copula distribution and density functions, the complex correlation between the multivariate random variables can be accurately described. Typical static copula distribution functions include N-, T-, Gumbel, Clayton, and Frank copula functions. Additionally, the dynamic N-, T-, Clayton, and SJC copula functions can be used to model dynamic copula distribution functions. Hence, the correlation characteristics of the random variables can be analyzed using different copula distribution functions, which can reflect different perspectives.

2.2.2. Correlation Analysis Based on Copula Functions

Generally, it is difficult to obtain a clear marginal distribution function because of the complexity of the multi-energy load time-series data. Hence, we applied the maximum likelihood estimation based on nonparametric kernel density (MLK) to estimate the correlation coefficient of the copula function [23,24]. The MLK method is not limited by the exact expression of the marginal distribution function. Instead, it uses the nonparametric kernel density estimation function of the analyzed variables.

For the normalized variable sequences Lp and Lc, the corresponding nonparametric kernel density estimation is conducted using the following equations:

$$f_P(x_1)=\frac{1}{T}{\displaystyle \sum _{t=1}^T{K}_w(x_1-L_P)},$$

(13)

$$f_C(x_2)=\frac{1}{T}{\displaystyle \sum _{t=1}^T{K}_w(x_2-L_C)},$$

(14)

where $f_P\left(x_1\right)$ and $f_C\left(x_2\right)$ represent the probability density functions of $L_P$ and $L_C$, respectively. $K_w$ is the kernel function and T the size of the sequence of variables.

Using the MLK method, the static correlation coefficients can be obtained by substituting the marginal distributions with the probability density functions $f_P\left(x_1\right)$ and $f_C\left(x_2\right)$ in the likelihood function expressed in equation (15). Further, its extreme points can be calculated using equation (16).

$$S(\theta )={\displaystyle \sum \ln c[f_P(x_1),f_C(x_2)]}$$

(15)

$$\widehat{\theta }=\arg \max S(\theta )$$

(16)

However, to calculate the dynamic time-varying correlation coefficient series, a dynamic copula function must be considered. Further, the likelihood function must be obtained using the parameters of the dynamic distributions associated to the variables and corresponding evolution equations. Hence, we define the dynamic N-copula function using the dynamic distribution parameter ${\rho }_{N,t}$ and dynamic T-copula function with parameters ${\rho }_{T,t}$ and $k_t$ degrees of freedom. The correlation coefficient matrices of the dynamic N- and T-copula functions based on the DCC(1,1) decomposition can be expressed as follows:

$${\mathit{R}}_t={({\mathit{Q}}_t^{*})}^{-1/2}\cdot {\mathit{Q}}_t\cdot {({\mathit{Q}}_t^{*})}^{-1/2},$$

(17)

where $Q_t^{*}=diag{Q}_t$ and its evolution equation is expressed as follows

$${\mathit{Q}}_t=\mathit{R}(1-\alpha -\beta )+\alpha ({\epsilon }_{t-1}{{\epsilon }^{\prime }}_{t-1})+\beta {\mathit{Q}}_{t-1},$$

(18)

where $\alpha$ and $\beta$ represent the estimated evolution parameters, which satisfy the constraints $0<\alpha <1$,$0<\beta <1$ and $0<\alpha +\beta <1$, respectively. ${\epsilon }_t$ is the pseudo-inverse of the threshold distribution function.

The evolution equation of the dynamic Clayton copula function is defined as follows:

$${\theta }_{C,t}=\Lambda (\omega +\beta {\theta }_{t-1}+\alpha \cdot \frac{1}{10}{\displaystyle \sum _{j=1}^{10}(|L_{P,t-j}-L_{C,t-j}|}),$$

(19)

where $\omega$, $\alpha$, and $\beta$ represent the estimated evolution parameters. $\Lambda (x)=\frac{(1-e^{-x})}{(1+e^{-x})}$ is the restriction function.

The evolution equation of the dynamic SJC copula function is expressed as follows:

$${\tau }_t^U=\tilde{\Lambda }({\omega }_U+{\beta }_U{\tau }_{t-1}^U+{\alpha }_U\cdot \frac{1}{10}{\displaystyle \sum _{j=1}^{10}|L_{P,t-j}-L_{C,t-j}|}),$$

(20)

$${\tau }_t^L=\tilde{\Lambda }({\omega }_L+{\beta }_L{\tau }_{t-1}^L+{\alpha }_L\cdot \frac{1}{10}{\displaystyle \sum _{j=1}^{10}|L_{P,t-j}-L_{C,t-j}|}),$$

(21)

where ${\tau }_t^U$ is the upper-tail dependence coefficient, ${\omega }_U$, ${\alpha }_U$, and ${\beta }_U$ the estimated parameters of the upper-tail evolution, ${\tau }_t^L$ the lower-tail dependence coefficient, and ${\omega }_L$, ${\alpha }_L$, and ${\beta }_L$ the estimated parameters of the lower-tail evolution. The restriction function $\tilde{\Lambda }\left(x\right)$ satisfies $\tilde{\Lambda }\left(x\right)=(1+e^{-x}{)}^{-1}$.

In this study, the Akaike information criterion (AIC) was used to evaluate the adaptability of the different copula models presented above. The AIC index is calculated as follows

$$\mathrm{AIC}=2k-2\ln S,$$

(22)

where k represents the number of parameters of the copula function and S the associated maximum likelihood estimation.

2.2.3. Analysis of Multi-energy Load Characteristics Based on Copula Methods

The steps of the analysis of the multi-energy loads using the copula method are summarized in Table 1.

Copula functions can accurately describe the complex non-linear coupling relationship characterizing the multi-energy loads, which can considerably contribute to the improvement of the overall performance of the proposed short-term multi-energy load prediction model. As shown in Figure 1, the copula function-based feature analysis process proposed in this study includes the following steps:

• Multi-energy load data are normalized in the [0,1] interval to ensure data uniformity.
• The kernel density estimation function is calculated using the MLK method, which determines the marginal density function of the variable sequence.

To obtain static correlation coefficients, the marginal density copula functions are used to calculate the extreme points of the likelihood function.

• To obtain dynamic correlation coefficients, the dynamic copula distributions are used to construct the likelihood function considering the corresponding evolution equation parameters ($\omega$, $\alpha$, and $\beta$).
• Once the maximum likelihood estimates and the corresponding evolution equation parameters are obtained, they are substituted into the evolution equation parameters to calculate the required time-varying cross-correlation coefficients.
• Simultaneously, the selected copula functions are optimized based on the maximum likelihood estimate, and later the optimal copula model is obtained by comparing the corresponding AIC indexes.

Figure 1. Schematic of the copula feature analysis process.

Figure 1. Schematic of the copula feature analysis process.

3. Results and Discussion

4. Conclusions

References

1. Nti, I.K.; Teimeh, M.; Nyarko-Boateng, O.; Adekoya, A.F. Electricity load forecasting: A systematic review. J Electr Syst Inf Technol. 2020, 7, 1–19. [CrossRef]
2. Hammad, M.A.; Jereb, B.; Rosi, B.; Dragan, D. Methods and models for electric load forecasting: A comprehensive review. Logist Supply Chain Sustain Glob Chall. 2020, 11, 51–76. [CrossRef]
3. Kim, J.H.; Seong, N.C.; Choi, W. Cooling load forecasting via predictive optimization of a non-linear autoregressive exogenous (NARX) neural network model. Sustainability. 2019, 11, 6535. [CrossRef]
4. Fu, G. Deep belief network based ensemble approach for cooling load forecasting of air-conditioning system. Energy. 2018, 148, 269–282. [CrossRef]
5. Li, F.; Zheng, H.; Li, X.; Yang, F. Day-ahead city natural gas load forecasting based on decomposition-fusion technique and diversified ensemble learning model. Appl Energy. 2021, 303, 117623. [CrossRef]
6. Lu, H.; Azimi, M.; Iseley, T. Short-term load forecasting of urban gas using a hybrid model based on improved fruit fly optimization algorithm and support vector machine. Energy Rep. 2019, 5, 666–677. [CrossRef]
7. Zhang, Z.; Liu, Y.; Cao, L.; Si, H. A forecasting method of district heat load based on improved wavelet neural network. J Energy Resour Technol. 2020, 142, 1–29. [CrossRef]
8. Kannari, L.; Kiljander, J.; Piira, K.; Piippo, J.; Koponen, P. Building heat demand forecasting by training a common machine learning model with physics-based simulator. Forecasting. 2021, 3, 290–302. [CrossRef]
9. Zhu, J.; Dong, H.; Zheng, W.; Li, S.; Huang, Y.; Xi, L. Review and prospect of data-driven techniques for load forecasting in integrated energy systems. Appl Energy. 2022, 321, 119269. [CrossRef]
10. Liu, X.; Zhang, Z.; Song, Z. A comparative study of the data-driven day-ahead hourly provincial load forecasting methods: From classical data mining to deep learning. Renew Sustain Energy Rev. 2020, 119, 109632. [CrossRef]
11. Chen, J.; Li, T.; Zou, Y.; Wang, G.; Ye, H.; Lv, F. An ensemble feature selection method for short-term electrical load forecasting. In Proceedings of the 2019 IEEE 3rd Conference on Energy Internet and Energy System Integration, 2019; Vol. EI2; pp. 1429–1432. [CrossRef]
12. Salami, M.; Sobhani, F.M.; Ghazizadeh, M.S. A hybrid short-term load forecasting model developed by factor and feature selection algorithms using improved grasshopper optimization algorithm and principal component analysis. Electr Eng. 2020, 102, 437–460. [CrossRef]
13. Yahaya, A.S.; Javaid, N.; Latif, K.; Rehman, A. An enhanced very short-term load forecasting scheme based on activation function. In International Conference on Computer and Information Sciences (ICCIS); IEEE Publications, 2019; pp. 1–6. [CrossRef]
14. Kalkstein, L.S.; Tan, G.; Skindlov, J.A. An evaluation of three clustering procedures for use in synoptic climatological classification. J Climate Appl Meteor. 1987, 26, 717–730. [CrossRef]
15. Douglas, D.H.; Peucker, T.K. Algorithms for the reduction of the number of points required to represent a digitized line or its caricature. Cartographica. 1973, 10, 112–122. [CrossRef]
16. Visvalingam, M.; Whyatt, J.D. The Douglas-Peucker algorithm for line simplification: Re-evaluation through visualization. In Comput Graph Forum. 1990, 9, 213–225. [CrossRef]
17. Hershberger, J.; Snoeyink, J. Speeding up the Douglas-Peucker line simplification algorithm. In Proceedings of the 5th International Syrnposium on Spatial Data Handling, 1992; pp. 134–143.
18. Hershberger, J.; Snoeyink, J. An O(n log n) implementation of the Douglas-Peucker algorithm for line simplification. In Proceeding of the 10th Annual Symposium on Computational Geometry, 1994; pp. 383–384.
19. Wang, J.; Liu, Z. Implementation and analysis of curve vector data compression algorithm. Geom Spat Inf Technol. 2006, 122–124.
20. van der Vlist, R.; Taal, C.; Heusdens, R. Tracking recurring patterns in time series using dynamic time warping. In 27th European Signal Processing Conference (EUSIPCO), 2019; pp. 1–5. [CrossRef]
21. Luo, M.; Tian, Y.; Li, C.; Chen, Y. The non-linear test and empirically study on the financial crisis contagion based on copula method. In 2nd International Conference on Artificial Intelligence, Management Science and Electronic Commerce (AIMSEC), 2011; pp. 2587–2591.
22. Yundai, X.; Yue, Y. Analysis of aggregated wind power dependence based on optimal vine copula. In IEEE Innovative Smart Grid Technologies, Asia, I., Ed.; Asia Publishing, 2019; pp. 1788–1792. [CrossRef]
23. Patton, A.J. Modelling asymmetric exchange rate dependence. Int Economic Rev. 2006, 47, 527–556. [CrossRef]
24. Patton, A.J. Estimation of multivariate models for time series of possibly different lengths. J Appl Econ. 2006, 21, 147–173. [CrossRef]
25. Leontaritis, I.J.; Billings, S.A. Input-output parametric models for non-linear systems part I: Deterministic non-linear systems. Int J Control. 1985, 41, 303–328. [CrossRef]
26. Leontaritis, I.J.; Billings, S.A. Input-output parametric models for non-linear systems part II: Stochastic non-linear systems. Int J Control. 1985, 41, 329–344. [CrossRef]
27. Wang, H.; Yeung, D.Y. Towards Bayesian deep learning: A framework and some existing methods. IEEE Trans Knowl Data Eng. 2016, 28, 3395–3408. [CrossRef]
28. Mackay, D.J.C. Bayesian interpolation. Neural Comput. 1992, 4, 415–447. [CrossRef]