Estimating Dynamic Load Parameters from Ambient PMU Measurements

TL;DR

This paper introduces a new method for estimating dynamic load parameters using ambient PMU data, avoiding the need for large disturbances and providing continuous, real-time updates with demonstrated accuracy and robustness.

Contribution

The paper presents a novel ambient measurement-based algorithm for dynamic load parameter estimation that operates continuously without large disturbances.

Findings

Method accurately estimates load parameters in simulations.

Algorithm is robust to measurement noise.

Provides real-time updates of load parameters.

Abstract

In this paper, a novel method to estimate dynamic load parameters via ambient PMU measurements is proposed. Unlike conventional parameter identification methods, the proposed algorithm does not require the existence of large disturbance to power systems, and is able to provide up-to-date dynamic load parameters consistently and continuously. The accuracy and robustness of the method are demonstrated through numerical simulations.

Tables

${\delta }_i$	generator rotor angle
${\omega }_i$	generator angular frequency
$M_i$	inertial constant
$P_{mi}$	mechanical power input
$P_{Gi}({\delta }_i,{\theta }_i,V_i)$	real power injection
$Q_{Gi}({\delta }_i,{\theta }_i,V_i)$	reactive power injection
$D_i$	damping coefficient
$N$	total number of buses
${\theta }_{ij}$	voltage angle difference between bus $i$ and $j$
$\left|V_i\right|$	voltage magnitude
$G_{ij}$	line conductance between bus $i$ and $j$
$B_{ij}$	line susceptance between bus $i$ and $j$

$g_k$	effective conductance of the load
$b_k$	effective susceptance of the load
${\tau }_{gk}$	active power time constant of the load
${\tau }_{bk}$	reactive power time constant of the load
$P_k$	real power demand of the load
$Q_k$	reactive power demand of the load
$P_k^s$	steady-state real power demand of the load
$Q_k^s$	steady-state reactive power demand of the load

$𝒈={[g_1,\dots ,g_m]}^T$,	$𝒃={[b_1,\dots ,b_m]}^T$,
$T_g=\text{diag}[{\tau }_{g1},\dots ,{\tau }_{gm}]$,	$T_b=\text{diag}[{\tau }_{b1},\dots ,{\tau }_{bm}]$,
$𝑷={[P_1,\dots ,P_m]}^T$,	$𝑸={[Q_1,\dots ,Q_m]}^T$,
$P^s=\text{diag}[P_1^s,\dots ,P_m^s]$,	$Q^s=\text{diag}[Q_1^s,\dots ,Q_m^s]$,
${\Sigma }^p=\text{diag}[{\sigma }_1^p,\dots ,{\sigma }_m^p]$,	${\Sigma }^q=\text{diag}[{\sigma }_1^q,\dots ,{\sigma }_m^q]$,
${𝝃}^{𝒑}={[{\xi }_1^p,\dots ,{\xi }_m^p]}^T$,	${𝝃}^{𝒒}={[{\xi }_1^q,\dots ,{\xi }_m^q]}^T$.

Table 4. TABLE I : A comparison between the actual and the estimated time constants in the 9-bus system

	actual value (s)	estimated value (s)	error
${\tau }_{g1}$	1.0000	0.9145	8.55%
${\tau }_{g2}$	3.0000	2.9867	0.44%
${\tau }_{g3}$	0.2000	0.2122	6.1%
${\tau }_{b1}$	5.0000	4.7974	4.05%
${\tau }_{b2}$	7.0000	6.9777	0.32%
${\tau }_{b3}$	0.8000	0.7462	6.72%

Table 5. TABLE II : A comparison between the actual and the estimated time constants in the 9-bus system with the measurement noises

	actual value (s)	estimated value (s)	error
${\tau }_{g1}$	1.0000	0.9144	8.56%
${\tau }_{g2}$	3.0000	2.9819	0.60%
${\tau }_{g3}$	0.2000	0.2121	6.06%
${\tau }_{b1}$	5.0000	4.7752	4.50%
${\tau }_{b2}$	7.0000	6.9426	0.82%
${\tau }_{b3}$	0.8000	0.7443	6.97%

Table 6. TABLE III : A comparison between the actual and the estimated time constants in the 39-bus system

	actual value (s)	estimated value (s)	error
