Recently, some novel researches regarding the state estimation for power distribution system have been focused on the dynamic state estimation algorithm due to trend abrupt changes in the distribution system such as actions of prosumers. The previous works studied the Kalman filter for nonlinear system such as extended Kalman filter, unscented Kalman filter. They assume that the state transient model follows smoothing model of the previous states. However, it is hard to set the smoothing coefficient for the estimated states. Thus, this paper proposes the mathematical model of state transient model of power distribution system. Based on the proposed mathematical model, extended Kalman filter algorithm is adopted for dynamic state estimation. All jacobian matrices for the extended Kalman filter is derived as a function of system states, and the proposed algorithm is verified by using IEEE 15 test bus.

#### The Transactions of the Korean Institute of Electrical Engineers

**ISO Journal Title**Trans. Korean. Inst. Elect. Eng.

- SCOPUS
- KCI Accredited Journal

### Journal Search

## 1. Introduction

Power systems are controlled to maintain the stage which is planned by power system
operators. To operating the power sys- tem efficiently, state estimation on power
system such as distribu- tion system and transmission system is crucial. Thus, many
papers on state estimation for power system have been published such as least squre
algorithm ^{(1)}, weighted least square estimation [2–4] and Gauss-Newton method ^{(5,}^{6)}. However, these approaches are static state estimation, which are hard to capture
the properties of system dynamics.

Nowadays, distribution system becomes larger and complex. It may have unexpected changes
such as abrupt disconnection of generator and electric vehicles. Increasing number
of photovoltaic panel and wind turbine also causes the dynamic behavior of power system.
Therefore, instead of static estimation, the researches on dynamic state estimation
algorithm has been published. Especially, the authors’ of ^{(7)} studied the extended Kalman filter for distri-algorithm requires state transient
model which describes the state of next time step by using current information. Because
the exact state transient model of distribution power system is unknown, the exist-
ing dynamic state estimation algorithm have used the state transient model by smoothing
the previous states. However, the smoothing model may not reflect the current power
system, which is the moti- vation behind this study.

This paper suggests mathematical distribution power system state transient model by analyzing a nonlinear function of powerw injection model. Measurement model is also described as nonlin- ear function of state. Then, extended Kalman filter is adopted for dynamic state estimation. In this step, jacobian matrices of nonlin- ear functions are derived as a function of system state. To verify proposed algorithm, IEEE 15 test bus with scheduled generator is used.

The notations in this paper is fairly standard. For vector $x, x(i)$ denotes $i-\mathrm{th}$ component of $x$. For matrix $A$, $A_{i j}$ and $[A]_{i j}$ denote $(i, j)-$element of matrix $A$. $\mathbf{0}_{a \times b} \in \mathcal{R}^{a \times b}$ is the zero matrix, $\mathbf{1}_{a \times b} \in$ $\mathcal{R}^{a \times b}$ is the matrix whose all components are one.

## 2. Problem Statements

To use dynamic state estimation algorithm, the state space representation for the distribution system is required. Namely, setting the following system equations is essential.

where $x_{k} \in \mathcal{R}^{n}$ is the state, $z_{k}$ is the measurement, $w_{k}$ and $v_{k}$ are modeling error of system model and measurement model, respec- tively. $f(\cdot)$ and $h(\cdot)$ are model functions. $n$ is the number of nodes in power system. We set the estimation target as complex voltage in the distribution node and use power flows/injections in the distri- bution system as measurement values through phasor measurement units (PMUs). Therefore, $x_{k}$ and $z_{k}$ are defined as follows:

##### (3)

$$ x_{k}=\left[\begin{array}{ccccccc}{\left|V_{1, k}\right|} & {\cdots} & {\left|V_{n, k}\right|} & {\theta_{1, k}} & {\cdots} & {\theta_{n, k}}\end{array}\right]^{T}, $$

##### (4)

