# Accurate Parameter Estimation for Unbalanced Three-Phase System

Authors: Yuan Chen; Hing Cheung So
Journal: The Scientific World Journal
 (2014)

---

## Abstract

Smart grid is an intelligent power generation and control console in modern electricity networks, where the unbalanced three-phase power system is the commonly used model. Here, parameter estimation for this system is addressed. After converting the three-phase waveforms into a pair of orthogonal signals via theα
β-transformation, the nonlinear least squares (NLS) estimator is developed for accurately finding the frequency, phase, and voltage parameters. The estimator is realized by the Newton-Raphson scheme, whose global convergence is studied in this paper. Computer simulations show that the mean square error performance of NLS method can attain the Cramér-Rao lower bound. Moreover, our proposal provides more accurate frequency estimation when compared with the complex least mean square (CLMS) and augmented CLMS.


---

## Body

## 1. Introduction

Due to the increasing demand for electricity and the finite supply of nonrenewable energy sources, electrical power generation systems have faced a huge challenge. In order to improve the efficiency and reliability and to reduce the cost of electricity network, the concept of smart grid [1] is proposed, which can utilize the renewable and sustainable resources such as wind and solar energies. Traditional grid is a star network including a central point (e.g., power generation station) and leaf nodes (e.g., user terminals), whereas the smart grid is a mesh network whose nodes can act as both users and generators. In smart grid, because of this role conversion and the operations in terminals [2], frequency and amplitude variations exist and they can cause many serious problems such as loss of synchronism, power system stabilization, and equipment connection [3]. As a result, it is important to monitor the variations via accurately estimating the corresponding parameters [4].It is common to use the unbalanced three-phase power system [5] for modeling in smart grid applications. Conventional estimators exhibit poor behavior applied for the three-phase system directly, because they work well only for the single-phase signal, which cannot truly characterize the unbalanced system [6]. Although we can perform estimation in each phase separately, accurate results may not be obtained because fixed phase displacement does not hold. Nevertheless, by making use of the α
β-transformation [7], the three-phase waveforms can be mapped into a pair of in-phase and quadrature signals. Based on this model, a number of approaches for frequency estimation have been proposed, including the complex least mean square (CLMS) [8] and augmented complex least mean square (ACLMS) [9] methods. Nevertheless, both estimators only focus on finding the frequency of the unbalanced system. In this work, we contribute to the development of an optimal estimator for the frequency, phase, and amplitudes from the orthogonal signals.The rest of this paper is organized as follows. In Section2, the problem is formulated and then nonlinear least squares (NLS) estimator is devised. We apply the Newton-Raphson scheme to solve the corresponding nonlinear optimization problem where algorithm initialization and global convergence are examined. Computer simulations are included in Section 3, which show that the mean square error (MSE) performance of NLS method can attain the optimum benchmark of the Cramér-Rao lower bound (CRLB) in the presence of white Gaussian disturbances and its superiority over the CLMS and ACLMS algorithms in frequency estimation is demonstrated. Finally, conclusions are drawn in Section 4.


## 2. Proposed Method

### 2.1. Development and Convergence Analysis