${\tau }_{g1}$	0.1000	0.1225	22.55%
${\tau }_{g2}$	0.6000	0.5860	2.33%
${\tau }_{g3}$	1.1000	1.1086	0.78%
${\tau }_{g4}$	1.6000	1.6924	5.78%
${\tau }_{g5}$	2.1000	2.1323	1.54%
${\tau }_{g6}$	2.6000	2.6130	0.50%
${\tau }_{g7}$	3.1000	3.0715	0.92%
${\tau }_{g8}$	3.6000	3.3004	8.32%
${\tau }_{g9}$	4.1000	4.5080	9.95%
${\tau }_{g10}$	4.6000	4.6636	1.38%
${\tau }_{b1}$	0.5000	0.5277	5.53%
${\tau }_{b2}$	1.0000	0.9831	1.69%
${\tau }_{b3}$	1.5000	1.5325	2.17%
${\tau }_{b4}$	2.0000	2.0559	2.79%
${\tau }_{b5}$	2.5000	2.5784	3.14%
${\tau }_{b6}$	3.0000	3.3522	11.74%
${\tau }_{b7}$	3.5000	3.1273	10.65%
${\tau }_{b8}$	4.0000	4.1105	2.76%
${\tau }_{b9}$	4.5000	4.3084	4.26%
${\tau }_{b10}$	5.0000	5.1723	3.45%

Equations

$${\dot{\delta }}_i$$

$$M_i\dot{{\omega }_i}$$

$$P_{Gi}({\delta }_i,{\theta }_i,V_i)$$

$$Q_{Gi}({\delta }_i,{\theta }_i,V_i)$$

$${\dot{g}}_k$$

$${\dot{b}}_k$$

$$P_k=g_k{V}_k^2$$

$$Q_k=b_k{V}_k^2$$

$${\dot{g}}_k$$

$${\dot{b}}_k$$

$$\displaystyle\left[\begin{array}[]{c}\dot{\bm{g}}\\ \dot{\bm{b}}\end{array}\right]$$

$$+$$

$$=$$

$$A{C}_{\mathbf{xx}}+C_{\mathbf{xx}}{A}^T=-B{B}^T$$

