Bayesian Estimation on Load Model Coefficients of ZIP and Induction Motor Model

Abstract

Parameter identification in load models is a critical factor for power system computation, simulation, and prediction, as well as stability and reliability analysis. Conventional point estimation based composite load modeling approaches suffer from disturbances and noises, and provide limited information of the system dynamics. In this work, a statistics (Bayesian Estimation) based distribution estimation approach is proposed for both static and dynamic load models. When dealing with multiple parameters, Gibbs sampling method is employed. The proposed method samples all parameters in each iteration and updates one parameter while others remain fixed. The proposed method provides a distribution estimation for load model coefficients and is robust for measuring errors. The proposed parameter identification approach is generic and can be applied to both transmission and distribution networks. Simulations using a 33-feeder system illustrated the efficiency and robustness of the proposal.

1. Introduction

1.1. Composite Load Modeling Methods

Load modeling is a traditional topic in power systems and has been widely studied over the years [1,2,3,4,5,6,7,8]. The importance of an accurate model for voltage stability studies is established in [9]. In [10,11], the authors discussed load modeling in the context of planning, operation, and control systems. The emergence of smart grid technologies, such as distributed generation (DG) and load/demand side management, has created the need for appropriate load models in order to further improve the accuracy and reliability of the power system [12,13].

Load models can be categorized into static models and dynamic models. Static models include the ZIP (constant impedance (Z), constant current (I), and constant power (P)) model, the exponential model, and the frequency dependent model. Examples of widely used dynamic models include the induction motor (IM) and the exponential recovery load model (ERL) [1]. The research focus to date has been on composite load modeling, which consists of static and dynamic models [1]. Specifically, in [1], it is shown that the ZIP+IM model can be used under various conditions and in different locations and compositions. A simplified composite model is shown in Figure 1.

Component-based and measurement-based approaches are two typical load model parameter identification methods that are used today [11,14]. The component-based approach requires information of load compositions. As introduced in [1], this approach requires models of individual components, component compositions (consumption of different load type in percentage), as well as class composition (residential, industrial or commercial load percentage). Measurement-based approaches calculate the parameters using algorithms such as least-squares (LS) and genetic algorithm (GA) based on selected model structure by minimizing the objective function as follows.

$$\min \frac{1}{n}\sum _{i=1}^n[{(P_i^{\mathrm{m}}-P_i^{\mathrm{e}})}^2+{(Q_i^{\mathrm{m}}-Q_i^{\mathrm{e}})}^2]$$

(1)

where superscript “m” indicates the measured value, and “e” presents the estimated value. Measurements should be obtained under different conditions and disturbances. Different techniques have been used in [14,15,16] to identify the parameters for the composite model. Additional methods for load model parameter identification include the Hammerstein model based on real-time system identification for dynamic load [17] and the artificial neural network (ANN), which is based on modeling proposed in [18].

1.2. Motivations for Using Bayesian Estimation

Traditional measurement based parameter estimation approaches do not consider the robustness of the model, and therefore the model sometimes may fail to present the characteristic of the load accurately due to its time-varying feature. Thus, in [19], the robust time-varying load parameter identification approach uses batch-model regression to obtain the updated system parameters. Moreover, the measurement based approach is highly dependent on measurements from meters, and so measurement anomalies may affect the robustness of the estimations in both time-varying and time-independent load models. Influence of measurement noises when implementing various parameter identification methods is an issue discussed in [20]. The challenges in component based approaches are that they rely on individual load components and the compositions of these load. In [4], for example, ZIP coefficients, commonly used in electrical appliances, have to be first determined in order to develop an accurate model and before critical values can be applied. Moreover, the process of determining different components involves such factors as characterization of consumers, clustering, and weather. Since the data are highly dependent on the type of day, weather, and service class, it is a computationally wasted approach to get ZIP coefficients for these appliances. Difficulties involved in obtaining the comprehensive load composition information is another fact that needs to be carefully considered when implementing component based approaches.