$$ z_{k}=\left[\begin{array}{cccc}{P_{I, k}^{T}} & {Q_{I, k}^{T}} & {P_{F, k}^{T}} & {Q_{F, k}^{T}}\end{array}\right]^{T} $$
where the notations $\left|V_{i, k}\right|$ and $\theta_{i, k}$ denote voltage magnitude,
voltage phase of bus $i$ at time $k$. Also, $P_{I, k}\left(Q_{I, k}\right)$ and $P_{F,
k}\left(Q_{F, k}\right)$ mean real part (imaginary part) of power injections and that
of power flows, respectively. While the measurement model $h(\cdot)$ can be derived
by using the topology of power system, the exact model $f(\cdot)$ in (1) is unknown.
Therefore, as an alternative model, the recent papers ^{(7,}^{8)} assumed the quisi steady-state behavior of distribution system and used the following
linear time-varying model.

### 2.1 Previous Works: Linear Time-varying State Transition Model

where $F_{k}$ and $g_{k}$ are the state transition process, and they are up- dated
online as the Holt-Winters double exponential smoothing method ^{(10)}:

##### (7)

$$ \begin{aligned} g_{k}=&\left(1+\beta_{k}\right)\left(1-\alpha_{k}\right) x_{k-1 | k-2}-\beta_{k} \mathbf{a}_{k-2} \\ &+\left(1-\beta_{k}\right) \mathbf{b}_{k-2}, \end{aligned} $$where $I$ is an identity matrix with appropriate dimension, $a_{k}$, $β_{k}$ $∈ (0, 1)$ are smoothing coefficient to be adjusted, $x_{k | k^{\prime}}$ denotes the esti- mation of $x_{k}$ given the measurements up to time $k^{\prime}$. The terms $m_{k}$, $l_{k}$ are recursively updated as

Remark 1: There are some reason why this model is not always robust. First, because the terms (8), (9) are composed of the estimation state $x_{k | k-1}$, the model accuracy may drop if the estimation state $x_{k | k-1}$ is not correct. Second, it is hard to adjust the adequate smoothing coefficient $\alpha_{k}, \beta_{k}$ for each time. Finally, quasi steady-state behavior of distribution system cannot be always assumed since power systems in nowadays can have some abrupt load changes. Therefore, we propose the following power system model which comes from mathematical analysis.

## 3. Distribution system modeling

In this section, the nonlinear models of measurement system and state transient system are provided via mathematical approach.

### 3.1 Measurement system model

We use both power injection and flow as measurement values. In distribution power system, they are represented as complex values and can be modeled by using power flow equation. First, power injection model is represented as follows:

##### (10)

$$ \begin{aligned} P_{I, k}(i)=& \sum_{j=1}^{2 n}\left|V_{i, k}\right|\left|V_{j, k}\right|\left\{G_{i j} \cos \left(\theta_{i, k}-\theta_{j, k}\right)\right.\\ &+B_{i j} \sin \left(\theta_{i, k}-\theta_{j, k}\right) \},\end{aligned} $$

##### (11)

$$ \begin{aligned} Q_{I, k}(i)=& \sum_{j=1}^{2 n}\left|V_{i, k}\right|\left|V_{j, k}\right|\left\{G_{i j} \sin \left(\theta_{i, k}-\theta_{j, k}\right)\right.\\ &-B_{i j} \cos \left(\theta_{i, k}-\theta_{j, k}\right) \}, \end{aligned} $$where $P_{I, k}(i)\left(Q_{I, k}(i)\right)$ is real(imaginary) part of power injection in $i-\mathrm{th}$ node at time $k$, and $G_{i, j}\left(B_{i, j}\right)$ is real(imaginary) part of $(i, j)-$component of $\boldsymbol{Y}_{b u s}$, which is an admittance matrix determined by distribution power system topology.

Next, power flow model is obtained such that:

##### (12)

$$ \begin{aligned} P_{i, j, k}=&\left|V_{i, k}\right|^{2} g_{i j}-\left|V_{i, k}\right|\left|V_{j, k}\right|\left\{g_{i j} \cos \left(\theta_{i, k}-\theta_{j, k}\right)\right.\\ &+b_{i j} \sin \left(\theta_{i, k}-\theta_{j, k}\right) \}, \end{aligned} $$

##### (13)

