A data-driven approach for online aggregated load modeling through intelligent terminals

Abstract

Load model has significant impact on power system simulation. Current load modeling approaches are inadequate on revealing the accuracy and time-variation of load compositions. The application of wireless sensors dispersed in power distribution networks provides further opportunities for load modeling. In this article, a data-driven online aggregated load modeling approach is proposed systematically. First, all the electricity consumers are clustered according to big data of power consumption behaviors. In each cluster, typical users are designated to stand for the characteristics of the cluster, and intrusive measurement is adapted to capture these typical users’ time-varying information by employing wireless intelligent terminals, which can identify the composition of static load and induction motor load online. Second, the load models of other users in each cluster are assumed identical to typical users, including static impedance–current–power models and induction motor models. Finally, the composite load model is achieved by hierarchical aggregation and bottom-to-up stepwise equivalence. Simulations demonstrate that the load model built by proposed approach reflects higher accuracy and adaptability in power system.

Introduction

With the large-scale construction of ultra high voltage direct current/alternating current (UHVAC/DC) transmission project in China, the risk of large-capacity power shortage is increasing correspondingly in power system. ¹ The bipolar locking of Jinping-Sunan DC transmission line resulted in 4900 MW power shortage in East China power grid on 19 September 2015, and the lowest frequency declined to 49.56 Hz, which threatened the stable operation of power system severely. ² After accident, simulations show that the frequency response trajectory could not be reproduced exactly while adapting existing model parameters, among which the load model is one of the most important factors. The impact of induction motor (IM) proportion in load model on dynamic frequency characteristic is analyzed in Figure 1. Based on the adopted load model and simulation boundary, the frequency response characteristic slightly deteriorates with the increase of IM load proportion.

Figure 1.

The dynamic frequency trajectory under different proportions of IM load in East China power grid.

Moreover, the procedure of DC commutation failure is affected by load model as well. By means of fault scan in East China power grid, it is found that the three-phase N-1 fault of 500 kV AC line is likely to trigger simultaneous commutation failure of three DC lines (Fulong-Fengxian, JinPing-Sunan, and Yibin-Jinhua). The number of these “perilous” 500 kV-AC lines will significantly decrease while altering part of constant-power (CP) load to IM load. The detailed simulation results are shown in Table 1.

Table 1.

The impact of load model on DC commutation failure.

Load model	N-1 fault of 500 kV AC line lead to simultaneous commutation failureof three DC lines
	Number of“perilous” lines	Ratio in all 500 kV AC lines
60% CP+40% CZ	366	94%
60% IM+40% CZ	48	12%
40% IM+40%CP+20% CI	109	28%

CZ, CI, and CP represent the constant-impedance, constant-current, and constant-power part in static ZIP load model; DC: direct current.

Obviously, load model has crucial impact on power system simulation. Applying inaccurate load model may lead to discrepancy between simulation results and actual system response, and obtain ‘conservative” or “aggressive” conclusions, which is bad for the control strategy decision.^3–6

Various model structures have been studied over the history of load modeling practice. Nowadays, the composite load model composed of static constant impedance–current–power (ZIP) component and dynamic IM component is widely applied for power system simulation due to its preferable generalization capability.^7–9 As expected, two kinds of approaches for specifying the parameters of this composite load model are developed, namely, the component-based approach ¹⁰ and the measurement-based approach.^11–13 The component-based approach captures the inventory of load components connected to feeders so as to synthesize composite load model. ¹⁴ However, the statistical work is time consuming and laborious, and dynamic performance of load model is poorly expressed. The measurement-based approach regards various loads as a black box model, and its parameters are identified to fit the real field measurement, which reflects real dynamic behaviors and allows for direct application to practical power system. ¹⁵ Whereas, the measurement devices used by measurement-based approach are generally deployed at high-voltage feeders, so that the detailed load information at user-side is unavailable. Identifying accurate load model is complicated due to the following issues: (1) load composition is time variant and stochastic; (2) a wide range of diverse load components is penetrating into power system.

In recent years, continuously expansive power system and penetration of new load components have intensified the load diversity and time-variation, brought severe challenges to load modeling practice.^16,17 More literatures have involved in further innovative works for load modeling, expect to improve the accuracy of load models. These studies can be approximately concluded as two aspects:

1. Innovation of load model modality. Yu et al. ¹⁸ propose a new simplified composite load model by dominant parameters selection, in which the IM model is converted as a second-ordered state equation. Kim et al. ¹⁹ develop a fast-computational technique for deriving the composite load model parameters based on analytical similarity of parameter sensitivity. These literatures mainly focus on reducing the number of load model parameters to be identified and facilitating the parameter estimation in optimization process.

2. Optimization of load modeling approach. Kontis et al. ²⁰ employ the vector fitting technique for measurement-based dynamic load modeling. Proposed method in this paper exploits multi-signal from phasor measurement unit (PMU) measurement and facilitate the development of robust dynamic load models. Ge et al. ²¹ present a new event-oriented method for online load modeling, which is based on an adjustable sliding window applied to load step changes. These approaches can be regarded as the evolution of conventional measurement-based method.