The discrete-time observations of the unbalanced three-phase power system are modeled as [5]:
(1)
v
a
[
n
]
=
V
a
cos
⁡
(
ω
n
+
ϕ
)
+
ξ
a
[
n
]
,
v
b
[
n
]
=
V
b
cos
⁡
(
ω
n
+
ϕ
-
2
π
3
)
+
ξ
b
[
n
]
,
v
c
[
n
]
=
V
c
cos
⁡
(
ω
n
+
ϕ
+
2
π
3
)
+
ξ
c
[
n
]
,
where V
a, V
b, and V
c are the inequivalent amplitudes of different phase components, ω
=
Ω
/
F
s is the discrete frequency with Ω and F
s being the voltage frequency in radian and sampling frequency in Hz, respectively, and ϕ is the initial phase. The nominal value of Ω is Ω
*
=
100
π (or 120
π) rads
-
1. According to [10], the noise terms ξ
a
[
n
], ξ
b
[
n
], and ξ
c
[
n
], are independent and identically distributed additive white Gaussian noise sequences with same variance σ
2. The task is to find the unknown parameters, namely, ω, ϕ, V
a, V
b, and V
c. In this study, we apply the α
β-transformation [7] on (1) to achieve accurate parameter estimation. The transformed signals, denoted by v
α
[
n
] and v
β
[
n
], are computed as
(2)
[
v
α
[
n
]
v
β
[
n
]
]
=
2
3
[
1
-
1
2
-
1
2
0
3
2
-
3
2
]
[
v
a
[
n
]
v
b
[
n
]
v
c
[
n
]
]
.
Based on (1)-(2), v
α
[
n
] and v
β
[
n
] can also be expressed as
(3)
v
α
[
n
]
=
s
α
[
n
]
+
q
α
[
n
]
,
v
β
[
n
]
=
s
β
[
n
]
+
q
β
[
n
]
,
where
(4)
s
α
[
n
]
=
(
A
cos
⁡
(
ϕ
)
-
B
sin
(
ϕ
)
)
cos
⁡
(
ω
n
)
-
(
A
sin
(
ϕ
)
+
B
cos
⁡
(
ϕ
)
)
sin
(
ω
n
)
,
s
β
[
n
]
=
(
-
B
cos
⁡
(
ϕ
)
+
C
sin
(
ϕ
)
)
cos
⁡
⁡
(
ω
n
)
+
(
B
sin
(
ϕ
)
+
C
cos
⁡
(
ϕ
)
)
sin
(
ω
n
)
,
with
(5)
A
=
6
12
(
4
V
a
+
V
b
+
V
c
)
,
B
=
2
4
(
V
b
-
V
c
)
,
C
=
6
4
(
V
b
+
V
c
)
,
q
α
[
n
]
=
2
3
(
ξ
a
[
n
]
-
1
2
ξ
b
[
n
]
-
1
2
ξ
c
[
n
]
)
,
q
β
[
n
]
=
2
2
(
ξ
b
[
n
]
-
ξ
c
[
n
]
)
.
Although both q
α
[
n
] and q
β
[
n
] contain ξ
b
[
n
] and ξ
c
[
n
], it is easy to show that the noise terms are uncorrelated; that is, E
{
q
α
[
n
]
q
β
[
n
]
}
=
0 where E denotes the expectation operator, and they have identical variance σ
2.Assuming that we haveN samples for each channel, (3) can be written in matrix form as follows:
(6)
v
=
s
+
q
,
where
(7)
v
=
[
v
α
[
1
]
v
α
[
2
]
⋯
v
α
[
N
]
v
β
[
1
]
v
β
[
2
]
⋯
v
β
[
N
]
]
T
,
s
=
[
s
α
[
1
]
s
α
[
2
]
⋯
s
α
[
N
]
s
β
[
1
]
s
β
[
2
]
⋯
s
β
[
N
]
]
T
,
s
=
G
x
,
G
=
[
H
0
0
H
]
,
H
=
[
cos
⁡
(
ω
)
sin
(
ω
)
cos
⁡
(
2
ω
)
sin
(
2
ω
)
⋮
⋮
cos
⁡
(
N
ω
)
sin
(
N
ω
)
]
,
x
=
[
A
cos
⁡
(
ϕ
)
-
B
sin
(
ϕ
)
-
A
sin
(
ϕ
)
-
B
cos
⁡
(
ϕ
)
-
B
cos
⁡
(
ϕ
)
+
C
sin
(
ϕ
)
B
sin
(
ϕ
)
+
C
cos
⁡
(
ϕ
)
]
,
q
=
[
q
α
[
1
]
q
α
[
2
]
⋯
q
α
[
N
]
q
β
[
1
]
q
β
[
2
]
⋯
q
β
[
N
]
]
T
.
Here, T denotes the transpose operator and0 is the N
×
2 zero matrix. We see that x corresponds to the linear unknowns, while ω is the nonlinear unknown in (6). Employing NLS [11, 12], the estimates of ω and x, denoted by ω
^ and x
^, are
(8)
{
ω
^
,
x
^
}
=
arg
min
⁡
ω
,
x
J
(
ω
,
x
)
,
J
(
ω
,
x
)
=
(
v
-
G
x
)
T
(
v
-
G
x
)
.Based on the Newton-Raphson procedure, the updating rule forω
^ is
(9)
ω
^
(
k
+
1
)
=
ω
^
(
k
)
-
∇
ω
J
(
ω
^
(
k
)
,
x
^
(
k
)
)
∇
ω
2
J
(
ω
^
(
k
)
,
x
^
(
k
)
)
,
where
(10)
∇
ω
J
(
ω
,
x
)
=
-
(
G
1
x
)
T
(
v
-
G
x
)
,
∇
ω
2
J
(
ω
,
x
)
=
x
T
(
G
1
T
G
1
)
x
+
(
G
2
x
)
T
(
v
-
G
x
)
,
with
(11)
G
1
=
[
H
1
0
0
H
1
]
,
H
1
=
[
-
sin
(
ω
)
cos
⁡
(
ω
)
-
2
sin
(
2
ω
)
2
cos
⁡
(
2
ω
)
⋮
⋮
-
N
sin
(
N
ω
)
N
cos
⁡
(
N
ω
)
]
,
G
2
=
[
H
2
0
0
H
2
]
,
H
2
=
[
cos
⁡
(
ω
)
sin
(
ω
)
4
cos
⁡
(
2
ω
)
4
sin
(
2
ω
)
⋮
⋮
N
2
cos
⁡
(
N
ω
)
N
2
sin
(
N
ω
)
]
.
Here, ω
^
(
k
) and x
^
(
k
) are the estimates of ω and x at the kth iteration. Once we have ω
^
(
k
), x
^
(
k
) is easily obtained from (8) as
(12)
x
^
(
k
)
=
(
G
T
G
)
-
1
G
T
v
|
ω
=
ω
^
(
k
)
,
where −1 denotes the matrix inverse. To start the algorithm of (9) and (12), we need ω
^
(
0
). Noting that ω should be around its nominal value ω
*
=
Ω
*
/
F
s, that is, ω
∈
[
ω
*
-
τ
,
ω
*
+
τ
] where τ is the maximum deviation from ω
*, ω
^
(
0
) is computed using grid search as follows. We assign K uniformly-spaced grid points in the range [
ω
*
-
τ
,
ω
*
+
τ
] where one of them is ω
^
(
0
). For each possible candidate ω
^
(
0
), we determine x
^
(
0
) according to (12). The pair {
ω
^
(
0
)
,
x
^
(
0
)
} which gives the minimum value of J
(
ω
,
x
) will be chosen as the initial guess for (9). In our study, the iterative algorithm is terminated when |
(
ω
^
(
k
+
1
)
-
ω
^
(
k
)
)
/
(
ω
^
(
k
+
1
)
)
|
<
ϵ, where ϵ
>
0 is a small tolerance constant, is reached. After obtaining ω
^, the NLS estimates of ϕ, V
a, V
b, and V
c are straightforwardly computed from x
^ as
(13)
ϕ
^
=
tan
⁡
-
1
(
x
^
3
-
x
^
2
x
^
1
+
x
^
4
)
,
V
^
a
=
(
6
2
x
^
1
-
6
6
x
^
3
)
cos
⁡
⁡
(
ϕ
^
)
-
(
6
2
x
^
2
+
6
6
x
^
3
)
sin
(
ϕ
^
)
,
V
^
b
=
(
x
^
4
-
2
x
^
2
)
cos
⁡
(
ϕ
^
)
+
(
6
3
x
^
3
-
2
x
^
1
)
sin
(
ϕ
^
)
,
V
^
c
=
(
x
^
4
+
2
x
^
2
)
cos
⁡
(
ϕ
^
)
+
(
6
3
x
^
3
+
2
x
^
1
)
sin
(
ϕ
^
)
,
where x
^
m (m
=
1,2
,
3,4) denotes the mth element of x
^.Finally, we examine the convergence of (9). According to [13], global convergence with quadratic rate is guaranteed when M
|
ω
-
ω
^
(
0
)
|
<
1 is satisfied, where M is as follows:
(14)
M
=
sup
⁡
ω
∈
{
-
π
,
π
}
|
∇
ω
3
J
(
ω
,
x
)
∇
ω
2
J
(
ω
,
x
)
|
.
To determine the value of M, we first relate G, G
1, and G
2 as follows:
(15)
G
1
T
G
1
+
G
2
T
G
=
c
E
,
G
1
=
N
1
G
F
,
G
2
=
N
2
G
,
where
(16)
E
=
[
1
0
0
0
1
0
0
2
0
]
,
F
=
[
0
-
1
0
0
1
0
0
0
0
0
0
-
1
0
0
1
0
]
,
c
=
N
(
N
+
1
)
(
2
N
+
1
)
6
,
N
1
=
diag
⁡
(
[
1
2
⋯
N
1
2
⋯
N
]
)
,
N
2
=
diag
⁡
(
[
1
2
2
⋯
N
2
1
2
2
⋯
N
2
]
)
.
Based on (10)-(11), ∇
ω
2
J
(
ω
,
x
) can be expressed as
(17)
∇
ω
2
J
(
ω
,
x
)
=
x
T
(
c
E
-
G
2
T
G
)
x
+
x
T
G
2
T
(
v
-
G
x
)
≥
x
T
G
T
N
2
v
-
2
x
T
G
T
N
2
G
x
≥
x
T
G
T
N
2
v
-
2
N
2
x
T
(
G
T
G
)
x
.
Equality holds if and only if x
1
2
+
x
2
2
+
2
x
2
x
3
=
0, ω
=
m
π with m
=
-
1
,
-
1
/
2
,
0
,
1
/
2
,
1 and
(18)
∇
ω
3
J
(
ω
,
x
)
=
x
T
G
1
T
N
2
v
=
x
T
F
T
G
T
N
3
v
,
where N
3
=
diag
⁡
(
[
1
2
3
⋯
N
3
1
2
3
⋯
N
3
]
). Substituting (12) into (17)-(18) yields
(19)
∇
ω
2
J
(
ω
,
x
)
≥
v
T
(
G
(
G
T
G
)
-
1
G
T
)
N
2
v
-
2
N
2
v
T
(
G
(
G
T
G
)
-
1
G
T
)
v
≥
v
T
W
1
v
,
∇
ω
3
J
(
ω
,
x
)
=
v
T
W
2
v
,
where
(20)
W
1
=
G
(
G
T
G
)
-
1
G
T
(
N
2
-
2
N
2
)
,
W
2
=
G
(
G
T
G
)
-
1
F
T
G
T
N
3
.
We can then write M as
(21)
M
=
sup
⁡
ω
∈
{
-
π
,
π
}
{
|
∇
ω
3
J
(
ω
,
x
)
v
T
v
|
|
1
(
∇
ω
2
J
(
ω
,
x
)
)
/
(
v
T
v
)
|
}
=
|
λ
2
λ
1
|
,
where λ
1 and λ
2 are the minimum eigenvalue of W
1 and the maximum eigenvalue of W
2, respectively. It is easily shown that λ
1
=
N
3 and λ
2
=
N
2. Hence M
=
N. As a result, if the initial estimate is chosen such that N
|
ω
-
ω
^
(
0
)
|
<
1 is satisfied, global solution will be obtained.


