# Model Predictive Control for Load Frequency Control with Wind Turbines

Authors: Yi Zhang; Xiangjie Liu; Yujia Yan
Journal: Journal of Control Science and Engineering
 (2015)

---

## Abstract

Reliable load frequency (LFC) control is crucial to the operation and design of modern electric power systems. Considering the LFC problem of a four-area interconnected power system with wind turbines, this paper presents a distributed model predictive control (DMPC) based on coordination scheme. The proposed algorithm solves a series of local optimization problems to minimize a performance objective for each control area. The scheme incorporates the two critical nonlinear constraints, for example, the generation rate constraint (GRC) and the valve limit, into convex optimization problems. Furthermore, the algorithm reduces the impact on the randomness and intermittence of wind turbine effectively. A performance comparison between the proposed controller with and that without the participation of the wind turbines is carried out. Good performance is obtained in the presence of power system nonlinearities due to the governors and turbines constraints and load change disturbances.


---

## Body

## 1. Introduction

Wind energy is considered as a promising and encoring renewable energy alternative for power generation owing to environmental and economical benefits. The world market of wind installation set a new record in the year of 2014 and reached a total size of 51 GW [1]. Nowadays, due to the interconnection of more distributed generators, especially wind turbines that are committed to grid operation, electric power system has become more complicated than ever.Power system LFC incorporating WTGs can be a quite challengeable issue. The output power of WTGs varies with wind speed fluctuation [2]. This wind power fluctuation imposed additional power imbalance to the power system and may cause frequency deviation from the nominal value [3]. Significant frequency deviations may lead to the disconnection of some loads and generations and even can lead to whole power system oscillations. Previous studies [4–6] provide extensive overviews of the primary and secondary frequency control strategies of power systems with wind power plants.LFC, secondary frequency control, has been performed by integrating the area control error (ACE), which acts on the load reference settings of the governors. LFC tasks are maintaining tie-line power flow and system frequency close to nominal value for the multiarea interconnected power system [7]. As a fundamental characteristic of electric power operations, frequency of the system deviates from its nominal value due to generation-demand imbalance. Conventional generators, in which the turbine rotational speed is nearly constant, provide inertia and governor response against frequency deviations; however, the speed of a wind turbine is not synchronous with the grid and is usually controlled to track the maximum power point. It implies that the wind turbines will have less time to react to the power imbalance, probably resulting in lager frequency deviations.Thus, it is thus necessary to establish the optimal profile of the WTGs power surge in coordination with the characteristics of conventional plants to achieve a more economical and reliable operation of power system. With the large amount of realistic constraints, for example, generation rate constraints (GRCs) in the conventional units, the pitch angle, and generator torque constraints in WTGs, the LFC becomes a large-scale, distributed, multiconstraints optimization problem.Recently, a few attempts studied the idea of wind turbines in the issue of LFC [8–10]. Two types of wind farm models are derived and demonstrated to portray the capability of set-point tracking under automatic generation control (AGC) [8]. This inference leads to the development of a simplified wind farm model that is specially designed for the set-point control in the power system study. However, the durability of inertia effect depends on the allowable rotor speed range. An adaptive fuzzy logic structure was used to propose a new LFC scheme in the interconnected large-scale power system in the presence of wind turbines [9]. The performance against sudden load change and wind power fluctuations in different wind power penetration rates is confirmed by simulation. A flatness-based method to control frequency and power flow for multiarea power system with wind turbine is presented in [10]. And, practical constraints such as generator ramping rates of wind turbine generator can be considered in designing the controllers. As abovementioned reference, the control schemes are designed for each area to maintain the frequency at nominal value and to keep power flows near scheduled values. However, local controller in each area does not work cooperatively towards satisfying systemwide control objectives. In addition, the control scheme [8–10] mentioned above could yield unsatisfactory performance since the effects of nonlinearities such as generation rate constraint and generation ramping rate were not considered.Model predictive control (MPC), also called receding horizon control, was originally developed to be an effective method for processing industrial control. It transforms the control problem into a finite horizon optimal control problem that can also satisfy multivariable constraints on the input, output, and state variables. In the power industry, MPC has been successfully used in controlling power plant steam-boiler generation processes [11–13]. In power system control, MPC was first developed to be an economic-oriented LFC [14], which generates the control action based on the open-loop optimization method over a finite horizon. MPC has subsequently been developed to realize the constrained optimal algorithm for LFC problem. In [15], the constraint handling ability of MPC is employed to effectively account for the generation rate constraints (GRCs) but without the analysis of closed-loop stability and robustness. Recently, MPC has been successfully used in LFC design of multiarea power system with wind turbines [16]. However, each area controller is designed independently and the communication between the local controllers is not considered. On the other hand, with the size and capacity of wind farms increasing in recent years, traditional centralized MPCs encounter many difficulties due to limitations in exchanging information with large-scale, geographically extensive control areas. In order to deal with these issues, advanced distributed control strategies have to be investigated and implemented.Developing decentralized/distributed LFC structures can be an effective way of solving this problem. The decentralized model predictive control scheme for the LFC of multiarea interconnected power system is presented in [17]. However, the local controller does not consider generation rate constraint that is only imposed on the turbine in the simulation. In the distributed MPC (DMPC), the benefits from using a decentralized structure are partially preserved, and the plant-wide performance and stability are improved through coordination [18, 19]. In [20], feasible cooperation-based MPC method is used in distributed LFC instead of centralized MPC. It is noted that the range of load change used in the cases is very large and inappropriate for the LFC issue.This paper studies the effect of merging the wind turbines on the system frequency of multiarea power system. The first control area includes an aggregated wind turbine model (which consists of 60 wind turbine units) beside the thermal power plant. According to the distributed LFC structure, the dynamics model of the four-area interconnected power system is established. In our scheme, the overall power system is decomposed into four areas and each area has its own local MPC controller. These areas-based MPCs exchange their measurement and predictions by communication and incorporate the information from other controllers into their local objective so as to coordinate with each other. The controllers calculate the optimal control signal while respecting constraints over the wind turbines output frequency deviation and the load change. Not only do the effects of the physical constraints conclude generating rate constraints (GRCs) and the limit of governor position in conventional power plant, but also the wind speed constraints in wind turbines are considered. Comparisons of response to step load change, computational burden, and robustness have been made between DMPC, centralized MPC, and decentralized MPC. The results confirm the superiority of the proposed DMPC technique.The remainder of the paper is organized as follows. Modeling of wind turbines participation in LFC is presented in Section2, and the proposed DMPC algorithm is presented in Section 3. Section 4 presents the application of the algorithm in a four-area interconnected power system. The conclusions are presented in Section 5.


## 2. Distributed Model of Hybrid Power System

Figure1 illustrates the interconnected power system consisting of four control areas connected by tie-lines, which consists of thermal power plant, variable speed wind turbines (VSWTs), and hydro power plant. In area 1, wind turbine is taken into consideration as it can provide a new solution to the contradiction between economic development and environment pollution. Area 4 is the thermal power plant, while area 2 and area 3 are hydro power plants.Figure 1
The four-area interconnected hybrid power system.Detailed compositions of each area are shown in Figures2–4. In addition, area 1 includes an aggregated wind turbine model which consists of 40 VSWT units while the capacity of thermal plant is 600 MW. The variables and parameters are listed in Table 1. In each control area, a change in local demand (load) alters the nominal frequency. The DMPC in each control area i manipulates the load reference set-point Δ
P
r
e
f
,
i to drive the frequency deviations Δ
f
i and tie-line power flow deviations Δ
P
t
i
e
,
i
j to zero.Table 1
Power system variables and parameter.
Parameter/variable
Description
Unit
ω
r
Angular velocity of rotor
rad/s
ω
g
Angular velocity of high speed shaft and generator
rad/s
T
g
Generator reaction torque
Nm
T
r
Aerodynamic torque
Nm
K
s
Total stiffness on low speed shaft
Nm/rad
J
r
Inertia of the rotor (low speed shaft and gearbox)
Kgm2
J
g
Inertia of the rotor (high speed shaft and gearbox)
Kgm2
n
gear
Exchange ratio
Null
η
gear
Efficiency of the gear box
%
τ
β
Actuator time constant
s
K
β
Actuator gain
Hz/p.u.MW
v
m
Wind speed
m/s
θ
ref
Pitch demand
rad
θ
Pitch angle
rad
P
e
ref
Power demand
p.u.MW
P
e
The output of wind turbine
p.u.MW
Δ
f
i
(
t
)
Frequency deviation
Hz
Δ
P
g
i
(
t
)
Generator output power deviation
p.u.MW
Δ
X
g
i
(
t
)
Governor valve position deviation
p.u.
Δ
X
g
h
i
(
t
)
Governor valve servomotor position deviation
p.u.
Δ
P
t
i
e
,
i
(
t
)
Tie-line active power deviation
p.u.MW
Δ
P
d
i
(
t
)
Load disturbance
p.u.MW
K
P
i
Power system gain
Hz/p.u.MW
K
r
i
Reheat turbine gain
Hz/p.u.MW
T
P
i
Power system time constant
s
T
r
i
Reheat turbine time constant
s
T
W
i
Water starting time
s
T
i
i, T
R
i
Hydro governor time constants
s
T
G
i
Thermal governor time constant
s
T
T
i
Turbine time constant
s
K
S
i
j
Interconnection gain between control areas
p.u.MW
K
B
i
Frequency bias factor
p.u.MW/Hz
R
i
Speed drop due to governor action
Hz/p.u.MW
ACE
i
Area control error
p.u.MWFigure 2
Block diagram of a thermal power plant and wind turbines (i
=
1).Figure 3
Block diagram of a hydro power plant (i
=
2,3).Figure 4
Block diagram of a thermal power plant (i
=
4).


