Bayesian estimation based parameter estimation for composite load model

Abstract

A method for managing a power load of a grid includes performing a statistic-based distribution estimation of a composite load model using static and dynamic models with Gibbs sampling; deriving a distribution estimation of load model coefficients; and controlling grid power based on a simulation, a prediction, a stability analysis or a reliability analysis with the load model coefficients.

Claims

1. A method for managing a power load of a grid, comprising: performing a statistic-based distribution estimation of a composite load model using static and dynamic models with Gibbs sampling, wherein the load comprises resistive and inductive loads; deriving a distribution estimation of load model coefficients; and controlling grid power based on a simulation, a prediction, a stability analysis or a reliability analysis with the load model coefficients by applying the simulation, prediction, stability analysis or reliability analysis to control power of electrical loads having impedance, current or power provided to one or more resistive and inductive loads.

2. The method of claim 1, wherein the statistic-based estimation comprises a Bayesian estimation.

3. The method of claim 1, wherein the static model comprises a ZIP model.

4. The method of claim 3, wherein the ZIP model describes load power changes as a voltage varies in a steady-state condition.

5. The method of claim 3, wherein the ZIP model comprises: $\begin{array}{cc}P=P_0\left({\alpha }_1{\overline{V}}^2+{\alpha }_2\overline{V}+{\alpha }_3\right)Q=Q_0\left({\alpha }_4{\overline{V}}^2+{\alpha }_5\overline{V}+{\alpha }_6\right)\underset{i=1}{\overset{3}{.Math.}}{\alpha }_i=\underset{i=4}{\overset{6}{.Math.}}{\alpha }_i=1P=P_0\left({\alpha }_1\overline{V^2}+{\alpha }_2\overline{V}+{\alpha }_3\right)Q=Q_0\left({\alpha }_4\overline{V^2}+{\alpha }_5\overline{V}+{\alpha }_6\right)\underset{i=1}{\overset{3}{.Math.}}{\alpha }_i=\underset{i=4}{\overset{6}{.Math.}}{\alpha }_i=1& \end{array}$ where Σ.sub.i=1.sup.3α.sub.i=λ.sub.p, Σ.sub.i=4.sup.6α.sub.i=λ.sub.q, V=V/V.sub.0.

6. The method of claim 1, wherein the dynamic model comprises an induction motor (IM) model.