### 2.2. Cramér-Rao Lower Bound

Let the unknown parameter vector beθ
=
[
ω
ϕ
V
a
V
b
V
c
]
T. Then, the CRLB of θ
^ is obtained from the diagonal elements of the inverse of the Fisher information matrix, denoted by F
(
θ
) [14]. The (
i
,
j
) entry of F
(
θ
) is written as
(22)
[
F
(
θ
)
]
i
,
j
=
(
∂
s
∂
θ
i
)
T
C
-
1
(
∂
s
∂
θ
j
)
,
i
,
j
=
1
,
…
,
5
,
where(23)
∂
s
∂
θ
i
=
[
∂
s
α
[
1
]
∂
θ
i
∂
s
α
[
2
]
∂
θ
i
⋯
∂
s
α
[
N
]
∂
θ
i
∂
s
β
[
1
]
∂
θ
i
∂
s
β
[
2
]
∂
θ
i
⋯
∂
s
β
[
N
]
∂
θ
i
]
T
,with
(24)
∂
s
α
[
n
]
∂
ω
=
-
n
(
A
cos
⁡
(
ϕ
)
-
B
sin
(
ϕ
)
)
sin
(
ω
n
)
-
n
(
A
sin
(
ϕ
)
+
B
cos
⁡
(
ϕ
)
)
cos
⁡
(
ω
n
)
,
∂
s
β
[
n
]
∂
ω
=
n
(
B
cos
⁡
(
ϕ
)
-
C
sin
(
ϕ
)
)
sin
⁡
(
ω
n
)
+
n
(
B
sin
(
ϕ
)
+
C
cos
⁡
(
ϕ
)
)
cos
⁡
(
ω
n
)
,
∂
s
α
[
n
]
∂
ϕ
=
(
-
A
sin
(
ϕ
)
-
B
cos
⁡
(
ϕ
)
)
cos
⁡
⁡
(
ω
n
)
-
(
A
cos
⁡
(
ϕ
)
-
B
sin
(
ϕ
)
)
sin
(
ω
n
)
,
∂
s
β
[
n
]
∂
ϕ
=
(
B
sin
(
ϕ
)
+
C
cos
⁡
(
ϕ
)
)
cos
⁡
⁡
(
ω
n
)
+
(
B
cos
⁡
(
ϕ
)
-
C
sin
(
ϕ
)
)
sin
(
ω
n
)
,
∂
s
α
[
n
]
∂
V
a
=
6
3
(
cos
⁡
(
ϕ
)
cos
⁡
(
ω
n
)
-
sin
(
ϕ
)
sin
(
ω
n
)
)
,
∂
s
β
[
n
]
∂
V
a
=
0
,
∂
s
α
[
n
]
∂
V
b
=
(
6
12
cos
⁡
⁡
(
ϕ
)
-
2
4
sin
⁡
(
ϕ
)
)
cos
⁡
⁡
(
ω
n
)
-
(
6
12
sin
(
ϕ
)
+
2
4
cos
⁡
(
ϕ
)
)
sin
(
ω
n
)
,
∂
s
β
[
n
]
∂
V
b
=
(
-
2
4
cos
⁡
⁡
(
ϕ
)
+
6
4
sin
⁡
(
ϕ
)
)
cos
⁡
⁡
(
ω
n
)
+
(
2
4
sin
(
ϕ
)
+
6
4
cos
⁡
(
ϕ
)
)
sin
(
ω
n
)
,
∂
s
α
[
n
]
∂
V
c
=
(
6
12
cos
⁡
⁡
(
ϕ
)
+
2
4
sin
⁡
(
ϕ
)
)
cos
⁡
⁡
(
ω
n
)
-
(
6
12
sin
(
ϕ
)
-
2
4
cos
⁡
(
ϕ
)
)
sin
(
ω
n
)
,
∂
s
β
[
n
]
∂
V
c
=
(
2
4
cos
⁡
⁡
(
ϕ
)
+
6
4
sin
⁡
(
ϕ
)
)
cos
⁡
⁡
(
ω
n
)
+
(
-
2
4
sin
(
ϕ
)
+
6
4
cos
⁡
(
ϕ
)
)
sin
(
ω
n
)
and C is the covariance matrix of the noise term which is C
=
σ
2
I with I denoting the 2
N
×
2
N identity matrix. The (
i
,
j
) element can be simplified as
(25)
[
F
(
θ
)
]
i
,
j
=
1
σ
2
∑
n
=
1
N
‍
(
∂
s
α
[
n
]
∂
θ
i
∂
s
α
[
n
]
∂
θ
j
+
∂
s
β
[
n
]
∂
θ
i
∂
s
β
[
n
]
∂
θ
j
)
.


