Mid-Term Electricity Market Clearing Price Forecasting with Sparse Data: A Case in Newly-Reformed Yunnan Electricity Market

Abstract

For the power systems, for which few data are available for mid-term electricity market clearing price (MCP) forecasting at the early stage of market reform, a novel grey prediction model (defined as interval GM(0, N) model) is proposed in this paper. Over the traditional GM(0, N) model, three major improvements of the proposed model are: (i) the lower and upper bounds are firstly identified to give an interval estimation of the forecasting value; (ii) a novel whitenization method is then established to determine the definite forecasting value from the forecasting interval; and (iii) the model parameters are identified by an improved particle swarm optimization (PSO) instead of the least square method (LSM) for the limitation of LSM. Finally, a newly-reformed electricity market in Yunnan province of China is studied, and input variables are contrapuntally selected. The accuracy of the proposed model is validated by observed data. Compared with the multiple linear regression (MLR) model, the traditional GM(0, N) model and the artificial neural network (ANN) model, the proposed model gives a better performance and its superiority is further ensured by the use of the modified Diebold–Mariano (MDM) test, suggesting that it is suitable for mid-term electricity MCP forecasting in a data-sparse electricity market.

1. Introduction

Electricity price is one of the most important parameters of deregulated and competitive electricity markets; wherein, the electricity market clearing price (MCP) is the final outcome of the market bid price and generally determined by the aggregated demand and supply bid curves, which exists when an electricity market is clear of shortage and surplus. On 15 March 2015, the Communist Party of China (CPC) Central Committee and State Council’s opinions on further deepening the reforms of the electric power system [1] marked the beginning of new electricity market reform in China, after Investment Decentralization in 1986, Unbundling of Government Admin and Business Operation in 1997 and Unbundling of Generation and Transmission in 2002 [2]. Unlike mature electricity markets in other countries, to guarantee the security and stability operation at the early stage of reform, the Chinese electricity market mainly comprises the mid- and long-term exchanges so far. Therefore, the mid-term electricity MCP forecasting has become one of the most concerning issues, as it is essential for resource reallocation, maintenance scheduling, financial risk reducing, bilateral contracting, budgeting and planning purposes [3].

Over the last 15–20 years, electricity price forecasting (EPF) has been carried out with various methods, bringing a variety of success. Till now, there is no consensus criterion on the classification of these models [4,5]. To our knowledge, they can be mainly classified into the following categories: (1) equilibrium and simulation models, such as the Nash–Cournot model [6,7,8], production-cost models (PCM) [9], strategic PCM (SPCM) [10], and so on; (2) statistical methods, such as multiple regression [11], exponential smoothing [12,13], time series [14,15,16], wavelet transform (WT) [17,18,19], AR-type models ( such as auto-regressive moving average (ARMA) [16], auto-regressive integrated moving average (ARIMA) [20,21], seasonal ARIMA (SARIMA) [22]), and so on; (3) artificial intelligence models, such as artificial neural network (ANN) [23,24], support vector machine (SVM) [25,26] and data mining models [27,28]; and (4) hybrid models, which combine two or more single methods or models described above [29,30,31]. According to these research works, most works are focused on short-term EPF, commonly known as day-ahead EPF. In contrast, not much has been performed in the mid-term horizon [3,32,33,34,35]. Compared to short-term forecasting, mid-term forecasting is much more challenging. This is because the prediction horizon is much longer, and the price may not be contiguous to the immediate past periods, which means the trend information of the immediate past is not as valuable as for short-term forecasting. In addition, the influence factors of mid-term EPF are various and, worse, some are unavailable. When it comes to the newly-reformed electricity market in China, being at the early stage, another bigger challenge is how to build accurate models with relatively few available samples, which is beyond the capability of most forecasting models above. Appropriate prediction models need to be studied.

In systems theory, if the description of a system is completely known, it is known as a white system. In contrast, a black-box system means that it is completely unknown. When in an intermediate state, it is called a grey system [36]. The grey system theory (GST) was firstly proposed by Deng in 1982 [37], and several grey prediction models were created later. For the grey model, it is characterized by its prominent capability in modeling with small samples and poor information, which commonly exist in the natural world. Thus, it has been successfully employed in various fields and demonstrated satisfactory results, such as social [38,39], economic [40,41], agricultural [42], industrial [43], etc. The electricity MCP is a synthetic variable of the market mechanism with many uncertainties, obviously in accordance with the characteristics of GST, so the electricity MCP forecasting is adapted to be solved through grey prediction models. Furthermore, as the mid-term electricity MCP is not closely related to the immediate past periods and affected by multiple factors, a widely-used static multi-variable grey prediction model in GST named GM(0, N) is concerned in this paper. However, to directly forecast the mid-term electricity MCP, there are some difficulties or disadvantages in the traditional GM(0, N) model, including:

• The input variables of the model would affect the forecasting result directly. Thus, the selection of appropriate influence factors of electricity MCP is extremely important, but also tough, as it varies from different electricity markets.
• The mid-term electricity MCP is obviously volatile in nature and out of a monotonous trend. The forecasting accuracy is unsatisfactory when the value of the next forecasting point is away from the fitting model constructed only with known data.
• The least square method (LSM) used in the traditional GM(0, N) model to identify parameters is effective on condition of the existence of the inverse of matrix $B^{T}B$, but this does not work when $B^{T}B$ is a singular matrix in some cases.

To overcome these problems and achieve a successful implementation for mid-term electricity MCP forecasting with few available data, a novel grey prediction model defined as the interval GM(0, N) model is proposed in this paper. The main contributions of this work can be summarized as follows:

• Based on depth analysis of the newly-reformed electricity market, the influence factors of electricity MCP are studied, and input variables are carefully selected according to three aspects: supply factors, demand factors and supplemented factors.
• In the proposed interval GM(0, N) model, two improved GM(0, N) models are included, respectively estimating the upper and lower bounds of the forecasting value. Firstly, to reduce randomness and increase smoothness, all input sequences (not containing values of the next forecasting point) including the MCP sequence and factor sequences are ranked in accordance with the ascending order of the MCP sequence. Then, one of the factor sequences is selected as the benchmark ranking sequence according to its ranking order and correlation with the MCP sequence. The position of the forecasting point in the ranked MCP sequence thus can be determined by sorting the benchmark ranking sequence, which contains the factor value of the forecasting point. Finally, two neighboring (upper and lower) points of the forecasting point in the ranked MCP sequence are regarded as two virtual values to construct two new MCP sequences, which are respectively used as characteristic sequences for two improved GM(0, N) models, obtaining the forecasting interval. In the two improved GM(0, N) models, the input sequences used for building the model include not only known data, but also the virtual MCP and predicted influence factors of the forecasting point.
• Based on the forecasting interval, a novel whitenization method considering correlations between electricity MCP and influence factors is established to determine the definite forecasting value.
• The parameters of the GM(0, N) model are identified by an improved particle swarm optimization (PSO) instead of LSM.
• The performance of the proposed model has been validated by applying it to the newly-reformed Yunnan electricity market. Further comparisons between the proposed model and other models, including the multiple linear regression (MLR) model, the traditional GM(0, N) model and the artificial neural network (ANN) model, are carefully discussed and also evaluated by using the modified Diebold–Mariano (MDM) test. The results indicate that the proposed model is an effective means of mid-term electricity MCP forecasting with sparse data.

The remainder of this paper is organized as follows. In Section 2, the method and theory are introduced. Section 3 is devoted to the detailed presentation of the proposed model. In Section 4, electricity market conditions are analyzed, and input variables are selected. In Section 5, the case study and forecasting results are carefully discussed. Finally, Section 6 outlines the main conclusions.

2. Method and Theory

2.1. Principle of the Traditional GM(0, N) Model