### 2.1. Wind Turbine Model

A wind turbine is an installation for converting kinetic energy extracted from wind to electrical energy. Figure5 illustrates the basic model structure of a wind turbine and the interactions between the different dynamic components in the model. The whole wind turbine can be divided into four subsystems: aerodynamics subsystem, mechanical subsystem, electrical, and actuator subsystem [21].Figure 5
Diagram of a variable speed wind turbine.The linearization model for the variable speed wind turbine in Figure6 can be represented by(1a)
Δ
φ
˙
ε
=
Δ
ω
r
-
1
n
g
e
a
r
Δ
ω
g
,
(1b)
Δ
ω
˙
r
=
-
K
s
J
r
Δ
φ
ε
+
1
J
r
Δ
T
r
,
(1c)
Δ
ω
˙
g
=
η
g
e
a
r
K
s
n
g
e
a
r
Δ
φ
ε
-
1
J
r
Δ
T
g
,
(1d)
Δ
θ
˙
=
-
1
τ
β
Δ
θ
+
K
β
τ
β
Δ
θ
r
e
f
.Figure 6
Diagram of wind power plant in area 1.The generator reaction torqueT
g and the reference pitch angle θ
r
e
f are used as indicator of the input of VSWT, as u
e
=
Δ
β
r
e
f
Δ
T
g
T
∈
R
2. Moreover, η is the efficiency of the generator and ω
g and T
g are used as indicator of the output power as P
e
=
η
ω
g
T
g
∈
R
1, where ω
g is the angular velocity of generator shaft. A generalized representation of the state-space model of the variable speed turbine can be described as(2a)
x
˙
e
t
=
A
v
m
x
e
t
+
B
1
v
m
ω
t
+
B
2
u
e
t
,
(2b)
z
e
t
=
C
x
e
t
+
D
1
ω
t
+
D
2
u
e
twith(3)
A
v
m
=
0
1
-
1
n
g
e
a
r
0
-
K
s
J
r
1
J
r
∂
T
r
∂
ω
r
o
p
0
1
J
r
∂
T
r
∂
θ
o
p
η
g
e
a
r
K
s
n
g
e
a
r
J
g
0
0
0
0
0
0
-
1
τ
θ
,
B
1
v
m
=
0
1
J
r
∂
T
r
∂
v
o
p
0
0
,
B
2
=
0
0
0
0
0
-
1
J
g
K
β
τ
β
0
,
C
=
0
0
1
0
0
0
0
0
,
D
1
=
0
0
,
D
2
=
0
0
0
1
,
x
e
=
Δ
φ
ε
Δ
ω
r
Δ
ω
g
Δ
θ
T
,
u
e
=
Δ
θ
r
e
f
Δ
T
g
,
r
e
f
T
,
z
e
=
Δ
ω
g
Δ
T
g
T
,
y
e
=
P
e
=
η
ω
g
T
g
.


### 2.2. Four-Area Power System with Wind Turbine

Denoting that the control areai
(
i
=
1,2
,
3,4
) is to be interconnected with the control area j, j
≠
i, through a tie-line, a linear continuous time-varying model of control area i can be written as(4)
x
˙
i
=
A
i
i
x
i
+
B
i
i
u
i
+
F
i
i
d
i
+
∑
i
≠
j
A
i
j
x
j
+
B
i
j
u
j
+
F
i
j
d
j
,
y
i
=
C
i
i
x
i
,where x
i
∈
R
n, u
i
∈
R
m, d
i
∈
R
k, and y
i
∈
R
l are the state vector, the control signal vector, the disturbance vector, and the vector of output of control area i, respectively. x
j
∈
R
p, u
j
∈
R
q, and d
j
∈
R
s are the state vector, the control signal vector, and the disturbance vector of neighbor control area, respectively. Matrices A
i
i, B
i
i
C
i
i, and F
i
i represent appropriate system matrices of control area i, and A
i
j, B
i
j, and F
i
j represent the matrices of interaction variables between area i and area j. Tie-line power for area i is represented by(5)
Δ
P
t
i
e
,
i
=
∑
j
=
1
j
≠
i
4
Δ
P
t
i
e
i
j
=
∑
j
=
1
j
≠
i
4
K
s
i
j
Δ
f
i
-
Δ
f
j
,
Δ
P
t
i
e
i
j
=
-
Δ
P
t
i
e
j
i
.The state, disturbance, and output vectors for areai are defined by(6)
x
i
=
Δ
f
i
Δ
P
t
i
e
,
i
Δ
P
g
i
Δ
X
g
i
Δ
φ
ε
Δ
ω
r
Δ
ω
g
Δ
θ
T
i
=
1
,
x
i
=
Δ
f
i
Δ
P
t
i
e
,
i
Δ
P
g
i
Δ
X
g
i
Δ
X
g
h
i
T
i
=
2,3
,
x
i
=
Δ
f
i
Δ
P
t
i
e
,
i
Δ
P
g
i
Δ
X
g
i
Δ
P
r
i
t
T
i
=
4
,
d
i
=
Δ
P
d
i
i
=
1,2
,
3,4
,
u
1
=
Δ
P
c
1
Δ
θ
r
e
f
Δ
T
g
T
,
y
i
=
A
C
E
i
=
K
B
i
Δ
f
i
+
Δ
P
t
i
e
,
i
i
=
1,2
,
3,4
.The state, control, and disturbance matrices for area 1 are as follows:(7)
A
11
=
-
1
T
P
1
-
K
P
1
T
P
1
K
P
1
T
P
1
0
0
0
0
0
∑
j
K
s
i
j
0
0
0
0
0
0
0
0
0
-
1
T
T
1
0
1
T
T
1
0
0
0
1
T
G
1
0
0
-
1
T
G
1
0
0
0
0
0
0
0
0
0
1
-
1
n
g
e
a
r
0
0
0
0
0
-
K
s
J
r
1
J
r
∂
T
r
∂
ω
r
o
p
0
1
J
r
∂
T
r
∂
θ
o
p
0
0
0
0
η
g
e
a
r
K
s
n
g
e
a
r
J
g
0
0
0
0
0
0
0
0
0
0
-
1
τ
θ
,
B
11
=
0
0
0
1
T
G
1
0
0
0
0
0
0
0
0
0
0
0
K
β
τ
θ
0
0
0
0
0
0
-
1
J
g
0
T
,
F
11
=
-
K
P
1
T
P
1
0
0
0
0
0
0
0
0
0
0
0
0
1
J
r
∂
T
r
∂
v
m
o
p
0
0
T
,
C
11
=
K
b
1
1
0
0
0
0
0
0
.However for hydro plants in areas 2 and 3 they are as follows:(8)
A
22
=
A
33
=
-
1
T
P
i
-
K
P
i
T
P
i
K
P
i
T
P
i
0
0
∑
j
K
S
i
j
0
0
0
0
2
α
0
-
2
T
W
i
2
κ
2
β
-
α
0
0
-
1
T
2
i
-
β
-
1
T
1
i
R
i
0
0
0
-
1
T
1
i
,
B
22
=
B
33
=
0
0
-
2
R
i
α
R
i
α
1
T
1
i
T
,
C
22
=
C
33
=
K
B
i
1
0
0
0
,
F
22
=
F
33
=
-
K
p
i
T
p
i
0
0
0
0
T
,where α
=
T
R
i
/
T
1
i
T
2
i
R
i, β
=
T
R
i
-
T
1
i
/
T
1
i
T
2
i, and κ
=
T
2
i
+
T
W
i
/
T
2
i
T
W
i.However for thermal power plants in area 4 they are as follows:(9)
A
44
=
-
1
T
P
i
-
K
P
i
T
P
i
K
P
i
T
P
i
0
0
∑
j
K
S
i
j
0
0
0
0
0
0
-
1
T
T
i
0
1
T
T
i
-
1
T
G
i
R
i
0
0
-
1
T
G
i
0
-
K
r
i
T
G
i
R
i
0
0
1
T
r
i
-
K
r
i
T
G
i
-
1
T
r
i
,
B
44
=
0
0
0
1
T
G
i
0
T
,
C
44
=
K
B
i
1
0
0
0
,
F
44
=
-
K
p
i
T
p
i
0
0
0
0
T
.The interaction matrices between the four control areas are as follows:(10)
A
i
j
=
0
0
0
0
0
0
0
0
-
K
S
i
j
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
i
=
1
,
j
=
2,3
,
4
,
A
i
j
=
0
0
0
0
0
0
0
0
-
K
S
i
j
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
i
=
1
,
j
=
2,3
,
i
≠
j
,
A
i
j
=
0
0
0
0
0
-
K
S
i
j
0
0
0
0
0
0
0
0
0
0
0
0
0
0
i
=
j
=
2,3
,
4
,
i
≠
j
,
B
i
j
=
0
8
×
4
,
F
i
j
=
0
8
×
2
i
=
1
,
j
=
2,3
,
4
,
i
≠
j
,
B
i
j
=
0
5
×
1
,
F
i
j
=
0
5
×
1
i
=
j
=
2,3
,
4
,
i
≠
j
.The GRCs for the thermal plants are Δ
P
˙
g
i
≤
0.0017p.u.MW/s and the hydro units are Δ
P
˙
g
i
≤
0.045p.u.MW/s. In addition, the load disturbance is constrained to Δ
P
˙
d
i
≤
0.3.