## 2.1. Development and Convergence Analysis

The discrete-time observations of the unbalanced three-phase power system are modeled as [5]:
(1)
v
a
[
n
]
=
V
a
cos
⁡
(
ω
n
+
ϕ
)
+
ξ
a
[
n
]
,
v
b
[
n
]
=
V
b
cos
⁡
(
ω
n
+
ϕ
-
2
π
3
)
+
ξ
b
[
n
]
,
v
c
[
n
]
=
V
c
cos
⁡
(
ω
n
+
ϕ
+
2
π
3
)
+
ξ
c
[
n
]
,
where V
a, V
b, and V
c are the inequivalent amplitudes of different phase components, ω
=
Ω
/
F
s is the discrete frequency with Ω and F
s being the voltage frequency in radian and sampling frequency in Hz, respectively, and ϕ is the initial phase. The nominal value of Ω is Ω
*
=
100
π (or 120
π) rads
-
1. According to [10], the noise terms ξ
a
[
n
], ξ
b
[
n
], and ξ
c
[
n
], are independent and identically distributed additive white Gaussian noise sequences with same variance σ
2. The task is to find the unknown parameters, namely, ω, ϕ, V
a, V
b, and V
c. In this study, we apply the α
β-transformation [7] on (1) to achieve accurate parameter estimation. The transformed signals, denoted by v
α
[
n
] and v
β
[
n
], are computed as
(2)
[
v
α
[
n
]
v
β
[
n
]
]
=
2
3
[
1
-
1
2
-
1
2
0
3
2
-
3
2
]
[
v
a
[
n
]
v
b
[
n
]
v
c
[
n
]
]
.
Based on (1)-(2), v
α
[
n
] and v
β
[
n
] can also be expressed as
(3)
v
α
[
n
]
=
s
α
[
n
]
+
q
α
[
n
]
,
v
β
[
n
]
=
s
β
[
n
]
+
q
β
[
n
]
,
where
(4)
s
α
[
n
]
=
(
A
cos
⁡
(
ϕ
)
-
B
sin
(
ϕ
)
)
cos
⁡
(
ω
n
)
-
(
A
sin
(
ϕ
)
+
B
cos
⁡
(
ϕ
)
)
sin
(
ω
n
)
,
s
β
[
n
]
=
(
-
B
cos
⁡
(
ϕ
)
+
C
sin
(
ϕ
)
)
cos
⁡
⁡
(
ω
n
)
+
(
B
sin
(
ϕ
)
+
C
cos
⁡
(
ϕ
)
)
sin
(
ω
n
)
,
with
(5)
A
=
6
12
(
4
V
a
+
V
b
+
V
c
)
,
B
=
2
4
(
V
b
-
V
c
)
,
C
=
6
4
(
V
b
+
V
c
)
,
q
α
[
n
]
=
2
3
(
ξ
a
[
n
]
-
1
2
ξ
b
[
n
]
-
1
2
ξ
c
[
n
]
)
,
q
β
[
n
]
=
2
2
(
ξ
b
[
n
]
-
ξ
c
[
n
]
)
.
Although both q
α
[
n
] and q
β
[
n
] contain ξ
b
[
n
] and ξ
c
[
n
], it is easy to show that the noise terms are uncorrelated; that is, E
{
q
α
[
n
]
q
β
[
n
]
}
=
0 where E denotes the expectation operator, and they have identical variance σ
2.Assuming that we haveN samples for each channel, (3) can be written in matrix form as follows:
(6)
v
=
s
+
q
,
where
(7)
v
=
[
v
α
[
1
]
v
α
[
2
]
⋯
v
α
[
N
]
v
β
[
1
]
v
β
[
2
]
⋯
v
β
[
N
]
]
T
,
s
=
[
s
α
[
1
]
s
α
[
2
]
⋯
s
α
[
N
]
s
β
[
1
]
s
β
[
2
]
⋯
s
β
[
N
]
]
T
,
s
=
G
x
,
G
=
[
H
0
0
H
]
,
H
=
[
cos
⁡
(
ω
)
sin
(
ω
)
cos
⁡
(
2
ω
)
sin
(
2
ω
)
⋮
⋮
cos
⁡
(
N
ω
)
sin
(
N
ω
)
]
,
x
=
[
A
cos
⁡
(
ϕ
)
-
B
sin
(
ϕ
)
-
A
sin
(
ϕ
)
-
B
cos
⁡
(
ϕ
)
-
B
cos
⁡
(
ϕ
)
+
C
sin
(
ϕ
)
B
sin
(
ϕ
)
+
C
cos
⁡
(
ϕ
)
]
,
q
=
[
q
α
[
1
]
q
α
[
2
]
⋯
q
α
[
N
]
q
β
[
1
]
q
β
[
2
]
⋯
q
β
[
N
]
]
T
.
Here, T denotes the transpose operator and0 is the N
×
2 zero matrix. We see that x corresponds to the linear unknowns, while ω is the nonlinear unknown in (6). Employing NLS [11, 12], the estimates of ω and x, denoted by ω
^ and x
^, are
(8)
{
ω
^
,
x
^
}
=
arg
min
⁡
ω
,
x
J
(
ω
,
x
)
,
J
(
ω
,
x
)
=
(
v
-
G
x
)
T
(
v
-
G
x
)
.Based on the Newton-Raphson procedure, the updating rule forω
^ is
(9)
ω
^
(
k
+
1
)
=
ω
^
(
k
)
-
∇
ω
J
(
ω
^
(
k
)
,
x
^
(
k
)
)
∇
ω
2
J
(
ω
^
(
k
)
,
x
^
(
k
)
)
,
where
(10)
∇
ω
J
(
ω
,
x
)
=
-
(
G
1
x
)
T
(
v
-
G
x
)
,
∇
ω
2
J
(
ω
,
x
)
=
x
T
(
G
1
T
G
1
)
x
+
(
G
2
x
)
T
(
v
-
G
x
)
,
with
(11)
G
1
=
[
H
1
0
0
H
1
]
,
H
1
=
[
-
sin
(
ω
)
cos
⁡
(
ω
)
-
2
sin
(
2
ω
)
2
cos
⁡
(
2
ω
)
⋮
⋮
-
N
sin
(
N
ω
)
N
cos
⁡
(
N
ω
)
]
,
G
2
=
[
H
2
0
0
H
2
]
,
H
2
=
[
cos
⁡
(
ω
)
sin
(
ω
)
4
cos
⁡
(
2
ω
)
4
sin
(
2
ω
)
⋮
⋮
N
2
cos
⁡
(
N
ω
)
N
2
sin
(
N
ω
)
]
.
Here, ω
^
(
k
) and x
^
(
k
) are the estimates of ω and x at the kth iteration. Once we have ω
^
(
k
), x
^
(
k
) is easily obtained from (8) as
(12)
x
^
(
k
)
=
(
G
T
G
)
-
1
G
T
v
|
ω
=
ω
^
(
k
)
,
where −1 denotes the matrix inverse. To start the algorithm of (9) and (12), we need ω
^
(
0
). Noting that ω should be around its nominal value ω
*
=
Ω
*
/
F
s, that is, ω
∈
[
ω
*
-
τ
,
ω
*
+
τ
] where τ is the maximum deviation from ω
*, ω
^
(
0
) is computed using grid search as follows. We assign K uniformly-spaced grid points in the range [
ω
*
-
τ
,
ω
*
+
τ
] where one of them is ω
^
(
0
). For each possible candidate ω
^
(
0
), we determine x
^
(
0
) according to (12). The pair {
ω
^
(
0
)
,
x
^
(
0
)
} which gives the minimum value of J
(
ω
,
x
) will be chosen as the initial guess for (9). In our study, the iterative algorithm is terminated when |
(
ω
^
(
k
+
1
)
-
ω
^
(
k
)
)
/
(
ω
^
(
k
+
1
)
)
|
<
ϵ, where ϵ
>
0 is a small tolerance constant, is reached. After obtaining ω
^, the NLS estimates of ϕ, V
a, V
b, and V
c are straightforwardly computed from x
^ as
(13)
ϕ
^
=
tan
⁡
-
1
(
x
^
3
-
x
^
2
x
^
1
+
x
^
4
)
,
V
^
a
=
(
6
2
x
^
1
-
6
6
x
^
3
)
cos
⁡
⁡
(
ϕ
^
)
-
(
6
2
x
^
2
+
6
6
x
^
3
)
sin
(
ϕ
^
)
,
V
^
b
=
(
x
^
4
-
2
x
^
2
)
cos
⁡
(
ϕ
^
)
+
(
6
3
x
^
3
-
2
x
^
1
)
sin
(
ϕ
^
)
,
V
^
c
=
(
x
^
4
+
2
x
^
2
)
cos
⁡
(
ϕ
^
)
+
(
6
3
x
^
3
+
2
x
^
1
)
sin
(
ϕ
^
)
,
where x
^
m (m
=
1,2
,
3,4) denotes the mth element of x
^.Finally, we examine the convergence of (9). According to [13], global convergence with quadratic rate is guaranteed when M
|
ω
-
ω
^
(
0
)
|
<
1 is satisfied, where M is as follows:
(14)
M
=
sup
⁡
ω
∈
{
-
π
,
π
}
|
∇
ω
3
J
(
ω
,
x
)
∇
ω
2
J
(
ω
,
x
)
|
.
To determine the value of M, we first relate G, G
1, and G
2 as follows:
(15)
G
1
T
G
1
+
G
2
T
G
=
c
E
,
G
1
=
N
1
G
F
,
G
2
=
N
2
G
,
where
(16)
E
=
[
1
0
0
0
1
0
0
2
0
]
,
F
=
[
0
-
1
0
0
1
0
0
0
0
0
0
-
1
0
0
1
0
]
,
c
=
N
(
N
+
1
)
(
2
N
+
1
)
6
,
N
1
=
diag
⁡
(
[
1
2
⋯
N
1
2
⋯
N
]
)
,
N
2
=
diag
⁡
(
[
1
2
2
⋯
N
2
1
2
2
⋯
N
2
]
)
.
Based on (10)-(11), ∇
ω
2
J
(
ω
,
x
) can be expressed as
(17)
∇
ω
2
J
(
ω
,
x
)
=
x
T
(
c
E
-
G
2
T
G
)
x
+
x
T
G
2
T
(
v
-
G
x
)
≥
x
T
G
T
N
2
v
-
2
x
T
G
T
N
2
G
x
≥
x
T
G
T
N
2
v
-
2
N
2
x
T
(
G
T
G
)
x
.
Equality holds if and only if x
1
2
+
x
2
2
+
2
x
2
x
3
=
0, ω
=
m
π with m
=
-
1
,
-
1
/
2
,
0
,
1
/
2
,
1 and
(18)
∇
ω
3
J
(
ω
,
x
)
=
x
T
G
1
T
N
2
v
=
x
T
F
T
G
T
N
3
v
,
where N
3
=
diag
⁡
(
[
1
2
3
⋯
N
3
1
2
3
⋯
N
3
]
). Substituting (12) into (17)-(18) yields
(19)
∇
ω
2
J
(
ω
,
x
)
≥
v
T
(
G
(
G
T
G
)
-
1
G
T
)
N
2
v
-
2
N
2
v
T
(
G
(
G
T
G
)
-
1
G
T
)
v
≥
v
T
W
1
v
,
∇
ω
3
J
(
ω
,
x
)
=
v
T
W
2
v
,
where
(20)
W
1
=
G
(
G
T
G
)
-
1
G
T
(
N
2
-
2
N
2
)
,
W
2
=
G
(
G
T
G
)
-
1
F
T
G
T
N
3
.
We can then write M as
(21)
M
=
sup
⁡
ω
∈
{
-
π
,
π
}
{
|
∇
ω
3
J
(
ω
,
x
)
v
T
v
|
|
1
(
∇
ω
2
J
(
ω
,
x
)
)
/
(
v
T
v
)
|
}
=
|
λ
2
λ
1
|
,
where λ
1 and λ
2 are the minimum eigenvalue of W
1 and the maximum eigenvalue of W
2, respectively. It is easily shown that λ
1
=
N
3 and λ
2
=
N
2. Hence M
=
N. As a result, if the initial estimate is chosen such that N
|
ω
-
ω
^
(
0
)
|
<
1 is satisfied, global solution will be obtained.