Set $X_1^{(0)}=\{x_1^{(0)}(1),x_1^{(0)}(2),\cdots ,x_1^{(0)}(n)\}$ as the raw characteristic sequence (here, the electricity MCP sequence); $n$ is the data length. Set $X_i^{(0)}=\{x_i^{(0)}(1),x_i^{(0)}(2),\cdots ,x_i^{(0)}(n)\},i\in I,I=\{2,3,\cdots ,N\}$ as the relevant factor sequences (here, the influence factors of mid-term electricity MCP); $N$ is the number of all sequences, and the total number of factor sequences is $N-1$. Therefore, the AGO (accumulated generation operation) sequence $X_i^{(1)}=\{x_i^{(1)}(1),x_i^{(1)}(2),\cdots ,x_i^{(1)}(n)\}$ is firstly generated by Equation (1).

$$x_i^{(1)}(k)={\displaystyle \sum _{m=1}^k{x}_i^{(0)}(m),(k=1,2,\cdots ,n;i=1,2,\cdots ,N)}$$

(1)

Then, call formula Equation (2):

$$x_1^{(1)}(k)={\displaystyle \sum _{i=2}^N{b}_i{x}_i^{(1)}(k)+a,(k=1,2,\cdots ,n)}$$

(2)

as the traditional GM(0, N) model. Where, $b_i$ and $a$ are identified parameters, which can be solved by the least square method (LSM), if the inverse of matrix $B^{T}B$ exists. That is:

$$P_{0N}={\left[b_2,b_3,\cdots ,b_N,a\right]}^T={\left(B^{T}B\right)}^{-1}{B}^T{y}_N$$

(3)

where, $B=\left[\begin{array}{ccccc}{x}_2^{(1)}(1)& x_3^{(1)}(1)& \cdots & x_N^{(1)}(1)& 1\\ x_2^{(1)}(2)& x_3^{(1)}(2)& \cdots & x_N^{(1)}(2)& 1\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ x_2^{(1)}(n)& x_3^{(1)}(n)& \cdots & x_N^{(1)}(n)& 1\end{array}\right]$, $y_N={[x_1^{(1)}(1),x_1^{(1)}(2),\cdots ,x_1^{(1)}(n),1]}^T$.

Thus, we can get the next forecasting value for the characteristic sequence by Equation (4).