7. The method of claim 6, wherein the IM model can be expressed as follows: $\{\begin{array}{c}{y}_{E_d}\left[t\right]={\beta }_1{E^{\prime }}_d\left[t\right]+{\beta }_2{I}_q\left[t\right]-\left(\omega -1\right){E^{\prime }}_q\left[t\right]+{\mathrm{.Math.}}_{E_d}\left[t\right]\\ y_{E_q}\left[t\right]={\beta }_1{E^{\prime }}_q\left[t\right]-{\beta }_2{I}_d\left[t\right]+\left(\omega -1\right){E^{\prime }}_d\left[t\right]+{\mathrm{.Math.}}_{E_q}\left[t\right]\\ y_{\omega }\left[t\right]={\beta }_3\left({\omega }^2-{E^{\prime }}_d\left[t\right]I_d\left[t\right]-{E^{\prime }}_q\left[t\right]I_q\left[t\right]\right)+{\mathrm{.Math.}}_{\omega }\left[t\right]\end{array}\{\begin{array}{c}{y}_{I_d}\left[t\right]={\alpha }_b\left(U_d\left[t\right]-{E^{\prime }}_d\left[t\right]\right)+{\alpha }_c\left(U_q\left[t\right]-{E^{\prime }}_q\left[t\right]\right)+{\mathrm{.Math.}}_{I_d}\left[t\right]\\ y_{I_q}\left[t\right]={\alpha }_b\left(U_q\left[t\right]-{E^{\prime }}_q\left[t\right]\right)+{\alpha }_c\left(U_d\left[t\right]-{E^{\prime }}_d\left[t\right]\right)+{\mathrm{.Math.}}_{I_q}\left[t\right]\end{array}$ where ε.sub.E.sub.d[t], ε.sub.E.sub.q[t], ε.sub.ω[t], ε.sub.I.sub.d[t], ε.sub.I.sub.q[t] are all normal distributions with mean equals 0, variance equals 1/τ.sub.E, 1/τ.sub.E, 1/τ.sub.Ω, 1/τ.sub.I, 1/τ.sub.I, respectively.

8. The method of claim 1, comprising determining: $\begin{array}{cc}\{\begin{array}{c}{P}_{ZIP}=P_{ZIP0}\left({\alpha }_1\overline{V^2}+{\alpha }_2\overline{V}+{\alpha }_3\right),\\ Q_{\mathrm{ZIP}}=Q_{\mathrm{ZIP}0}\left({\alpha }_4\overline{V^2}+{\alpha }_5\overline{V}+{\alpha }_6\right),\end{array}& \\ \{\begin{array}{c}\frac{{{\mathrm{dE}}^{\prime }}_d}{\mathrm{dt}}=-\frac{1}{T\text{'}}\left[{E^{\prime }}_d+\left(X-X^{\prime }\right)I_q\right]-\left(\omega -1\right){E^{\prime }}_q\\ \frac{{{\mathrm{dE}}^{\prime }}_q}{\mathrm{dt}}=-\frac{1}{T\text{'}}\left[{E^{\prime }}_q-\left(X-X^{\prime }\right)I_d\right]+\left(\omega -1\right){E^{\prime }}_d\\ \frac{d\omega }{\mathrm{dt}}=-\frac{1}{2H}\left[\left(A{\omega }^2+B\omega +C\right)T_0-\left({E^{\prime }}_d{I}_d+{E^{\prime }}_q{I}_q\right)\right]\end{array}\{\begin{array}{c}{I}_d=\frac{1}{\left(R_s^2+{X^{\prime }}^2\right)}\left[R_s\left(U_d-{E^{\prime }}_d\right)+X^{\prime }\left(U_q-{E^{\prime }}_q\right)\right]\\ I_q=\frac{1}{\left(R_s^2+{X^{\prime }}^2\right)}\left[R_s\left(U_q-{E^{\prime }}_q\right)-X^{\prime }\left(U_d-{\left(E^{\prime }\right)}_d\right)\right]\end{array}& \end{array}$ where parameters α.sub.1, . . . α.sub.6, X.sub.m, X.sub.r, X.sub.s, R.sub.r, A, B, C, H and wherein:
X′X.sub.s+X.sub.mX.sub.r/(X.sub.m+X.sub.r),
XX.sub.s+X.sub.m,
T′(X.sub.r+X.sub.m)/R.sub.r.

9. The method of claim 1, comprising determining a prior.

10. The method of claim 9, comprising applying gradient descent (GD) to update the prior.

11. The method of claim 10, wherein the GD comprises applying an Adaptive Moment Estimation (Adam) to determine adaptive learning rates for each parameter.

12. The method of claim 1, comprising determining a total active power and reactive power as: $\{\begin{array}{c}P={\lambda }_p{P}_{ZIP}+\left(1-{\lambda }_p\right)P_{IM}\\ Q={\lambda }_q{Q}_{ZIP}+\left(1-{\lambda }_q\right)Q_{IM}\end{array}\hspace{1em}$ where λ.sub.p+λ.sub.q=1.

13. The method of claim 1, comprising determining real power and reactive power of the load.

14. The method of claim 1, wherein a mean value is selected as an estimated value of each parameter.

15. The method of claim 1, comprising performing Gibbs sampling of the IM model.

16. The method of claim 15, comprising determining τ.sub.I.sup.(k+1)|E′.sub.d,E′.sub.q,U.sub.d,U.sub.q,I.sub.d,I.sub.q,α.sub.b.sup.(k+1),α.sub.c.sup.(k+1)˜G(α.sub.I.sup.(k+1),β.sub.I.sup.(k+1)).

Description

BRIEF DESCRIPTIONS OF FIGURES

(1) The following Figures are for illustration purposes only and are not drawn to scale. The exemplary embodiments, both as to organization and method of operation, may best be understood by reference to the detailed description which follows taken in conjunction with the accompanying drawings in which:

(2) FIG. 1 shows an exemplary composite model of ZIP+IM.

(3) FIG. 2 shows a flowchart of an exemplary algorithm.

(4) FIG. 3 illustrates the estimated parameters of the ZIP model.

(5) FIG. 4 shows the estimated parameters of the ZIP model.

(6) FIG. 5 illustrates the comparison of system active response between the original system and the estimated system.

(7) FIG. 6 shows the active power response comparison by using LS and GA.

(8) FIG. 7 shows the error of active power using GS, GA, LS with 10% measurement error.

(9) FIG. 8 shows the error of active power using GS, GA, LS with measurement anomalies.

(10) FIG. 9 shows estimated and real value of the parameters for ZIP+IM mode with 0 measurement error.

(11) FIG. 10 shows estimated and real value of the parameters for ZIP+IM mode with 10% measurement error and anomalies.

(12) FIG. 11 shows composite model parameter estimations without measurement error by different methods.

(13) FIG. 12 shows composite model parameter estimation with different methods considering 10% measurement errors.

(14) FIG. 13 shows composite model parameter estimation with different methods considering outliers.

DETAILED DESCRIPTION OF THE ILLUSTRATED EMBODIMENTS

(15) Nomenclature

(16) P.sub.ZIP, Q.sub.ZIP: Real/reactive power measurement of a ZIP load

(17) α.sub.1, . . . , α.sub.6: Coefficients to estimate in the ZIP model

(18) V: voltage magnitude at the load terminal

(19) P.sub.ZIP0, Q.sub.ZIP0, V.sub.0: corresponding measurement of a ZIP load at the initial operating condition.

(20) V=V/V.sub.0

(21) P.sub.IM, Q.sub.IM: Real/reactive power measurement of an IM load

(22) {circumflex over (P)}, {circumflex over (Q)}: estimate of real/reactive power

(23) θ: the coefficient vector to estimate

(24) q(θ): the prior distribution of the coefficients

(25) R.sub.S: motor stator winding resistance

(26) X.sub.S: motor stator leakage reactance

(27) X.sub.m: motor magnetizing reactance

(28) R.sub.r: rotor resistance

(29) X.sub.r: rotor leakage reactance

(30) H: rotor inertia constant

(31) ω: rotor speed

(32) I.sub.d, I.sub.q: stator current in d-axis and q-axis

(33) U.sub.d, U.sub.q: bus voltage in d-axis and q-axis

(34) E.sub.d′, E.sub.q′: stator voltage in d-axis and q-axis

(35) T.sub.0: the initial load torque

(36) X.sub.r, X.sub.m, X.sub.s, R.sub.r, R.sub.T, A, B, C, H: Coefficients to estimate in the IM model

(37) (x, y){(x.sub.ZIP[1], y.sub.ZIP[1]), . . . , (x.sub.ZIP[n], y.sub.ZIP[n])}: time series measurement of the ZIP model

(38) μ, μ.sub.1, μ.sub.2, μ.sub.3, μ.sub.β.sub.1, μ.sub.β.sub.2, μ.sub.β.sub.3, μ.sub.α.sub.b, μ.sub.α.sub.c: mean of the corresponding distribution.

(39) τ, τ.sub.1, τ.sub.2, τ.sub.3, τ.sub.β.sub.1, τ.sub.β.sub.2, τ.sub.β.sub.3, τ.sub.α.sub.b, τ.sub.α.sub.c, τ.sub.E, τ.sub.ω, τ.sub.I: reciprocal of the standard deviation of the corresponding distribution

(40) ε[t], ε.sub.E.sub.d [t], ε.sub.E.sub.q[t], ε.sub.ω[t], ε.sub.I.sub.d[t], ε.sub.I.sub.q[t]: noise of the corresponding model

(41) λ.sub.p, λ.sub.q: real/reactive power fraction of the ZIP model

(42) I. Approaches

(43) 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. The mathematical representations of the ZIP and IM model can be found in 1 and 2, respectively.

(44) $\begin{array}{cc}\{\begin{array}{c}{P}_{ZIP}=P_{ZIP0}\left({\alpha }_1\stackrel{\_}{V^2}+{\alpha }_2\stackrel{\_}{V}+{\alpha }_3\right),\\ Q_{\mathrm{ZIP}}=Q_{\mathrm{ZIP0}}\left({\alpha }_4\stackrel{\_}{V^2}+{\alpha }_5\stackrel{\_}{V}+{\alpha }_6\right),\end{array}& \left(1\right)\\ \{\begin{array}{c}\frac{{\mathrm{dE}}_d^{\prime }}{\mathrm{dt}}=-\frac{1}{T^{\prime }}\left[E_d^{\prime }+\left(X-X^{\prime }\right)I_q\right]-\left(\omega -1\right)E_q^{\prime }\\ \frac{{\mathrm{dE}}_q^{\prime }}{\mathrm{dt}}=-\frac{1}{T^{\prime }}\left[E_q^{\prime }\left(X-X^{\prime }\right)I_d\right]+\left(\omega -1\right)E_d^{\prime }\\ \frac{d\omega }{\mathrm{dt}}=-\frac{1}{2H}\left[\left(A{\omega }^2+\mathrm{B\omega }+C\right)T_0-\left(E_d^{\prime }{I}_d+E_q^{\prime }{I}_q\right)\right]\end{array}\{\begin{array}{c}{I}_d=\frac{1}{\left(R_s^2+X^{\mathrm{\prime 2}}\right)}\left[R_s\left(U_d-E_d^{\prime }\right)+X^{\prime }\left(U_q-E_q^{\prime }\right)\right]\\ I_q=\frac{1}{\left(R_s^2+X^{\mathrm{\prime 2}}\right)}\left[R_s\left(U_q-E_q^{\prime }\right)-X^{\prime }\left(U_d-E_d^{\prime }\right)\right]\end{array}& \left(2\right)\end{array}$

(45) The parameters need to be identified in this composite model are α.sub.1, . . . α.sub.6, X.sub.m, X.sub.r, X.sub.s, R.sub.r, A, B, C, H. Some notations in equations (1) and (2) are defined as follows:
X′X.sub.s+X.sub.mX.sub.r/(X.sub.m+X.sub.r),
XX.sub.s+X.sub.m,
T′(X.sub.r+X.sub.m)/R.sub.r  (3)

(46) A. ZIP Model

(47) 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.

(48) $\begin{array}{cc}{P}_{ZIP}=P_{ZIP0}\left({\alpha }_1\stackrel{\_}{V^2}+{\alpha }_2\stackrel{\_}{V}+{\alpha }_3\right)Q_{\mathrm{ZIP}}=Q_{\mathrm{ZIP0}}\left({\alpha }_4\stackrel{\_}{V^2}+{\alpha }_5\stackrel{\_}{V}+{\alpha }_6\right)\underset{i=1}{\overset{3}{.Math.}}{\alpha }_i=\underset{i=4}{\overset{6}{.Math.}}{\alpha }_i=1& \left(4\right)\end{array}$
where Σ.sub.i=1.sup.3α.sub.i=λ.sub.p, Σ.sub.i=4.sup.6α.sub.i=λ.sub.q, V=V/V.sub.0. When voltage V deviates from V.sub.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.

(49) One implementation defines y.sub.ZIP[t]=P.sub.ZIP[t]/P.sub.ZIP0, x.sub.ZIP[t]=V[t]/V.sub.0, where P.sub.ZIP[t] and V[t] are the t.sup.th measurements of the real power and voltage magnitude of the ZIP load. The following assumptions are made.

(50) The measurement noise follows a normal distribution, i.e., y.sub.ZIP[t]=α.sub.1x.sub.ZIP[t].sup.2+α.sub.2 x.sub.ZIP[t]+α.sub.3+ε[t], where ε[t]˜(0,1/τ), and 1/τ 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 experiments is large enough. 2) Since normal distribution is a conjugate distribution, it is easier for model parameter updating when implementing Gibbs sampling.

(51) Total number of n independent and identically distributed (i.i.d.) samples were drawn, namely (x, y){(x.sub.ZIP[1], y.sub.ZIP [1]), . . . , (x.sub.ZIP[n], y.sub.ZIP[n])}. Thus, the likelihood

(52) $\begin{array}{cc}p\left(x,y|{\alpha }_1,{\alpha }_2,\tau \right)\propto \underset{t=1}{\overset{n}{.Math.}}\exp \left[{\left(y_{ZIP}\left[t\right]-\mu \left[t\right]\right)}^2\tau /2\right],& \left(5\right)\end{array}$
where μ[t]α.sub.1x[t].sup.2+α.sub.2x[t]+α.sub.3 is the mean.

(53) B. IM Model

(54) The parameters to be identified in an IM model are X.sub.r, X.sub.m, X.sub.s, R.sub.r, R.sub.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. Consequently B and C are both equal to 0. For simplicity, define Y.sub.E.sub.ddE.sub.q/dt, y.sub.E.sub.qdE.sub.d/dt, y.sub.ωdω/dt, y.sub.I.sub.dI.sub.d, y.sub.I.sub.qI.sub.q, β.sub.11/T, β.sub.2−(X−X′)/T′, β.sub.3−½H, α.sub.bR.sub.s/R.sub.s.sup.2+X.sup.2, α.sub.cX′/R.sub.s.sup.2+X.sup.2. Assuming that the measurement noise is i.i.d. and follows a normal distribution, the IM model can be expressed as follows:

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

(56) where ε.sub.E.sub.d[t], ε.sub.E.sub.q[t], ε.sub.ω[t], ε.sub.I.sub.d[t], ε.sub.I.sub.q[t] are all normal distributions with mean equals 0, variance equals 1/τ.sub.E.sub.d, 1/τ.sub.E.sub.q, 1/τ.sub.ω, 1/τ.sub.I, 1/τ.sub.I, respectively. The detailed simulation results with different measurement error will be introduced in the Case Studies section below to illustrate the robustness of the instant method.

(57) For a composite model with ZIP+IM model, the total active power and reactive power can be respectively computed as:

(58) $\{\begin{array}{c}P={\lambda }_p{P}_{ZIP}+\left(1-{\lambda }_p\right)P_{IM}\\ Q={\lambda }_q{Q}_{ZIP}+\left(1-{\lambda }_q\right)Q_{IM}\end{array}$
where λ.sub.p+λ.sub.q=1. Since the available data can be used for identification of aforementioned parameters is limited by the fact that the composite model proposed in FIG. 1 is an equivalent model, some of the parameters are not able to be direct measured, such as E.sub.d, E.sub.q, ω, etc. Thus, it is challenging to identify all the parameters only use P, r, Q, and U, especially for parameters in equation (2). Therefore, in literature such problems were considered as optimization problems. Typical values of some parameters, e.g. R.sub.s, X.sub.m, X.sub.r, were selected to help solving the differential equations in (2), and the calculated values (E.sub.d, E.sub.q, ω) will be substituted into the constraints, for instance, in (5), to solve optimization problems. If (5) is smaller than the predefined threshold value, the identification procedure will be terminated. Otherwise, a new typical value will be selected to follow the same steps until the threshold is met.

(59) $\left[{\left(-\right)}^2+{\left(-\right)}^2\right]$

(60) In one embodiment, a BE based method is used without optimizing the error between the measured value and the estimated value. The algorithm in FIG. 2 is inspired by the traditional method of: selecting typical values to solve the differential equations in (2) as well as using priors to calculate power consumed ZIP load. The results of calculated E.sub.d, E.sub.q, ω) are used in GS to sample the parameters in 2. The sampled values are used to calculate the updated value of PIM by using 6, and PZIP by using 1, which will be used in next iteration as the updated values of the parameters. The detailed steps of this algorithm can be found in FIG. 2.

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

(62) TABLE-US-00001 Algorithm 1 Gibbs Sampling Draw initial samples θ.sup.(0) ~ q(θ), where q(θ) is the prior For iteration i = 1,2, ... , M do Calculate p (θ.sub.1|θ.sub.2.sup.(i−1), θ.sub.3.sup.(i−1), ... , θ.sub.n.sup.(i−1)) and sample θ.sub.1.sup.(i) ~ p (θ.sub.1|θ.sub.2.sup.(i−1), θ.sub.3.sup.(i−1), ... , θ.sub.n.sup.(i−1)) Calculate p (θ.sub.2|θ.sub.1.sup.(i), θ.sub.3.sup.(i−1), ... , θ.sub.n.sup.(i−1)) and sample θ.sub.2.sup.(i) ~ p (θ.sub.2|θ.sub.1.sup.(i), θ.sub.3.sup.(i−1), ... , θ.sub.n.sup.(i−1)) Calculate p (θ.sub.n|θ.sub.1.sup.(i), θ.sub.3.sup.(i), ... , θ.sub.n−1.sup.(i)) and Sample θ.sub.n.sup.(i) ~ p (θ.sub.n|θ.sub.1.sup.(i), θ.sub.3.sup.(i), ... , θ.sub.n−1.sup.(i)) End for The distribution estimate is the histogram of θ.sup.i, i = m, ... , M. Others are burn-in data and discarded

(63) 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.

(64) C. Gibbs in ZIP Model

(65) Start with the prior guess as follows:
α.sub.1˜(μ.sub.1.sup.(0),1/τ.sub.1.sup.(0))  (7)
α.sub.2˜(μ.sub.2.sup.(0),1/τ.sub.2.sup.(0))  (8)
α.sub.3˜(μ.sub.3.sup.(0),1/τ.sub.3.sup.(0))  (9)
τ˜(a.sup.(0),b.sup.(0))  (10)
where the distribution of τ is a gamma distribution follows (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 α.sub.1.sup.(0), α.sub.2.sup.(0), α.sub.3.sup.(0), τ.sup.(0) are initialized. After some algebraic yields:
p(α.sub.1|α.sub.2.sup.(0),α.sub.3.sup.(0),τ.sup.(0),x,y)∝p(x,y|α.sub.1,α.sub.2.sup.(0),α.sub.3.sup.(0),τ.sup.(0))p(α.sub.1)  (11)
where p(α.sub.1|α.sub.2.sup.(0), τ.sup.(0), x, y) is the posterior probability given the samples of x, y, α.sub.2, and τ, p(x, y|α.sub.1, α.sub.2.sup.(0), τ.sup.(0)) is the likelihood in (3), and p(α.sub.1) is the prior estimation following (7). Taking the log form on both sides of (11) yields

(66) $\begin{array}{cc}\log p\left({\alpha }_1|{\alpha }_2^{\left(0\right)},{\alpha }_3^{\left(0\right)},{\tau }^{\left(0\right)},x,y\right)\propto -\frac{{\tau }_1^{\left(0\right)}}{2}{\left({\alpha }_1-{\mu }_1^{\left(0\right)}\right)}^2-\frac{{\tau }^{\left(0\right)}}{2}\underset{t=1}{\overset{n}{.Math.}}{\left(y\left[t\right]-\left({\alpha }_1{x\left[t\right]}^2+{\alpha }_2^{\left(0\right)}x\left[t\right]+{\alpha }_3^{\left(0\right)}\right)\right)}^2& \left(12\right)\end{array}$

(67) Taking log form helps convert the multiplications of the probabilities to summations, which can significantly simplify calculations when updating the distributions of the posteriors. For a normal distribution y˜(μ, 1/τ), the log dependence on y is

(68) $-\frac{\tau }{2}{\left(y-\mu \right)}^2\propto -\frac{\tau }{2}{y}^2+\tau \mu y$

(69) The right-hand side of equation (12) can be further written as the following if the terms not related to α.sub.1 are omitted.

(70) $\begin{array}{cc}-\left({\tau }_1^{\left(0\right)}+{\tau }^{\left(0\right)}\underset{t=1}{\overset{n}{.Math.}}{\left({x\left[t\right]}^2-1\right)}^2\right){\alpha }_1^2/2+\left({\tau }_1^{\left(0\right)}{\mu }_1^{\left(0\right)}-{\tau }^{\left(0\right)}\underset{t=1}{\overset{n}{.Math.}}\left({\alpha }_3^{\left(0\right)}\left(y\left[t\right]-1\right){x\left[t\right]}^2-{\alpha }_2^{\left(0\right)}{x\left[t\right]}^3\right)\right){\alpha }_1\mathrm{Define}{\mu }_1^{\left(1\right)}\frac{\begin{array}{c}{\tau }_1^{\left(0\right)}{\mu }_1^{\left(0\right)}-{\tau }^{\left(0\right)}{\Sigma }_{t=1}^n\\ \left(\left({\alpha }_3^{\left(0\right)}\left(y\left[t\right]-1\right){x\left[t\right]}^2-{\alpha }_2^{\text{(0)}}{x\left[t\right]}^3\right)\right)\end{array}}{{\tau }_1^{\left(0\right)}+{\tau }^{\left(0\right)}{{\Sigma }_{t=1}^n\left({x\left[t\right]}^2-1\right)}^2},\mathrm{and}{\tau }_1^{\left(1\right)}{\tau }_1^{\left(0\right)}+{\tau }^{\left(0\right)}{{\Sigma }_{t=1}^n\left({x\left[t\right]}^2-1\right)}^2,\mathrm{yields}{\alpha }_1.Math.{\alpha }_2^{\left(0\right)},{\alpha }_3^{\left(0\right)},{\tau }^{\left(0\right)},{{\mu }_1}^{\left(0\right)},{{\tau }_1}^{\left(0\right)},x,y\sim \mathrm{��}\left({{\mu }_1}^{\left(1\right)},1/{{\tau }_1}^{\left(1\right)}\right)& \left(13\right)\end{array}$

(71) Then sample a new α.sub.1 from the estimated distribution (μ.sub.1.sup.(1)), τ.sub.1.sup.(1)) as α.sub.1.sup.(1). Following the similar procedures α.sub.2 and α.sub.3 can be derived. Define

(72) 0${\mu }_2^{\left(1\right)}\stackrel{\Delta }{=}\frac{{\tau }_2^{\left(0\right)}{\mu }_2^{\left(0\right)}-{\tau }^{\left(0\right)}{\Sigma }_{t=1}^n\left(\left({\alpha }_3^{\left(0\right)}\left(y\left[t\right]-1\right)x\left[t\right]-{\alpha }_1^{\left(1\right)}{x\left[t\right]}^3\right)\right)}{{\tau }_2^{\left(0\right)}+{\tau }^{\left(0\right)}{{\Sigma }_{t=1}^n\left(x\left[t\right]-1\right)}^2}$
and τ.sub.2.sup.(1)τ.sub.2.sup.(0)+τ.sup.(0)Σ.sub.i=1.sup.n(x[t]−1).sup.2. Therefore, the following can be derived: α.sub.2|α.sub.1.sup.(1), α.sub.2.sup.(0), τ.sup.(0), τ.sub.2.sup.(0), μ.sub.2.sup.(0), x, y˜(μ.sub.2.sup.(1), 1/τ.sub.2.sup.(1))

(73) Define

(74) ${\mu }_3^{\left(1\right)}\stackrel{\Delta }{=}\frac{{\tau }_3^{\left(0\right)}{\mu }_3^{\left(0\right)}-{\tau }^{\left(0\right)}{\Sigma }_{t=1}^n\left(y\left[t\right]-{\alpha }_1^{\left(1\right)}{x\left[t\right]}^2-{\alpha }_2^{\left(1\right)}x\left[t\right]\right)}{{\tau }_3^{\left(0\right)}+{\tau }^{\left(0\right)}{{\Sigma }_{t=1}^n\left(x\left[t\right]-1\right)}^2}\mathrm{and}$${\tau }_3^{\left(1\right)}\stackrel{\Delta }{=}{\tau }_3^{\left(0\right)}+{\tau }^{\left(0\right)}{.Math.}_{i=1}^n{\left(x\left[t\right]-1\right)}^2.$

(75) Therefore, the following can be derived:
α.sub.3|α.sub.1.sup.(1),α.sub.2.sup.(1),τ.sup.(0),τ.sub.3.sup.(0),μ.sub.4.sup.(0),x,y˜(μ.sub.3.sup.(1),1/τ.sub.3.sup.(1)).

(76) For τ, the posterior given new samples of α.sub.1.sup.(1), and α.sub.2.sup.(1) can be written as p(τ|α.sub.1.sup.(1), α.sub.2.sup.(1), x, y)∝p(x, y|α.sub.1.sup.(1), α.sub.2.sup.(1), τ)p(τ). Taking the log form of both sides of the posterior and yields

(77) $\log p\left(\tau .Math.{\alpha }_1^{\left(1\right)},{\alpha }_2^{\left(1\right)},\varkappa ,y\right)\propto \frac{n}{2}\log \tau -\frac{\tau }{2}\underset{t=1}{\overset{n}{.Math.}}\left(y\left[t\right]-{\alpha }_1^{\left(1\right)}{\varkappa \left[t\right]}^2-{\alpha }_2^{\left(1\right)}\varkappa \left[t\right]-{\alpha }_3^{\left(1\right)}+\left(a^{\left(0\right)}-1\right)\log \tau -b^{\left(0\right)}\tau \right).$

(78) Define a.sup.(1)=a.sup.(0)+n/2, b.sup.(1)=.sup.(0)+Σ.sub.t=1.sup.n(y[t]−α.sub.1.sup.(1)x[t].sup.2−α.sub.2.sup.(1)x[t]−α.sub.3.sup.(1)).sup.2/2. The posterior τ|α.sub.1.sup.(1), α.sub.2.sup.(1), x, y, ˜(a.sup.(1), b.sup.(1)) can be obtained.

(79) D. Gibbs in IM Models

(80) GS in IM model will follow the same steps as that in ZIP model. Start with the priors as follows.
β.sub.1˜(μ.sub.β1.sup.(0),1/τ.sub.β1.sup.(0)),β.sub.2˜(μ.sub.β2.sup.(0),1/τ.sub.β2.sup.(0)),
β.sub.3˜(μ.sub.β.sub.3.sup.(0),1/τ.sub.β.sub.3.sup.(0)),α.sub.b˜(μ.sub.α.sub.b.sup.(0),1/τ.sub.α.sub.b.sup.(0))
α.sub.c˜(μ.sub.αc.sup.(0),1/τ.sub.αc.sup.(0)),τ.sub.E˜(a.sub.E.sup.(0),b.sub.E.sup.(0))
τ.sub.ω˜(a.sub.ω.sup.(0),b.sub.ω.sup.(0)),τ.sub.I˜(a.sub.I.sup.(0),b.sub.I.sup.(0)).

(81) Given the measurement of y.sub.E.sub.d, y.sub.E.sub.q, y.sub.ω, E.sub.q, E.sub.d, I.sub.d, I.sub.q, U.sub.d, U.sub.q, Gibbs sampling of the IM model is stated as follows.

(82) Define

(83) ${\mu }_{{\beta }_1}^{\left(k+1\right)}=\frac{\begin{array}{c}{\tau }_{{\beta }_1}^{\left(k\right)}{\mu }_{{\beta }_1}^{\left(k\right)}+{\tau }_E^{\left(k\right)}{.Math.}_{t=1}^n({E^{\prime }}_d\left[t\right]y_{E_d}\left[t\right]+{E^{\prime }}_q\left[t\right]y_{E_q}\left[t\right]-\\ {\beta }_2^{\left(k\right)}\left({E^{\prime }}_d\left[t\right]I_q\left[t\right]-{E^{\prime }}_q\left[t\right]I_d\left[t\right]\right))\end{array}}{{\tau }_{{\beta }_1}^{\left(k\right)}+{\tau }_E^{\left(k\right)}{.Math.}_{t=1}^n\left({{E^{\prime }}_d\left[t\right]}^2+{{E^{\prime }}_q\left[t\right]}^2\right)},{\tau }_{{\beta }_1}^{\left(k+1\right)}=\frac{1}{{\tau }_{{\beta }_1}^{\left(k\right)}+{\tau }_E^{\left(k\right)}{.Math.}_{t=1}^n\left({{E^{\prime }}_d\left[t\right]}^2+{{E^{\prime }}_q\left[t\right]}^2\right)}$
We have, β.sub.1|β.sub.2.sup.(k), y.sub.E.sub.d, y.sub.E.sub.q, E.sub.d, E.sub.q, U.sub.d, U.sub.q, I.sub.q, I.sub.d, ω, τ.sub.E.sup.(k), τ.sub.β.sub.1.sup.(k), μ.sub.β.sub.1.sup.(k)˜(μ.sub.β.sub.1.sup.(k+1), τ.sub.β1.sup.(k+1))
Define:

(84) ${\mu }_{{\beta }_2}^{\left(k+1\right)}=\left[{\tau }_{{\beta }_2}^{\left(k\right)}{\mu }_{{\beta }_2}^{\left(k\right)}+{\tau }_E^{\left(k\right)}\stackrel{n}{\underset{t=1}{.Math.}}\left(y_{E_d}\left[t\right]I_q\left[t\right]-y_{E_q}\left[t\right]I_d\left[t\right]+\hspace{1em}\hspace{1em}{\beta }_1^{\left(k+1\right)}\left({E^{\prime }}_d\left[t\right]I_q\left[t\right]-{E^{\prime }}_q\left[t\right]I_d\left[t\right]\right)+\left(\omega \left[t\right]-1\right)\left({E^{\prime }}_d\left[t\right]I_d\left[t\right]+E_q\left[t\right]I_q\left[t\right]\right)\right)\right]/[{\tau }_{{\beta }_2}^{\left(k\right)}+\hspace{1em}{\tau }_E^{\left(k\right)}\underset{t=1}{\overset{n}{.Math.}}\left({I_d\left[t\right]}^2+{I_q\left[t\right]}^2\right)],{\tau }_{{\beta }_2}^{\left(k+1\right)}=\frac{1}{{\tau }_{{\beta }_2}^{\left(k\right)}+{\tau }_E^{\left(k\right)}{.Math.}_{i=1}^n\left({I_d\left[t\right]}^2+{I_q\left[t\right]}^2\right)},$
We have β.sub.2|β.sub.1.sup.(k+1), y.sub.E.sub.d, y.sub.E.sub.q, E.sub.d, E.sub.q, U.sub.d, U.sub.q, I.sub.d, I.sub.q, ω, τ.sub.E.sup.(k), τ.sub.β.sub.2.sup.(k), μ.sub.β.sub.2.sup.(k)˜(μ.sub.β.sub.2.sup.(k+1), τ.sub.β.sub.2.sup.(k+1))
Define

(85) ${\mu }_{{\beta }_3}^{\left(k+1\right)}=[{\tau }_{{\beta }_3}^{\left(k\right)}{\mu }_{{\beta }_3}^{\left(k\right)}-\hspace{1em}{\tau }_{\omega }^{\left(k\right)}\stackrel{n}{\underset{t=1}{.Math.}}\left(y_{\omega }\left[t\right]\left(-{\omega \left[t\right]}^2{T}_0+{E^{\prime }}_d\left[t\right]I_d\left[t\right]+E_q\left[t\right]I_q\left[t\right]\right)\right)][{\tau }_{{\beta }_3}^{\left(k\right)}+\hspace{1em}{\tau }_{\omega }^{\left(k\right)}\stackrel{n}{\underset{t=1}{.Math.}}\left({\omega \left[t\right]}^4{T}_0^2+{\left({E^{\prime }}_d\left[t\right]I_d\left[t\right]+E_q\left[t\right]I_q\left[t\right]\right)}^2-2{\omega \left[t\right]}^2{T}_0\left({E^{\prime }}_d\left[t\right]I_d\left[t\right]+E_q\left[t\right]I_q\left[t\right]\right)\right)],{\mu }_{{\beta }_3}^{\left(k+1\right)}=1/\left[{\tau }_{{\beta }_3}^{\left(k\right)}+{\tau }_{\omega }^{\left(k\right)}\stackrel{n}{\underset{t=1}{.Math.}}\left({\omega \left[t\right]}^4{T}_0^2+{\left({E^{\prime }}_d\left[t\right]I_d\left[t\right]+E_q\left[t\right]I_q\left[t\right]\right)}^2-2{\omega \left[t\right]}^2{T}_0\left({E^{\prime }}_d\left[t\right]I_d\left[t\right]+E_q\left[t\right]I_q\left[t\right]\right)\right)\right],$
We have, β.sub.3|E.sub.d, E.sub.q, U.sub.d, U.sub.q, I.sub.d, I.sub.q, ω, τ.sub.E.sup.(k), τ.sub.β.sub.3.sup.(k), μ.sub.β.sub.3.sup.(k)˜(μ.sub.β.sub.3.sup.(k+1), τ.sub.β.sub.3.sup.(k+1))
Define

(86) ${\mu }_{{\alpha }_b}^{\left(k+1\right)}=\frac{\begin{array}{c}{\tau }_{{\alpha }_b}^{\left(k\right)}{\mu }_{{\alpha }_b}^{\left(k\right)}+{\tau }_I^{\left(k\right)}{.Math.}_{t=1}^n(I_d\left[t\right]\left(U_d\left[t\right]-{E^{\prime }}_d\left[t\right]\right)+\\ I_q\left[t\right]\left(U_q\left[t\right]-{E^{\prime }}_q\left[t\right]\right))\end{array}}{{\tau }_{{\alpha }_b}^{\left(k\right)}+{\tau }_I^{\left(k\right)}{.Math.}_{t=1}^n\left({\left(U_d\left[t\right]-E_d\left[t\right]\right)}^2+{\left(U_q\left[t\right]-E_q\left[t\right]\right)}^2\right)},{\tau }_{{\alpha }_b}^{\left(k+1\right)}=\frac{1}{{\tau }_{{\alpha }_b}^{\left(k\right)}+{\tau }_I^{\left(k\right)}{.Math.}_{t=1}^n\left({\left(U_d\left[t\right]-E_d\left[t\right]\right)}^2+{\left(U_q\left[t\right]-E_q\left[t\right]\right)}^2\right)}$
We have, α.sub.b|E.sub.d, E.sub.q, U.sub.d, U.sub.q, I.sub.d, I.sub.q, ω, τ.sub.I.sup.(k), τ.sub.α.sub.b.sup.(k), μ.sub.α.sub.b.sup.(k)˜(μ.sub.α.sub.b.sup.(k+1), τ.sub.α.sub.b.sup.(k+1))
Define

(87) ${\mu }_{{\alpha }_c}^{\left(k+1\right)}=\frac{\begin{array}{c}{\tau }_{{\alpha }_c}^{\left(k\right)}{\mu }_{{\alpha }_c}^{\left(k\right)}+{\tau }_I^{\left(k\right)}{.Math.}_{t=1}^n(I_d\left[t\right]\left(U_d\left[t\right]-{E^{\prime }}_d\left[t\right]\right)-\\ I_q\left[t\right]\left(U_q\left[t\right]-{E^{\prime }}_q\left[t\right]\right))\end{array}}{{\tau }_{{\alpha }_c}^{\left(k\right)}+{\tau }_I^{\left(k\right)}{.Math.}_{t=1}^n\left({\left(U_d\left[t\right]-E_d\left[t\right]\right)}^2+{\left(U_q\left[t\right]-E_q\left[t\right]\right)}^2\right)},{\tau }_{{\alpha }_c}^{\left(k+1\right)}=\frac{1}{{\tau }_{{\alpha }_c}^{\left(k\right)}+{\tau }_I^{\left(k\right)}{.Math.}_{t=1}\left({\left(U_d\left[t\right]-E_d\left[t\right]\right)}^2+{\left(U_q\left[t\right]-E_q\left[t\right]\right)}^2\right)}$
We have, α.sub.c|E.sub.d, E.sub.q, U.sub.d, U.sub.q, I.sub.d, I.sub.q, ω, τ.sub.I.sup.(k), τ.sub.α.sub.c.sup.(k), μ.sub.α.sub.c.sup.(k)˜(μ.sub.α.sub.c.sup.(k+1), τ.sub.α.sub.c.sup.(k+1))
Define

(88) $a_E^{\left(k+1\right)}=a_E^{\left(k\right)}+n/2,b_E^{\left(k+1\right)}=b_E^{\left(k\right)}+\stackrel{n}{\underset{t=1}{.Math.}}\left[{\left(y_{E_d}\left[t\right]-{\beta }_1^{\left(k+1\right)}{E^{\prime }}_q\left[t\right]-{\beta }_2^{\left(k+1\right)}{I}_d\left[t\right]+\left(\omega \left[t\right]-1\right){E^{\prime }}_q\left[t\right]\right)}^2+{\left(y_{E_q}\left[t\right]-{\beta }_1^{\left(k+1\right)}{E^{\prime }}_d\left[t\right]+{\beta }_2^{\left(k+1\right)}{I}_q\left[t\right]-\left(\omega \left[t\right]-1\right){E^{\prime }}_d\left[t\right]\right)}^2\right]/2$
We have, τ.sub.E.sup.(k+1)|I.sub.d, I.sub.q, y.sub.E.sub.d, y.sub.E.sub.q, E.sub.d, E.sub.q, ω, β.sub.1.sup.(k+1), β.sub.2.sup.(k+1)˜(a.sub.E.sup.(k+1), b.sub.E.sup.(k+1))
Define

(89) $a_{\omega }^{\left(k+1\right)}=a_{\omega }^{\left(k\right)}+n/2,b_{\omega }^{\left(k+1\right)}=b_{\omega }^{\left(k\right)}+\frac{1}{2}\underset{t=1}{\overset{n}{.Math.}}{\left[y_{\omega }\left[t\right]-{\beta }_3^{\left(k+1\right)}\left({\omega \left[t\right]}^2{T}_0-{E^{\prime }}_d\left[t\right]I_d\left[t\right]-{E^{\prime }}_q\left[t\right]I_q\left[t\right]\right)\right]}^2,$
We have, τ.sub.ω.sup.(k+1)|I.sub.d, I.sub.q, E.sub.d, E.sub.q, y.sub.ω, ω, U.sub.d, U.sub.q, β.sub.3.sup.(k+1)˜(a.sub.ω.sup.(k+1), b.sub.ω.sup.(k+1))
Define

(90) 0$a_I^{\left(k+1\right)}=a_I^{\left(k\right)}+n/2,b_I^{\left(k+1\right)}=b_I^{\left(k\right)}+\underset{t=1}{\overset{n}{.Math.}}\left[{\left(I_d\left[t\right]-{\alpha }_b^{\left(k+1\right)}\left(U_d\left[t\right]-{E^{\prime }}_d\left[t\right]\right)-{\alpha }_c^{\left(k+1\right)}\left(U_q\left[t\right]-{E^{\prime }}_q\left[t\right]\right)\right)}^2+{\left(I_q\left[t\right]-{\alpha }_b^{\left(k+1\right)}\left(U_q\left[t\right]-{E^{\prime }}_q\left[t\right]\right)+{\alpha }_c^{\left(k+1\right)}\left(U_d\left[t\right]-{E^{\prime }}_d\left[t\right]\right)\setminus big\right)}^2\right]/2,$
We have, τ.sub.I.sup.(k+1)|E.sub.d, E.sub.q, U.sub.d, U.sub.q, I.sub.d, I.sub.q, α.sub.b.sup.(k+1), α.sub.c.sup.(k+1)˜(a.sub.I.sup.(k+1), b.sub.I.sup.(k+1))

(91) E. Updates of the Prior

(92) Gradient descent is a first-order iterative optimization algorithm for finding the minimum of a function. To find a local minimum of a function using gradient descent, one takes steps proportional to the negative of the gradient (or approximate gradient) of the function at the current point. In one embodiment, ADAM method is used for parameter updating as stated in Algorithm 2.

(93) TABLE-US-00002 Algorithm2 For a slightly more efficient (but less clear) order of computation. g.sub.t .sup.2 indicates the elementwise square g.sub.t ⊙ g.sub.t. Good default settings for the tested machine learning problems are s = 0. 001, β.sub.1 = 0. 9, β.sub.2 = 0. 999 and ϵ = 10.sup.−8. All operations on vectors are element-wise. With β.sub.1.sup.t and β.sub.2.sup.t we denote β.sub.1 and β.sub.2 to the power t. Require: Δ: Step size Require: β.sub.1, β.sub.2 ∈ [0, 1): Exponential decay rates for the moment estimates Require: f(θ) Stochastic objective function with parameters θ Require: θ.sub.0: Initial parameter vector m.sub.0 ← 0 (Initialize 1st moment vector) v.sub.0 ← 0 (Initialize 2nd moment vector) t ← 0 (Initialize timestep) while θ.sub.t not converged do t ← t + 1 g.sub.t ← ∇.sub.θf.sub.t(θ.sub.t−1) (Get gradients w.r.t. stochastic objective at timestep t) m.sub.t ← β.sub.1 .Math. m.sub.t−1 + (1 − β.sub.1) .Math. g.sub.t (Update biased first moment estimate) v.sub.t ← β.sub.2 .Math. v.sub.t−1 + (1 − β.sub.2) .Math. g.sub.t.sup.2 (Update biased second raw moment estimate)  ← m.sub.t/(1 − β.sub.1 .sup.t) (Compute bias-corrected first moment estimate)  ← v.sub.t /(1 − β.sub.2.sup.t) (Compute bias-corrected second raw moment estimate) θ.sub.t ← θ.sub.t−1 − Δ .Math.   /(   + ϵ) (Update parameters) end while return (Resulting parameters)

(94) II. Case Studies

(95) Simulation studies are carried out in this section for composite models with ZIP+IM considering measurement errors and outliers. Benchmark test using least square and genetic algorithm methods are also proposed. Only active power is shown as an example, the reactive power will follow the same steps. In GS, the length of burn-in process m is chosen as 5000.

(96) A. Composite Model Parameters Identification

(97) The 33-bus test feeder, as shown in FIG. 3, is simulated to generate the testing data. Load at bus 18 is connected by a ZIP+IM model as an example. The ZIP factors of the load at node 18 were assigned as α.sub.1=0.25, α.sub.2=0.25, α.sub.3=0.5 with λ.sub.p=0.5. Therefore, the expected values of the parameters are all scaled λ. By implementing the proposed GS method, the coefficients for ZIP model and IM model are estimated as shown in FIG. 4. The discrete distributions for each parameter can be derived by implementing the instant method. Distributions describe the probabilities of different values falling into different ranges. In this study, the mean value is selected as the estimated value of each parameter as shown in FIG. 9. The estimation error for each parameter compared with the real values are also included in FIG. 9. According to FIG. 9, the estimation errors are less than 6% with 0 measurement error. In most cases, the error is even smaller. The active power response by using the identified parameters is shown in FIG. 5. Comparing with the original active response, only tiny difference can be observed. It is worth to mention that the estimation of parameters highly relies on the prior distributions. A good prior can significantly increase the estimation accuracy and shorten the burn-in period.

(98) GS is able to handle certain measurement error and anomalies. In FIG. 10, the estimated values by considering 10% measurement error and measurement anomalies by using the proposed method are presented. The measurement anomalies are simulated by replacing the results of 2500 measurements by 0 (both active power and voltage) among total 250,000 readings, which can be considered as 1% measurement anomaly. The results shown in FIG. 10 indicate the identifications are influenced by the error and anomalies, but the accuracy is still within an acceptable range.

(99) B. Comparison with Benchmarks

(100) The ZIP and IM model parameters derived by the proposed GS method are compared against least square (LS) and Genetic Algorithm (GA) methods. The “Isqcurvefit” and “ga” function in Matlab was used to derive the coefficients in LS and GA, respectively. Parameters identified by three different approaches are listed in FIG. 11. Comparisons of the active power responses using parameters identified by LS, and GA techniques are presented in FIG. 6. According to FIGS. 6 and 10, it is obvious that GS approach has the best performance, while GA and LS methods larger error than GS method. The scenarios with measurement errors and anomalies are also considered in LS and GA approaches. Same measurement error and anomaly is used, and the comparisons are shown in FIGS. 12 and 13, respectively. The system responses using the aforementioned parameters by GS, LS, GA approaches with 10% measurement error and 1% measurement anomaly are shown in FIGS. 7 and 8, respectively. It can be shown that under the same circumstance, GS method is able to provide the best estimation, while LS and GA methods are significantly impacted by measurement error and anomalies.

(101) With the integration of renewable resources, electrical vehicles and other uncertain resources into power system, the accuracy of the instant load model plays an important role in simulation, prediction, as well as stability and reliability analysis. The statistic (Bayesian Estimation) based distribution estimation approach generates composite load models, including static (ZIP) and dynamic (Induction Motor) parts, by implementing Gibbs sampling. The system provides a distribution estimation of load models coefficients and is robust to measurement errors.