$$ \begin{aligned} Q_{i, j, k}=&-\left|V_{i, k}\right|^{2}\left(b_{i j}+b_{i i}\right)-\left|V_{i, k}\right|\left|V_{j, k}\right|\left\{g_{i j} \sin \left(\theta_{i, k}-\theta_{j, k}\right)\right.\\ &-b_{i j} \cos \left(\theta_{i, k}-\theta_{j, k}\right) \}, \end{aligned} $$where $P_{i, j, k}\left(Q_{i, j, k}\right)$ is real(imaginary) part of power flow from node $i$ to node $j$ at time $k$, and $g_{i j}\left(b_{i j}\right)$ is real(imaginary) part of the ad-mittance between node $i$ and node $j$. Therefore, the following mea-surement model is used in this paper.

where $i \in[1,2, \cdots, n], t \in m_{i}$ and $m_{l}$ denotes index that is linked to node $l$ and larger than $l$, such that $m_{l} \in\left\{i | l<i \leq n, Y_{b u s}(l, i) \neq 0\right\}$. The modeling error wk is obtained as offline, and it is assumed that the white Gaussian noise vector with known covariance $R$. Additionally, the dimension of $z_{k}$ is determined by distribution system topology.

### 3.2 State transient system

As mentioned in previous section, power injection model is rep- resented as a function of system state, and it has same dimension as system state.

##### (16)

$$ g\left(x_{k}\right)=\left[\begin{array}{c}{P_{I, k}} \\ {Q_{I, k}}\end{array}\right] \triangleq y_{k}, $$where nonlinear function $g(\cdot)$ is constructed by using the equations (10)-(11). Then taking partial derivative of $x_{k}$ to both sides earns following relations:

##### (17)

$$ \frac{\partial y_{k}}{\partial x_{k}}=\frac{\partial g(x)}{\partial x} \left.\right|_{x=x_{k}} \triangleq J_{k}\left(x_{k}\right)‧ $$By setting $\partial y_{k} \approx y_{k+1}-y_{k}$ and $\partial x_{k} \approx x_{k+1}-x_{k}, x_{k+1}$ is calculated as follows:

(18), additional assumption is used such that

which is valid under the $PMU$. The matrix $J_{k}\left(x_{k}\right)$, which is a jaco-bian matrix of $g\left(x_{k}\right)$, is defined as

##### (20)

$$ J_{k}\left(x_{k}\right)=\left[\begin{array}{cc}{\frac{\partial P_{I, k}}{\partial\left|V_{k}\right|}} & {\frac{\partial P_{I, k}}{\partial \theta_{k}}} \\ {\frac{\partial Q_{L, k}}{\partial\left|V_{k}\right|}} & {\frac{\partial Q_{I, k}}{\partial \theta_{k}}}\end{array}\right]‧ $$
It is known that active power injection $P_{I, k}$ and reactive power injection $Q_{I,
k}$ sensitivity of voltage magnitude $\left|V_{k}\right|$ and voltage phase $\theta_{k}$
is very small ^{(11)}. Therefore, by considering the computation efficiency, we suppose that $\frac{\partial
P_{I, k}}{\partial\left|V_{k}\right|}=0$ and $\frac{\partial Q_{I, k}}{\partial \theta_{k}}=0$.
This supposition reduces the computational complexity of $J_{k}^{-1}\left(x_{k}\right)$
from $O\left(8 n^{3}\right)$ to 2$O\left(n^{3}\right)$. Then state transient model
of distribution system is determined as follows:

##### (22)

$$ \begin{aligned} \triangleq x_{k}+\left[\begin{array}{cc}{\mathbf{0}_{n \times n}} & {J_{2}^{-1}\left(x_{k}\right)} \\ {J_{1}^{-1}\left(x_{k}\right)} & {\mathbf{0}_{n \times n}}\end{array}\right]\left[\begin{array}{c}{I_{2 n} | \mathbf{0}_{2 n \times 2 m}}\end{array}\right] \\ \times\left\{z_{k}-z_{k-1}\right\}+w_{k} \end{aligned} $$where $I_{a} \in \mathcal{R}^{a \times a}$ is the identity matrix, and $w_{k}$ is state transition modeling error which is assumed to be white Gaussian noise vector with known covariance matrix Q. The components of matrices $J_{1}\left(x_{k}\right), J_{2}\left(x_{k}\right)$ are described in (23), (24). The detailed proof is omitted due to lack of page limit.