Bayesian estimation (BE) based composite load parameter identification approach, on the other hand, can successfully overcome the aforementioned disadvantages in measurement-based and component-based approaches. First, BE is a time-varying estimation method. Parameters can be updated according to the requirements. Second, BE does not require the complicated procedure used in component based approaches. The information of the composition of the loads and the coefficients of the appliances are not necessarily needed. Third, BE is a robust estimation method due to its statistic characteristics: measurement anomalies will not significantly affect the results if the number of samples is large enough, and the measurement error will not affect the estimation results since in each step of sampling the estimated parameter is drawn from a posterior distribution calculated according to the data and eventually forms a distribution of the target parameter.

1.3. Bayesian Probability and Gibbs Sampling

Bayesian probability is one of several interpretations of probability theory. To evaluate the probability of a hypothesis, Bayesian probability specifies some prior probability, which is then updated to a posterior probability in the light of new, relevant data (likelihood). The formula of Bayesian estimation can be written as:

$$p(\theta |x)=\frac{p(\theta )·p(x|\theta )}{\sum p(\theta )·p(x|\theta )}$$

where $p(\theta )$ is the prior, $p(\theta |x)$ is the posterior given a observation x, and $p(x|\theta )$ is the likelihood. The denominator is called the normalization factor. Bayesian estimation gives the probability density function (PDF) of a parameter, which is given by some already known observations. It is the product of the combination of prior information as well as observations of happened events (likelihood). Different from a specific number given by ML, a PDF is able to provide more information and can be used in such functions as forecasting and optimization.

Gibbs sampling is an extension of the Monte Carlo Markov Chain method [21], which performs well when there are multiple parameters to identify. Following the principle of Bayesian probability, Gibbs sampling calculates the posterior distribution of one parameter while fixing the value of others. Then, a new sample of the target parameter is drawn from the posterior, and the next new parameter is sampled with the updated sample. The detailed sampling algorithm is shown in Algorithm 1.

Algorithm 1 Gibbs sampling.

Get measurement of $\mathbf{x}$. Draw initial samples ${\mathit{\theta }}^{(0)}\sim q(\mathit{\theta })$, where $q(\mathit{\theta })$ is the prior and $\mathit{\theta }=\{{\theta }_1,{\theta }_2,\cdots ,{\theta }_M\}$. For iteration $k=1,2,\dots ,K$ Calculate $p({\theta }_1|{\theta }_2^{(k-1)},{\theta }_3^{(k-1)},\dots ,{\theta }_M^{(k-1)},\mathbf{x})$ and draw sample ${\theta }_1^{(k)}$, Calculate $p({\theta }_2|{\theta }_1^{(k)},{\theta }_3^{(k-1)},\dots ,{\theta }_M^{(k-1)},\mathbf{x})$ and draw sample ${\theta }_2^{(k)}$,  ⋮ Calculate $p({\theta }_M|{\theta }_1^{(k)},{\theta }_2^{(k)},\dots ,{\theta }_{M-1}^{(k)},\mathbf{x})$ and draw sample ${\theta }_M^{(k)}$. End For The distribution estimate is the histogram of ${\mathit{\theta }}^k$, $k=\underline{K},\cdots ,K$, where K is the preset maximum iteration number, and $\underline{K}$ is the burn-in number. The results ${\mathit{\theta }}^k$, $k=1,\cdots ,\underline{K}-1,$ are burn-in data and discarded.

1.4. Summary of Contributions

The following contributions are made in this work.

• To the best knowledge of the authors, the proposed Bayesian approach has not been investigated in the context of load model identification. Bayesian estimation-based approaches have been used in wind forecasting and state monitoring [1,2,22], as well as load forecasting [3,4,5]. The promising performance has shown plentiful possibility in power systems.
• The proposal is a generic approach. Detailed formulation and derivation of both static (ZIP) and dynamic (IM) models are presented. The proposed algorithm can be applied in both transmission and distribution network context.
• Numerical experiments illustrated that the proposal provides robust distribution estimation despite the existence of noise. Furthermore, the proposal estimation is more accurate in point estimation than conventional algorithms.

The rest of the paper is organized as follows. In Section 2, a ZIP model for steady state is presented, and the Gibbs sampling is derived. In Section 3, an induction motor model in the context of dynamic state is presented. A modified version of Gibbs sampling is proposed. Numerical experiments are presented in Section 4 to illustrate the effectiveness and robustness of the proposed algorithms. Section 5 discusses existing difficulties in applying the proposal to the real-time operation. Section 6 concludes the paper.