## 2.2. Cramér-Rao Lower Bound

Let the unknown parameter vector beθ
=
[
ω
ϕ
V
a
V
b
V
c
]
T. Then, the CRLB of θ
^ is obtained from the diagonal elements of the inverse of the Fisher information matrix, denoted by F
(
θ
) [14]. The (
i
,
j
) entry of F
(
θ
) is written as
(22)
[
F
(
θ
)
]
i
,
j
=
(
∂
s
∂
θ
i
)
T
C
-
1
(
∂
s
∂
θ
j
)
,
i
,
j
=
1
,
…
,
5
,
where(23)
∂
s
∂
θ
i
=
[
∂
s
α
[
1
]
∂
θ
i
∂
s
α
[
2
]
∂
θ
i
⋯
∂
s
α
[
N
]
∂
θ
i
∂
s
β
[
1
]
∂
θ
i
∂
s
β
[
2
]
∂
θ
i
⋯
∂
s
β
[
N
]
∂
θ
i
]
T
,with
(24)
∂
s
α
[
n
]
∂
ω
=
-
n
(
A
cos
⁡
(
ϕ
)
-
B
sin
(
ϕ
)
)
sin
(
ω
n
)
-
n
(
A
sin
(
ϕ
)
+
B
cos
⁡
(
ϕ
)
)
cos
⁡
(
ω
n
)
,
∂
s
β
[
n
]
∂
ω
=
n
(
B
cos
⁡
(
ϕ
)
-
C
sin
(
ϕ
)
)
sin
⁡
(
ω
n
)
+
n
(
B
sin
(
ϕ
)
+
C
cos
⁡
(
ϕ
)
)
cos
⁡
(
ω
n
)
,
∂
s
α
[
n
]
∂
ϕ
=
(
-
A
sin
(
ϕ
)
-
B
cos
⁡
(
ϕ
)
)
cos
⁡
⁡
(
ω
n
)
-
(
A
cos
⁡
(
ϕ
)
-
B
sin
(
ϕ
)
)
sin
(
ω
n
)
,
∂
s
β
[
n
]
∂
ϕ
=
(
B
sin
(
ϕ
)
+
C
cos
⁡
(
ϕ
)
)
cos
⁡
⁡
(
ω
n
)
+
(
B
cos
⁡
(
ϕ
)
-
C
sin
(
ϕ
)
)
sin
(
ω
n
)
,
∂
s
α
[
n
]
∂
V
a
=
6
3
(
cos
⁡
(
ϕ
)
cos
⁡
(
ω
n
)
-
sin
(
ϕ
)
sin
(
ω
n
)
)
,
∂
s
β
[
n
]
∂
V
a
=
0
,
∂
s
α
[
n
]
∂
V
b
=
(
6
12
cos
⁡
⁡
(
ϕ
)
-
2
4
sin
⁡
(
ϕ
)
)
cos
⁡
⁡
(
ω
n
)
-
(
6
12
sin
(
ϕ
)
+
2
4
cos
⁡
(
ϕ
)
)
sin
(
ω
n
)
,
∂
s
β
[
n
]
∂
V
b
=
(
-
2
4
cos
⁡
⁡
(
ϕ
)
+
6
4
sin
⁡
(
ϕ
)
)
cos
⁡
⁡
(
ω
n
)
+
(
2
4
sin
(
ϕ
)
+
6
4
cos
⁡
(
ϕ
)
)
sin
(
ω
n
)
,
∂
s
α
[
n
]
∂
V
c
=
(
6
12
cos
⁡
⁡
(
ϕ
)
+
2
4
sin
⁡
(
ϕ
)
)
cos
⁡
⁡
(
ω
n
)
-
(
6
12
sin
(
ϕ
)
-
2
4
cos
⁡
(
ϕ
)
)
sin
(
ω
n
)
,
∂
s
β
[
n
]
∂
V
c
=
(
2
4
cos
⁡
⁡
(
ϕ
)
+
6
4
sin
⁡
(
ϕ
)
)
cos
⁡
⁡
(
ω
n
)
+
(
-
2
4
sin
(
ϕ
)
+
6
4
cos
⁡
(
ϕ
)
)
sin
(
ω
n
)
and C is the covariance matrix of the noise term which is C
=
σ
2
I with I denoting the 2
N
×
2
N identity matrix. The (
i
,
j
) element can be simplified as
(25)
[
F
(
θ
)
]
i
,
j
=
1
σ
2
∑
n
=
1
N
‍
(
∂
s
α
[
n
]
∂
θ
i
∂
s
α
[
n
]
∂
θ
j
+
∂
s
β
[
n
]
∂
θ
i
∂
s
β
[
n
]
∂
θ
j
)
.


