Bayesian Estimation Based Parameter Estimation for Composite Load

TL;DR

This paper introduces a Bayesian estimation approach for composite load modeling in power systems, providing robust parameter distribution estimates for static and dynamic loads, improving over traditional point estimation methods.

Contribution

It proposes a Bayesian estimation framework using Gibbs sampling for static and dynamic load models, offering a more robust and informative parameter estimation method.

Findings

Provides distribution estimates of load model coefficients

Robust to measurement errors

Applicable to both static and dynamic load models

Abstract

Accurate identification of parameters of load models is essential in power system computations, including simulation, prediction, and 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 statistic (Bayesian Estimation) based distribution estimation approach is proposed for both static (ZIP) and dynamic (Induction Motor) load modeling. When dealing with multiple parameters, Gibbs sampling method is employed. In each iteration, the proposal samples each parameter while keeps others fixed. The proposed method provides a distribution estimation of load models coefficients and is robust to measurement errors.

Figures

Tables

Table 1. TABLE I: 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

Table 2. TABLE II: Active power & 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

Table 3. TABLE III: 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

Equations

$$\begin{array}[]{l}P=P_{0}(\alpha_{1}\bar{V}^{2}+\alpha_{2}\bar{V}+\alpha_{3}),Q=Q_{0}(\alpha_{4}\bar{V}^{2}+\alpha_{5}\bar{V}+\alpha_{6}),\end{array}$$

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