## 2.1. Wind Turbine Model

A wind turbine is an installation for converting kinetic energy extracted from wind to electrical energy. Figure5 illustrates the basic model structure of a wind turbine and the interactions between the different dynamic components in the model. The whole wind turbine can be divided into four subsystems: aerodynamics subsystem, mechanical subsystem, electrical, and actuator subsystem [21].Figure 5
Diagram of a variable speed wind turbine.The linearization model for the variable speed wind turbine in Figure6 can be represented by(1a)
Δ
φ
˙
ε
=
Δ
ω
r
-
1
n
g
e
a
r
Δ
ω
g
,
(1b)
Δ
ω
˙
r
=
-
K
s
J
r
Δ
φ
ε
+
1
J
r
Δ
T
r
,
(1c)
Δ
ω
˙
g
=
η
g
e
a
r
K
s
n
g
e
a
r
Δ
φ
ε
-
1
J
r
Δ
T
g
,
(1d)
Δ
θ
˙
=
-
1
τ
β
Δ
θ
+
K
β
τ
β
Δ
θ
r
e
f
.Figure 6
Diagram of wind power plant in area 1.The generator reaction torqueT
g and the reference pitch angle θ
r
e
f are used as indicator of the input of VSWT, as u
e
=
Δ
β
r
e
f
Δ
T
g
T
∈
R
2. Moreover, η is the efficiency of the generator and ω
g and T
g are used as indicator of the output power as P
e
=
η
ω
g
T
g
∈
R
1, where ω
g is the angular velocity of generator shaft. A generalized representation of the state-space model of the variable speed turbine can be described as(2a)
x
˙
e
t
=
A
v
m
x
e
t
+
B
1
v
m
ω
t
+
B
2
u
e
t
,
(2b)
z
e
t
=
C
x
e
t
+
D
1
ω
t
+
D
2
u
e
twith(3)
A
v
m
=
0
1
-
1
n
g
e
a
r
0
-
K
s
J
r
1
J
r
∂
T
r
∂
ω
r
o
p
0
1
J
r
∂
T
r
∂
θ
o
p
η
g
e
a
r
K
s
n
g
e
a
r
J
g
0
0
0
0
0
0
-
1
τ
θ
,
B
1
v
m
=
0
1
J
r
∂
T
r
∂
v
o
p
0
0
,
B
2
=
0
0
0
0
0
-
1
J
g
K
β
τ
β
0
,
C
=
0
0
1
0
0
0
0
0
,
D
1
=
0
0
,
D
2
=
0
0
0
1
,
x
e
=
Δ
φ
ε
Δ
ω
r
Δ
ω
g
Δ
θ
T
,
u
e
=
Δ
θ
r
e
f
Δ
T
g
,
r
e
f
T
,
z
e
=
Δ
ω
g
Δ
T
g
T
,
y
e
=
P
e
=
η
ω
g
T
g
.


## 2.2. Four-Area Power System with Wind Turbine

Denoting that the control areai
(
i
=
1,2
,
3,4
) is to be interconnected with the control area j, j
≠
i, through a tie-line, a linear continuous time-varying model of control area i can be written as(4)
x
˙
i
=
A
i
i
x
i
+
B
i
i
u
i
+
F
i
i
d
i
+
∑
i
≠
j
A
i
j
x
j
+
B
i
j
u
j
+
F
i
j
d
j
,
y
i
=
C
i
i
x
i
,where x
i
∈
R
n, u
i
∈
R
m, d
i
∈
R
k, and y
i
∈
R
l are the state vector, the control signal vector, the disturbance vector, and the vector of output of control area i, respectively. x
j
∈
R
p, u
j
∈
R
q, and d
j
∈
R
s are the state vector, the control signal vector, and the disturbance vector of neighbor control area, respectively. Matrices A
i
i, B
i
i
C
i
i, and F
i
i represent appropriate system matrices of control area i, and A
i
j, B
i
j, and F
i
j represent the matrices of interaction variables between area i and area j. Tie-line power for area i is represented by(5)
Δ
P
t
i
e
,
i
=
∑
j
=
1
j
≠
i
4
Δ
P
t
i
e
i
j
=
∑
j
=
1
j
≠
i
4
K
s
i
j
Δ
f
i
-
Δ
f
j
,
Δ
P
t
i
e
i
j
=
-
Δ
P
t
i
e
j
i
.The state, disturbance, and output vectors for areai are defined by(6)
x
i
=
Δ
f
i
Δ
P
t
i
e
,
i
Δ
P
g
i
Δ
X
g
i
Δ
φ
ε
Δ
ω
r
Δ
ω
g
Δ
θ
T
i
=
1
,
x
i
=
Δ
f
i
Δ
P
t
i
e
,
i
Δ
P
g
i
Δ
X
g
i
Δ
X
g
h
i
T
i
=
2,3
,
x
i
=
Δ
f
i
Δ
P
t
i
e
,
i
Δ
P
g
i
Δ
X
g
i
Δ
P
r
i
t
T
i
=
4
,
d
i
=
Δ
P
d
i
i
=
1,2
,
3,4
,
u
1
=
Δ
P
c
1
Δ
θ
r
e
f
Δ
T
g
T
,
y
i
=
A
C
E
i
=
K
B
i
Δ
f
i
+
Δ
P
t
i
e
,
i
i
=
1,2
,
3,4
.The state, control, and disturbance matrices for area 1 are as follows:(7)
A
11
=
-
1
T
P
1
-
K
P
1
T
P
1
K
P
1
T
P
1
0
0
0
0
0
∑
j
K
s
i
j
0
0
0
0
0
0
0
0
0
-
1
T
T
1
0
1
T
T
1
0
0
0
1
T
G
1
0
0
-
1
T
G
1
0
0
0
0
0
0
0
0
0
1
-
1
n
g
e
a
r
0
0
0
0
0
-
K
s
J
r
1
J
r
∂
T
r
∂
ω
r
o
p
0
1
J
r
∂
T
r
∂
θ
o
p
0
0
0
0
η
g
e
a
r
K
s
n
g
e
a
r
J
g
0
0
0
0
0
0
0
0
0
0
-
1
τ
θ
,
B
11
=
0
0
0
1
T
G
1
0
0
0
0
0
0
0
0
0
0
0
K
β
τ
θ
0
0
0
0
0
0
-
1
J
g
0
T
,
F
11
=
-
K
P
1
T
P
1
0
0
0
0
0
0
0
0
0
0
0
0
1
J
r
∂
T
r
∂
v
m
o
p
0
0
T
,
C
11
=
K
b
1
1
0
0
0
0
0
0
.However for hydro plants in areas 2 and 3 they are as follows:(8)
A
22
=
A
33
=
-
1
T
P
i
-
K
P
i
T
P
i
K
P
i
T
P
i
0
0
∑
j
K
S
i
j
0
0
0
0
2
α
0
-
2
T
W
i
2
κ
2
β
-
α
0
0
-
1
T
2
i
-
β
-
1
T
1
i
R
i
0
0
0
-
1
T
1
i
,
B
22
=
B
33
=
0
0
-
2
R
i
α
R
i
α
1
T
1
i
T
,
C
22
=
C
33
=
K
B
i
1
0
0
0
,
F
22
=
F
33
=
-
K
p
i
T
p
i
0
0
0
0
T
,where α
=
T
R
i
/
T
1
i
T
2
i
R
i, β
=
T
R
i
-
T
1
i
/
T
1
i
T
2
i, and κ
=
T
2
i
+
T
W
i
/
T
2
i
T
W
i.However for thermal power plants in area 4 they are as follows:(9)
A
44
=
-
1
T
P
i
-
K
P
i
T
P
i
K
P
i
T
P
i
0
0
∑
j
K
S
i
j
0
0
0
0
0
0
-
1
T
T
i
0
1
T
T
i
-
1
T
G
i
R
i
0
0
-
1
T
G
i
0
-
K
r
i
T
G
i
R
i
0
0
1
T
r
i
-
K
r
i
T
G
i
-
1
T
r
i
,
B
44
=
0
0
0
1
T
G
i
0
T
,
C
44
=
K
B
i
1
0
0
0
,
F
44
=
-
K
p
i
T
p
i
0
0
0
0
T
.The interaction matrices between the four control areas are as follows:(10)
A
i
j
=
0
0
0
0
0
0
0
0
-
K
S
i
j
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
i
=
1
,
j
=
2,3
,
4
,
A
i
j
=
0
0
0
0
0
0
0
0
-
K
S
i
j
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
i
=
1
,
j
=
2,3
,
i
≠
j
,
A
i
j
=
0
0
0
0
0
-
K
S
i
j
0
0
0
0
0
0
0
0
0
0
0
0
0
0
i
=
j
=
2,3
,
4
,
i
≠
j
,
B
i
j
=
0
8
×
4
,
F
i
j
=
0
8
×
2
i
=
1
,
j
=
2,3
,
4
,
i
≠
j
,
B
i
j
=
0
5
×
1
,
F
i
j
=
0
5
×
1
i
=
j
=
2,3
,
4
,
i
≠
j
.The GRCs for the thermal plants are Δ
P
˙
g
i
≤
0.0017p.u.MW/s and the hydro units are Δ
P
˙
g
i
≤
0.045p.u.MW/s. In addition, the load disturbance is constrained to Δ
P
˙
d
i
≤
0.3.