$$\{\begin{array}{l}{\widehat{x}}_1^{(0)}(1)=x_1^{(0)}(1)\\ {\widehat{x}}_1^{(0)}(k+1)={\widehat{x}}_1^{(1)}(k+1)-{\widehat{x}}_1^{(1)}(k),k=1,2,\cdots ,n\end{array}$$

(4)

where ${\widehat{x}}_1^{(1)}\left(k+1\right)$ is the forecasting value of $x_1^{(1)}(k+1)$; ${\widehat{x}}_1^{(0)}\left(k+1\right)$ is the forecasting value of $x_1^{(0)}(k+1)$.

2.2. Limitation and Requirement for Implementation

2.2.1. Limitation of Least Square Method (LSM) in the Traditional GM(0, N) Model

In the traditional GM(0, N) model, parameters $b_i$ and $a$ are identified by LSM, shown in Equation (3). However, this is conditioned on the existence of the inverse of matrix $B^{T}B$.

Proof: The determinant of $B^{T}B$ is $Det(B^{T}B)$; the inverse matrix of $B^{T}B$ exists only when:

$$Det(B^{T}B)\ne 0$$

(5)

that is when the $B$ is a column full rank matrix. That is when, for any constant $u\in (-\infty ,+\infty )$,

$$X_i^{(1)}+u{X}_j^{(1)}\ne \overrightarrow{0},(\forall i,j\in I,i\ne j)$$

(6)

otherwise, the inverse of matrix $B^{T}B$ does not exist.

From the equations above, it can be observed that the limitation of LSM depends on the input dataset. However, Equation (6) does not hold for all cases, especially with the increasing length of available data. Therefore, the LSM cannot be adapted to identify parameters $b_i$ and $a$ for the GM(0, N) model in some cases. To avoid any mistakes caused by this limitation and to ensure the general applicability of the model, the parameter identification method should be improved.

2.2.2. Requirements for Input Variables

Although the great feature of the GM(0, N) model is allowing the randomness of data sequences and modeling with small samples and poor conditions, the volatile nature of the electricity price has a negative effect on forecasting results. Despite the fact that AGO has been applied to convert the sequences to a monotonic increasing trend, the modeling accuracy is still limited by directly using raw known data. For better accuracy, the input data need to be pre-processed. Furthermore, the related influence factors of electricity MCP are various, and input variables should be carefully selected, as they also have an impact on the forecasting result.

2.3. Performance Evaluation

2.3.1. Checking Method of Grey Prediction Models

To evaluate the fitting precision of the proposed grey prediction model, two main parameters in the after-test residue checking method, posterior-error ($C$) and micro-error-probability ($P$), are adopted, respectively defined as:

$$C=\frac{S_2}{S_1}$$

(7)

$$P=probability\{|\Delta x_1^{(0)}(k)-\Delta {\overline{x}}_1^{(0)}(k)|<0.6745S_1\}$$

(8)

where:

$$\{\begin{array}{l}{S}_1=\sqrt{\frac{1}{n-1}{\displaystyle \sum _{k=1}^n{\left[x_1^{(0)}(k)-{\overline{x}}_1^{(0)}(k)\right]}^2}}\hfill \\ S_2=\sqrt{\frac{1}{n-1}{\displaystyle \sum _{k=1}^n{\left[\Delta x_1^{(0)}(k)-\Delta {\overline{x}}_1^{(0)}(k)\right]}^2}}\hfill \\ {\overline{x}}_1^{(0)}(k)=\frac{1}{n}{\displaystyle \sum _{k=1}^n{x}_1^{(0)}(k)}\hfill \\ \Delta x_1^{(0)}(k)=x_1^{(0)}(k)-{\widehat{x}}_1^{(0)}(k)\hfill \\ \Delta {\overline{x}}_1^{(0)}(k)=\frac{1}{n}{\displaystyle \sum _{k=1}^n\Delta x_1^{(0)}(k)}\hfill \end{array}$$

(9)

The fitting precision grade is shown in

Table 1. Reference for the fitting precision grade.

Parameters	Fitting Precision Grade
	Good	Qualified	Just	Unqualified
$C$	<0.35	0.35–0.50	0.50–0.65	≥0.65
$P$	>0.95	0.80–0.95	0.70–0.80	≤0.70

.

2.3.2. Performance Evaluation of Forecasting Models

For evaluating electricity MCP forecasting values, mean absolute error (MAE), mean squares error (MSE) and mean absolute percentage error (MAPE) are the most widely-used measurements, respectively shown in Equations (10)–(12). The smaller the MAE, MSE and MAPE, the better forecasting performance the model shows.

$$MAE=\frac{1}{n}{\displaystyle \sum _{k=1}^n|x_1^{(0)}(k)-{\widehat{x}}_1^{(0)}(k)|}$$

(10)

$$MSE=\frac{1}{n}{\displaystyle \sum _{k=1}^n{\left[x_1^{(0)}(k)-{\widehat{x}}_1^{(0)}(k)\right]}^2}$$

(11)

$$MAPE=\frac{1}{n}{\displaystyle \sum _{k=1}^n\left|\frac{x_1^{(0)}(k)-{\widehat{x}}_1^{(0)}(k)}{x_1^{(0)}(k)}\right|}\times 100\%$$

(12)

To further judge the forecasting performances of different forecasting models from a statistical point of view, a statistical evaluation method named modified Diebold–Mariano (MDM) test [44] is also proposed. The MDM test is an extension of the Diebold–Mariano test [45] and has been widely used in forecasting research. The MDM test statistic for the h-step ahead forecast is specified as:

$$MDM={\left[\frac{n+1-2h+n^{-1}h(h-1)}{n}\right]}^{\frac{1}{2}}{\left[V(\overline{d})\right]}^{-\frac{1}{2}}\left[\overline{d}\right]$$

(13)

where $V(\overline{d})=\left[n^{-1}\left({\gamma }_0+2{\displaystyle \sum _{k=1}^{h-1}{\gamma }_k}\right)\right]$ is the variance of $\overline{d}$; ${\gamma }_k=n^{-1}{\displaystyle \sum _{t=k+1}^n\left(d_t-\overline{d}\right)}(d_{t-k}-\overline{d})$ is the $k$-th autocovariance of $d_t$; $\overline{d}=n^{-1}{\displaystyle \sum _{t=1}^n{d}_t}$ is the mean of $d_t$; $d_t=L\left(e_{1,t}\right)-L\left(e_{2,t}\right),t=1,2,3,\cdots ,n$ is the loss differentials; $L\left(e_{i,t}\right)$ is the loss function of forecast error $e_{i,t},i=1,2$; $e_{i,t}=x_t-{\widehat{x}}_{i,t}$, $x_t$ denotes the actual data series, and ${\widehat{x}}_{i,t}$ denotes the forecasting series of forecasting model $i$. In our study, the popular loss functions used in power systems, i.e., the MAE and MSE loss functions, are adopted. Under the null hypothesis of the equal forecasting performances of forecasting models, the MDM test statistic follows a t-distribution with (n − 1) degrees of freedom, so that tests can be carried out. As the detailed principle of the MDM test can be easily found in related papers, we will not repeat it here.

3. Novel Interval GM(0, N) Model

Sorting is a common way to get a descending (or ascending) sequence. However, the position of the forecasting point in the MCP sequence cannot be obtained, as its actual value is unknown. For this, one of the factor sequences is selected as the benchmark ranking sequence to determine the position of the forecasting point in the MCP sequence. Then, the actual value of the forecasting point is respectively replaced by its two neighboring (upper and lower) values in the MCP sequence. Two new MCP sequences thus are constructed, which are respectively used as the characteristic sequence for two improved GM(0, N) models, obtaining the upper and lower bounds of the forecasting value (i.e., forecasting interval). Furthermore, the model parameters are identified by an improved PSO instead of LSM. Finally, the definite forecasting value is believed to be related to this forecasting interval and further determined by a novel whitenization method. In this section, the implementation of this improved grey prediction model (defined as the interval GM(0, N) model) is given in detail.

3.1. Selection of the Benchmark Ranking Sequence

The factor sequence that shares the most similarities with the characteristic sequence is selected to help with determining the position of the forecasting point in the MCP sequence. The chosen factor sequence is called the benchmark ranking sequence, and the steps are given below.

Step 1. Get factor sequences that have the same (or reverse) ranking order as the characteristic sequence, defined as $X_i^{sub(0)}$. Two sequences $X=(x_1,\cdots ,x_k,\cdots ,x_n)$ and $Y=(y_1,\cdots ,y_k,\cdots ,y_n)$ are regarded as having the same (or reverse) ranking order, if:

(1)

The lengths of $X$ and $Y$ are equal.

(2)

A new sequence $\tilde{X}=\{x_1,\cdots ,x_k,\cdots ,x_n\}$ is obtained by ranking $X$ in ascending order of value. Corresponding with the same index of the subscript, $\tilde{Y}=\{y_1,\cdots ,y_k,\cdots ,y_n\}$ is formed. For any $k$, there is $y_k-y_{k-1}\ge 0$ (or $<0$).

Similarly, the factor sequences $X_i^{sub(0)}$ can be obtained, $i\in I^{sub},I^{sub}\subset I$.

Step 2. The correlation coefficient $r_{1,i}$ between $X_1^{(0)}$ and $X_i^{(0)}$ is calculated by correlation analysis, shown in Equation (14).

$$r_{1,i}=r(X_1^{(0)},X_i^{(0)})=\frac{{\displaystyle \sum _{k=1}^n(x_1^{(0)}(k)-{\overline{x}}_1^{(0)})(x_i^{(0)}(k)-{\overline{x}}_i^{(0)})}}{\sqrt{{\displaystyle \sum _{k=1}^n(x_1^{(0)}(k)-{\overline{x}}_1^{(0)}){\displaystyle \sum _{k=1}^n(x_i^{(0)}(k)-{\overline{x}}_i^{(0)})}}}},(i\in I)$$

(14)

where ${\overline{x}}_i^{(0)}=\frac{1}{n}{\displaystyle \sum _{k=1}^n{x}_i^{(0)}(k)}$ is the mean value of $X_i^{(0)}$. Then, the benchmark ranking sequence (define its index in factor sequences as $\mathsf{\tau }$) is determined by the following rules: if $I^{sub}\ne \varnothing$, $\mathsf{\tau }$ is the index when $r_{1,\tau }=\underset{i\in I^{sub}}{\mathrm{MAX}}\{r_{1,i}\}$, else if $I^{sub}=\varnothing$, $\mathsf{\tau }$ is the index when $r_{1,\tau }=\underset{i\in I}{\mathrm{MAX}}\{r_{1,i}\}$. Therefore, the characteristic sequence and factor sequences are ranked in accordance with the ascending or descending order of the benchmark ranking sequence. In what follows, sequences mean ranked sequences.

3.2. Calculation for Forecasting Interval

By comparing the corresponding factor value of the forecasting point in the benchmark ranking sequence, the position of the forecasting point can be obtained, which is also regarded as the position in the ranked MCP sequence. Suppose the position is $k$; the ranked MCP sequence that contains the forecasting point is:

$$X_1^{(0)}=\{x_1^{(0)}(1),x_1^{(0)}(2),\cdots ,x_1^{(0)}(k-1),\varnothing (k),x_1^{(0)}(k+1),\cdots ,x_1^{(0)}(n+1)\}$$

(15)

where the $\varnothing (k)$ is the point to be forecasted and the data length is changed from $n$ to $n+1$. Defining the lower- and upper-bound virtual values of $x_1^{(0)}(k)$ as $x_1^{low(0)}(k)$ and $x_1^{up(0)}(k)$, respectively. Therefore, the lower- and upper-bound MCP sequences can be constructed, respectively named as $X_1^{low(0)}$ and $X_1^{up(0)}$, and expressed as:

$$\{\begin{array}{l}{X}_1^{low(0)}=\{x_1^{(0)}(1),x_1^{(0)}(2),\cdots ,x_1^{(0)}(k-1),x_1^{low(0)}(k),x_1^{(0)}(k+1),\cdots ,x_1^{(0)}(n+1)\}\\ X_1^{up(0)}=\{x_1^{(0)}(1),x_1^{(0)}(2),\cdots ,x_1^{(0)}(k-1),x_1^{up(0)}(k),x_1^{(0)}(k+1),\cdots ,x_1^{(0)}(n+1)\}\end{array}$$

(16)

two neighboring (lower and upper) points of $x_1^{(0)}(k)$ in the ranked MCP sequence, i.e., $x_1^{(0)}(k-1)$ and $x_1^{(0)}(k+1)$, are used as the virtual values. Thus, Equation (16) is:

$$\{\begin{array}{l}{X}_1^{low(0)}=\{x_1^{(0)}(1),x_1^{(0)}(2),\cdots ,x_1^{(0)}(k-1),x_1^{(0)}(k-1),x_1^{(0)}(k+1),\cdots ,x_1^{(0)}(n+1)\}\\ X_1^{up(0)}=\{x_1^{(0)}(1),x_1^{(0)}(2),\cdots ,x_1^{(0)}(k-1),x_1^{(0)}(k+1),x_1^{(0)}(k+1),\cdots ,x_1^{(0)}(n+1)\}\end{array}$$

(17)

Then, the $X_1^{low(0)}$ and $X_1^{up(0)}$ are respectively used as characteristic sequences for two improved GM(0, N) models, correspondingly offering the lower- and upper-bound forecasting values, defined as ${\widehat{x}}_1^{low(0)}(k)$ and ${\widehat{x}}_1^{up(0)}(k)$. Especially, when the position is at the beginning (or ending) of the MCP sequence, there only exists an upper-bound (or lower-bound) MCP sequence. In this case, a one-sided interval is obtained. In the improved GM(0, N) model, all model parameters are identified by an improved PSO (detailed in Section 3.4), and input sequences include not only known data, but also the virtual MCP and predicted factors of the forecasting point.

3.3. Definite Forecasting Value Determination

When the forecasting interval $\left[{\widehat{x}}_1^{low(0)}(k),{\widehat{x}}_1^{up(0)}(k)\right]$ is gained, the purpose of the whitenization method in GST is to obtain the definite forecasting value. Generally, the expression is:

$${\widehat{x}}_1^{(0)}(k)=\mathsf{\alpha }{\widehat{x}}_1^{low(0)}(k)+(1-\mathsf{\alpha }){\widehat{x}}_1^{up(0)}(k)$$

(18)

where $\mathsf{\alpha }$ is a coefficient that represents the proportion of ${\widehat{x}}_1^{low(0)}(k)$ to ${\widehat{x}}_1^{(0)}(k)$. This means that the definite forecasting value is a linear combination of the two boundary values. The most used is the equal proportion method, namely $\mathsf{\alpha }=0.5$. However, this may result in additional error because the gap between ${\widehat{x}}_1^{(0)}(k)$ and ${\widehat{x}}_1^{low(0)}(k)$ is not equal to that between ${\widehat{x}}_1^{(0)}(k)$ and ${\widehat{x}}_1^{up(0)}(k)$ in most cases. To give a more detailed explanation to coefficient $\mathsf{\alpha }$ and to calculate it accurately, a novel whitenization method is needed. Defining the gap between $x_i^{(0)}(k)$ and $x_i^{(0)}(k-1)$ as $\Delta x_i^{(0)}(k)$, the absolute error between ${\widehat{x}}_1^{(0)}(k)$ and ${\widehat{x}}_1^{low(0)}(k)$ as $\Delta {\widehat{x}}_1^{low(0)}(k)$, that is:

$$\{\begin{array}{l}\Delta x_i^{(0)}(k)=x_i^{(0)}(k)-x_i^{(0)}(k-1),i=\{1,2,\cdots ,N\}\\ \Delta {\widehat{x}}_1^{low(0)}(k)=|{\widehat{x}}_1^{(0)}(k)-{\widehat{x}}_1^{low(0)}(k)|\end{array}$$

(19)

as the ${\widehat{x}}_1^{low(0)}(k)$ is simulated from the sequence $X_1^{low(0)}$ in which the lower-bound value $x_1^{(0)}(k-1)$ is used as the replacement value of point $k$; the gap $\Delta x_1^{(0)}(k)=x_1^{(0)}(k)-x_1^{(0)}(k-1)=x_1^{(0)}(k)-x_1^{low(0)}(k)$ is the root of the existence of $\Delta {\widehat{x}}_1^{low(0)}(k)$. The amount of $\Delta x_1^{(0)}(k)$ certainly reflects the impact of ${\widehat{x}}_1^{low(0)}(k)$ on ${\widehat{x}}_1^0(k)$. Similarly, $\Delta x_1^{(0)}(k+1)$ reflects the impact of ${\widehat{x}}_1^{up(0)}(k)$ on ${\widehat{x}}_1^0(k)$. Thus, the coefficient $\mathsf{\alpha }$, which represents the proportion of ${\widehat{x}}_1^{low(0)}(k)$ to ${\widehat{x}}_1^{(0)}(k)$, is related with $\Delta x_1^{(0)}(k)$ and $\Delta x_1^{(0)}(k+1)$ and can be calculated by Equation (20). A schematic diagram is also shown in Figure 1. Especially with the one-sided interval as explained in Section 3.2, the definite forecasting value is only related with $\Delta x_1^{(0)}(k)$ or $\Delta x_1^{(0)}(k+1)$, and the value of $\mathsf{\alpha }$ respectively equals zero and one when the position is the beginning and ending.

$$\mathsf{\alpha }=\frac{\Delta x_1^{(0)}(k+1)}{\Delta x_1^{(0)}(k)+\Delta x_1^{(0)}(k+1)}$$

(20)

However, $x_1^{(0)}(k)$ is unknown; thus, $\Delta x_1^{(0)}(k)$ and $\Delta x_1^{(0)}(k+1)$ are also unknown. Therefore, coefficient $\mathsf{\alpha }$ cannot be calculated by Equation (20) directly. In this paper, considering the correlations between electricity MCP and influence factors, a novel whitenization method is proposed. The value of coefficient $\mathsf{\alpha }$ is calculated by correlation coefficient weighted gaps of factor sequences, expressed in Equation (21).

$$\{\begin{array}{l}{\mathsf{\alpha }}_i=\frac{\Delta x_i^{(0)}(k+1)}{\Delta x_i^{(0)}(k)+\Delta x_i^{(0)}(k+1)}\\ \mathsf{\alpha }={\displaystyle \sum _{i=2}^N\left(\frac{\left|r_{1,i}\right|}{{\displaystyle \sum _{i=2}^N\left|r_{1,i}\right|}}\times {\mathsf{\alpha }}_i\right)}\end{array},i=\{2,3,\dots ,N\}$$

(21)

where ${\mathrm{r}}_{1,i}$ is the correlation coefficient between $X_1^{(0)}$ and $X_i^{(0)}$ expressed in Equation (14).

3.4. Parameters Identification by Improved Particle Swarm Optimization (PSO)

As mentioned in Section 2.1, model parameters to be identified are $b_i,i\in \{2,3,\cdots ,N\}$ and $a$. Under given input data including the characteristic sequence and factor sequences, the goal is to seek a vector of these parameters, which provides the best forecasting results. The optimal model can be formulated as follows:

$$\begin{array}{l}\min f(b_i,a)=MAPE=\frac{1}{n+1}{\displaystyle \sum _{k=1}^{n+1}\left|\frac{x_1^{(0)}(k)-{\widehat{x}}_1^{(0)}(k)}{x_1^{(0)}(k)}\right|}\times 100\%\\ \begin{array}{ll}\mathrm{subject}\mathrm{to}\hfill & b_i\in (-\infty ,+\infty ),i\in \{2,3,\cdots ,N\}\hfill \\ & a\in (-\infty ,+\infty )\hfill \end{array}\end{array}$$

(22)

where $n+1$ is the data length of the characteristic sequence (containing the forecasting point).

Particle swarm optimization (PSO) is an evolutionary algorithm introduced by Kennedy and Eberhart in 1995 [46]. Due to its simplicity of implementation and ability to quickly converge to a reasonably good solution, PSO has been very widely used in many fields. However, for multimodal functions or high-dimensional problems, the canonical PSO tends to be trapped in premature convergence and provides a poor solution. Therefore, a PSO based on time-varying parameters is used in this paper to identify parameters for the proposed interval GM(0, N) model. It has been proven that the time-varying PSO can significantly improve algorithm performance in terms of the global search ability and the convergence speed [47,48,49,50]. The procedure of time-varying PSO is described in detail as follows.

Step 1. Initialize particles. The number of particles is set to 60 after repeated testing. The position vector of the $t$-th particle is defined as $X_t=\{x_{t1},x_{t2},\cdots ,x_{tN}\}$ and the velocity vector as $V_t=\{v_{t1},v_{t2},\cdots ,v_{tN}\}$, where $t=1,2,3,\cdots ,60$. The positions of particles correspond to the solution of $b_i,i\in \{2,3,\cdots ,N\}$ and $a$. For faster convergence, the initial particles are constructed with the assistance of the multiple linear regression (MLR) model. The MLR model takes the form:

$$y=\mathsf{\beta }X+\mathsf{\epsilon }$$

(23)

where $y$, $X$ are respectively the AGO characteristic sequence $X_1^{(1)}$ and factor sequences $X_i^{(1)},i=\{2,3,\cdots ,N\}$. Suppose the solution is $X^{ini}=\{\mathsf{\beta },\mathsf{\epsilon }\}=\{x_1^{ini},x_2^{ini},\cdots ,x_N^{ini}\}$, then:

$$x_{td}=x_d^{ini}+r(x_{max}-x_{min})+x_{min},t=1,2,3,\cdots ,60,d=1,2,3,\cdots ,D$$

(24)

where $x_{max}=\underset{1\le d\le N}{\mathrm{MAX}}\{|x_d^{ini}|\}$, $x_{min}=\underset{1\le d\le N}{-\mathrm{MAX}}\{|x_d^{ini}|\}$, $r$ is a random number in interval [0, 1], $d$ is the index of parameters and $D$ is the dimension of particle, $D=N$.

Step 2. Evaluate fitness for each particle and further identify the best position (defined as $P_t^k=\{p_{t1}^k,p_{t2}^k,\cdots ,p_{tD}^k\}$) of each particle and the best position (defined as $P_g^k=\{p_{g1}^k,p_{g2}^k,\cdots ,p_{gD}^k\}$) of the swarm. Up to the $k$-th evolution generation, let $p_{tbest}$ denote the best fitness of the $t$-th particle. For each particle, if its current fitness value is better than $p_{tbest}$, then set this value as $p_{tbest}$ and its position as $P_t^k$. Similarly, the particle with the best fitness of the swarm is identified, and its position is defined as $P_g^k$.

Step 3. Update the position and velocity for each particle, shown in Equation (25).

$$\{\begin{array}{l}{v}_{td}^{k+1}=\mathsf{\omega }{v}_{td}^k+c_1{r}_1(p_{td}^k-x_{td}^k)+c_2{r}_2(p_{gd}^k-x_{td}^k)\\ x_{td}^{k+1}=x_{td}^k+v_{td}^{k+1}\end{array}$$

(25)

where $\mathsf{\omega }$ is the inertia weight; $d=1,2,3,\cdots ,D$; $t=1,2,3,\cdots ,60$; $k$ is the current evolution generation; $c_1,c_2$ are two acceleration coefficients in which $c_1$ is named the cognitive learning rate and $c_2$ is named the social learning rate; $r_1,r_2$ are two random numbers in the interval [0, 1].

$\mathsf{\omega }$, $c_1$ and $c_2$ are three main parameters in PSO, which can significantly affect the performance of PSO. It has been demonstrated that a relatively larger value of $\mathsf{\omega }$ is beneficial for global search, and a smaller one is good for local search [51]. Therefore, Shi [47] introduced a linear decreasing inertia weight, and later, Chatterjee [48] developed the trajectory of $\mathsf{\omega }$ according to a nonlinear function, which is also used in this article, shown in Equation (26).

$$\mathsf{\omega }={\mathsf{\omega }}_{start}-({\mathsf{\omega }}_{start}-{\mathsf{\omega }}_{end}){(\frac{g}{G_{max}})}^2$$

(26)

where ${\mathsf{\omega }}_{start}$ and ${\mathsf{\omega }}_{end}$ are the initial value and end value, respectively. $g$ and $G_{max}$ represent the present and maximum evolution generation, respectively. Generally, ${\mathsf{\omega }}_{start}=0.9$, ${\mathsf{\omega }}_{end}=0.4$. $G_{max}$ is set as constant 500 in this paper. In addition, Eberhart [46] described that a relatively high cognitive learning rate $c_1$ is beneficial for particles to wander through the entire search space, and a large social learning rate $c_2$ performs better in local search. Therefore, Ratnaweera [49] and Wang [50] proposed time-varying acceleration coefficients, which can be mathematically represented as follows:

$$c_1=c_{1,max}+(c_{1,min}-c_{1,max})(\frac{g}{G_{max}})$$

(27)

$$c_2=c_{2,min}+(c_{2,max}-c_{2,min})(\frac{g}{G_{max}})$$

(28)

where, $c_{1,min}$, $c_{1,max}$, $c_{2,min}$ and $c_{2,max}$ are constant factors. To find out the best ranges of $c_1$ and $c_2$, series numerical simulations are also carried out (detailed in Section 5.1).

Step 4. Repeat Steps 2–4 until the iteration reaches its maximum value. Therefore, the best fitness value can be confirmed, and its corresponding particle is obtained. The position of this particle is the optimal solution for parameters $b_i$ and $a$.

3.5. Summary of Calculation Process

As described above, the calculation process of our proposed interval GM(0, N) model can be summarized as below in Figure 2.

4. Study Area and Influence Factors

4.1. Newly-Reformed Yunnan Electricity Market

The Chinese government is starting its implementation of further electricity market reform in several piloting provinces, and Yunnan province is one of them [52]. The insufficient electricity demand and great surplus of clean energy are the main reasons for Yunnan to promote the reform of its electricity sector. At the early stage of electricity market reform, to guarantee the stability and security operation of the power grid, some units are treated as “must run units”, which do not participate in the market competition. Meanwhile, to meet the requirements of environmental policies, electricity generated by clean and renewable energy sources should be fully acquired in a competitive market, such as wind, solar and small hydro. The generation contracts are made in monthly electricity market. The process and key times of Yunnan’s monthly electricity market are shown in Figure 3. At the closure of bidding, the Yunnan Power Exchange Center (YPEC) aggregates the submitted seller and buyer bids. Then, the market supply and demand curves are generated, and the electricity MCP is calculated.

4.2. Identification of Input Variables

The biggest difference between electricity and other commodities is that it cannot be stored in large quantities, which leads to a high volatility of electricity price. Many factors may determine or affect electricity price, such as historical price, historical/forecasted load, weather conditions, macroeconomic policies, competitors’ irrational behaviors, and so on [5,10,53]. Till now, there are no consensus factors for price forecasting. Generally, the researchers tend to utilize past experience, also taking into consideration the unique economic environments of different markets, in selecting the input variables for their respective models. Therefore, we need to understand Yunnan’s supply and demand conditions and select appropriate input variables. In this paper, factors that cannot be quantified and expressed in numerical values are not considered, such as macroeconomic policies and competitors’ irrational behaviors.

Yunnan has extremely rich hydropower resources, and three of the 13 hydropower bases in China are located here. About 74.2% of the total installed capacity is hydropower, with 14.0% of small hydropower and 60.2% of medium/large hydropower. Besides hydropower, thermal, wind and solar power are also included, accounting for 16.2%, 8.0% and 1.5% of the total installed capacity respectively. Over 77% of the yearly generation in Yunnan is medium/large hydropower based; around 10% is generated by thermal power; and the remaining 13% is a mixture of production types including small hydro, wind and solar power. By the end of 2015, Yunnan’s installed capacity and electricity generation of each energy resource are shown in Figure 4 and Figure 5. Yunnan is also an important electricity supplier for China’s west-to-east electricity transmission project, Laos and Vietnam. In 2015, about 44% of the total energy production has been exported. Therefore, its electricity demand can be classified into two types: one is provincial electricity demand, and the other is export electricity demand.

Based on the experience, electricity supply and demand conditions of the Yunnan market, three types of influence factors are considered.

4.2.1. Supply Factors

The medium/large hydropower and thermal power are the main energy sources of Yunnan, which would directly affect the electricity MCP. As clean and renewable energy sources, the wind power, small hydropower and solar power all have the priority right to be purchased, which would increase the competition between other generators without the priority right and affect the electricity MCP, as well. Therefore, the available energy of each energy resource is concerned, and the energy production is considered as the input variable on the supply factors’ side. It is noted that the factors that can further influence energy production are not included. The reasons for this are: (1) compared with energy production, the factors that can influence energy production are having relatively indirect influences on the market price; to a certain degree, their influences on the market price have already been reflected in that of energy production; and (2) the factors that influence energy production can be various, and some of them are unavailable. Take hydropower as an example: its energy production is strongly correlated with runoffs, reservoir’s capacity and level, unit maintenance, flood prevention requirements, operation rules, and so on. It seems extremely difficult to consider all of these factors; the adopted energy production is a synthetic effect of these on electricity MCP. Besides, the portion that does not participate in the market competition (such as must run units) has been excluded from the available energies.

4.2.2. Demand Factors

The electricity MCP is obviously related to the market demand. Currently, in Yunnan, more than 40% of generated electricity is transmitted to outside areas through mid- or long-term contracts. The amount of export electricity has an influence on the competition of provincial electricity consumers and should affect the electricity MCP. For the provincial electricity demand, about 90% is provided in the monthly electricity market. The export electricity demand and provincial electricity demand are the main demand factors and selected as input variables.

4.2.3. Supplemented Factors

In Yunnan’s developing electricity market, there are only two types of participants so far: the generating companies (GENCOs) and consumption companies (CONCOs). The number of GENCOs and CONCOs reflects the competition in the supply and demand sides. With more market participants, the market tends to be a perfectly competitive market, and the electricity MCP changes. Here, the number of GENCOs and CONCOs in the electricity market is used as supplemented input variables.

The input variables for electricity MCP forecasting considered in this article are listed in the following

Table 2. The input variables for electricity market clearing price (MCP) forecasting.

Index	Types	Input Variables	Symbols
1	Supply factors	Energy production of medium/large hydro power	$F_1$
2		Energy production of thermal power	$F_2$
3		Energy production of small hydro power	$F_3$
4		Energy production of wind power	$F_4$
5		Energy production of solar power	$F_5$
6	Demand factors	Export electricity demand	$F_6$
7		Provincial electricity demand	$F_7$
8	Supplemented factors	Number of GENCOs	$F_8$
9		Number of CONCOs	$F_9$

. Meanwhile, we should notice that the selection of appropriate input variables plays an important role in ensuring accuracy. For different electricity markets, the energy source compositions and economic environments are obviously different, resulting in specific input variables. The proposed model is promising to be applied to other electricity markets, but the input variables should be carefully re-identified before modeling.

4.3. Data Collection and Normalization

For the case study, actual monthly data including MCP sequence and factor sequences in the Yunnan electricity market, from April 2015–April 2016 are used. The data shown in

Table 3. Monthly MCP and influence factors of the Yunnan electricity market from April 2015–April 2016.

Month	MCP	${\mathit{F}}_1$	${\mathit{F}}_2$	${\mathit{F}}_3$	${\mathit{F}}_4$	${\mathit{F}}_5$	${\mathit{F}}_6$	${\mathit{F}}_7$	${\mathit{F}}_8$	${\mathit{F}}_9$
April 2015	0.9174	0.2015	0.8689	0.0374	0.3354	0.0945	0.2972	0.7269	0.3462	0.0044
May 2015	0.4128	0.3630	0.8592	0.0290	0.2925	0.0199	0.5315	0.7343	0.2308	0.0000
June 2015	0.0494	0.5985	0.3532	0.3125	0.1701	0.0000	0.6965	0.8876	0.3846	0.1027
July 2015	0.0000	1.0000	0.3647	0.8816	0.0000	0.0995	0.9156	0.8774	0.6154	0.1754
August 2015	0.0105	0.9298	0.2992	1.0000	0.0241	0.0995	1.0000	0.9368	0.5000	0.1789
September 2015	0.0296	0.7841	0.0000	0.8854	0.0377	0.0945	0.9657	1.0000	0.2308	0.2170
October 2015	0.0303	0.7627	0.7197	0.5264	0.2543	0.1592	0.8949	0.9885	0.2692	0.2365
November 2015	0.1692	0.3667	0.7306	0.3354	0.3747	0.2985	0.3693	0.5218	0.2692	0.2294
December 2015	0.4726	0.3381	1.0000	0.2895	0.4741	0.3483	0.3101	0.5220	0.0000	0.1798
January 2016	0.9259	0.1890	0.6754	0.1360	0.6410	0.5323	0.2268	0.8476	1.0000	0.8335
February 2016	1.0000	0.0000	0.4339	0.0000	0.8095	0.6169	0.0000	0.0000	0.9231	0.6740
March 2016	0.7791	0.2034	0.5140	0.0053	1.0000	1.0000	0.2783	0.4350	0.9231	0.7741
April 2016	0.7330	0.2547	0.5546	0.0031	0.6599	0.7761	0.2992	0.7086	1.0000	1.0000

are normalized by the min-max normalization ($Nor.$), which is formulated as:

$$X^{\prime }=\frac{X-X_{min}}{X_{max}-X_{min}}$$

(29)

where $X_{max}$ and $X_{min}$ are the maximum and minimum values of the $X$ sequence, respectively. Correspondingly, the inverse normalization transformation ($Inv.$) is:

$$X=X^{\prime }(X_{max}-X_{min})+X_{min}$$

(30)

5. Case Study

5.1. An Example for Forecasting

To detail the forecasting procedure of the proposed interval GM(0, N) model, a forecasting example for April 2016 is given. In this case study, the proposed model is trained by data from April 2015–March 2016. The procedure is described as follows.

Step 1. According to Section 3.1, the influence factor $F_1$ is selected as the benchmark ranking sequence because it has nearly the reverse ranking order and the highest correlation coefficient with the MCP sequence. The correlation coefficient is 0.9059. Then, all data sequences that contain the forecasting point are ranked in accordance with the descending order of $F_1$, partly shown in

Table 4. Ranked monthly MCP and the influence factors that contain the forecasting point.

Month	MCP	${\mathit{F}}_1$	${\mathit{F}}_2$	${\mathit{F}}_3$	${\mathit{F}}_4$	${\mathit{F}}_5$	${\mathit{F}}_6$	${\mathit{F}}_7$	${\mathit{F}}_8$	${\mathit{F}}_9$
…	-	-	-	…		-	-	-	…
May 2015	0.4128	0.3630	0.8592	0.0290	0.2925	0.0199	0.5315	0.7343	0.2308	0.0000
December 2015	0.4726	0.3381	1.0000	0.2895	0.4741	0.3483	0.3101	0.5220	0.0000	0.1798
April 2016	$\varnothing (k)$	0.2547	0.5546	0.0031	0.6599	0.7761	0.2992	0.7086	1.0000	1.0000
March 2016	0.7791	0.2034	0.5140	0.0053	1.0000	1.0000	0.2783	0.4350	0.9231	0.7741
April 2015	0.9174	0.2015	0.8689	0.0374	0.3354	0.0945	0.2972	0.7269	0.3462	0.0044
…	-	-	-	…		-	-	-	…

.

Step 2. The position of forecasting point is determined by comparing the $F_1$ value of April 2016 with other months, i.e., between December 2015 and March 2016. Two neighboring values, respectively 0.7791 and 0.4726, are used as two virtual values for April 2016. Thus, the upper- and lower-bound MCP sequences are constructed. They are:

$$X_1^{up(0)}=\{0.0000,0.0105,\cdots ,0.4726,0.7791,0.7791,\cdots ,1.0000\}$$

$$X_1^{low(0)}=\{0.0000,0.0105,\cdots ,0.4726,0.4726,0.7791,\cdots ,1.0000\}$$

Step 3. Two improved GM(0, N) models are respectively built from the upper- and lower-bound MCP sequences combined with ranked factor sequences. Additionally, the model parameters are identified by the time-varying PSO as described in Section 3.4. To find out the best ranges of $c_1$ and $c_2$, numerical simulations are carried out. Based on the values used in [49,50,54,55,56,57,58], the test range of $c_1+c_2$ is set to [3.0, 5.0]. The average and standard deviation of the optimum solutions for 100 trials of different values of $c_1+c_2$ are shown in

Table 5. The average optimum value (AOV) and standard deviation (SD) for 100 trials of different values of $c_1+c_2$.

	Different Values and Ranges of ${\mathit{c}}_{\mathbf{1}}$ and ${\mathit{c}}_{\mathbf{2}}$
Statistical Index	$c_1+c_2$ = 3.1	$c_1+c_2$ = 3.2	$c_1+c_2$ = 3.3	$c_1+c_2$ = 3.4	$c_1+c_2$ = 3.5
	$c_1$ = 2.6~0.5	$c_1$ = 2.7~0.5	$c_1$ = 2.8~0.5	$c_1$ = 2.9~0.5	$c_1$ = 3.0~0.5
	$c_2$ = 0.5~2.6	$c_2$ = 0.5~2.7	$c_2$ = 0.5~2.8	$c_2$ = 0.5~2.9	$c_2$ = 0.5~3.0
AOV (%)	4.41	4.16	4.06	4.13	4.10
SD (×10−2)	1.67	0.81	0.60	0.64	1.04
Statistical Index	$c_1+c_2$ = 3.6	$c_1+c_2$ = 3.7	$c_1+c_2$ = 3.8	$c_1+c_2$ = 3.9	$c_1+c_2$ = 4.0
	$c_1$ = 3.1~0.5	$c_1$ = 3.2~0.5	$c_1$ = 3.3~0.5	$c_1$ = 3.4~0.5	$c_1$ = 3.5~0.5
	$c_2$ = 0.5~3.1	$c_2$ = 0.5~3.2	$c_2$ = 0.5~3.3	$c_2$ = 0.5~3.4	$c_2$ = 0.5~3.5
AOV (%)	4.12	4.26	3.99	3.76	3.74
SD (×10−2)	0.76	1.46	0.48	0.73	0.49
Statistical Index	$c_1+c_2$ = 4.1	$c_1+c_2$ = 4.2	$c_1+c_2$ = 4.3	$c_1+c_2$ = 4.4	$c_1+c_2$ = 4.5
	$c_1$ = 3.6~0.5	$c_1$ = 3.7~0.5	$c_1$ = 3.8~0.5	$c_1$ = 3.9~0.5	$c_1$ = 4.0~0.5
	$c_2$ = 0.5~3.6	$c_2$ = 0.5~3.7	$c_2$ = 0.5~3.8	$c_2$ = 0.5~3.9	$c_2$ = 0.5~4.0
AOV (%)	3.80	3.67	3.87	3.90	3.93
SD (×10−2)	1.01	0.35	0.72	0.91	1.59
Statistical Index	$c_1+c_2$ = 4.6	$c_1+c_2$ = 4.7	$c_1+c_2$ = 4.8	$c_1+c_2$ = 4.9	$c_1+c_2$ = 5.0
	$c_1$ = 4.1~0.5	$c_1$ = 4.2~0.5	$c_1$ = 4.3~0.5	$c_1$ = 4.4~0.5	$c_1$ = 4.5~0.5
	$c_2$ = 0.5~4.1	$c_2$ = 0.5~4.2	$c_2$ = 0.5~4.3	$c_2$ = 0.5~4.4	$c_2$ = 0.5~4.5
AOV (%)	4.19	4.17	4.16	4.10	4.28
SD (×10−2)	0.69	0.57	0.71	0.56	0.74

Notes: $c_1$ = 2.6 ~ 0.5 means that $c_1$ is changing from 2.6 to 0.5, and similar to the others.

, and curves are presented in Figure 6 accordingly. From the result, we can see that when $c_1+c_2=4.2$, the performance is best. Thus, in our study, with the increasing of iterations, $c_1$ is decreasing from 3.7 down to 0.5, and $c_2$ is increasing from 0.5–3.7. In this way, the resulting MCP forecasting interval is [0.4787, 0.7910]. The evolution process of best fitness when calculating upper- and lower-bound forecasting values by the time-varying PSO is shown below in Figure 7.

Step 4. As described in Section 3.3, the proportion coefficient $\alpha$ is calculated by Equation (21). Thus, the definite forecasting value of April 2016 is:

$${\widehat{x}}_1^{(0)}=\alpha \cdot {\widehat{x}}_1^{low(0)}+(1-\alpha )\cdot {\widehat{x}}_1^{up(0)}=0.4700\times 0.4787+(1-0.4700)\times 0.7910=0.6442$$

Therefore, the forecasting value of April 2016 is 0.2803 after inverse normalization transformation. As the observed value is 0.2893, the MAE and MAPE are respectively $|0.2893-0.2803|=0.0090$ and $|0.2893-0.2803|÷0.2893\times 100\%=3.10\%$. The forecasting result is acceptable. It is important to note that the forecasting result is based on the assumption that the input factor sequences of the forecasting point are accurately predicted.

5.2. Sensitivity Analysis of Position

To validate the applicability of the proposed interval GM(0, N) model, a sensitivity analysis of the position of the forecasting point is discussed. Here, also take the MCP forecasting of April 2016 as an example. The true position of April 2016 in the ranked MCP sequence is Position 1. Suppose that the position is somehow misjudged; define the left position as Position 2 and the right as Position 3, as shown in Figure 8.

The proposed interval GM(0, N) model is applied to three different positions respectively, and the results are shown in

Table 6. Forecasting results of April 2016 in Positions 1–3.

Position	Observed Value		Virtual Values		Forecasting Interval		Forecasting Value			MAPE (%)
	Nor.	Inv.	Lower	Upper	Lower	Upper	α	Nor.	Inv.
Position 1	0.7330	0.2893	0.4726	0.7791	0.4787	0.7910	0.4700	0.6442	0.2803	3.10
Position 2			0.4128	0.4726	0.4186	0.4702	−0.5047	0.4962	0.2653	8.28
Position 3			0.7791	0.9174	0.7735	0.9098	2.3012	0.5961	0.2755	4.79
Average	0.7330	0.2893	0.5548	0.7230	0.5569	0.7237	0.7555	0.5789	0.2737	5.39

. The MAPE for Position 1 is 3.10%, which is the best among the three positions. It obviously shows that the correct position estimation can provide a more accurate forecasting interval, which has a direct influence on the forecasting result. For Positions 2 and 3, the MAPEs are 8.28% and 4.79%, respectively. Compared to Position 1, although their forecasting performances are deteriorating, they are also acceptable in mid-term forecasting. The MAPE of the three positions is 5.39%. By taking into account the correlation between electricity MCP and influence factors, the coefficient $\mathsf{\alpha }$ given by the developed novel whitenization method has a certain degree of correction to the forecasting interval. As shown in Figure 9, whether the position of the forecasting point is correct or not, the coefficient $\mathsf{\alpha }$ will guide the forecasting interval towards the actual value, which also shows a better result than the equal proportion method ($\mathsf{\alpha }=0.5$). For Positions 1 and 2, the forecasting interval is on the left of the observed value, so the direction of $\mathsf{\alpha }$ is from left to right. On the other hand, the forecasting interval of Position 3 is on the right; the direction is accordingly from right to left. Although the position of the forecasting point is not exactly judged in some cases, if the deviation of the position is not too big, the proposed model can obtain an acceptable forecasting result.

5.3. Comparison with Other Models

To further illustrate the forecasting performance of the proposed interval GM(0, N) model (Model 4), three other forecasting models are established for comparison, i.e., the multiple linear regression (MLR) model (Model 1), the traditional GM(0, N) model (Model 2) and the artificial neural network (ANN) model (Model 3). The MLR model is a typical static model for price forecasting by investigating the relationship between a dependent variable and several independent variables, which is also shaped like the proposed model. The used ANN model is a typical three-layered back propagation neural network model, including one input layer, one hidden layer and one output layer, which is usually preferred in practical engineering applications [23,59,60]. The sigmoid transfer function is used in both the neurons of the hidden and output layers. Nine inputs (influence factors) and one output (market price) are applied to the ANN model. At the training stage, the Levenberg–Marquardt algorithm is used to train the ANN model. Various numbers of neurons in the hidden layer are tested, and all statistical criteria of these network structures are recorded and compared. The best results are produced with 11 hidden units. The maximum iterations used in this paper is 1000. All four models are applied to forecast electricity MCP for each month in turn from April 2015–April 2016 under the same experimental conditions. For each month, the data of other months are used as input data for building these models respectively.

Forecasting results of these models are shown in

Table 7. Forecasting results of the proposed model compared with the other three models.

Month	Observed Value	Forecasting Value				Absolute Percentage Error (%)
		Model 1	Model 2	Model 3	Model 4	Model 1	Model 2	Model 3	Model 4
April 2015	0.4471	0.3394	0.0000	0.3627	0.3929	13.08	54.32	10.25	6.58
May 2015	0.3107	0.3253	0.1370	0.2841	0.2778	2.14	25.30	3.87	4.79
June 2015	0.2124	0.3101	0.1957	0.2815	0.1920	16.62	2.83	11.75	3.46
July 2015	0.1990	0.1023	0.1901	0.1794	0.2334	16.81	1.55	3.42	5.98
August 2015	0.2018	0.2403	0.1369	0.2595	0.2005	6.65	11.24	9.98	0.23
September 2015	0.2070	0.0733	0.3286	0.1949	0.2056	22.94	20.85	2.09	0.25
October 2015	0.2072	0.1795	0.3583	0.2202	0.2033	4.75	25.90	2.23	0.68
November 2015	0.2448	0.3530	0.1696	0.3759	0.3111	17.43	12.11	21.12	10.68
December 2015	0.3268	0.1926	0.0917	0.3543	0.3069	19.10	33.45	3.91	2.83
January 2016	0.4494	0.6114	1.0000	0.3468	0.4981	19.62	66.71	12.43	5.89
February 2016	0.4695	0.0991	0.3321	0.3601	0.4307	43.80	16.24	12.93	4.59
March 2016	0.4097	0.9014	0.5564	0.4701	0.4125	62.58	18.66	7.68	0.35
April 2016	0.3972	0.1793	0.2502	0.4157	0.3732	28.18	19.01	2.39	3.10
MAPE	-	-	-	-	-	22.06	23.70	8.00	3.80

and Figure 10. We can see that the proposed interval GM(0, N) model shows a higher precision than the other compared models. The MAPE of the proposed model is 3.80%, and that of the MLR model is 22.06%; the traditional GM(0, N) model is 23.70%; the ANN model is 8.00%. Although the MAPE of the ANN model does not differ very greatly from that of the proposed model, the fluctuation of its results is much higher. The maximum absolute percentage error of the ANN model is 21.12%, which is not acceptable. The error variance of the proposed model is 0.0297, which is also much less than the MLR model (0.1581), the traditional GM(0, N) model (0.1801) and the ANN model (0.0550). For fitting precision checking, the $C$ and $P$ of the proposed model are 0.32 and 100.0%, respectively, showing a good fitting precision grade. While these of the traditional GM(0, N) model are 2.22 and 30.78%, respectively, the forecasting result is unqualified. This truly shows that our improvements on the traditional GM(0, N) model are effective in mid-term electricity MCP forecasting, and good forecasting results are obtained.

5.4. Forecasting Evaluation Based on the Modified Diebold–Mariano (MDM) Test

In this part, the forecasting performance is compared via the MDM test. The MDM test is carried out respectively between our proposed model and the other three forecasting models. The evaluation results are summarized in

Table 8. Forecasting evaluation results based on the modified Diebold–Mariano (MDM) test.

Indicators	MDM Test Results between Different Models
	Models 1 and 4	Models 2 and 4	Models 3 and 4
MDM-MAE	3.2614 ***	3.5409 ***	3.4856 ***
MDM-MSE	2.0354 *	2.0850 *	2.8789 **

Notes: MDM-MAE denotes the MDM test based on the MAE loss function; MDM-MSE denotes the MDM test based on the MSE loss function; ***, ** and * denote significance at the 1%, 5% and 10% levels, respectively.

.

From Table 8, conclusions can be drawn by comparison of Model 1 and Model 4:

• According to the MDM test based on the MAE loss function, the null hypothesis is rejected at the 1% level of significance. In other words, the observed differences are pretty significant, and the forecasting performance of Model 4 is better than Model 1.
• According to the MDM test based on the MSE loss function, the null hypothesis is rejected at the 10% level of significance. That is to say, the observed differences are also significant, and the forecasting performance of Model 4 is better than Model 1.

Similarly, according to the comparisons of Models 2 and 4 and Models 3 and 4, both the MDM test by MAE and MSE loss functions evaluate that the forecasting performance of Model 4 is better than Models 2 and 3. Therefore, the superiority of the proposed model is further validated statistically.

6. Conclusions

Accurate mid-term electricity MCP forecasting is essential for all market participants. In this paper, a novel interval GM(0, N) model for mid-term electricity MCP forecasting is proposed, aiming to present a feasible method for the electricity market with few available data. In the proposed model, two improved GM(0, N) models are included to respectively estimate the upper- and lower-bound forecasting values. Firstly, all input sequences (not containing further values) including the MCP sequence and factor sequences are ranked in accordance with the ascending order of the MCP sequence. Then, one of the factor sequences is selected as the benchmark ranking sequence according to its ranking order and correlation with the MCP sequence. The position of forecasting point in the MCP sequence thus can be determined by the benchmark ranking sequence. Finally, the upper and lower points of the forecasting point in the MCP sequence are regarded as two virtual values to construct two new MCP sequences, which are respectively used as input for two improved GM(0, N) models, obtaining the forecasting interval. In the two improved GM(0, N) models, the input sequences used for modeling include not only known data, but also the virtual MCP and predicted factors of the forecasting point, and model parameters are identified by an improved PSO instead of LSM. Based on the forecasting interval, a novel whitenization method, taking correlation coefficient weighted gaps of factor sequences into consideration, is also established to determine the definite forecasting value.

The proposed model has been applied to the newly-reformed Yunnan electricity market. Being at the early stage of reform, this market shares many common features with restructuring practices in other countries, while simultaneously exhibiting many unique characteristics, as well. After careful analysis of its market conditions, the input variables are appropriately determined, among which the energy production of medium/large hydro power is selected as the benchmark ranking sequence. For April 2016, the absolute percentage error of the forecasting value is 3.10%. Additionally, the MAPE of the proposed model for MCP forecasting from April 2015–April 2016 is 3.80%, which shows a better forecasting performance compared with the MLR model, the traditional GM(0, N) model and the ANN model, and its superiority is also confirmed by the modified Diebold–Mariano test. The good forecasting results indicate that the proposed model is suitable for mid-term electricity MCP forecasting in a data-sparse electricity market.