$$\begin{array}[]{l}\left\{\begin{array}[]{l}dE^{\prime}_{d}/dt=-[E^{\prime}_{d}+(X-X^{\prime})I_{q}]/T^{\prime}-(\omega-1)E^{\prime}_{q}\\ dE^{\prime}_{q}/dt=-[E^{\prime}_{q}-(X-X^{\prime})I_{d}]/T^{\prime}+(\omega-1)E^{\prime}_{d}\\ d\omega/dt=-[(A^{2}\omega+B\omega+C)T_{0}-(E^{\prime}_{d}I_{d}+E^{\prime}_{q}I_{q})]/2H\end{array}\right.\\ \left\{\begin{array}[]{l}I_{d}=[R_{s}(U_{d}-E^{\prime}_{d})+X^{\prime}(U_{q}-E^{\prime}_{q})]/(R_{s}^{2}+X^{\prime 2})\\ I_{q}=[R_{s}(U_{q}-E^{\prime}_{q})-X^{\prime}(U_{d}-E^{\prime}_{d})]/(R_{s}^{2}+X^{\prime 2})\end{array}\right.\end{array}$$

$$\begin{array}[]{l}\left\{\begin{array}[]{l}y_{E_{d}}[t]=\beta_{1}E^{\prime}_{d}[t]+\beta_{2}I_{q}[t]-(\omega-1)E^{\prime}_{q}[t]+\varepsilon_{E_{d}}[t]\\ y_{E_{q}}[t]=\beta_{1}E^{\prime}_{q}[t]-\beta_{2}I_{d}[t]+(\omega-1)E^{\prime}_{d}[t]+\varepsilon_{E_{q}}[t]\\ y_{\omega}[t]=\beta_{3}(\omega^{2}-E^{\prime}_{d}[t]I_{d}[t]-E^{\prime}_{q}[t]I_{q}[t])+\varepsilon_{\omega}[t]\\ \end{array}\right.\\ \left\{\begin{array}[]{l}y_{I_{d}}[t]=\alpha_{b}(U_{d}[t]-E^{\prime}_{d}[t])+\alpha_{c}(U_{q}[t]-E^{\prime}_{q}[t])+\varepsilon_{I_{d}}[t]\\ y_{I_{q}}[t]=\alpha_{b}(U_{q}[t]-E^{\prime}_{q}[t])+\alpha_{c}(U_{d}[t]-E^{\prime}_{d}[t])+\varepsilon_{I_{q}}[t]\end{array}\right.\end{array}$$

$${\alpha }_1$$

$${\alpha }_2$$

$$\tau$$

$$p({\alpha }_1|{\alpha }_2^{(0)},{\tau }^{(0)},\text{x},\text{y})\propto p(\text{x},\text{y}|{\alpha }_1,{\alpha }_2^{(0)},{\tau }^{(0)})p({\alpha }_1)$$

$$\begin{array}[]{l}\log p(\alpha_{1}|\alpha_{2}^{(0)},\tau^{(0)},\textbf{x},\textbf{y})\propto\\ ~{}~{}-\dfrac{\tau_{1}^{(0)}}{2}(\alpha_{1}-\mu_{1}^{(0)})^{2}-\\ ~{}~{}\dfrac{\tau^{(0)}}{2}\sum_{t=1}^{n}\big{(}y[t]-(\alpha_{1}x[t]^{2}+\alpha_{2}^{(0)}x[t]+1-\alpha_{1}-\alpha_{2}^{(0)})\big{)}^{2}\end{array}$$

$$\begin{array}[]{l}-\big{(}\tau_{1}^{(0)}+\tau^{(0)}\sum_{t=1}^{n}(x[t]^{2}-1)^{2}\big{)}\alpha_{1}^{2}/2+\big{(}\tau_{1}^{(0)}\mu_{1}^{(0)}\\ ~{}~{}-\tau^{(0)}\sum_{t=1}^{n}(\alpha_{2}^{(0)}-1-\alpha_{2}^{(0)}x[t]+y[t])(1-x[t]^{2})\big{)}\alpha_{1}\end{array}$$

$$\begin{array}[]{l}\mu_{1}^{(1)}\triangleq\\ \dfrac{\tau_{1}^{(0)}\mu_{1}^{(0)}-\tau^{(0)}\sum_{t=1}^{n}\big{(}(\alpha_{2}^{(0)}-1-\alpha_{2}^{(0)}x[t]+y[t])(1-x[t]^{2})\big{)}}{\tau_{1}^{(0)}+\tau^{(0)}\sum_{t=1}^{n}(x[t]^{2}-1)^{2}},\\ \end{array}$$

$$\begin{array}[]{l}\tau_{1}^{(1)}\triangleq\tau_{1}^{(0)}+\tau^{(0)}\sum_{t=1}^{n}(x[t]^{2}-1)^{2},\end{array}$$

$${\alpha }_1|{\alpha }_2^{(0)},{\tau }^{(0)},{\tau }_1^{(0)},{\mu }_1^{(0)},\text{x},\text{y}\sim \mathcal{N}({\mu }_1^{(1)},1/{\tau }_1^{(1)}).$$

$$\begin{array}[]{l}\mu_{2}^{(1)}\triangleq\\ \dfrac{\tau_{2}^{(0)}\mu_{2}^{(0)}-\tau^{(0)}\sum_{t=1}^{n}\big{(}(\alpha_{1}^{(1)}-1-\alpha_{1}^{(1)}x[t]+y[t])(1-x[t])\big{)}}{\tau_{2}^{(0)}+\tau^{(0)}\sum_{t=1}^{n}(x[t]-1)^{2}},\\ \tau_{2}^{(1)}\triangleq\tau_{2}^{(0)}+\tau^{(0)}\sum_{i=1}^{n}(x[t]-1)^{2}.\end{array}$$

$${\alpha }_2|{\alpha }_1^{(1)},{\tau }^{(0)},{\tau }_1^{(0)},{\mu }_1^{(0)},\text{x},\text{y}\sim \mathcal{N}({\mu }_2^{(1)},1/{\tau }_2^{(1)}).$$

$$\begin{array}[]{l}\log p(\tau|\alpha_{1}^{(1)},\alpha_{2}^{(1)},\textbf{x},\textbf{y})\propto\\ ~{}~{}\dfrac{n}{2}\log\tau-\dfrac{\tau}{2}\sum_{t=1}^{n}\big{(}y[t]-\alpha_{1}^{(1)}x[t]^{2}-\alpha_{2}^{(1)}x[t]-\\ ~{}~{}1+\alpha_{1}^{(1)}+\alpha_{2}^{(1)}+(a^{(0)}-1)\log\tau-b^{(0)}\tau\big{)}.\end{array}$$

$$\begin{array}[]{l}a^{(1)}=a^{(0)}+n/2,\\ b^{(1)}=b^{(0)}+\\ \sum_{t=1}^{n}\big{(}y[t]-\alpha_{1}^{(1)}x[t]^{2}-\alpha_{2}^{(1)}x[t]-(1-\alpha_{1}^{(1)}-\alpha_{2}^{(1)})\big{)}^{2}/2.\end{array}$$

$${\beta }_1$$

$${\beta }_3$$

$${\alpha }_c$$

$${\tau }_{\omega }$$

Full text

I Introduction

Load models can be categorized into static and dynamic models. Conventional static models include ZIP model, exponential model, frequency dependent model, etc., and the induction motor (IM), exponential recovery load model (ERL) [1] are the major dynamic load models used in research studies recently. It is reported that a composite load under ZIP+IM model can be used with various conditions, locations and composition, and it has been widely used in industry and studied in voltage stability as well as planning, operation and control of power systems[1, 2, 3, 4].

Load models are typically identified by component and measurement based approaches in state of the art [1, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. Component based approaches highly rely on characteristics of components and compositions of load. In [7], ZIP coefficients for widely used electrical appliances are determined by experiments, and the overall ZIP model was established with respect to the predetermined appliances. However, the computational cost is high, and the accuracy of individual consumers is greatly affected by non-electric factors such as data and weather [7, 1]. Difficulties in obtaining load composition information is another factor to be considered when implementing this type of methods.

Measurement-based approaches, such as least-squares (LS) and genetic algorithm (GA), are another mainstream for load parameter identifications. Different techniques are reported in [9, 15] to identify the parameters for composite load models. In [16], a robust time-varying load parameters identification approach is proposed by using batch-model regression in order to obtain the updated system parameters. One disadvantage of measurement based approach is the dependence on data quality. Measurement anomalies may affect the robustness of estimation in both time-varying and time-independent load models. Moreover, these approaches estimate the expected value of coefficients, which provides limited information.

Bayesian estimation (BE) [17] based composite load parameter identification approach can successfully overcome the aforementioned disadvantages in measurement and component based approaches. First, BE is a distribution estimation rather than point estimation technique, which provides the likelihood of each parameter. In this case, this method is able to provide accurate estimation of the parameters in both time constant and time-varying cases. Second, BE does not require information of load compositions, nor coefficients of appliances. Third, BE is a robust estimation method due to its statistical characteristics: measurement anomalies will not significantly affect the results if the number of samples is large enough, and the effect of measurement error is also not significant either because the distribution interval contains the real value inside the distribution with expected probabilities.

The rest of the paper is organized as follows. The problem of ZIP and IM model estimation is formulated in Section II. The Bayesian estimation technology and the detailed parameter identification method (Gibbs Sampling) is introduced in Section III. The advantage of the proposal is illustrated in Section IV by simulations. Section V draws the conclusion.

II Formulation

In this section, a composite load model consisting of ZIP and IM model is introduced [2], followed by the proposed parameter estimation method, Gibbs Sampling (GS).

A composite load with ZIP+IM model is shown in Fig. 1, in which the ZIP model describes the steady-state behavior and the IM model corresponds to the dynamic process.

II-A ZIP Model

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

[TABLE]

where ${\sum_{i=1}^{3}{\alpha_{i}}=\sum_{i=4}^{6}{\alpha_{i}}=1}$, $\bar{V}=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 at the initial operating condition. When voltage $V$ deviates from $V_{0}$, real and reactive power of the load is assumed to follow a quadratic model. As the real power and reactive power follow the same form of model, we will discuss only the real power as an example in the following. The reactive power follows the same procedure.

Define ${y[t]=P[t]/P_{0}}$, ${x[t]=V[t]/V_{0}}$, where $P[t]$ and $V[t]$ are the $t$th 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}+\varepsilon[t]}$, where ${\varepsilon[t]\sim\mathcal{N}(0,1/\tau)}$, and $1/\tau$ is the variance111The 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 were drawn, namely ${(\textbf{x},\textbf{y})\triangleq\{(x[1],y[1]),\cdots,(x[n],y[n])\}}$. Thus the likelihood

[TABLE]

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

II-B * IM Model*

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

[TABLE]

where ${X^{\prime}\triangleq X_{s}+X_{m}X_{r}/(X_{m}+X_{r})}$, ${X\triangleq X_{s}+X_{m}}$, ${T^{\prime}\triangleq(X_{r}+X_{m})/R_{r}}$, ${A+B+C=1}$. Here $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^{\prime}_{d}$ and $E^{\prime}_{q}$ 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 [10]. Consequently B and C are both equal to 0. For simplicity, define $y_{Ed}\triangleq dE^{\prime}_{q}/dt$, $y_{Eq}\triangleq dE^{\prime}_{d}/dt$, $y_{\omega}\triangleq d\omega/dt$, $y_{Id}\triangleq I_{d}$, $y_{Iq}\triangleq I_{q}$, $\beta_{1}\triangleq-1/T^{\prime}$, $\beta_{2}\triangleq-(X-X^{\prime})/T^{\prime}$, $\beta_{3}\triangleq-{1}/{2H}$, $\alpha_{b}\triangleq{R_{s}}/{R_{s}^{2}+X^{\prime 2}}$, $\alpha_{c}\triangleq{X^{\prime}}/{R_{s}^{2}+X^{\prime 2}}$. Assuming that the measurement noise is i.i.d. and follows a normal distribution, we can rewrite the IM model as follows:

[TABLE]

where ${\varepsilon_{E_{d}}[t]\sim\mathcal{N}(0,1/\tau_{E})}$, ${\varepsilon_{E_{q}}[t]\sim\mathcal{N}(0,1/\tau_{E})}$, ${\varepsilon_{\omega}[t]\sim\mathcal{N}(0,1/\tau_{\omega})}$, ${\varepsilon_{I_{d}}[t]\sim\mathcal{N}(0,1/\tau_{I})}$, and ${\varepsilon_{I_{q}}[t]\sim\mathcal{N}(0,1/\tau_{I})}$.

III Bayesian Estimation in Composite Load Parameter Identification

III-A Gibbs sampling

Gibbs sampling is an extension of Monte Carlo Markov Chain method [17], which performs well when there are multiple parameters to identify. The detailed sampling algorithm is shown in Algorithm 1.

Starting with priors, Gibbs sampling estimates the posterior of one parameter while fixing others’ values as samples from previous estimated posteriors. This process repeats for all parameters in one iteration.

III-B Gibbs in ZIP Model

Start with the prior guess as follows:

[TABLE]

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 remains the same form. Only the first iteration will be shown here as an example. First, samples are drawn from the prior and $\alpha_{1}^{(0)},\alpha_{2}^{(0)},\tau^{(0)}$ are initialized. After some algebraic yields:

[TABLE]

where $p(\alpha_{1}|\alpha_{2}^{(0)},\tau^{(0)},\textbf{x},\textbf{y})$ is the posterior probability given the samples of x, y, $\alpha_{2}$, and $\tau$, $p(\textbf{x},\textbf{y}|\alpha_{1},\alpha_{2}^{(0)},\tau^{(0)})$ is the likelihood in (1), and $p(\alpha_{1})$ is the prior estimation following (3). Taking the log form on both sides of (6) yields

[TABLE]

Taking log form helps convert the multiplications of the probabilies to summations, which can significantly simplify calculations when updating the distributions of the posteriors. For a normal distribution $y\sim\mathcal{N}(\mu,1/\tau)$, the log dependence on $y$ is $-\frac{\tau}{2}(y-\mu)^{2}\propto-\frac{\tau}{2}y^{2}+\tau\mu y$.

The right hand side of equation (7) can be further written as the following if the terms not related to $\alpha_{1}$ are omitted.

[TABLE]

Define

[TABLE]

[TABLE]

yields

[TABLE]

Then sample a new $\alpha_{1}$ from the estimated distribution $\mathcal{N}(\mu_{1}^{(1)},\tau_{1}^{(1)})$ as $\alpha_{1}^{(1)}$. Following the similar procedures $\alpha_{2}$ can be derived. Define

[TABLE]

Therefore, the following can be derived:

[TABLE]

For $\tau$, the posterior given new samples of $\alpha_{1}^{(1)}$ and $\alpha_{2}^{(1)}$ can be written as ${p(\tau|\alpha_{1}^{(1)},\alpha_{2}^{(1)},\textbf{x},\textbf{y})\propto p(\textbf{x},\textbf{y}|\alpha_{1}^{(1)},\alpha_{2}^{(1)},\tau)p(\tau)}$. Taking the log form of both sides of the posterior and yields

[TABLE]

Define

[TABLE]

the posterior $\tau|\alpha_{1}^{(1)},\alpha_{2}^{(1)},\textbf{x},\textbf{y}\sim\mathcal{G}(a^{(1)},b^{(1)})$ can be obtained.

III-C Gibbs in IM models

Start with the priors as follows.

[TABLE]

The parameter identification in the IM model follows a similar procedure as the ZIP model and is omitted here for space consideration. Details can be found in [18].

IV Case Study

Simulation studies are carried out in this section for ZIP and IM models, separately. In this work, the maximum Monte Carlo runs $M$ is set to 40000 and the burn-in length $m$ is chosen as 5000 [17].

IV-A ZIP Model Identification

The 33-bus test feeder [19] is used to generate the testing data. A detailed description of the test system can be found in [18]. Load at bus 18 is replaced by a ZIP model. Randomness is added to loads at other buses by multiplying a random factor drawing from a uniform distribution following $\mathcal{U}[0.1,4.5]$. The random load is stated as $P_{w,i}=P_{i}\cdot\mathcal{U}[0.1,4.5]$, 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 simulates the disturbances and uncertainties in the system, which significantly impact 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 converges in each iteration with different $P_{w,i}$, and the corresponding voltage $V_{w,i}$ is recorded. The measurement noise $\varepsilon$ follows $\mathcal{N}(0,0.1)$, which means a 10% measurement error.

By implementing the proposed GS method, the coefficients of the ZIP model are estimated as shown in Fig. 2.

Different from other estimation approaches, the GS approach can generate a distribution which describes the probability of the real value falling into a certain range. As shown in Fig. 2, the mean values of the distributions of $\alpha_{1}$ and $\alpha_{2}$ are 0.25 and 0.249, respectively, and $\alpha_{3}=1-\alpha_{1}-\alpha_{2}$.

Figs. 3 and 4 show the voltage and voltage/real power differences, respectively, using the estimated and real model in 100 i.i.d. experiments. The randomness comes from the factors multiplied to loads at other buses. The dash line with triangle in Fig. 3 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 Fig. 4 is the distributions of voltage and active power differences, $\Delta V$ and $\Delta P$. It can be seen that the differences are in the range of 10*-4* and 10*-3*, respectively. Fig. 5 shows the burn-in process observed. It shows that the process is very fast, and there is no significant burn-in process that can be observed from this figure.

IV-B IM model Identification

The estimation of the parameters in the IM model follows the same procedure. A dynamic model is built in Matlab/Simulink to generate the data that is used in (2). The sampling results are shown in Fig. 6 and compared with the real data in Table I. According to Table I, the estimation error are less than 6% with 5% measurements error. It is worth to mention that estimation of parameters highly relies on the prior distributions. A good prior can significantly increase the estimation accuracy and shorten the burn-in period.

IV-C Benchmarks

The ZIP and IM model parameters derived by the proposed GS method are compared with least square (LS) [20] and Kalman Filter (KF) methods [21] with 10% noise. The “fit” function in Matlab was used to derive the coefficients in LS. The parameter tuning in KF was introduced in [21], and these parameters were directly used in this study. The comparison results of the voltage and active power using GS, LS, and KF in ZIP model are shown in Fig 7. The average absolute mean of voltage and active power errors in p.u. are listed in Table II. According to Fig 7, it is obvious that GS approach has the best performance, while KF method only falls slightly behind. In contrast, LS gives the largest error among all three approaches when there is large measurement error. However, KF only works well with time-invariant parameters estimation [22]. 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 III. The estimation errors from GS are the smallest. The performance of KF is worse than itself in the ZIP model, but still better than the LS method.

V Conclusion

In this paper, a novel Bayesian estimation based load model parameter identification method is proposed. The proposed method can accurately identify the parameter distributions for both ZIP and IM models. Compared with other measurement based algorithms, the proposed method provides a distribution estimation and is robust to measurement errors. The accuracy and robustness of the proposed method compared with conventional load modeling method are also demonstrated by numerical experiments. The identified distributions can be further used in load prediction, stability and reliability analysis, as well as other related areas. Further study may include how to estimate the ZIP and IM model jointly and extending the process in an online manner.