## 3. Distributed Model Predictive Controller

### 3.1. Distributed Model Predictive Controller

The block diagram of the DMPC scheme for a four-area interconnected power system is illustrated in Figure7. Though there exists large amount of variables in the interconnected power system, the 30 state variables expressed in (1a), (1b), (1c), and (1d) concerning the frequency, the generator output power, the governor valve (servomotor) position, the tie-line active power, the wind power, and the 4 load disturbance Δ
P
d
i are crucial to LFC problem. They can be measured or estimated directly by the local controller. The DMPC in each area exchange control information through the power line communication, which is a sole networking technology with high reliability that can provide high speed communication to power grids applications [22].Figure 7
Block diagram of DMPC for power system with wind turbines.Distributed MPC. The partitioned discrete-time model for control area i of the continuous-time four-area interconnected power system ((1a), (1b), (1c), and (1d)) can be expressed as follows:(11)
x
i
k
+
1
=
A
-
i
i
x
i
k
+
B
-
i
i
u
i
k
+
F
-
i
i
d
i
k
+
∑
i
≠
j
A
-
i
j
x
j
k
+
B
-
i
j
u
j
k
+
F
-
i
j
d
j
k
,
y
i
k
=
C
-
i
i
x
i
k
,where A
-
i
i, B
-
i
i, C
-
i
i, F
-
i
i, A
-
i
j, B
-
i
j, and F
-
i
j represent the discrete new matrices obtained from original matrices in (4) based on the Zero-Order Hold (ZOH) method.Assume that the state variablesx
i
k and the disturbance D
i can be measured or estimated directly by the controller in area i at sampling time k. Optimizations and exchange of variables are termed iterate. The iteration number is denoted by p.For DMPC, the optimal state-input trajectory(
x
i
,
u
i
) for each area i
,
i
=
1,2
,
3,4 at iterate p is obtained as the solution to the optimization problem: (12)
m
i
n
u
i
k
+
n
∣
k
J
i
k
(13)
J
i
k
=
∑
n
=
0
N
x
i
T
k
+
n
∣
k
Q
i
x
i
k
+
n
∣
k
+
u
i
T
k
+
n
∣
k
R
i
u
i
k
+
n
∣
k
(14a)
Subject to
x
i
3
k
+
n
∣
k
2
≤
0.0017
,
i
=
1,4
(14b)
x
i
3
k
+
n
∣
k
2
≤
0.0045
,
i
=
2,3
(14c)
x
i
4
k
+
n
∣
k
2
≤
σ
i
,
i
=
1,2
,
3,4
.For notational convenience, we drop thek dependence of x
i
k, u
i
k, i
=
1,2
,
3,4. It is shown in [20] that each x
-
i can be expressed as(15)
x
-
i
=
E
-
i
i
u
-
i
+
f
-
i
i
x
i
k
+
β
-
i
i
d
i
k
+
∑
i
≠
j
E
-
i
j
u
-
j
+
g
-
i
j
x
-
j
+
f
-
i
j
x
j
k
+
β
-
i
j
d
j
kwith(16)
x
-
i
=
x
i
k
+
1
∣
k
T
x
i
k
+
2
∣
k
T
⋯
x
i
k
+
N
p
∣
k
T
T
,
u
-
i
=
u
i
k
∣
k
T
u
i
k
+
1
∣
k
T
⋯
u
i
k
+
N
c
-
1
∣
k
T
T
.LetN
c denote the control horizon and let N
p denote the predictive horizon. x
-
i is no more a vector but a matrix after iteration obtained from original equation (4). The matrices in (15) have detailed expressions as follows:(17)
E
-
i
i
=
B
-
i
i
0
⋯
0
A
-
i
i
B
-
i
i
B
-
i
i
⋯
0
⋮
⋮
⋮
A
-
i
i
N
-
1
B
-
i
i
A
-
i
i
N
-
2
⋯
0
,
E
-
i
j
=
B
-
i
j
0
⋯
0
A
-
i
i
B
-
i
j
B
-
i
j
⋯
0
⋮
⋮
⋮
A
-
i
i
N
-
1
B
-
i
j
A
-
i
i
N
-
2
⋯
0
,
f
-
i
i
=
A
-
i
i
A
-
i
i
A
-
i
i
⋮
A
-
i
i
N
-
1
A
-
i
i
,
f
-
i
j
=
A
-
i
j
A
-
i
i
A
-
i
j
⋮
A
-
i
i
N
-
1
A
-
i
j
,
β
-
i
i
=
F
-
i
i
A
-
i
i
F
-
i
i
⋮
A
-
i
i
N
-
1
F
-
i
i
,
β
-
i
j
=
F
-
i
j
A
-
i
i
F
-
i
j
⋮
A
-
i
i
N
-
1
F
-
i
j
,
g
-
i
j
=
0
0
0
⋯
0
A
-
i
j
0
0
⋯
0
A
-
i
i
A
-
i
j
A
-
i
j
0
⋯
0
⋮
⋮
⋮
⋯
⋮
A
-
i
i
N
-
2
A
-
i
j
A
-
i
i
N
-
3
A
-
i
j
⋯
A
-
i
j
0
,where E
-
i
i, f
-
i
i, β
-
i
i, E
-
i
j, f
-
i
j, β
-
i
j, and g
-
i
j are the new matrices obtained from A
-
i
i, B
-
i
i, C
-
i
i, F
-
i
i, A
-
i
j, B
-
i
j, and F
-
i
j after iteration.Combining the models in (15) gives the following system of equations:(18)
Λ
x
~
=
ε
u
~
+
μ
x
~
k
+
ϕ
d
kwith(19)
Λ
=
I
-
g
-
12
-
g
-
13
-
g
-
14
-
g
-
21
I
-
g
-
23
-
g
-
24
-
g
-
31
-
g
-
32
I
-
g
-
34
-
g
-
41
-
g
-
42
-
g
-
43
I
,
ε
=
E
-
11
E
-
12
E
-
13
E
-
14
E
-
21
E
-
22
E
-
23
E
-
24
E
-
31
E
-
32
E
-
33
E
-
34
E
-
41
E
-
42
E
-
43
E
-
44
,
μ
=
f
-
11
f
-
12
f
-
13
f
-
14
f
-
21
f
-
22
f
-
23
f
-
24
f
-
31
f
-
32
f
-
33
f
-
34
f
-
41
f
-
42
f
-
43
f
-
44
,
ϕ
=
β
-
11
β
-
12
β
-
13
β
-
14
β
-
21
β
-
22
β
-
23
β
-
24
β
-
31
β
-
32
β
-
33
β
-
34
β
-
41
β
-
42
β
-
43
β
-
44
,
x
~
=
x
-
1
x
-
2
x
-
3
x
-
4
T
,
u
~
=
u
-
1
u
-
2
u
-
3
u
-
4
T
.Since matrix Λ is invertible, we can write it as(20)
x
-
i
=
E
-
i
i
u
-
i
+
f
-
i
i
x
i
k
+
β
-
i
i
d
i
k
+
∑
i
≠
j
E
-
i
j
u
-
j
+
f
-
i
j
x
j
k
+
β
-
i
j
d
j
kin which(21)
E
-
i
j
=
Λ
-
1
ε
,
f
-
i
j
=
Λ
-
1
μ
,
β
-
i
j
=
Λ
-
1
ϕ
.To do so, we eliminate the unknown matrixx
-
j because we have knowledge of x
j
k since it is just a vector at time k.In the distributed MPC algorithm, for subsystemi, the control signal U
-
i is designed at each time interval k
≥
0. By solving the following optimization problem denoted by J
i, it is usually defined as(22)
J
i
=
min
u
-
i
⁡
1
2
u
-
i
T
Φ
i
u
-
i
T
+
γ
i
+
Γ
i
+
∑
i
≠
j
H
i
j
u
-
j
T
u
-
iin which(23)
Q
i
=
diag
⁡
ω
i
Q
i
,
…
,
ω
i
Q
i
︷
N
p
,
R
i
=
diag
⁡
ω
i
R
i
,
…
,
ω
i
R
i
︷
N
c
,
Φ
i
=
R
i
+
E
-
i
i
T
Q
i
E
-
i
i
+
∑
j
=
1
j
≠
i
4
E
-
j
i
T
Q
j
E
-
j
i
,
γ
i
=
E
-
i
i
T
Q
i
g
-
i
i
+
∑
j
=
1
j
≠
i
4
E
-
j
i
T
Q
j
g
-
j
i
,
g
i
i
=
f
i
i
x
i
k
+
∑
j
=
1
4
f
i
j
x
j
k
,
Γ
i
=
E
i
i
T
Q
i
ρ
i
+
∑
j
=
1
4
E
j
i
T
Q
j
ρ
j
,
ρ
i
=
β
i
i
d
i
k
+
∑
j
=
1
4
β
i
j
d
j
k
,
H
i
j
=
E
-
i
i
T
Q
i
E
-
i
j
+
∑
j
=
1
j
≠
i
4
E
-
j
i
T
Q
j
E
-
j
i
.At time intervalk, (22) is implemented based on the future states and manipulated variables. The first input in the optimal sequence is injected into the processes, and the procedure is repeated at subsequent time intervals.Q
i
≥
0
,
R
i
≥
0 are symmetric weighting matrices and ω
i
>
0
,
∑
i
=
1
4
ω
i
=
1.Defineη
i
=
γ
i
+
Γ
i
+
∑
j
≠
i
H
i
j
u
-
j.Then (22) is rewritten as(24)
J
i
=
min
u
-
i
⁡
1
2
u
-
i
T
Φ
i
u
i
T
+
η
i
T
u
-
i
.