Actually, recent studies have provided some preferable program for load modeling. Whereas, the inherent issue is unsettled, detailed information of bottomed load components is unavailable and its real-time characteristics cannot be tracked. Only by mining users’ power consumption data and commanding the dynamic characteristics of load components, the accurate time-varying load model can be built.

Data-driven online aggregated load modeling strategy

With the application of industrial wireless sensor technology, a large number of dispersed sensors are penetrating into electricity users, and more power data are available and allow real-time access.^22–25 These messages can be collected within a short time interval and imposed for online load modeling practices. Load monitoring approaches can be concluded as intrusive load monitoring (ILM) and non-intrusive load monitoring (NILM). NILM approach aims to identify appliances according to the power characteristics, which are disaggregated from user’s integrated power consumption.^26–28 High-frequency smart meter is usually employed as its monitoring device. ²⁹ The limitations of NILM are mainly reflected in two aspects: (1) high data sampling frequency is needed, generally required up to 1 kHz; (2) load components may not be accurately identified. ILM approach captures the electrical information by means of connecting sensors with load components directly inside users.

Due to the high identification accuracy of ILM approach, it is adopted for online load measurement and identification in our research. Specially, a kind of wireless intelligent terminal is studied and developed, which is used for capturing detailed load information through ILM approach. In practice, for a large number of users in distribution network, the installation cost of intelligent terminals is excessive. Consequently, a “point-to-area” concept is developed by selecting typical users based on user clustering results.

Combined with load identification theory and intelligent terminal devices, a data-driven online load modeling approach is presented. Detailed processes are stated as follows:

1. For numerous users in distribution network, the daily load data are adapted for user clustering, which can be provided by intelligent meter reading system. The dimension of daily load data is reduced and all users are clustered consequently. The result of user clustering is updated for several days interval.

2. In each cluster, several users that are closest to cluster center are designated as “typical users” to stand for the power consumption characteristics of the cluster. Specially, the amount of selected “typical users” can be determined according to practical scenarios.

3. The intelligent terminals are adopted for intrusive measurement to capture the real-time power information of load components in “typical users,” which is wirelessly uploaded to control center for a short time interval, such as 10 min.

4. In control center, load information is conducted for sample training and feature matching to achieve load type identification, and corresponding model parameters are matched for different types of load components. Consequently, the load models of “typical users” are attained.

5. The load models of other users in each cluster are substituted by “typical users,” which consist of static ZIP models and IM models generally. Above process is summarized as “point-to-area” concept.

6. The aggregation models of static ZIP loads and IM loads in all users are derived respectively by model parameters equivalence.

7. Considering the impact of distribution transformer and lines between different voltage grades, the composite load model is built by bottom-to-up stepwise equivalence, which is updated for 10-min interval.

The flow chart of proposed load modeling approach is depicted in Figure 2.

Figure 2.

Diagram of user clustering and online load modeling.

Based on above analysis, the contributions of this article are mainly reflected in three aspects:

1. A data-driven approach is proposed and applied for user clustering, and based on which a “point-to-area” concept is developed for mining the power consumption characteristics from typical users.

2. A bottom-to-up stepwise equivalence approach is proposed, by converting static ZIP model paralleled with IM model to high-voltage side from low-voltage side.

3. The time-variation of load composition can be more accurately reflected by proposed online aggregated load modeling approach.

The following contents of this article are arranged as follows: section “Data-driven user clustering approach” proposes a data-driven user clustering approach based on dimension-reduced load pattern indexes. Section “Online load measurement and aggregated load modeling” introduces the online load measurement and aggregated load modeling procedure. Section “Case study” validates the rationality and effectiveness of proposed online load modeling approach in case study. Concluding remarks are addressed in section “Conclusion.”

Data-driven user clustering approach

User’s daily load curve comprises sequential power consumption data retrieved from all appliances; its pattern is determined by operating state of appliances. Users owning similar daily load curves generally comprise nearly identical appliances and reflect similar power consumption behavior. ³⁰ A new user clustering approach originated from data mining is proposed, by reducing the dimension of user’s daily load data. First, six-dimensional load pattern indexes are applied for dimension reduction. Then, based on the expert weights of load pattern indexes, k-means algorithm is utilized for user clustering. Finally, the expert weights are revised according to clustering validity test, so as to achieve the optional clustering results.

Dimension reduction of daily load data

The load pattern describes the consumption behavior and characteristics of load components, which can be expressed by load pattern index. Six-dimensional load pattern indexes are appointed from four different periods: whole day, peak period, smooth period, and valley period, so as to reflect user’s daily power consumption characteristics. The definitions of load pattern indexes are shown in Table 2.

Table 2.

Daily load pattern indexes.

Period		Index	Definition
Whole day	00:00–24:00	Load rate	a 1 = Pav/Pmax
		Maximum-hours rate	a 2 = Psum/(24Pmax)
		Peak-valley rate	a 3 = (Pmax – Pmin)/Pmax