2. ZIP Model and Coefficient Estimation

2.1. ZIP Model

The ZIP model describes how the power of load changes as voltage varies in the steady-state condition. The ZIP model is presented as follows.

$$\begin{array}{c}P=P_0({\alpha }_1{\overline{V}}^2+{\alpha }_2\overline{V}+{\alpha }_3),Q=Q_0({\alpha }_4{\overline{V}}^2+{\alpha }_5\overline{V}+{\alpha }_6),\hfill \\ \sum _{i=1}^3{\alpha }_i=\sum _{i=4}^6{\alpha }_i=1,\hfill \end{array}$$

where $\overline{V}={\displaystyle \frac{V}{V_0}}$, P and Q are real and reactive power, and V is voltage magnitude at the load terminal. Variables $P_0$, $Q_0$ and $V_0$ represent the values of the respective variables under the initial operating condition. When voltage V deviates from $V_0$, real and reactive power of the load are assumed to follow a quadratic model. As the real power and reactive power follow the same form of model, we discuss only the real power as an example in the following. The identification of reactive power parameters follows the same procedure.

Define $y[t]={\displaystyle \frac{P[t]}{P_0}}$, $x[t]={\displaystyle \frac{V[t]}{V_0}}$, where $P[t]$ and $V[t]$ are the tth measurements of the real power and voltage magnitude of the ZIP load. The following assumptions are made.

• The measurement noise follows a normal distribution, i.e., $y[t]={\alpha }_{1}x{[t]}^2+{\alpha }_{2}x[t]+{\alpha }_3+\epsilon [t]$, where $\epsilon [t]\sim \mathcal{N}(0,1/\tau )$, and $1/\tau$ is the variance The reasons for making this assumption are: (1) According to the law of large numbers, a normal distribution would be the best one to represent the characteristics of the noise if the number of experiment is large enough. (2) Since normal distribution is a conjugate distribution, it is easier for model parameter updating when implementing Gibbs sampling.
• Total number of n independent and identically distributed (i.i.d.) samples are drawn, namely $(\mathbf{x},\mathbf{y})\triangleq \{(x[1],y[1]),\cdots ,(x[n],y[n])\}$. Thus, the likelihood is

  $$p(\mathbf{x},\mathbf{y}|{\alpha }_1,{\alpha }_2,\tau )\propto \prod _{t=1}^{n}exp[{(y[t]-\mu [t])}^2\tau /2],$$

  (2)

  where $\mu [t]\triangleq {\alpha }_{1}x{[t]}^2+{\alpha }_{2}x[t]+1-{\alpha }_1-{\alpha }_2$ is the mean.

2.2. Gibbs in the ZIP Model

Gibbs sampling can be directly applied in the ZIP model estimation. Assume the prior distributions are as follows.

$$\begin{array}{ccc}\hfill {\alpha }_1& \sim & \mathcal{N}({\mu }_1^{(0)},1/{\tau }_1^{(0)}),\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill {\alpha }_2& \sim & \mathcal{N}({\mu }_2^{(0)},1/{\tau }_2^{(0)}),\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill \tau & \sim & \mathcal{G}(a^{(0)},b^{(0)}).\hfill \end{array}$$

(5)

where the distribution of $\tau$ is a gamma distribution follows $\mathcal{G}(a,b)$. It can be shown that, after each iteration of Gibbs sampling, the post distributions of these three parameters remain in the same form. Gibbs sampling of the ZIP model is stated in Algorithm 2. The detailed derivation of the algorithm can be found in Appendix A.

Algorithm 2 Gibbs Sampling in the ZIP model.