### 3.2. Constraint Handling

The two crucial nonlinearities, for example, the GRCs and the valve position limits of the governor, have been considered as the state constraints in the designed DMPC, as shown in Figures8 and 9.Figure 8
Thermal power plant with GRC.Figure 9
Hydro power plant with GRC.In power system, the GRC can be expressed asΔ
P
˙
g
k
m
i
n
≤
Δ
P
˙
g
k
≤
Δ
P
˙
g
k
m
a
x, and then the constraints on Δ
P
g can be expressed as follows: (25)
T
Δ
P
˙
g
k
m
i
n
+
Δ
P
g
k
-
1
≤
Δ
P
g
k
≤
T
Δ
P
˙
g
k
m
a
x
+
Δ
P
g
k
-
1
,
(26)
Δ
P
-
g
=
Δ
P
g
k
+
1
∣
k
Δ
P
g
k
+
2
∣
k
⋯
Δ
P
g
k
+
N
p
∣
k
T
.Since Δ
P
g
i
=
X
i
3, the state form can be expressed as(27)
Δ
P
-
g
=
S
i
x
-
i
,where S
i
=
d
i
a
g
ω
i
S
i
i
,
…
,
ω
i
S
i
i
︷
N
p.Wheni
=
1,4, S
i
i
=
0
0
1
0
0, and when i
=
2,3, S
i
i
=
0
0
1
0
0, with (25) and (27), the constraints on Δ
P
g
k are expressed as N
-
i
≤
S
i
x
-
i
≤
M
-
i.Define(28)
N
-
i
=
N
i
N
i
⋯
N
i
︷
N
p
T
,
M
-
i
=
M
i
M
i
⋯
M
i
︷
N
p
T
,where N
i and M
i are obtained from (15).Consider the constraints onΔ
P
g
k:(29)
S
i
E
i
i
-
S
i
E
i
i
u
-
i
≤
M
-
i
-
S
i
f
i
i
x
i
k
+
β
i
i
d
i
k
+
∑
j
≠
i
E
i
j
u
-
j
+
f
i
j
x
j
k
+
β
i
j
d
j
k
-
N
-
i
+
S
i
f
i
i
x
i
k
+
β
i
i
d
i
k
+
∑
i
≠
j
E
i
j
u
-
j
+
f
i
j
x
j
k
+
β
i
j
d
j
k
.Define(30)
Ψ
i
=
S
i
E
i
i
-
S
i
E
i
i
,
Π
i
=
M
-
i
-
S
i
f
i
i
x
i
k
+
β
i
i
d
i
k
+
∑
j
≠
i
E
i
j
u
-
j
+
f
i
j
x
j
k
+
β
i
j
d
j
k
-
N
-
i
+
S
i
f
i
i
x
i
k
+
β
i
i
d
i
k
+
∑
i
≠
j
E
i
j
u
-
j
+
f
i
j
x
j
k
+
β
i
j
d
j
k
.Then, distributed MPC algorithm (24) for multiple-interconnected system can be transformed into the following optimization problem with GRC constraints:(31)
J
i
=
m
i
n
u
-
i
1
2
u
~
i
T
Φ
i
u
-
i
T
+
η
i
T
u
-
i
Subject to
Ψ
i
u
-
i
≤
Π
i
.


### 3.3. The DMPC Algorithm

Step 1 (initialization).
The constant matricesR
i, R
j and Q
i, Q
j, at control interval k
=
0, are given. Choose the specified error tolerance ε
i. Set iteration p
=
0.Step 2 (communication).
The controller in each subsystemi exchanges its previous predictions x
-
i
(
k
), x
-
j
k, set u
-
i
0
(
k
), and u
-
j
0
(
k
) at initial instant.Step 3 (optimization and iteration).
Whilep
<
p
m
a
x,
u
-
i
∗
(
p
) is solved by the optimal problem (31)
If u
i
p
-
u
i
p
-
1
≤
ε
i
∀
i
∈
1,2
,
3,4
Break
End if
Exchange the solutionsu
-
i
p and u
-
j
p and set p
=
p
+
1
If ε
i
=
0
∀
i
∈
1,2
,
3,4
Break
End if
End whileStep 4 (assignment and prediction).
Send outu
i
(
k
)
=
u
-
i
(
k
). Otherwise, u
i
(
k
)
=
u
-
i
(
k
-
1
). Predict the future states.Step 5 (implementation).
Setk
=
k
+
1, and repeat Step 1.


## 3.1. Distributed Model Predictive Controller