##### (23)

$$ \left[J_{1}\left(x_{k}\right)\right]_{i, j}=\left\{\begin{array}{c}{-x_{k}(i) x_{k}(j)\left\{G_{i j} \sin \left(x_{k}(n+i)-x_{k}(n+j)\right)\right.} \\ {-B_{i j} \cos \left(x_{k}(n+i)-x_{k}(n+j)\right) \} \text { if } i \neq j}, \\ {-x_{k}(i) \sum_{j=1, j \neq i}^{n} x_{k}(j)\left\{G_{i j} \sin \left(x_{k}(n+i)-x_{k}(n+j)\right)\right.} \\ {-B_{i j} \cos \left(x_{k}(n+i)-x_{k}(n+j)\right) \} \text { otherwise }},\end{array}\right. $$

##### (24)

$$ \left[J_{2}\left(x_{k}\right)\right]_{i, j}=\left\{\begin{array}{c}{-x_{k}(i)\left\{G_{i j} \sin \left(x_{k}(n+i)-x_{k}(n+j)\right)\right.} \\ {-B_{i j} \cos \left(x_{k}(n+i)-x_{k}(n+j)\right) \} \text { if } i \neq j}, \\ {\sum_{j=1, j \neq i}^{n} x_{k}(j)\left\{G_{i j} \sin \left(x_{k}(n+i)-x_{k}(n+j)\right)\right.} \\ {-B_{i j} \cos \left(x_{k}(n+i)-x_{k}(n+j)\right) \}} \\ {-2 x_{k}(i) B_{i i}} & {\text { otherwise }},\end{array}\right. $$

##### (25)

$$ \left[H_{1}\left(x_{k}\right)\right]_{i, j}=\left\{\begin{array}{c}{x_{k}\left(i_{x}\right) x_{k}\left(i_{y}\right)\left\{g_{i_{i} i_{y}} \sin \left(x_{k}\left(n+i_{x}\right)-x_{k}\left(n+i_{y}\right)\right)\right.} \\ {-b_{i_{i} i_{y}} \cos \left(x_{k}\left(n+i_{x}\right)-x_{k}\left(n+i_{y}\right)\right) \} \text { if } j=i_{x}}, \\ {-x_{k}\left(i_{x}\right) x_{k}\left(i_{y}\right)\left\{g_{i_{i} i_{y}} \sin \left(x_{k}\left(n+i_{x}\right)-x_{k}\left(n+i_{y}\right)\right)\right.} \\ {-b_{i_{i} i_{y}} \cos \left(x_{k}\left(n+i_{x}\right)-x_{k}\left(n+i_{y}\right)\right) \} \text { if } j=i_{y}}, \\ {0 \text { otherwise }},\end{array}\right. $$

##### (26)

$$ \left[H_{2}\left(x_{k}\right)\right]_{i, j}=\left\{\begin{array}{c}{-2 x_{k}\left(i_{x}\right)\left\{b_{i_{i} i_{x}}+b_{i_{i} i_{y}}\right\}} \\ {-x_{k}\left(i_{y}\right)\left\{g_{i_{x} i_{y}} \sin \left(x_{k}\left(n+i_{x}\right)-x_{k}\left(n+i_{y}\right)\right)\right.} \\ {\cdots} \\ {-b_{i_{t} i_{y}} \cos \left(x_{k}\left(n+i_{x}\right)-x_{k}\left(n+i_{y}\right)\right) \} \text { if } j=i_{x}} \\ {-x_{k}\left(i_{x}\right)\left\{g_{i_{x} i_{y}} \sin \left(x_{k}\left(n+i_{x}\right)-x_{k}\left(n+i_{y}\right)\right)\right.} \\ {-b_{i_{x} i_{y}} \cos \left(x_{k}\left(n+i_{x}\right)-x_{k}\left(n+i_{y}\right)\right) \} \text { if } j=i_{y}} \\ {0 \quad \text { otherwise }},\end{array}\right. $$where $\mathbf{m}=\left\{(x, y) | x \in[1,2, \cdots, n], y \in m_{i}\right\}$ and $\mathbf{m}_{i} \triangleq\left(i_{x}, i_{y}\right)$.

### 3.3 Modeling result

By using a simulation data from Matlab Simulink of IEEE 15 radial test bus, we verify the effect of the proposed models. Fig. 1 shows the actual states from Matlab simulink (red-solid) and cal- culated states from the proposed model in (18) (blue-dashed). The used data in (18) is sampled as 30Hz. Fig. 2 illustrates the actual measurement from Matlab simulink (red-solid) and calculated val- ues by using states in Fig. 1 (blue-dashed). It looks like that the measurement error which is the difference between actual value of measurement and calculated measurement value follows white Gaussian distribution. Thus, this model can be used to extended kalman filter.

## 4. Dynamic state estimation algorithm for power system

represented as nonlinear function and model errors are assumed as white Gaussian
noise vector, extended Kalman filter (EKF) can be adopted ^{(12)}.

### 4.1 Jacobian Calculation

Jacobians of nonlinear functions $f(\cdot), h(\cdot)$ are defined such that

##### (27)

$$ F_{k}=\frac{\partial f(x)}{\partial x} | x=x_{k-1 | k-1}, H_{k}=\frac{\partial h}{\partial x} | x=x_{k-1 | k-1}‧ $$First, $H_{k}$ is derived as follows:

##### (28)

$$ H_{k} \triangleq\left[\begin{array}{cc}{\mathbf{0}_{n \times n}} & {J_{1}\left(x_{k}\right)} \\ {J_{2}\left(x_{k}\right)} & {\mathbf{0}_{n \times n}} \\ {\mathbf{0}_{n \times n}} & {H_{1}\left(x_{k}\right)} \\ {H_{2}\left(x_{k}\right)} & {\mathbf{0}_{n \times n}}\end{array}\right], $$where each component of $J_{1}\left(x_{k}\right), J_{2}\left(x_{k}\right), H_{1}\left(x_{k}\right), H_{2}\left(x_{k}\right)$ is given in (23),(24),(25),(26). Here, the assumption that $\frac{\partial P_{F, k}}{\partial\left|V_{k}\right|}=0$ and $\frac{\partial Q_{F, k}}{\partial \theta_{k}}=0$ is still valid. Next, it is hard to calculate $F_{k}$ because $f\left(x_{k}, z_{k}, z_{k}-1\right)$ consists of inverse of $J_{1}\left(x_{k}\right), J_{2}\left(x_{k}\right)$. Therefore, instead of directly calculating partial derivative of $f\left(x_{k}, z_{k}, z_{k}-1\right)$, the complex step differentiation algorithm is used.

Lemma 1 ^{(13)}: Let $f(z)$ be a function with complex variable $z=r_{0}+i h$, where $r_{0}$ is a
point on the real axis and $h$ is a real parameter. By taking Taylor series expansion
on $f(z)$ off the real axis, then

where $f^{(n)}(\cdot)$ is $n-\mathrm{th}$ order partial differential function. After that, taking the imaginary part of both sides, dividing by $h$, and setting $h$ as very small scalar holds the following relation:

where $\mathfrak{I}\{Z\}$ denotes the imaginary part of $Z$.

### 4.2 State Prediction

This stage is to propagate the state into the next time step by using state transition model:

where $P_{a | b}$ is the covariance matrix of predicted state between time step $a$ and $b$, and $Q$ is covariance matrix of transition model noise. $x_{k | k-1}$ is predicted state estimation at time step $k − 1$ and also called as state forecast value.

### 4.3 State Update

After acquiring measurement $z_{k}$, state estimation can be updated:

##### (33)

$$ \begin{aligned} x_{k | k} &=x_{k | k-1}+K_{k}\left\{z_{k}-h\left(x_{k | k-1}\right)\right\} \end{aligned}, $$

##### (35)

$$ \begin{aligned} K_{k} &=P_{k | k-1} H_{k}^{T}\left(H_{k} P_{k | k-1} H_{k}^{T}+R\right)^{-1} \end{aligned}, $$where $K_{k}$ is Kalman filter gain and $R$ is covariance matrix of mea- surement noise, which is estimated as off-line.

## 5. Simulation Results

To verify the validity of the proposed dynamic model (21) and the effectiveness of dynamic state estimation of distribution power system, IEEE 15-bus radial test system is considered. In Fig. 3, bus 1 is substation, where its voltage phase can act a role of reference, and there are 3 generators (G1,G2,G3) which can be scheduled by system operator. To show the robustness of the model under the suddenly changes in power system, power generation from G1-G3 is scheduled as (Fig. 4). The measurements are obtained as MATLAB simulink model, and the sampling frequency of measuring is 30Hz, which is same as PMU data.

By conducting several simulations under difference generator scheduling data, and taking its average values, we also verify that the modeling error has a pattern of normal distribution, such that $w_{k} \mathcal{N}(\overline{\mathbf{m}}, \sigma)$, where $\overline{\mathbf{m}}=\left[\begin{array}{ll}{\mathbf{0}_{1 \times 15}} & {3 \cdot 10^{-3} \cdot \mathbf{1}_{1 \times 15}}\end{array}\right]^{T} \approx\left[\mathbf{0}_{30 \times 1}\right]$, $\sigma=$ $\left[4 \cdot 10^{-4} \cdot 1_{1 \times 15} \quad 4.5 \cdot 10^{-3} \cdot 1_{1 \times 15}\right]^{T}$, and the unit is p.u for voltage magnitude and degree for voltage phase.

Next, state estimation result by using proposed model and EKF algorithm is shown in Fig. 5. They show the voltage magnitude and phase of node 13. The actual signal is acquired from Matlab simulink.

## 6. Conclusion

This paper suggested the distribution system state transient model by using mathematical analysis of nonlinear functions, whereas the existing papers have been used only smoothing model. Also, setting the measurement model as a set of nonlinear function of system state was done for dynamic state estimation for power system. Jacobian matrices for nonlinear state transient model and measurement model were calculated to apply extended Kalman fil- ter algorithm. The validity of proposed modeling and dynamic es- timation algorithm were verified by using simulation data of IEEE 15 bus radial test system. Our future work is extending our result under the bad data problem and cyber attack.

### References

## 저자소개

She received her B.S. degree from School of Electrical En- gineering and Computer Science at Ulsan National Insti- tute of Science and Technology (UNIST), Korea, in 2014, and the M.S. degree in the Division of Electrical and Com- puter Engineering at Pohang University of Science and Technology (POSTECH), Korea, in 2016.

She is currently a Ph. D candidate in Electric and Electrical Engineering at POSTECH. Her research interest includes stochastic system, singular system and robust control theory.

E-mail: chane2727@postech.ac.kr

He received his B.S. degree from School of Electrical En- gineering and Computer Science at Kyungpook National University (KNU), Korea, in 2014, and the M.S. degree in the Division of Electrical and Computer Engineering at Po- hang University of Science and Technology (POSTECH), Korea, in 2016.

He is currently a Ph. D candidate in Elec- tric and Electrical Engineering at POSTECH.

His research interest includes fuzzy system, stochastic system and ro- bust control theory.

E-mail: ispark90@postech.ac.kr

He received his B.S. degree and M.S. degree in Control and Instrumentation Engineering from Seoul National Uni- versity, Korea, in 1988 and 1990, respectively, and the Ph.D degree in Electrical Engineering from Stanford Uni- versity, U.S.A., in 1995.

Since 1996, he has been affili- ated with the Division of Electrical and Computer Engi- neering at Pohang University of Science and Technology (POSTECH), where he is currently a professor.

His cur- rent research interests include robust, LPV, and network- related control theories, delayed systems, fuzzy systems, and signal processing.

E-mail: ppg@postech.ac.kr