Peak period	08:00–11:00	Load rate in peak period	a 4 = Pav.peak/Pav
	18:00–21:00
Smooth period	06:00–08:00	Load rate in smooth period	a 5 = Pav.sh/Pav
	11:00–18:00
	21:00–22:00
Valley period	22:00–24:00	Load rate in valley period	a 6 = Pav.val/Pav
	00:00–06:00

The indexes a₁, a₂, and a₃ reflect the power consumption characteristics in whole day from three aspects of average load rate, time utilization efficiency, and peak-shaving capacity. The indexes a₄, a₅, and a₆ express the load rate in different periods of peak, smooth, and valley period. Proposed 6 kinds daily load pattern indexes can be used for reducing dimensions of daily load data matrix. The significance and importance of six-dimensional load pattern indexes are different in revealing load characteristics. Its weight vector is denoted as W  =[w₁, w₂, w₃, w₄, w₅, w₆], which satisfies $\sum _{j=1}^6{w}_j=1$ .

User clustering algorithm

Several approaches are available for users clustering, among which k-means, hierarchical clustering, and fuzzy C-means algorithm are maturely applied.^31–33 Particularly, k-means algorithm has become the most widely applied clustering approach in science and industry field, since its advantage of fast computing and adaptation for big data mining. Consequently, k-means algorithm is employed for user clustering in our research. The procedure of user clustering is as follows:

1. For N users’ pattern indexes matrix, in which K users are randomly arranged as the initial clustering center C_j , where C_j is the center of cluster-j, by denoted as C_j  =[C_j1, C_j2, C_j3, C_j4, C_j5, C_j6] (j = 1, 2,…, K).

2. For each user y_i  =[y_i1, y_i2, y_i3, y_i4, y_i5, y_i6], the weighted Euclidean distances from y_i to C_j are figured as

   $$d_{\mathrm{ij}}=\sqrt{w_1{|y_{i1}-C_{j1}|}^2+w_2{|y_{i2}-C_{j2}|}^2+\cdots +w_6{|y_{i6}-C_{j6}|}^2}$$

   (1)

3. All users are assigned to the nearest clustering center according to weighted Euclidean distances d_ij, and K clusters are formed correspondingly.

4. The number of users in cluster-j is denoted as m_j, and user-i in cluster-j is marked as y_{i, j} . The clustering center of each cluster is calculated again

   $${\mathit{C}}_j=\frac{1}{m_j}\sum _{i=1}^{m_j}{\mathit{y}}_{i,j}(j=1,2,\dots ,K)$$

   (2)

5. The summation of squared errors for all users’ pattern index matrix is figured out

   $$\epsilon =\sum _{j=1}^K\sum _{i=1}^{m_j}{\Vert {\mathit{y}}_{i,j}-{\mathit{C}}_j\Vert }^2$$

   (3)

6. For iteration times t, if $\epsilon (t)<\epsilon (t-1)$ , turn to step 2. Otherwise, the clustering process is finished.

Clustering validity test

Clustering validity test is an overall assessment of merits and demerits for user clustering, developed for determining the optimal clustering scheme. The Silhouette is a clustering validity test index which is suitable for evaluating the similarity of intra-cluster samples and the diversity of inter-cluster samples. It can be applied to determine the optimal cluster number and evaluate the clustering efficiency.

The pattern index matrix of N users is assigned as K clusters. For user-i in a cluster, D_in(i) is defined as the average distance of pattern indexes between user-i and all other users in this cluster. It is applied for evaluating the similarity of intra-cluster users. Contrary, the average distance from user-i to all other users in different clusters is depicted as D_ex(i). Which is used to describe the dissimilarity of inter-cluster users. The Silhouette index I_sil(i) of user-i is defined as

$$I_{\mathrm{sil}}(i)=\frac{D_{\mathrm{ex}}(i)-D_{\mathrm{in}}(i)}{max(D_{\mathrm{in}}(i),D_{\mathrm{ex}}(i))}$$

(4)

which satisfies −1 ≤ I_sil(i) ≤ 1. The smaller the D_in(i) and the larger the D_ex(i), the value of I_sil(i) is closer to 1. The larger value of I_sil(i) indicates that the cluster-i reflects greater intra-similarity and inter-dissimilarity, so that the cluster-i is regard as an effective cluster. If D_in(i) > D_ex(i), I_sil(i) < 0, the procedure of user clustering is suspected.

The average value of all users’ Silhouette index, denoted as I_silmean, is defined to evaluate the clustering efficiency. The larger value of I_silmean index indicates more beneficial clustering result. The number of clusters K is optimum when I_silmean index reaches its maximum value

$$I_{\mathrm{silmean}}=\frac{1}{N}\sum _{i=1}^N{I}_{\mathrm{sil}}(i)$$

(5)

Weight allocation of load pattern indexes

Considering the different significances of six-dimensional load pattern indexes, it is unreasonable to equally weight for them. Adopting the weight from multiple experts’ experience and averagely integrating them is regarded as the most commonly used method. However, the subjectivity of experts will affect clustering effect in some degree. A new user clustering approach is proposed by revising expert weights based on the clustering validity test index. The detailed procedures are as follows:

• 1.n experts are respectively invited to empower weights for six-dimensional load pattern indexes. The weight of expert-m is denoted as W _m  =[w_m1, w_m2, w_m3, w_m4, w_m5, w_m6], which satisfies $\sum _{j=1}^6{w}_{\mathrm{mj}}=1$ .

• 2.Applying the weight of expert-m for load pattern indexes, and k-means algorithm is utilized for user clustering. Correspondingly, the optimal cluster number $K_m^{*}$ and clustering validity test index I_silmean are obtained.

• 3.Figuring out the mean value of optimal cluster numbers from n experts and marked as ${\overline{K}}^{*}$ . For each expert weight-m, the corresponding optimal cluster number $k_{*}^m$ is compared with ${\overline{K}}^{*}$ . As a result, α group weight schemes are invalid if it reflects obvious deviation ( $|{\overline{K}}^{*}-K_m^{*}|>1$ ). The weight matrix of remaining n – α experts is depicted as

  $$\mathit{W}=\left[\begin{array}{c}{\mathit{W}}_1\hfill \\ \multicolumn{1}{c}{{\mathit{W}}_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{{\mathit{W}}_{n-\alpha }}\end{array}\right]$$

  (6)

• 4.Calculating the impact factor of n – α experts based on I_silmean index. Consequently, the impact factor matrix is formed as β  =[β₁, β₁,…, β_n – α]. Where the impact factor β_l of expert weight-l is figured as

$${\beta }_l=\frac{I_{\mathrm{silmean}}(l)}{\sum _{i=1}^{n-\alpha }{I}_{\mathrm{silmean}}(i)}(l=1,2,\dots ,n-\alpha )$$

(7)

Apparently, the larger value of I_silmean, the larger the impact factor β, and which is corresponding to the better clustering result.

• 5.Considering the impact factor of expert’s weight, the revised weight matrix of six-dimensional pattern index is obtained as

  $${\mathit{W}}^{*}=\mathbf{\beta }\mathit{W}\mathbf{=}\mathbf{[}{\mathbf{w}}_1^{\mathbf{*}}\mathbf{,}{\mathbf{w}}_2^{\mathbf{*}}\mathbf{,}{\mathbf{w}}_3^{\mathbf{*}}\mathbf{,}{\mathbf{w}}_4^{\mathbf{*}}\mathbf{,}{\mathbf{w}}_5^{\mathbf{*}}\mathbf{,}{\mathbf{w}}_6^{\mathbf{*}}\mathbf{]}$$

  (8)

• 6.Applying the revised weights to k-means clustering algorithm, the meanings of load pattern indexes can be reflect to a greater extent. Consequently, more reasonable clustering result is realized.

Online load measurement and aggregated load modeling

Online load measurement through intelligent terminal

Served as the online intrusive load measurement device in our research, the dispersed intelligent terminals are connected to Routers for wireless data transmission by In our research, numerous dispersed intelligent terminals are connected to Routers by wireless WiFi or ZigBee optionally, so as to capture and transmit load data. Consequently, the electricity information of load components is uploaded to control center via ISP (Internet service provider) by wired or wireless broadband, which can be displayed through Web interface. The network topology diagram of load information capturing and uploading is depicted in Figure 3.

Figure3.

The network topology diagram of load information capturing and uploading.

The intelligent terminal is composed of five parts including sampling and filter module, data processing module, communication module, execution module, and power module. The relationship of five parts is depicted in Figure 4.

Figure 4.

Diagram of intelligent terminal composition.

The sampling and filter module capture the voltage, current, and frequency from appliances, which is converter to low-voltage signal and provided to data processing module. As the center of intelligent terminal, the data processing module controls the whole system operation in hardware framework. Meanwhile, its software functions involve parameter calculation, data communication, advanced protection, and instruction execution. The communication module is one of interaction channel for intelligent terminal and control server, which supports WiFi and ZigBee signal optionally. The execution module is responsible for executing the control instructions from data processing module and feedback execution result. The power module is the energy source of entire intelligent terminal module, which converts 220 V AC voltage to 5 V DC voltage for system operation. Moreover, environmental parameters acquisition interface is reserved according to application requirements, which contains two analog input and one serial communication.

The electricity information in control center is conducted for sample training and feature matching so as to identify the load type, and typical model parameters are matched for different types of load components. The daily power curves of several typical appliances in households collected by intelligent terminal are shown in Figure 5. Obviously, the daily power curve of different appliance reflects distinct-different power consumption behavior, for example, the fridge reveals periodic characteristic and the microwave shows pulsing characteristic.

Figure 5.

Daily power curves of several typical appliances.

Aggregation of static model and IM model

Static load model is suitable for depicting static components such as resistive loads and light bulbs, which demonstrates the relations between load power with voltage and frequency. The frequency term is usually ignored since that frequency variation has little impact on load power. The polynomial static load model recommended by IEEE Task Force is adopted in this article as follows

$$\{\begin{array}{c}P=P_0[a{(U/U_0)}^2+b(U/U_0)+c]\hfill \\ \multicolumn{1}{c}{Q=Q_0[\alpha {(U/U_0)}^2+\beta (U/U_0)+\gamma ]}\end{array}$$

(9)

where a, b, c are the active power coefficients of constant-impedance, constant-current and CP in static load, which satisfy a + b + c = 1. Similarly, α, β, γ are the reactive power coefficients and α + β + γ = 1. U₀ is the rated load voltage; P₀ and Q₀ represent the active and reactive load power in the rated state.

Aggregating static load can be generalized as the equivalence of load model coefficient P₀, a, b, c and Q₀, α, β, γ. The active power and reactive power are weighted according to the coefficients of constant-impedance, constant-current, and CP part, respectively.

IM load occupies a large proportion in industrial users. The compressor load like air conditioning and refrigerators are typical IM components in commercial users and residential users. There are many methods to aggregate the IM load, almost all of which can be summarized as weighting its parameters based on rated capacity in some degree. The schematic structure of induction motor load model is shown in Figure 6, where R_s is the stator resistance, X_s is the stator leakage reactance, X_m is the excitation reactance, R_r is the rotor resistance, X_r is the rotor leakage reactance, s is the rotor slip.

Figure 6.

Equivalent circuit of induction motor load model.

Generally, the transient process of stator winding is much faster than rotor winding in IM. Therefore, the third-order electromechanical transient model of IM can be retrieved from fifth-order electromagnetic transient model by ignoring the transient characteristics of stator winding, which is shown as follows

$$\{\begin{array}{c}\frac{\mathrm{d}\stackrel{·}{E}\prime }{\mathrm{d}t}=-\frac{1}{{T\prime }_0}[\stackrel{·}{E}\prime -\mathrm{j}(X-X\prime )\stackrel{·}{I}]+\mathrm{j}{\omega }_0(\omega -1)\stackrel{·}{E}\prime \hfill \\ \multicolumn{1}{c}{\frac{\mathrm{d}\omega }{\mathrm{d}t}=\frac{1}{2H}(T_{\mathrm{E}}-T_{\mathrm{M}})}\\ \multicolumn{1}{c}{\stackrel{·}{U}=\stackrel{·}{E}\prime +(R_{\mathrm{s}}+\mathrm{j}X\prime )\stackrel{·}{I}}\end{array}$$

(10)

where $T{\prime }_0=(X_{\mathrm{r}}+X_{\mathrm{m}})/{\omega }_0{R}_{\mathrm{r}}$ is the transient open-circuit time constant; $X=X_{\mathrm{s}}+X_{\mathrm{m}}$ is the rotor open-circuit reactance; $X\prime =X_{\mathrm{s}}+X_{\mathrm{m}}{X}_{\mathrm{r}}/(X_{\mathrm{m}}+X_{\mathrm{r}})$ is the transient reactance and is equal to the blocked-rotor reactance; $T_{\mathrm{E}}=E{\prime }_{\mathrm{d}}{I}_{\mathrm{d}}+E{\prime }_{\mathrm{q}}{I}_{\mathrm{q}}$ is the electromagnetic torque; I_d and I_q are d-axis and q-axis stator current; $E{\prime }_{\mathrm{d}}$ and $E{\prime }_{\mathrm{q}}$ are respectively d-axis and q-axis transient electromotive force (EMF); $T_{\mathrm{M}}=T_0(A{\omega }^2+B\omega +C)$ is the mechanical torque; T₀ is the steady-state mechanical torque; A, B, C are mechanical torque coefficients; ω is the rotor rotation speed; and H is the rotor inertia constant.

The rated capacity of aggregated IM load is the summation of all individuals

$$S_{\mathrm{NM}}=\sum _{i=1}^k{S}_{\mathrm{N}i}$$

(11)

where k is the number of IM loads. By setting ρ_i = S_Ni/S_NM, the equivalent circuit parameters of aggregated IM load are the weighted average value of admittance in branches. The schematic of aggregating process is shown in Figure 7.

$$\frac{1}{Z_{\mathrm{M}}}=\sum _{i=1}^k{\rho }_i\frac{1}{Z_i}$$

(12)

where Z_M is the electrical branch impedance of aggregated IM load. For stator branch, Z = R_s + jX_s; for excitation branch, Z = jX_m; for rotor branch, Z = R_r/s + jX_r. The aggregated inertia time constant is

$$H_{\mathrm{M}}=\sum _{i=1}^k{\rho }_i{H}_i$$

(13)

Figure 7.

The aggregating process schematic of induction motor equivalent circuit.

Stepwise equivalence of composite load model

The load model is generally located on 220 kV/110 kV feeders in power system simulation. Nevertheless, the practical power load is impossibly supplied by 220 kV/110 kV bus directly. Therefore, the impact of distribution network should be further considered in load modeling process. The procedure of converting parallel IM load and static load model to high-voltage side is shown in Figure 8.

Figure 8.

Diagram of converting load model to high-voltage side through distribution network.

In Figure 8, IM and ZIP represent the induction motor load and static load; $\stackrel{·}{I}$ is the current on distribution line; ${\stackrel{·}{I}}_{\mathrm{ML}}$ and ${\stackrel{·}{I}}_{\mathrm{SL}}$ are the current of IM load and static load at low-voltage side; ${\stackrel{·}{U}}_{\mathrm{L}}$ and ${\stackrel{·}{U}}_{\mathrm{H}}$ are the voltage of low-voltage bus and high-voltage bus respectively; P_L and Q_L are the active power and reactive power inflowing to low-voltage bus; Z_D is the equivalent impedance of distribution network. The equivalent equations of IM load model connected to high-voltage side are depicted as

$$\{\begin{array}{c}\frac{\mathrm{d}{\stackrel{·}{E}}_{\mathrm{H}}\prime }{\mathrm{d}t}=-\frac{1}{{T\prime }_{0\mathrm{H}}}[{\stackrel{·}{E}}_{\mathrm{H}}\prime -\mathrm{j}(X_{\mathrm{H}}-X_{\mathrm{H}}\prime ){\stackrel{·}{I}}_{\mathrm{MH}}]+\mathrm{j}{\omega }_0(\omega -1){\stackrel{·}{E}}_{\mathrm{H}}\prime \hfill \\ \multicolumn{1}{c}{\frac{\mathrm{d}\omega }{\mathrm{d}t}=\frac{1}{2H_{\mathrm{H}}}(T_{\mathrm{EH}}-T_{\mathrm{MH}})}\\ \multicolumn{1}{c}{{\stackrel{·}{U}}_{\mathrm{H}}={\stackrel{·}{E}}_{\mathrm{H}}\prime +(R_{\mathrm{sH}}+\mathrm{j}{X}_{\mathrm{H}}\prime ){\stackrel{·}{I}}_{\mathrm{MH}}}\end{array}$$

(14)

The relevant model parameters in equation (14) are derived from equations (15) and (16). Detail computation process is deduced at Appendix 1.

$$\{\begin{array}{c}{R}_{\mathrm{sH}}+\mathrm{j}{X}_{\mathrm{sH}}=(R_{\mathrm{s}}+{\lambda }_1)+\mathrm{j}(X_{\mathrm{s}}+{\lambda }_2)\hfill \\ \multicolumn{1}{c}{R_{\mathrm{rH}}+\mathrm{j}{X}_{\mathrm{rH}}=R_{\mathrm{r}}+\mathrm{j}{X}_{\mathrm{r}}}\\ \multicolumn{1}{c}{X_{\mathrm{mH}}=X_{\mathrm{m}},{T\prime }_{0\mathrm{H}}={T\prime }_0}\\ \multicolumn{1}{c}{H_{\mathrm{H}}={\lambda }_3^{2}H,T_{\mathrm{EH}}={\lambda }_3^2{T}_{\mathrm{E}},T_{\mathrm{MH}}={\lambda }_3^2{T}_{\mathrm{M}}}\end{array}$$

(15)

$$\{\begin{array}{c}\frac{Z_{\mathrm{D}}}{Z_{\mathrm{D}}{Y}_{\mathrm{S}}+1}={\lambda }_1+\mathrm{j}{\lambda }_2\hfill \\ \multicolumn{1}{c}{Z_{\mathrm{D}}{Y}_{\mathrm{S}}+1={\stackrel{·}{\lambda }}_3}\end{array}$$

(16)

where Y_S = Y_Z + Y_I + Y_P, Y_Z is the equivalent admittance of constant-impedance load; Y_I and Y_P are the equivalent admittance of constant-current load and CP load respectively, which satisfy $Y_{\mathrm{I}}={\stackrel{·}{I}}_{\mathrm{const}}/{\stackrel{·}{U}}_{\mathrm{L}}$ , $Y_{\mathrm{P}}=(P_{\mathrm{const}}-\mathrm{j}{Q}_{\mathrm{const}})/U_{\mathrm{L}}^2$ . The value of which dynamically varied with the bus voltage perturbation. ${\stackrel{·}{I}}_{\mathrm{const}}$ is the current of constant-current load; P_const and Q_const are the active and reactive power consumed by CP load.

The equivalent impedance of static load at high-voltage side converted from low-voltage side can be described as

$$\begin{array}{c}\multicolumn{1}{c}{Y_{\mathrm{H}}=\frac{\stackrel{·}{I}-{\stackrel{·}{I}}_{\mathrm{MH}}}{{\stackrel{·}{U}}_{\mathrm{H}}}=Y_{\mathrm{S}}(Z_{\mathrm{D}}{Y}_{\mathrm{S}}+1)}\\ \multicolumn{1}{c}{-Z_{\mathrm{D}}{Y}_{\mathrm{S}}(Z_{\mathrm{D}}{Y}_{\mathrm{S}}+2)\frac{P-\mathrm{j}Q}{U_{\mathrm{H}}^2}}\end{array}$$

(17)

Y _H is related to the impedance of distribution system, initial voltage and power at high-voltage side additionally. The equivalent admittance of constant-impedance, constant-current, and constant-power static load can be separated from Y_H according to its ratio in static load model of low-voltage side.

Case study

User clustering

The data processing and user clustering procedure are developed in MATLAB R2017a. It comprises four major parts: data cleaning, data dimension reduction, data clustering, and clustering validity test. Specifically, the fuzzy clustering and data analysis toolbox in MATLAB R2017a are employed in simulation.

The daily load data of 1800 users in a regional power grid are derived in our research. After data cleaning, the remaining 1776 users are adopted as clustering object. And proposed six-dimensional load pattern indexes are weighted by four experts as shown in Table 3.

Table 3.

Expert weights of six-dimensional load pattern indexes.

Expertnumber	w 1	w 2	w 3	w 4	w 5	w 6
1	0.167	0.167	0.167	0.167	0.167	0.167
2	0.16	0.15	0.15	0.18	0.18	0.18
3	0.20	0.15	0.15	0.20	0.10	0.20
4	0.25	0.10	0.20	0.20	0.10	0.15

Based on four groups of expert weight schemes, all users are clustered by k-means algorithm. The value of optimal cluster number and corresponding I_silmean index under each expert weight are figured and shown in Table 4. Four groups of weights satisfy requirement and the optimal cluster numbers do not reflect major deviation. Based on the impact factor of four groups’ expert weight, the revised weight of six-dimensional load pattern indexes is figured by equation (8), shown in Table 5.

Table 4.

K-means clustering results under each expert weight.

Weights	Optimalcluster number	I silmeanindex	Impact factor β
Expert weight 1	3	0.377	0.206
Expert weight 2	3	0.428	0.234
Expert weight 3	3	0.529	0.289
Expert weight 4	3	0.496	0.271

Table 5.

Revised weight of six-dimensional load pattern indexes.

Revisedweight	w 1	w 2	w 3	w 4	w 5	w 6
Value	0.197	0.140	0.167	0.189	0.132	0.175

Revised weights of load pattern indexes are applied for k-means clustering. The index I_silmean is employed for clustering validity test, so as to achieve the most suitable user clustering scheme. Values of I_silmean indexes are figured out while setting different cluster numbers, as shown in Figure 9. It is apparent that I_silmean value can reach its maximum while all users are clustered for three clusters, in which case the clustering efficiency is optimal.

Figure 9.

The values of I_silmean index under different cluster numbers.

There are respectively 763, 635, and 378 users in aforementioned three clusters through k-means clustering algorithm; the power curve sets in three clusters are shown in Figure 10.

Figure 10.

Power curve sets in three clusters.

Power curves in same cluster reflect similar characteristics. In each cluster, typical users that are closest to cluster center can be designated for intrusive load measurement and identification by deployed intelligent terminals. Consequently, the load model of all users can be retrieved real-timely based on proposed “point-to-area” concept, the procedure of which is depicted in Figure 2.

Online aggregated load modeling

In order to test the accuracy performance of proposed approach for load modeling, a simulation system is designed and built in PSCAD/EMTDC X4.5, an effective power system simulation package. As shown in Figure 11, a distribution network is connected to node 6 of Western System Coordinating Council (WSCC) 9.

Figure 11.

Simulation system of distribution network connected with WSCC 9.

In Figure 11, the composite loads under 220V bus are assumed to denote the aggregated model of three type users, marked as load 1, load 2, and load 3. The load under 10 kV bus is a large-capacity industrial load, marked as load 4. Moreover, X_T1 and X_T2 are the reactance of 110 kV/10 kV and 10 kV/220 V transformers respectively; X_T1 = 0.03 p.u. and X_T2 = 0.02 p.u.

In Table 6, R_s, X_s, X_m, R_r, X_r, A, B, H are the basic parameters of IM load; Z%, I%, and P% represent the ratio of constant-impedance, constant-current, and CP in static load respectively; η_s is the power factor of static load; Z_D = R_D + jX_D is the impedance of distribution line; k_M is the ratio of IM load; above parameters are per-unit values. S_N(MW) is the rated capacity of composite load.

Table 6.

Parameters of load models.

Loads	R s	X s	X m	R r	X r	A	B	H	Z%	I%	P%	η s	R D	X D	k M	S N
Load 1	0.023	0.126	3.39	0.014	0.126	0.85	0.15	1.07	0.33	0.32	0.35	0.85	0.002	0.042	0.35	12
Load 2	0.032	0.096	2.69	0.032	0.096	1	0	0.50	0.20	0.50	0.30	0.85	0.001	0.04	0.20	10
Load 3	0.083	0.095	2.1	0.046	0.095	1	0	0.47	0.20	0.55	0.25	0.85	0.001	0.04	0.45	8
Load 4	0.018	0.117	3.6	0.009	0.117	1	0	1.40	0.10	0.85	0.05	0.85	0.003	0.04	0.60	20

The following two load modeling approaches are imposed as a comparison:

1. Building the composite load model on 110 kV bus by proposed bottom-to-up aggregation and stepwise equivalent approach. Detailed processes are (1) converting load 1, load 2, and load 3 to 220 V bus via distribution lines; (2) aggregating three type loads on 220 V bus; (3) converting the aggregated model on 220 V bus to 10 kV bus via transformer; (4) converting load 4–10 kV bus via distribution line; (5) aggregating all loads on 10 kV bus; (6) converting the aggregated model on 10 kV bus to 110 kV bus via transformer, and obtaining the composite load model on 110 kV bus finally.

2. Under system disturbance, capturing the voltage and power data on 110 kV bus directly and identifying the composite load model parameters by measurement-based approach, which can be summarized as integral identification.

The parameters of composite load model derived from above two approaches are shown in Table 7.

Table 7.

Parameters of built composite load model.

Approaches	R s	X s	X m	R r	X r	H	A	B	S M	k M	P 0	a	b	c	Q 0	α	β	γ
Proposed modeling approach	0.034	0.152	3.166	0.020	0.113	1.05	0.96	0.04	20.4	0.43	25.15	0.22	0.52	0.26	16.45	0.22	0.52	0.26
Integral identification approach	0.027	0.116	3.302	0.019	0.116	1.25	0.87	0.06	21.1	0.41	24.62	0.24	0.47	0.29	16.01	0.25	0.46	0.29

To validate the derived load models, a short-time grounding fault is imposed at bus 5. Subsequently, by considering the randomness and time-variation, an IM load in aggregated IM load 1 is switched off at 8 min. The real-time update can be achieved by proposed online load modeling, set as 10-min interval. In practical power system, the integral identification approach updates the load model for several months interval generally, and the variation of load components cannot be monitored real-timely. Table 8 shows the parameters of composite load model built by two approaches after load composition varied. Acted as the most important parameter in composite load model, the proportions of IM are compared, as shown in Table 9.

Table 8.

Parameters of composite load model after load composition varied.

Approaches	R s	X s	X m	R r	X r	H	A	B	S M	k M	P 0	a	b	c	Q 0	α	β	γ
Proposed modeling approach	0.034	0.150	3.174	0.021	0.116	1.09	0.98	0.02	19.2	0.414	24.72	0.24	0.50	0.26	16.21	0.24	0.50	0.26
Integral identification approach	0.027	0.116	3.302	0.019	0.116	1.25	0.87	0.06	21.1	0.41	24.62	0.24	0.47	0.29	16.01	0.25	0.46	0.29

Table 9.

Proportions of IM load at different simulation stages.

	Proportions of IM load
	0–8 min	10–15 min
Real load model	43.6%	41.9%
Proposed modeling approach	43.0%	41.4%
Integral identification approach	41.0%	41.0%

IM: induction motor.

The composite load model built by two approaches are applied for simulation and compared with time-domain simulation. The voltage and active power trajectories of 110 kV bus under different simulation conditions are shown in Figures 12 and 13.

Figure 12.

Voltage trajectories of 110 kV bus under different simulation conditions.

Figure 13.

Active power trajectories of 110 kV bus under different simulation conditions.

Based on above simulation results, the average errors between time-domain simulation approach and two load modeling approaches are calculated to evaluate the accuracy of derived load model, as shown in Table 10.

Table 10.

The average errors of two load modeling approaches.

Approaches	Average errors
	Voltage	Active power
Proposed modeling approach	1.88%	2.17%
Integral identification approach	3.16%	3.54%

Results indicated that the composite load model derived from proposed approach reflects lower errors. The errors in voltage and active power simulations are 1.88% and 2.17% respectively. Moreover, more accurate load dynamic characteristic is demonstrated; for example, the minimum values of active power under system disturbance are 34.9, 35.2, and 35.8 MW, for time-domain simulation, proposed modeling approach, and integral identification approach. Specially, the variation of load composition triggers a disturbance at 8 min, which can be rapidly responded at 10 min by proposed approach.

Conclusion

The complexity and time-variation have brought great challenges for load modeling practice. In this article, a data-driven online aggregated load modeling approach is proposed and studied. Based on the big data originated from electricity consumers, a data-driven user clustering algorithm is developed by dimension reduction of load pattern indexes. Detailed information of load components are captured by dispersed intelligent terminals in typical users, so as to identify the load composition. Specially, a “point-to-area” concept is innovatively proposed for deriving all users’ load models from regional typical users. The composite load model is achieved by hierarchical aggregation and stepwise equivalence. Crucially, the time-variation of detailed load composition can be detected real-timely by proposed load modeling approach, which is unavailable for traditional load modeling approach.

Actually, proposed online load modeling approach is still at the stage of theoretical study. For future work, it is expected to be applied in practical power system, to fit actual power consumption features and validate the effectiveness. Moreover, future research will involve further exploration on demand response based on big data collected by intelligent terminals. In particular, these data are employed to aggregate users and each cluster denotes different power consumption behaviors, which can be regarded as a virtual agent to response power variation proactively.