The block diagram of the DMPC scheme for a four-area interconnected power system is illustrated in Figure7. Though there exists large amount of variables in the interconnected power system, the 30 state variables expressed in (1a), (1b), (1c), and (1d) concerning the frequency, the generator output power, the governor valve (servomotor) position, the tie-line active power, the wind power, and the 4 load disturbance Δ
P
d
i are crucial to LFC problem. They can be measured or estimated directly by the local controller. The DMPC in each area exchange control information through the power line communication, which is a sole networking technology with high reliability that can provide high speed communication to power grids applications [22].Figure 7
Block diagram of DMPC for power system with wind turbines.Distributed MPC. The partitioned discrete-time model for control area i of the continuous-time four-area interconnected power system ((1a), (1b), (1c), and (1d)) can be expressed as follows:(11)
x
i
k
+
1
=
A
-
i
i
x
i
k
+
B
-
i
i
u
i
k
+
F
-
i
i
d
i
k
+
∑
i
≠
j
A
-
i
j
x
j
k
+
B
-
i
j
u
j
k
+
F
-
i
j
d
j
k
,
y
i
k
=
C
-
i
i
x
i
k
,where A
-
i
i, B
-
i
i, C
-
i
i, F
-
i
i, A
-
i
j, B
-
i
j, and F
-
i
j represent the discrete new matrices obtained from original matrices in (4) based on the Zero-Order Hold (ZOH) method.Assume that the state variablesx
i
k and the disturbance D
i can be measured or estimated directly by the controller in area i at sampling time k. Optimizations and exchange of variables are termed iterate. The iteration number is denoted by p.For DMPC, the optimal state-input trajectory(
x
i
,
u
i
) for each area i
,
i
=
1,2
,
3,4 at iterate p is obtained as the solution to the optimization problem: (12)
m
i
n
u
i
k
+
n
∣
k
J
i
k
(13)
J
i
k
=
∑
n
=
0
N
x
i
T
k
+
n
∣
k
Q
i
x
i
k
+
n
∣
k
+
u
i
T
k
+
n
∣
k
R
i
u
i
k
+
n
∣
k
(14a)
Subject to
x
i
3
k
+
n
∣
k
2
≤
0.0017
,
i
=
1,4
(14b)
x
i
3
k
+
n
∣
k
2
≤
0.0045
,
i
=
2,3
(14c)
x
i
4
k
+
n
∣
k
2
≤
σ
i
,
i
=
1,2
,
3,4
.For notational convenience, we drop thek dependence of x
i
k, u
i
k, i
=
1,2
,
3,4. It is shown in [20] that each x
-
i can be expressed as(15)
x
-
i
=
E
-
i
i
u
-
i
+
f
-
i
i
x
i
k
+
β
-
i
i
d
i
k
+
∑
i
≠
j
E
-
i
j
u
-
j
+
g
-
i
j
x
-
j
+
f
-
i
j
x
j
k
+
β
-
i
j
d
j
kwith(16)
x
-
i
=
x
i
k
+
1
∣
k
T
x
i
k
+
2
∣
k
T
⋯
x
i
k
+
N
p
∣
k
T
T
,
u
-
i
=
u
i
k
∣
k
T
u
i
k
+
1
∣
k
T
⋯
u
i
k
+
N
c
-
1
∣
k
T
T
.LetN
c denote the control horizon and let N
p denote the predictive horizon. x
-
i is no more a vector but a matrix after iteration obtained from original equation (4). The matrices in (15) have detailed expressions as follows:(17)
E
-
i
i
=
B
-
i
i
0
⋯
0
A
-
i
i
B
-
i
i
B
-
i
i
⋯
0
⋮
⋮
⋮
A
-
i
i
N
-
1
B
-
i
i
A
-
i
i
N
-
2
⋯
0
,
E
-
i
j
=
B
-
i
j
0
⋯
0
A
-
i
i
B
-
i
j
B
-
i
j
⋯
0
⋮
⋮
⋮
A
-
i
i
N
-
1
B
-
i
j
A
-
i
i
N
-
2
⋯
0
,
f
-
i
i
=
A
-
i
i
A
-
i
i
A
-
i
i
⋮
A
-
i
i
N
-
1
A
-
i
i
,
f
-
i
j
=
A
-
i
j
A
-
i
i
A
-
i
j
⋮
A
-
i
i
N
-
1
A
-
i
j
,
β
-
i
i
=
F
-
i
i
A
-
i
i
F
-
i
i
⋮
A
-
i
i
N
-
1
F
-
i
i
,
β
-
i
j
=
F
-
i
j
A
-
i
i
F
-
i
j
⋮
A
-
i
i
N
-
1
F
-
i
j
,
g
-
i
j
=
0
0
0
⋯
0
A
-
i
j
0
0
⋯
0
A
-
i
i
A
-
i
j
A
-
i
j
0
⋯
0
⋮
⋮
⋮
⋯
⋮
A
-
i
i
N
-
2
A
-
i
j
A
-
i
i
N
-
3
A
-
i
j
⋯
A
-
i
j
0
,where E
-
i
i, f
-
i
i, β
-
i
i, E
-
i
j, f
-
i
j, β
-
i
j, and g
-
i
j are the new matrices obtained from A
-
i
i, B
-
i
i, C
-
i
i, F
-
i
i, A
-
i
j, B
-
i
j, and F
-
i
j after iteration.Combining the models in (15) gives the following system of equations:(18)
Λ
x
~
=
ε
u
~
+
μ
x
~
k
+
ϕ
d
kwith(19)
Λ
=
I
-
g
-
12
-
g
-
13
-
g
-
14
-
g
-
21
I
-
g
-
23
-
g
-
24
-
g
-
31
-
g
-
32
I
-
g
-
34
-
g
-
41
-
g
-
42
-
g
-
43
I
,
ε
=
E
-
11
E
-
12
E
-
13
E
-
14
E
-
21
E
-
22
E
-
23
E
-
24
E
-
31
E
-
32
E
-
33
E
-
34
E
-
41
E
-
42
E
-
43
E
-
44
,
μ
=
f
-
11
f
-
12
f
-
13
f
-
14
f
-
21
f
-
22
f
-
23
f
-
24
f
-
31
f
-
32
f
-
33
f
-
34
f
-
41
f
-
42
f
-
43
f
-
44
,
ϕ
=
β
-
11
β
-
12
β
-
13
β
-
14
β
-
21
β
-
22
β
-
23
β
-
24
β
-
31
β
-
32
β
-
33
β
-
34
β
-
41
β
-
42
β
-
43
β
-
44
,
x
~
=
x
-
1
x
-
2
x
-
3
x
-
4
T
,
u
~
=
u
-
1
u
-
2
u
-
3
u
-
4
T
.Since matrix Λ is invertible, we can write it as(20)
x
-
i
=
E
-
i
i
u
-
i
+
f
-
i
i
x
i
k
+
β
-
i
i
d
i
k
+
∑
i
≠
j
E
-
i
j
u
-
j
+
f
-
i
j
x
j
k
+
β
-
i
j
d
j
kin which(21)
E
-
i
j
=
Λ
-
1
ε
,
f
-
i
j
=
Λ
-
1
μ
,
β
-
i
j
=
Λ
-
1
ϕ
.To do so, we eliminate the unknown matrixx
-
j because we have knowledge of x
j
k since it is just a vector at time k.In the distributed MPC algorithm, for subsystemi, the control signal U
-
i is designed at each time interval k
≥
0. By solving the following optimization problem denoted by J
i, it is usually defined as(22)
J
i
=
min
u
-
i
⁡
1
2
u
-
i
T
Φ
i
u
-
i
T
+
γ
i
+
Γ
i
+
∑
i
≠
j
H
i
j
u
-
j
T
u
-
iin which(23)
Q
i
=
diag
⁡
ω
i
Q
i
,
…
,
ω
i
Q
i
︷
N
p
,
R
i
=
diag
⁡
ω
i
R
i
,
…
,
ω
i
R
i
︷
N
c
,
Φ
i
=
R
i
+
E
-
i
i
T
Q
i
E
-
i
i
+
∑
j
=
1
j
≠
i
4
E
-
j
i
T
Q
j
E
-
j
i
,
γ
i
=
E
-
i
i
T
Q
i
g
-
i
i
+
∑
j
=
1
j
≠
i
4
E
-
j
i
T
Q
j
g
-
j
i
,
g
i
i
=
f
i
i
x
i
k
+
∑
j
=
1
4
f
i
j
x
j
k
,
Γ
i
=
E
i
i
T
Q
i
ρ
i
+
∑
j
=
1
4
E
j
i
T
Q
j
ρ
j
,
ρ
i
=
β
i
i
d
i
k
+
∑
j
=
1
4
β
i
j
d
j
k
,
H
i
j
=
E
-
i
i
T
Q
i
E
-
i
j
+
∑
j
=
1
j
≠
i
4
E
-
j
i
T
Q
j
E
-
j
i
.At time intervalk, (22) is implemented based on the future states and manipulated variables. The first input in the optimal sequence is injected into the processes, and the procedure is repeated at subsequent time intervals.Q
i
≥
0
,
R
i
≥
0 are symmetric weighting matrices and ω
i
>
0
,
∑
i
=
1
4
ω
i
=
1.Defineη
i
=
γ
i
+
Γ
i
+
∑
j
≠
i
H
i
j
u
-
j.Then (22) is rewritten as(24)
J
i
=
min
u
-
i
⁡
1
2
u
-
i
T
Φ
i
u
i
T
+
η
i
T
u
-
i
.