## 3. Simulation Results

To assess the proposed estimator for the unbalanced three-phase power system, computer simulations have been conducted. The MSEs,E
{
(
ω
^
-
ω
)
2
}, E
{
(
ϕ
^
-
ϕ
)
2
}, E
{
(
V
^
a
-
V
a
)
2
}, E
{
(
V
^
b
-
V
b
)
2
}, and E
{
(
V
^
c
-
V
c
)
2
}, and the mean frequency estimate, are employed as the performance measures. Comparisons with the CLMS and ACLMS algorithms as well as the CRLB are also made. We choose Ω
=
101
π rad/s and F
s
=
5000 Hz, and hence ω
=
0.0202
π with ω
*
=
0.02
π. The remaining parameters for (1) are assigned as ϕ
=
π
/
4, V
a
=
1, V
b
=
0.7, and V
c
=
0.6. The maximum frequency deviation is τ
=
π
/
2500 which corresponds to ±
2
% difference from the nominal value, while the number of grid points is chosen as K
=
30. When the condition M
<
1
/
τ is satisfied, global convergence is ensured. The tolerance parameter for the frequency estimate update is ϵ
=
1
0
-
8. All results are based on 1000 independent Monte Carlo runs.First, we study the MSEs forω, ϕ, V
a, V
b, and V
c versus σ
2 when the data length is assigned as N
=
40, and the results are plotted in Figure 1. It is seen that, when σ
2 is sufficiently small, the MSE performance aligns with the CRLB, indicating the optimality of NLS estimator. Figure 2 shows the MSEs versus N when the noise power is fixed at σ
2
=
-
10 dB. The high performance of the proposed scheme is again illustrated. Note that when N is small enough, the MSEs can be lower than the CRLB. It is because we have prior knowledge regarding the range of ω, but this information is not utilized in the CRLB derivation.Figure 1
Mean square errors versusσ
2.Figure 2
Mean square errors versusN.Next, we compare the performance of the proposed estimator with the CLMS and ACLMS methods. The proposed method is a batch-mode method, so we make it adaptive by utilizing a sliding window with a length ofL
=
60 on the observations. That means if we receive N samples, the frequency is estimated by the former (
N
-
L
+
1
) data. In this test, the noise power σ
2 is set to -
20 dB. The frequency and phase parameters are the same with the former test, while the voltages are shown in Figure 3. When N
<
250, V
a
=
V
b
=
V
c
=
1. We add 0.05 to V
a with 0.1 to V
b and V
c from N
=
250. Subsequently, V
c
=
0 after N
=
750. Figure 4 shows the mean frequency estimates under the time-varying case and it is seen that the proposed method is superior to the CLMS and ACLMS algorithms in both the performance of estimating frequency and the robustness to abruptly change of voltage.Figure 3
Three channel signals in unbalanced system.Figure 4
Mean estimated frequency versusN in unbalanced system.Figure5 shows the results when the observed data are contaminated by harmonics. We add a balanced 10
% third harmonic and a balanced 5
% fifth harmonic of the fundamental frequency ω to the system at N
=
250. It can be seen that although three methods give fluctuating performance, our estimate always oscillates around the true frequency value of ω
=
0.0202
π. Finally, Figure 6 addresses the impact of amplitude oscillation. In this test, the voltages are set as V
a
=
1
+
0.05
sin
(
2
π
n
/
F
s
), V
b
=
1
+
0.1
sin
(
2
π
n
/
F
s
) and V
c
=
1
+
0.15
sin
(
2
π
n
/
F
s
) at N
=
250. It is observed that our proposed method provides an estimate around the true value of ω
=
0.0202
π even when the amplitudes variation exist.Figure 5
Mean estimated frequency versusN under harmonic contamination.Figure 6
Mean estimated frequency versusN under amplitude variation.


## 4. Conclusion

An accurate estimator for the unbalanced three-phase power system in the presence of additive Gaussian noise has been developed. Theα
β-transformation is exploited to produce a pair of in-phase and quadrature signals from the three-phase waveforms, and then NLS cost function is constructed, where the frequency is the only nonlinear parameter. The Newton-Raphson scheme is employed to find NLS solution and its initialization and global convergence are studied. It is demonstrated that the MSE performance of the frequency, phase, and voltage estimates can achieve the CRLB and its mean frequency estimation accuracy is higher than that of the CLMS and ACLMS algorithms. A future work is to evaluate the developed algorithm using real three-phase power system measurements.



---

*Source: 539420-2014-08-04.xml*