$$\frac{\partial{P_{k}}}{\partial{g_{j}}}=\left\{\begin{array}[]{cc}V_{k}^{2}+2g_{j}V_{k}\frac{\partial{V_{k}}}{\partial{g_{j}}}\approx{V_{k}^{2}}&\mbox{if }j=k\\ 2g_{j}V_{k}\frac{\partial{V_{k}}}{\partial{g_{j}}}\approx 0&\mbox{if }j\not=k\end{array}\right.$$

$$\displaystyle A\approx\left[\begin{array}[]{cc}-{T_{g}}^{-1}V^{2}&\bm{0}\\ \bm{0}&-{T_{b}}^{-1}V^{2}\end{array}\right]$$

$$C_{\mathbf{gg}}=\frac{1}{2}{T}_g^{-1}(P^s{)}^2({\Sigma }^p{)}^2{V}^{-2}$$

$$C_{\mathbf{bb}}=\frac{1}{2}{T}_b^{-1}(Q^s{)}^2({\Sigma }^q{)}^2{V}^{-2}$$

$$C_{\mathbf{gb}}=C_{\mathbf{bg}}=0$$

$$g_k$$

$$b_k$$

$$Q_{g_k{g}_j}$$

$$Q_{b_k{b}_j}$$

$$T_g$$

$$T_b$$

$$\displaystyle Q_{\bm{gg}}=\left[\begin{array}[]{ccc}1.41\times 10^{-3}&9.84\times 10^{-5}&2.24\times 10^{-4}\\ 9.84\times 10^{-5}&4.16\times 10^{-4}&6.57\times 10^{-5}\\ 2.24\times 10^{-4}&6.57\times 10^{-5}&5.75\times 10^{-3}\end{array}\right]$$

$$\displaystyle Q_{\bm{bb}}=\left[\begin{array}[]{ccc}2.63\times 10^{-4}&-1.47\times 10^{-6}&6.46\times 10^{-6}\\ -1.47\times 10^{-6}&1.75\times 10^{-4}&-7.57\times 10^{-6}\\ 6.46\times 10^{-6}&-7.57\times 10^{-6}&1.62\times 10^{-3}\end{array}\right]$$

Full text

I Introduction

Load modelling and identification are of great importance to the security and stability of power systems. While the accurate models are available for generators, transmission lines and control devices, load modeling is still a challenging and open subject due to the fact that electric load at each substation is an aggregation of numerous individual loads with different behaviors [1]-[3]. In addition, the poor measurements, modeling, exchange information, as well as the uncertainties in customers behaviors/devices further result in load uncertainties [2]. Indeed, load uncertainty is one of the main factors that affect the accuracy of the power dynamic models implemented by system operators over the world[3].

Generally speaking, the load uncertainty comes from both model structure and parameter values. It has been shown in previous literature [4]-[6] that the use of different load models leads to different and even contradictory results for dynamic stability studies. Even though the applied model structure is verified, different parameter values may also yield different damping performances in small signal stability [2][7][8]. For instance, different time constants of loads may lead to either asymptotically stable system or systems experiencing oscillations (i.e., Hopf bifurcation occurs)[2]. Both load modelling and parameter identifications are essential in studying the dynamic behaviors of power systems. This paper mainly focuses on parameter identification for a generic dynamic load model that is suitable for small signal stability analysis.

Different methods for dynamic load parameter identification have been proposed, which can be classified into two categories: component-based approach [9] and measurement-based approach [10]-[16]. The latter approach is more commonly applied because real-time load variations and dynamic characteristics can be taken into account[14]. Measurement-based model identification is typically solved through optimization methods that minimize the error between the measured output variables and the simulated ones. In particular, the nonlinear least-square curve fitting method has been implemented in [1][10]-[12]. Genetic algorithms, neural network-based methods and other heuristic techniques have been applied in [13]-[16]. However, the optimization-based methods are time consuming and thus can not be implemented online[4]. More importantly, all those methods require measurement data from dynamic behaviors of system under big disturbances (e.g., during faults), which is not always available [17]. Indeed, the variation of load parameters may be much faster than the occurrence rate of natural disturbances [1].

In this paper, we propose a novel measurement-based method for dynamic load identification in ambient conditions, which does not require the existence of large disturbance. Particularly, the method combines the statistical properties extracted from PMU measurements and the inherent model knowledge, and is able to provide fairly accurate estimations for parameter values in near real-time. Note that a generic dynamic load model is implemented in this paper which is suitable for the purposes of small signal stability analysis and damping performance [2][7][8][18]. The proposed method can be implemented in online security analysis to provide up-to-date dynamic load parameters accurately.

The rest of the paper is organized as follows. Section II introduces the power system stochastic dynamic model. Particularly, the generic dynamic load model used in small signal stability is presented. Section III elaborates the proposed method for estimating parameters of dynamic loads. Section IV presents the validation of the proposed method through numerical simulations. The impact of measurement noise is also investigated. Conclusions and perspectives are given in Section V.

II power system stochastic dynamic model

Although we focus on load models, generator models are also incorporated to provide more realistic simulations. Specifically, the classical generator model which can reasonably represent the dynamics of generator in ambient conditions is implemented. The power system buses are numbered as follows: load buses: $k=1,2,...,m$, and generators: $i=m+1,...,N$. Particularly, to include the effects of the loads, the structure preserving model [19][20] is used:

[TABLE]

where

The detailed expressions of $P_{Gi}(\delta_{i},\theta_{i},V_{i})$ and $Q_{Gi}(\delta_{i},\theta_{i},V_{i})$ are neglected here for simplicity and can be found in many books (e.g., [19]).

Regarding dynamic loads, we use the following first-order load model proposed in [18] that can represent the common types of loads (e.g., induction motors, thermostatically controlled loads) in ambient conditions:

[TABLE]

where

The values $P_{k}^{s}$ and $Q_{k}^{s}$ describe the static (steady-state) power characteristics of the loads achieved in equilibrium. The instant real power and reactive power consumption can be characterized by the effective conductance $P_{k}=g_{k}V_{k}^{2}$ and susceptance $Q_{k}=b_{k}V_{k}^{2}$ at any time. The time constants $\tau_{gk}$ and $\tau_{bk}$ that typically depend on voltage and frequency represent the instant relaxation rate of the load.

To incorporate load variation, we apply a similar approach used in [21][22] and modify the set of load equations (5)-(6) as follows:

[TABLE]

where the steady-state real and reactive load demands are perturbed with independent Gaussian noise from their initial values. Specifically, $\xi_{k}^{p}$ and $\xi_{k}^{q}$ are standard Gaussian noise, and $\sigma_{k}^{p}$ and $\sigma_{k}^{q}$ represent the noise intensities for static real and reactive power, respectively.

As discussed in [2][18], this dynamic load model can naturally represent the most common types of loads in ambient conditions such as thermostatic load, induction motor, power electronic converter, aggregate effects of distribution load tap changer (LTC) transformers, etc. However, the range of time constants is considerably large ranging from cycles to several minutes, and even hours for different types of loads. For industrial plants, such as aluminum smelters, the time constants are in the range of $0.1$s to $0.5$s; for tap changers and other control devices, they are in the range of minutes; for heating load, they may range up to hours [8]. As a result, the uncertainty of composition of different types of loads can be aggregated in time constants $\tau_{g}$ and $\tau_{b}$[2]. This is reasonable in the situations when the network characteristics are known, generator models are validated and static load characteristics are understood better than their dynamic response which is the case in practical situations. In addition to a wide range of time constants, the variation of $\tau_{g}$ and $\tau_{b}$ can also be fast. For example, $\tau_{b}$ may change from $0.1$s to $24.1$s in one day (see Table I, II in [1]).

Because of wide range and fast variation of time constants $\tau_{g}$ and $\tau_{b}$, they need to be updated frequently to ensure the accuracy of dynamic load models used in online security and stability analysis. Conventionally, $\tau_{g}$ and $\tau_{b}$ are estimated from dynamic data by perturbing the system, for example, through changing the transformer tap[23]. However it is impractical to perturb the system frequently for estimating parameter values of loads. In this paper, we propose a novel method to estimate $\tau_{g}$ and $\tau_{b}$ for the loads of interests from ambient PMU measurements in daily operation. In particular, the estimation process does not require the existence of disturbance to the system.

III Methodology

In ambient conditions, the stochastic dynamic load equations (9)-(10) can be linearized as below:

[TABLE]

where

It is observed that $[\bm{g},\bm{b}]^{T}$ is a vector Ornstein-Uhlenbeck process that is stationary, Gaussian and Markovian [24][25]. Particularly, if the state matrix $A$ is stable, the stationary covariance matrix C_{\bm{xx}}=\left[\begin{array}[]{cc}C_{\bm{gg}}&C_{\bm{gb}}\\ C_{\bm{bg}}&C_{\bm{bb}}\end{array}\right] can be shown to satisfy the following Lyapunov equation[24][26]:

[TABLE]

which nicely combines the model knowledge and the statistical properties of state variables.

Since $P_{k}=g_{k}V_{k}^{2}$, we have:

[TABLE]

under the assumption that $\triangle V_{k}\approx 0$ in ambient conditions. Similar relation can be obtained for $\frac{\partial{Q_{k}}}{\partial{b_{j}}}$. As a result, the Jacobian matrix $A$ satisfies:

[TABLE]

where $V=\mbox{diag}[V_{1},...,V_{m}]$. Substituting (32) and the detailed expression of $C_{\bm{xx}}$ and $B$ into (28), and performing algebraic simplification, we have:

[TABLE]

Particularly, we utilize the relations (33)-(34) that link the measurements of stochastic load variation to the physical model, and provide an ingenious way to estimate the dynamic parameters $T_{g}$ and $T_{b}$ from measurements.

In practical applications, $V$, $C_{\bm{gg}}$, and $C_{\bm{bb}}$ need to be acquired or estimated from limited PMU measurements. A window size of $1000$s is used in the examples of this paper where time constants are up to several seconds. Note that the larger the time constants, the longer the sample window is needed to ensure accuracy. First, the sample mean $\bar{V}$ can be used as an estimation of $V$, then $\bm{g}$ and $\bm{b}$ can be estimated from PMU measurements (i.e., phasors $V_{k}$ and $I_{k}$) as follows:

[TABLE]

Regarding the covariance matrix $C_{\bm{g}\bm{g}}=\mathrm{E}[(\bm{g}-\mathrm{E}[\bm{g}])(\bm{g}-\mathrm{E}[\bm{g}])^{T}]$ and $C_{\bm{b}\bm{b}}=\mathrm{E}[(\bm{b}-\mathrm{E}[\bm{b}])(\bm{b}-E[\bm{b}])^{T}]$, we use their unbiased estimators—sample covariance matrixes $Q_{\bm{g}\bm{g}}$ and $Q_{\bm{b}\bm{b}}$ in practice, each entry of which is calculated as below:

[TABLE]

where $\bar{{g}}_{k}$ and $\bar{{b}}_{k}$ denote the sample mean of ${g_{k}}$ and ${b_{k}}$, respectively, and $n$ is the sample size.

Therefore, the proposed algorithm can be summarized as follows. We assume that PMUs are installed at the substations that the (aggregated) loads of interests are connected to. We also assume that the static characteristics of loads are well understood such that $P^{s}$, $Q^{s}$, $\Sigma^{p}$ and $\Sigma^{q}$ are prior known, which is reasonable as shown in[21][22]. Then the following algorithm provides an estimation of $T_{g}$ and $T_{b}$ for the dynamic loads from ambient PMU measurements:

Step 1.

Compute the sample mean $\bar{V}$ and estimate $\bm{g}$ and $\bm{b}$ from PMU measurements by (36)-(37).

Step 2.

Calculate the sample covariance matrix $Q_{\bm{g}\bm{g}}$ and $Q_{\bm{b}\bm{b}}$ by (38)-(39).

Step 3.

Approximate $T_{g}$ and $T_{b}$ as blow:

[TABLE]

Note that (40)-(41) are acquired by a simple algebraic manipulation of (33)-(34).

IV case studies

In this section, the proposed algorithm to estimate time constants of dynamic loads are validated through numerical simulations. Furthermore, the robustness of the proposed method to measurement noise is also demonstrated via simulation. All case studies were done in PSAT-2.1.9 [27].

IV-A Validation of the Method

We consider the standard WSCC 3-generator, 9-bus system model (see, e.g. [19]). The classical generator models (1)-(2) and the stochastic dynamic load models (9)-(10) are implemented in the structure preserving framework. The system parameters are available online: https://github.com/xiaozhew/PES-load-parameter-estimation. Particularly, there are three dynamic loads at buses 1, 2 and 3, the time constants of which are $\tau_{g}=1,3,0.2\mbox{s}$ and $\tau_{b}=5,7,0.8\mbox{s}$, respectively. The trajectories of some state variables and algebraic variables are shown in Fig. 1, from which we see that the state variables are fluctuating around their nominal values in ambient conditions, yet larger time constants lead to slower variations as expected (e.g., the variations of $g_{2}$ and $b_{2}$ are slower than $g_{3}$ and $b_{3}$).

By the proposed algorithm, we firstly compute the sample mean $\bar{V}=\mbox{diag}[0.9952,1.0126,1.0155]$. Then we estimate the dynamic conductance, susceptances and their corresponding sample covariance matrixes:

[TABLE]

[TABLE]

It is expected that both $Q_{\bm{gg}}$ and $Q_{\bm{bb}}$ are nearly diagonal as the stochastic perturbations are independent.

Since each entry of $(P^{s})^{2}(\Sigma^{p})^{2}$ and $(Q^{s})^{2}(\Sigma^{q})^{2}$ is set to be $0.0025$, $T_{g}$ and $T_{b}$ can be readily estimated from (40)-(41). A comparison between the estimated $\tau_{g}$, $\tau_{b}$ and their actual values are shown in Table. I. It’s observed that the proposed algorithm provides fairly accurate estimation for time constants of each load.

IV-B Impact of Measurement Noise

Like other measurement-based methods, the performance of the proposed algorithm may be affected by PMU measurement noise. In order to investigate the potential influence, measurement noises with standard deviation of $10^{-3}$ have been added to $\bm{g}$, $\bm{b}$ and ${V}$ in the 9-bus example shown in Section IV-A according to the IEEE Standards [28][29]. A comparison between the actual and the estimated time constants are presented in Table. II. It is observed that the proposed method provides similar accuracy to the case without the measurement noises, which indicates that the method is relatively robust under measurement noise.

IV-C Further Validation

For further validation, we apply the method to a larger system—the IEEE 39-bus 10-generator test system, the parameters of which are available online: https://github.com/xiaozhew/PES-load-parameter-estimation. In particular, 10 dynamic loads have been added to buses 1-10, and their corresponding time constants range from $0.1$s to $5$s. A comparison between the actual and the estimated time constants are presented in Table. III. The simulation results further demonstrate that the proposed method is able to provide good estimations for time constants of the dynamic loads.

V conclusions and perspectives

In this paper, we have proposed a novel method to estimate parameter values of dynamic load from ambient PMU measurements. The accuracy and robustness of the method have been demonstrated through numerical studies. Unlike conventional methods, the proposed technique does not require the existence of large disturbance to systems, and thus can be implemented continuously in daily operation to provide up-to-date dynamic load parameter values.

In the future, we plan to further validate the method by using real PMU data and extend the method to estimate dynamic load parameters without knowing their static characteristics.