## 3.2. Constraint Handling

The two crucial nonlinearities, for example, the GRCs and the valve position limits of the governor, have been considered as the state constraints in the designed DMPC, as shown in Figures8 and 9.Figure 8
Thermal power plant with GRC.Figure 9
Hydro power plant with GRC.In power system, the GRC can be expressed asΔ
P
˙
g
k
m
i
n
≤
Δ
P
˙
g
k
≤
Δ
P
˙
g
k
m
a
x, and then the constraints on Δ
P
g can be expressed as follows: (25)
T
Δ
P
˙
g
k
m
i
n
+
Δ
P
g
k
-
1
≤
Δ
P
g
k
≤
T
Δ
P
˙
g
k
m
a
x
+
Δ
P
g
k
-
1
,
(26)
Δ
P
-
g
=
Δ
P
g
k
+
1
∣
k
Δ
P
g
k
+
2
∣
k
⋯
Δ
P
g
k
+
N
p
∣
k
T
.Since Δ
P
g
i
=
X
i
3, the state form can be expressed as(27)
Δ
P
-
g
=
S
i
x
-
i
,where S
i
=
d
i
a
g
ω
i
S
i
i
,
…
,
ω
i
S
i
i
︷
N
p.Wheni
=
1,4, S
i
i
=
0
0
1
0
0, and when i
=
2,3, S
i
i
=
0
0
1
0
0, with (25) and (27), the constraints on Δ
P
g
k are expressed as N
-
i
≤
S
i
x
-
i
≤
M
-
i.Define(28)
N
-
i
=
N
i
N
i
⋯
N
i
︷
N
p
T
,
M
-
i
=
M
i
M
i
⋯
M
i
︷
N
p
T
,where N
i and M
i are obtained from (15).Consider the constraints onΔ
P
g
k:(29)
S
i
E
i
i
-
S
i
E
i
i
u
-
i
≤
M
-
i
-
S
i
f
i
i
x
i
k
+
β
i
i
d
i
k
+
∑
j
≠
i
E
i
j
u
-
j
+
f
i
j
x
j
k
+
β
i
j
d
j
k
-
N
-
i
+
S
i
f
i
i
x
i
k
+
β
i
i
d
i
k
+
∑
i
≠
j
E
i
j
u
-
j
+
f
i
j
x
j
k
+
β
i
j
d
j
k
.Define(30)
Ψ
i
=
S
i
E
i
i
-
S
i
E
i
i
,
Π
i
=
M
-
i
-
S
i
f
i
i
x
i
k
+
β
i
i
d
i
k
+
∑
j
≠
i
E
i
j
u
-
j
+
f
i
j
x
j
k
+
β
i
j
d
j
k
-
N
-
i
+
S
i
f
i
i
x
i
k
+
β
i
i
d
i
k
+
∑
i
≠
j
E
i
j
u
-
j
+
f
i
j
x
j
k
+
β
i
j
d
j
k
.Then, distributed MPC algorithm (24) for multiple-interconnected system can be transformed into the following optimization problem with GRC constraints:(31)
J
i
=
m
i
n
u
-
i
1
2
u
~
i
T
Φ
i
u
-
i
T
+
η
i
T
u
-
i
Subject to
Ψ
i
u
-
i
≤
Π
i
.


## 3.3. The DMPC Algorithm

Step 1 (initialization).
The constant matricesR
i, R
j and Q
i, Q
j, at control interval k
=
0, are given. Choose the specified error tolerance ε
i. Set iteration p
=
0.Step 2 (communication).
The controller in each subsystemi exchanges its previous predictions x
-
i
(
k
), x
-
j
k, set u
-
i
0
(
k
), and u
-
j
0
(
k
) at initial instant.Step 3 (optimization and iteration).
Whilep
<
p
m
a
x,
u
-
i
∗
(
p
) is solved by the optimal problem (31)
If u
i
p
-
u
i
p
-
1
≤
ε
i
∀
i
∈
1,2
,
3,4
Break
End if
Exchange the solutionsu
-
i
p and u
-
j
p and set p
=
p
+
1
If ε
i
=
0
∀
i
∈
1,2
,
3,4
Break
End if
End whileStep 4 (assignment and prediction).
Send outu
i
(
k
)
=
u
-
i
(
k
). Otherwise, u
i
(
k
)
=
u
-
i
(
k
-
1
). Predict the future states.Step 5 (implementation).
Setk
=
k
+
1, and repeat Step 1.


## 4. Simulation Results