Get measurement of $\mathbf{x},\mathbf{y}$. Draw initial samples ${\alpha }_1\sim \mathcal{N}({\mu }_1^{(0)},1/{\tau }_1^{(0)}),{\alpha }_2\sim \mathcal{N}({\mu }_2^{(0)},1/{\tau }_2^{(0)}),\tau \sim \mathcal{G}(a^{(0)},b^{(0)})$. For iteration $k=0,1,2,\dots ,K$ Draw a sample ${\alpha }_1^{(k+1)}\sim \mathcal{N}({\mu }_1^{(k+1)},1/{\tau }_1^{(k+1)})$, where $$\begin{array}{c}{\mu }_1^{(k+1)}\triangleq {\displaystyle \frac{{\tau }_1^{(k)}{\mu }_1^{(k)}-{\tau }^{(k)}{\sum }_{t=1}^n\left(({\alpha }_2^{(k)}-1-{\alpha }_2^{(k)}x[t]+y[t])(1-x{[t]}^2)\right)}{{\tau }_1^{(k)}+{\tau }^{(k)}{\sum }_{t=1}^n{(x{[t]}^2-1)}^2}},{\tau }_1^{(k+1)}\triangleq {\tau }_1^{(k)}+{\tau }^{(k)}\sum _{t=1}^n{(x{[t]}^2-1)}^2,\hfill \end{array}$$ Draw a sample ${\alpha }_2^{(k+1)}\sim \mathcal{N}({\mu }_2^{(k+1)},1/{\tau }_2^{(k+1)})$, where $$\begin{array}{c}{\mu }_2^{(k+1)}\triangleq {\displaystyle \frac{{\tau }_2^{(k)}{\mu }_2^{(k)}-{\tau }^{(k)}{\sum }_{t=1}^n\left(({\alpha }_1^{(k+1)}-1-{\alpha }_1^{(k+1)}x[t]+y[t])(1-x[t])\right)}{{\tau }_2^{(k)}+{\tau }^{(k)}{\sum }_{t=1}^n{(x[t]-1)}^2}},{\tau }_2^{(k+1)}\triangleq {\tau }_2^{(k)}+{\tau }^{(k)}\sum _{i=1}^n{(x[t]-1)}^2.\hfill \end{array}$$ Draw a sample ${\tau }^{(k+1)}\sim \mathcal{G}(a^{(k+1)},b^{(k+1)})$, where $$\begin{array}{c}{a}^{(k+1)}=a^{(k)}+n/2,b^{(k+1)}=b^{(k)}+\sum _{t=1}^n{\left(y[t]-{\alpha }_1^{(k+1)}x{[t]}^2-{\alpha }_2^{(k+1)}x[t]-(1-{\alpha }_1^{(k+1)}-{\alpha }_2^{(k+1)})\right)}^2/2.\hfill \end{array}$$ End For Discard burn-in data, and the distribution estimate is the histogram samples.

3. Induction Motor Model and Coefficient Estimation

3.1. Induction Motor Model

The dynamic part of load is usually represented by an IM model, which is discussed below [5].

$$\begin{array}{c}\left\{\begin{array}{c}{\displaystyle \frac{d{E}_d^{\prime }}{dt}}=-{\displaystyle \frac{1}{T^{\prime }}}[E_d^{\prime }+(X-X^{\prime })I_q]-(\omega -1)E_q^{\prime }\hfill \\ {\displaystyle \frac{d{E}_q^{\prime }}{dt}}=-{\displaystyle \frac{1}{T^{\prime }}}[E_q^{\prime }-(X-X^{\prime })I_d]+(\omega -1)E_d^{\prime }\hfill \\ {\displaystyle \frac{d\omega }{dt}}=-{\displaystyle \frac{1}{2H}}[(A{\omega }^2+B\omega +C)T_0-(E_d^{\prime }{I}_d+E_q^{\prime }{I}_q)]\hfill \end{array}\right.\hfill \\ \left\{\begin{array}{c}{I}_d={\displaystyle \frac{1}{R_s^2+X^{\prime }{}^2}}[R_s(U_d-E_d^{\prime })+X^{\prime }(U_q-E_q^{\prime })]\hfill \\ I_q={\displaystyle \frac{1}{R_s^2+X^{\prime }{}^2}}[R_s(U_q-E_q^{\prime })-X^{\prime }(U_d-E_d^{\prime })]\hfill \end{array}\right.\hfill \end{array}$$

(6)

where $T^{\prime }\triangleq \frac{X_r+X_m}{R_r}$, $X\triangleq X_s+X_m$, $X^{\prime }\triangleq X_s+\frac{X_m{X}_r}{X_m+X_r}$, $A+B+C=1$. In Equation (6), $R_s$ is motor stator winding resistance, $X_s$ is motor stator leakage reactance, $X_m$ is motor magnetizing reactance, $R_r$ is rotor resistance, $X_r$ is rotor leakage reactance, H is rotor inertia constant, $\omega$ is rotor speed, $I_d$ and $I_q$ are stator current in d-axis and q-axis, $U_d$ and $U_q$ are bus voltage in d-axis and q-axis, and $E_d^{\prime }$ and $E_q^{\prime }$ are stator voltage in d-axis and q-axis. $T_0$ is the initial load torque. The parameters to be identified in an IM model are: $X_r,X_m,X_s,R_r,R_s,A,B,C,H$. Typically, A is assumed to be 1 as the mechanical torque is assumed to be proportional to the square of the rotation speed of the motor [11]. Consequently, B and C are both equal to 0.

The differential equation (Equation (7)) is non-linear and the differentiations of $E_d^{\prime }$, $E_q^{\prime }$, $\omega$ cannot be measured directly. These features make the procedure in ZIP models unsuitable for IM model parameter identification and complicates the estimation process. In this work, the state of IM model is assumed to follow the transitions below.

$$\begin{array}{c}\left\{\begin{array}{c}{y}_{E_d}[t]={\beta }_1{E}_d^{\prime }[t]+{\beta }_2{I}_q[t]-(\omega -1)E_q^{\prime }[t]+{\epsilon }_{E_d}[t]\hfill \\ y_{E_q}[t]={\beta }_1{E}_q^{\prime }[t]-{\beta }_2{I}_d[t]+(\omega -1)E_d^{\prime }[t]+{\epsilon }_{E_q}[t]\hfill \\ y_{\omega }[t]={\beta }_3({\omega }^2{T}_0-E_d^{\prime }[t]I_d[t]-E_q^{\prime }[t]I_q[t])+{\epsilon }_{\omega }[t]\hfill \end{array}\right.\hfill \\ \left\{\begin{array}{c}{I}_d[t]={\alpha }_b(U_d[t]-E_d^{\prime }[t])+{\alpha }_c(U_q[t]-E_q^{\prime }[t])+{\epsilon }_{I_d}[t]\hfill \\ I_q[t]={\alpha }_b(U_q[t]-E_q^{\prime }[t])-{\alpha }_c(U_d[t]-E_d^{\prime }[t])+{\epsilon }_{I_q}[t]\hfill \end{array}\right.\hfill \end{array}$$

(7)

where $y_{Ed}\triangleq {\displaystyle \frac{d{E}_q^{\prime }}{dt}}$, $y_{Eq}\triangleq {\displaystyle \frac{d{E}_d^{\prime }}{dt}}$$y_{\omega }\triangleq {\displaystyle \frac{d\omega }{dt}}$, ${\beta }_1\triangleq -{\displaystyle \frac{1}{T^{\prime }}}$, ${\beta }_2\triangleq -{\displaystyle \frac{X-X^{\prime }}{T^{\prime }}}$, ${\beta }_3\triangleq -{\displaystyle \frac{1}{2H}}$, ${\alpha }_b\triangleq {\displaystyle \frac{R_s}{R_s^2+X^{\prime }{}^2}}$, ${\alpha }_c\triangleq {\displaystyle \frac{X^{\prime }}{R_s^2+X^{\prime }{}^2}}$, ${\epsilon }_{E_d}[t]\sim \mathcal{N}(0,1/{\tau }_E)$, ${\epsilon }_{E_q}[t]\sim \mathcal{N}(0,1/{\tau }_E)$, ${\epsilon }_{\omega }[t]\sim \mathcal{N}(0,1/{\tau }_{\omega })$, ${\epsilon }_{I_d}[t]\sim \mathcal{N}(0,1/{\tau }_I)$, and ${\epsilon }_{I_q}[t]\sim \mathcal{N}(0,1/{\tau }_I)$.

The prior distribution of these parameters is assumed to be as follows.

$$\begin{array}{cc}{\beta }_1\sim \mathcal{N}({\mu }_{{\beta }_1}^{(0)},1/{\tau }_{{\beta }_1}^{(0)}),\hfill & {\beta }_2\sim \mathcal{N}({\mu }_{{\beta }_2}^{(0)},1/{\tau }_{{\beta }_2}^{(0)})\hfill \\ {\beta }_3\sim \mathcal{N}({\mu }_{{\beta }_3}^{(0)},1/{\tau }_{{\beta }_3}^{(0)}),\hfill & {\alpha }_b\sim \mathcal{N}({\mu }_{{\alpha }_b}^{(0)},1/{\tau }_{{\alpha }_b}^{(0)})\hfill \\ {\alpha }_c\sim \mathcal{N}({\mu }_{{\alpha }_c}^{(0)},1/{\tau }_{{\alpha }_c}^{(0)}),\hfill & {\tau }_E\sim \mathcal{G}({\alpha }_E^{(0)},{\beta }_E^{(0)})\hfill \\ {\tau }_{\omega }\sim \mathcal{G}({\alpha }_{\omega }^{(0)},{\beta }_{\omega }^{(0)}),\hfill & {\tau }_I\sim \mathcal{G}({\alpha }_I^{(0)},{\beta }_I^{(0)}).\hfill \end{array}$$

That is, the transition of the IM model is assumed to be disturbed by some Gaussian noise and the prior distribution of the parameters follows Gaussian distributions. In this case, if the measurements of ${\mathbf{y}}_{E_d},{\mathbf{y}}_{E_q},{\mathbf{y}}_{\omega },{\mathbf{E}}_q^{\prime },{\mathbf{E}}_d^{\prime },{\mathbf{I}}_d,{\mathbf{I}}_q,{\mathbf{U}}_d,{\mathbf{U}}_q$ are available, the same procedure as the one in the ZIP model section can be followed and similar posterior distribution can be obtained after some algebra. In practice, however, not all of the measurements are available. Thus, we provide a modified Gibbs sampling of the IM model.

3.2. Gibbs Sampling in IM Model

If measurements of ${\mathbf{y}}_{E_d},{\mathbf{y}}_{E_q},{\mathbf{y}}_{\omega },{\mathbf{E}}_q^{\prime },{\mathbf{E}}_d^{\prime },{\mathbf{I}}_d,{\mathbf{I}}_q,{\mathbf{U}}_d,{\mathbf{U}}_q$ are available, Gibbs sampling of the IM model is similar to the procedure of the ZIP model. The Gibbs sampling of the IM model can be found in Algorithm 3.

However, measurements of the derivative equations, i.e., the values of $y_{E_d},y_{E_q},E_d^{\prime },E_q^{\prime }$, are not available. Thus, an iterative algorithm following the intuition of Gibbs sampling algorithm is proposed in this paper. Starting by using the typical value of corresponding parameters, Equation (6) is solved first, and a set of sampled value of $y_{E_d},y_{E_q},E_d^{\prime },E_q^{\prime }$ is obtained. Then, the Gibbs sampling is implemented to derive posterior distributions of the target parameters. The means of the posterior distributions are used to solve Equation (6) again and generate new samples of $y_{E_d},y_{E_q},E_d^{\prime },E_q^{\prime }$. This iteration continues until the predefined maximum iteration number is reached. The flow chart of the proposed method is shown in Figure 2.

Algorithm 3 Gibbs Sampling in the ZIP model.

Get measurements of ${\mathbf{y}}_{E_d},{\mathbf{y}}_{E_q},{\mathbf{y}}_{\omega },{\mathbf{E}}_q^{\prime },{\mathbf{E}}_d^{\prime },{\mathbf{I}}_d,{\mathbf{I}}_q,{\mathbf{U}}_d,{\mathbf{U}}_q$. Draw initial samples ${\beta }_1\sim \mathcal{N}({\mu }_{{\beta }_1}^{(0)},1/{\tau }_{{\beta }_1}^{(0)})$, ${\beta }_2\sim \mathcal{N}({\mu }_{{\beta }_2}^{(0)},1/{\tau }_{{\beta }_2}^{(0)})$, ${\beta }_3\sim \mathcal{N}({\mu }_{{\beta }_3}^{(0)},1/{\tau }_{{\beta }_3}^{(0)})$, ${\alpha }_b\sim \mathcal{N}({\mu }_{{\alpha }_b}^{(0)},1/{\tau }_{{\alpha }_b}^{(0)})$, ${\alpha }_c\sim \mathcal{N}({\mu }_{{\alpha }_c}^{(0)},1/{\tau }_{{\alpha }_c}^{(0)})$, ${\tau }_E\sim \mathcal{G}({\alpha }_E^{(0)},{\beta }_E^{(0)})$, ${\tau }_{\omega }\sim \mathcal{G}({\alpha }_{\omega }^{(0)},{\beta }_{\omega }^{(0)})$, ${\tau }_I\sim \mathcal{G}({\alpha }_I^{(0)},{\beta }_I^{(0)})$ from typical distributions. For iteration $k=0,1,2,\dots ,K$ Draw a sample ${\beta }_1^{(k+1)}\sim \mathcal{N}({\mu }_{{\beta }_1}^{(k+1)},1/{\tau }_{{\beta }_1}^{(k+1)})$ according to Equation (A6). Draw a sample ${\beta }_2^{(k+1)}\sim \mathcal{N}({\mu }_{{\beta }_2}^{(k+1)},1/{\tau }_{{\beta }_2}^{(k+1)})$ according to Equation (A7). Draw a sample ${\beta }_3^{(k+1)}\sim \mathcal{N}({\mu }_{{\beta }_3}^{(k+1)},1/{\tau }_{{\beta }_3}^{(k+1)})$ according to Equation (A8). Draw a sample ${\alpha }_b^{k+1}\sim \mathcal{N}({\mu }_{{\alpha }_b}^{(k+1)},1/{\tau }_{{\alpha }_b}^{(k+1)})$ according to Equation (A9). Draw a sample ${\alpha }_c^{k+1}\sim \mathcal{N}({\mu }_{{\alpha }_c}^{(k+1)},1/{\tau }_{{\alpha }_c}^{(k+1)})$ according to Equation (A10). Draw a sample ${\tau }_E^{k+1}\sim \mathcal{G}({\alpha }_E^{(k+1)},{\beta }_E^{(k+1)})$ according to Equation (A11). Draw a sample ${\tau }_{\omega }^{k+1}\sim \mathcal{G}({\alpha }_{\omega }^{(k+1)},{\beta }_{\omega }^{(k+1)})$ according to Equation (A12). Draw a sample ${\tau }_I^{k+1}\sim \mathcal{G}({\alpha }_I^{(k+1)},{\beta }_I^{(k+!)})$ according to Equation (A13). End For Discard burn-in data, and the distribution estimate is the histogram samples.

4. Numerical Results

Simulation studies were carried out for ZIP and IM models, separately. In this study, the maximum Gibbs sampling runs K was set to 40,000 and the burn-in length $\underline{K}$ was chosen as 5000 [21].

4.1. ZIP Model Identification

A 33-bus test feeder was used to generate the testing data. Load at bus 18 was replaced by a ZIP model. Randomness was added to loads at other buses by multiplying a random factor drawing from a uniform distribution following $\mathcal{U}[0.1,4.5]$ as shown in Equation (8).

$$P_{w,i}=P_i·\mathcal{U}[0.1,4.5]$$

(8)

where $P_i$ is the original load in the 33-bus test feeder, i indicates the node number, and $P_{w,i}$ is the weighted load at the corresponding node after multiplying a random weight. The changes at each load in every experiment simulated the disturbances and uncertainties in the system, which significantly impacted the voltage and power at the bus of interest.

The ZIP factors of the load at node 18 were assigned as ${\alpha }_1=0.25$, ${\alpha }_2=0.25$, ${\alpha }_3=0.5$. Power flow converged in each iteration with different $P_{w,i}$, and the corresponding voltage $V_{w,i}$ was recorded. The measurement noise $\epsilon$ followed $\mathcal{N}(0,0.1)$, which means a 10% measurement error.

By implementing the proposed GS method, the coefficients of the ZIP model were estimated, as shown in Figure 3.

Different from other estimation approaches, the GS approach can generate a distribution that describes the probability of the real value falling into a certain range. As shown in Figure 3, the mean values of the distributions of ${\alpha }_1$ and ${\alpha }_2$ were found to be 0.25 and 0.249, respectively, and ${\alpha }_3=1-{\alpha }_1-{\alpha }_2$.

Figure 4 and Figure 5 show the voltage and voltage/real power differences, respectively, using the estimated and real model in 100 i.i.d. experiments. The randomness came from the factors multiplied to loads at other buses. The dash line in Figure 4 indicates the voltage at bus 18 using the real ZIP coefficients, and the solid line is the voltage measured at the same bus using estimated coefficients. The comparison shown in Figure 5 is the distributions of voltage and active power differences, $\Delta V$ and $\Delta P$. It can be seen that the differences were in the range of 10⁻⁴ and 10⁻³ p.u., respectively. Figure 6 shows the burn-in process observed. It shows that the process was very quick, and no significant burn-in process that could be observed.

4.2. IM Model Identification

The estimation of the parameters in the IM model followed the same procedure. A dynamic model was built in Matlab/Simulink to generate the data that were used in Equation (7) in a 33 bus test feeder, as shown in Figure 7. The sampling results are shown in Figure 8 and compared with the real data in

Table 1. Estimated and real value of the parameters for IM mode with 5% measurement error.

Para.	Real Value	Est. Value (Mean)	Error (%)
${\beta }_1$	0.0077	0.007683	0.22
${\beta }_2$	0.018	0.01824	1.33
${\beta }_3$	25	24.8	0.8
${\alpha }_b$	0.20	0.211	5.5
${\alpha }_c$	0.80	0.813	1.63

. According to Table 1, the estimation error was less than 6%, with 5% measurement error. It is worth noting that estimation of parameters highly relied on the prior distributions. A good prior could significantly increase the estimation accuracy and shorten the burn-in period. A comparison of the active power response using test parameter and identified parameteris shown in Figure 9. The simulation results also indicate the differences in active power was small and the identified parameters were accurate enough to be used in other analyses.

4.3. Benchmarks

The ZIP and IM model parameters derived by the proposed GS method were compared with least square (LS) [23] and Kalman Filter (KF) methods [24] with 10% noise. The comparison results of the voltage and active power using GS, LS, and KF in ZIP model are shown in Figure 10. The average absolute mean errors of voltage and active power in p.u. are listed in

Table 2. Active power and voltage errors with different methods.

Para. Error	GS (%)	LS (%)	KF (%)
Voltage	0.007	0.062	0.024
Active Power	1.12	4.65	2.64

. As shown in Figure 10, it was obvious that the GS approach had the best performance, while KF method only fell slightly behind. In contrast, LS have the largest error among all three approaches when there was large measurement error. However, KF only worked well with time-invariant parameter estimation [25]. In practice, the ZIP component varies with time due to stochastic consumer behaviors.

The parameters estimation comparisons in the IM model are listed in

Table 3. IM model parameter estimation with different methods.

Para.	GS	LS	KF	Real Value
${\beta }_1$	0.007683	0.0076	0.0077	0.0077
${\beta }_2$	0.01824	0.1375	0.0185	0.018
${\beta }_3$	24.8	18	34	25
${\alpha }_b$	0.211	2.06	1.8	0.2
${\alpha }_c$	0.813	4.04	3	0.8

. The estimation errors from GS were the smallest. The performance of KF was worse than itself in the ZIP model, but still better than the LS method.

5. Discussion

Some difficulties during the identification are discussed.

1

It is assumed that the noise in ZIP and IM models follows a normal distribution. In practice, the true randomness of parameters of load models is not available. By the law of large numbers, a normal distribution would be the best approximation of the randomness if the sample size is large enough. However, the robustness of the proposal needs further investigation when the Gaussian assumption is violated.

2

We are only estimating the parameters of ZIP or IM model. It is necessary to consider multiple ZIP and MI models connected to different locations in the grid. There might be inferences to the estimation.

3

How to identify a composition model composed by ZIP and IM is another interesting topic, which is a future work of this project. Most of the work that has been done in composite load modeling needs to specify the percentage of ZIP by giving the active power load $P_{\mathrm{ZIP}}$.

6. Conclusions

Bayesian estimation based dynamic load modeling is an executable and reliable composite load modeling approach. This method can overcome some disadvantages in the traditional composite load modeling process. It can provide a robust time-varying dynamic load parameter estimation with higher accuracy and better performance. It is worth noting that this method requires a good estimation of each parameter, which can be achieved by using traditional identification methods, such as Least-square, Kalman Filter, Genetic Algorithm, etc. Future work will be focused on implementing those aforementioned approaches to better estimate the prior, parameter identification in composite load model, and model clustering using machine learning approaches. Another important study is evaluating the robustness of the proposal when assumptions are violated via experiments.