In this section, the four-area power system stability is analyzed, and the performances of the proposed DMPC have been tested in case of wind turbines participation at nominal parameters. The simulation of the proposed DMPC scheme is also verified by two cases. The performance and the implementation of the proposed DMPC are compared with other two types of typical LFC scheme.As comparison, we design the centralized MPC and decentralized MPC controller for four-area interconnected power system, respectively. The four-area interconnected power system can be described as(32)
x
k
+
1
=
A
x
k
+
B
u
k
+
F
d
k
,
y
k
+
1
=
C
x
k
,where(33)
A
=
A
11
A
12
A
13
A
14
A
21
A
22
A
23
A
24
A
31
A
32
A
33
A
34
A
41
A
42
A
43
A
44
,
B
=
B
11
B
12
B
13
B
14
B
21
B
22
B
23
B
24
B
31
B
32
B
33
B
34
B
41
B
42
B
43
B
44
,
C
=
C
11
0
0
0
0
C
22
0
0
0
0
C
33
0
0
0
0
C
44
,
F
=
F
11
0
0
0
0
F
22
0
0
0
0
F
33
0
0
0
0
F
44
,
x
=
x
1
T
x
2
T
x
3
T
x
4
T
T
,
u
=
u
1
T
u
2
T
u
3
T
u
4
T
T
,
y
=
y
1
T
y
2
T
y
3
T
y
4
T
T
,
d
=
d
1
T
d
2
T
d
3
T
d
4
T
Twith constraints (12), (13), (14a), (14b), and (14c) for each control area. In centralized MPC framework, the MPC for overall system (32) solves the following optimization problem:(34)
min
u
k
+
n
∣
k
⁡
J
k
(35)
J
k
=
∑
n
=
0
N
x
T
k
+
n
∣
k
Q
x
k
+
n
∣
k
+
u
T
k
+
n
∣
k
R
u
k
+
n
∣
ksubject to (14a), (14b), and (14c).The weighting matricesQ and R in objective function (35) are chosen as R
=
d
i
a
g
1,1
,
1,1 and(36)
Q
=
d
i
a
g
1000,0
,
0,1000,1000,0
,
0,1000,1000
,
0,0
,
1000,1000,0
,
0,1000
.In the decentralized modeling framework, it is assumed that the interaction between the control areas is negligible. Subsequently, the decentralized model for each control area is(37)
x
i
k
+
1
=
A
i
i
x
i
k
+
B
i
i
u
i
k
+
F
i
i
d
i
k
,
y
i
k
+
1
=
C
i
i
x
i
kwith the system matrices and constraints (12), (13), (14a), (14b), and (14c) for each control area denoted as in Section 2. In decentralized MPC framework, each control area based MPC solves the following optimization problem:(38)
min
u
i
k
+
n
∣
k
⁡
J
i
k
(39)
J
i
k
=
∑
n
=
0
N
x
i
T
k
+
n
∣
k
Q
i
x
i
k
+
n
∣
k
+
u
i
T
k
+
n
∣
k
R
i
u
i
k
+
n
∣
ksubject to (14a), (14b), and (14c).The weighting matricesQ
i and R
i in objective function (39) are chosen as R
1
=
R
2
=
R
3
=
R
4
=
1 and(40)
Q
1
=
Q
2
=
Q
3
=
Q
4
=
diag
⁡
1000,0
,
0,1000
.Choose the prediction horizon of the centralized MPC, decentralized MPC, and RDMPC to beN
=
15, choose the control horizon to be N
c
=
10, and choose the sample time T
s
=
0.1 and λ
=
0.1. Consider GRC for the thermal power plants in area 1 and area 4 to be Δ
P
˙
g
i
≤
r
=
0.1
p
.
u
.
M
W
/
m
i
n
=
0.0017
p
.
u
.
M
W
/
s and GRC for the hydro power plants in area 2 and area 3 to be Δ
P
˙
g
i
≤
r
=
2.7
p
.
u
.
M
W
/
m
i
n
=
0.045
p
.
u
.
M
W
/
s. In addition, area 1 includes an aggregated wind turbine model which consists of 30 wind turbine units of 2 MW rated VSWTs while the capacity of thermal plant is 600 MW. The wind turbine parameters and operating points [23] are indicated as follows:Operating point: 80 MW;wind speed: 12 m/s.
T
g
=
3781.9 Nm; ω
g
=
105 rad/s; ω
r
=
2.6869 rad/s.
K
s
=
7.87
e
6 Nm/rad; n
g
e
a
r
=
1
:
28.7; η
g
e
a
r
=
97.5
%.
J
r
=
2.8675 kgm2; J
g
=
54.5432 kgm2.
R
3
=
3.3Hz/p.u.MW; R
4
=
3Hz/p.u.MW.The parameters for the thermal and hydro plants used in the simulation are listed as follows:K
P
1
=
120Hz/p.u.MW; K
P
2
=
115Hz/p.u.MW.
K
P
3
=
80Hz/p.u.MW; K
P
4
=
75Hz/p.u.MW.
T
P
1
=
20
s; T
P
2
=
20
s; T
P
3
=
13
s; T
P
4
=
15
s.
R
1
=
2.4Hz/p.u.MW; R
2
=
2.5Hz/p.u.MW.
R
3
=
3.3Hz/p.u.MW; R
4
=
3Hz/p.u.MW.
K
B
1
=
0.425p.u.MW/Hz; K
B
2
=
0.409p.u.MW/Hz.
K
B
3
=
0.316p.u.MW/Hz; K
B
4
=
0.347p.u.MW/Hz.
T
G
1
=
0.08
s; T
G
2
=
0.1
s; T
G
3
=
0.08
s; T
G
4
=
0.2
s.
T
T
1
=
T
T
4
=
0.3
s; T
r
1
=
T
r
4
=
10
s; T
R
2
=
0.6
s.
T
R
3
=
0.513
s; T
22
=
5
s; T
23
=
10
s; T
W
2
=
1
s; T
W
3
=
2
s.
K
S
12
=
-
K
S
21
=
0.545p.u.MW.
K
S
23
=
-
K
S
32
=
0.444p.u.MW.
K
S
13
=
-
K
S
31
=
0.545p.u.MW.
K
S
14
=
-
K
S
41
=
0.5p.u.MW.
K
S
24
=
-
K
S
42
=
0.545p.u.MW.
K
S
34
=
-
K
S
43
=
0.545p.u.MW.Case 1 (response to step load change without wind turbines participation).
Wind turbine is present but it does not provide any power support in the event of grid frequency deviation. An event is simulated in which a system shown in Figure1 is subjected to step load disturbances as give in (41) at t
=
10
s. Consider(41)
Δ
P
d
1
=
Δ
P
d
2
=
Δ
P
d
3
=
Δ
P
d
4
=
0.1
.
Figure10 shows the simulation results of distributed MPC, centralized MPC, and decentralized MPC without wind turbine participation and only conventional integrator systems. The relative performance of distributed MPC, centralized MPC, and decentralized MPC rejecting the load disturbance in each area in Figure 10 is denoted by solid, dotted, and dashed lines, respectively. It has been noticed that the closed-loop trajectory of distributed MPC obtained by algorithm is little fast and almost indistinguishable from the closed-loop trajectory of centralized MPC. It successfully improves the dynamic response of area frequencies compared with decentralized MPC.
The control costs defined by [20] for different strategies are listed in Table 2. It is obviously seen that the DMPC controller needs nearly as much CPU time as decentralized MPC controller and significantly less CPU time than centralized MPC controllers. The proposed DMPC algorithm has significant computational advantages when compared to centralized MPC while achieving the best performance.Table 2
Cost of the different strategies.
Strategy
Cost [20]
Centralized MPC
0.10
Decentralized MPC
0.083
Distributed MPC
0.078Figure 10
Response of frequency deviation to step load disturbance in Case1: distributed MPC (solid line), centralized MPC (dotted line), and decentralized MPC (dashed line).Case 2 (response to step load change with wind turbines participation).
Wind turbine is present and it will provide active power support in the event of grid frequency deviation. An event is simulated in which a system shown in Figure1 is subjected to step load disturbances as give in (41) at t
=
10
s. Mean wind speed is assumed to be 17 m/s in area 1.
In Figures11 and 12, the behavior for the frequency is presented for Case 2 where the wind turbines are participating in load frequency control. The results from top to the bottom in Figure 11 are the frequency deviations for area 1 to area 4 and in Figure 12 are six tie-lines power change. In simulation, it is obvious that both the DMPC and the centralized MPC converge rapidly and drive the local frequency changes and tie-line power deviation to zero. The wind turbines that have participated in the interconnected power system do not affect the performance of the power system under distributed MPC and centralized MPC while satisfying all the physical constraints, for example, the GRC, the limit of the governors, and load step change constraints. However, with decentralized MPC, the rapid convergence cannot be guaranteed in the presence of wind turbines in area 1. This confirms the performance advantage of the proposed distributed model predictive control algorithm.
Figure13 shows the dynamic response of active power deviation Δ
P
e and rotor speed ω
g of wind turbine while participating in the load frequency control. When the control is activated, the frequency deviation becomes zero which consequently eliminated the additional active power deviation Δ
P
e and wind turbine is driven to operate again at the optimal rotor speed ω
g. It may be noted here that an increase in power step on top of the converter further reduces the rotor speed, thereby transferring more kinetic power to reduce the frequency dip. As shown in this figure, the distributed MPC in the presence of wind turbine has desirable performance in comparison to centralized MPC and decentralized MPC.
The distributed MPC control actions as shown in Figure14, Δ
θ
r
e
f, Δ
P
c
i, and Δ
T
g in four areas are depicted as solid, dotted, and dashed line, respectively. Δ
θ
r
e
f and Δ
T
g are the control signals of wind turbine in area 1, and Δ
P
c
i is the control signal of traditional power plants in the four areas. Figure 15 shows the generating outputs of traditional plants.Figure 11
Response of frequency deviation to step load disturbance in Case2: distributed MPC (solid line), centralized MPC (dotted line), and decentralized MPC (dashed line).Figure 12
Response of tie-line active power deviation in Case2: distributed MPC (solid line), centralized MPC (dotted line), and decentralized MPC (dashed line).Figure 13
Wind turbine response of electrical power and rotor speed in Case2: distributed MPC (solid line), centralized MPC (dotted line), and decentralized MPC (dashed line).Figure 14
Control signal of distributed MPC in Case2: Δ
θ
r
e
f in area 1 (solid line), Δ
P
c
i in four areas (dotted line), and Δ
T
g in area 1 (dashed line).Figure 15
Response of generated power deviation in Case2: distributed MPC (solid line), centralized MPC (dotted line), and decentralized MPC (dashed line).


## 5. Conclusions

In this paper, a DMPC scheme is presented for the LFC of a four-area interconnected power system with wind turbines. The state and input constraints including the valve position limit on the governor and the GRCs were incorporated into the system design. In our scheme, each control area has a local MPC controller, in which the four controllers coordinated with each other by exchanging their information. Comparisons of response to step load change and computational burden have been made between DMPC, centralized MPC, and decentralized MPC. The simulation results verified the reliability of the DMPC for achieving a performance that has advantages over the centralized MPC and distributed MPC in the presence of load changes. Moreover, the proposed DMPC scheme can guarantee a good performance under the wind turbines participation in LFC. Future work will be the extension of the proposed DMPC to different renewable energy contained LFC, since the greater utilization of intermittent renewable resources will induce greater power flow fluctuations.



---

*Source: 282740-2015-11-08.xml*