IEEE TRANSACTIONS ON SMART GRID
Abstract—The applicability of fractional order (FO) automatic
generation control (AGC) for power system frequency oscillation
damping is investigated in this paper, employing distributed
energy generation. The hybrid power system employs various
autonomous generation systems like wind turbine, solar
photovoltaic, diesel engine, fuel-cell and aqua electrolyzer along
with other energy storage devices like the battery and flywheel.
The controller is placed in a remote location while receiving and
sending signals over an unreliable communication network with
stochastic delay. The controller parameters are tuned using
robust optimization techniques employing different variants of
Particle Swarm Optimization (PSO) and are compared with the
corresponding optimal solutions. An archival based strategy is
used for reducing the number of function evaluations for the
robust optimization methods. The solutions obtained through the
robust optimization are able to handle higher variation in the
controller gains and orders without significant decrease in the
system performance. This is desirable from the FO controller
implementation point of view, as the design is able to
accommodate variations in the system parameter which may
result due to the approximation of FO operators, using different
realization methods and order of accuracy. Also a comparison is
made between the FO and the integer order (IO) controllers to
highlight the merits and demerits of each scheme.
Index Terms—distributed energy system; fractional order PID
controller; robust optimization; automatic generation control
I. INTRODUCTION
UE to the deregulation in the energy markets, the
environmental emission concerns and the rising costs of
electricity transmission and distribution, there is an increasing
trend of shifting from centralized power generation and
distribution to a more decentralized mode [1]. This has given
rise to distributed energy resources (DERs) with integration of
renewable energy technologies like wind and solar, along with
energy storage devices like flywheels, batteries etc. and
combined heat and power generation technologies [1]. Control
and communication play an important role in the efficient
operation of these distributed power systems [2]. To
effectively meet the challenges of control in DERs, the paper
looks at a novel controller design strategy for a hybrid power
system [3], where the sensor measurements and the control
Manuscript received July 09, 2014, revised May 11, 2015, Accepted.
I. Pan is with the Department of Earth Science and Engineering, Imperial
College London, Exhibition Road, SW7 2AZ, UK (email:
i.pan11@imperial.ac.uk).
S. Das is with the School of Electronics and Computer Science, University
of Southampton, Southampton SO17 1BJ, UK (e-mail:
s.das@soton.ac.uk).
signals are sent over an unreliable communication network
introducing random time delays in the control loop [4].
Fractional calculus [5] is a 300 year old mathematical
concept. However, in the last decade, it has found applications
in control systems and is gaining increasing interest from the
research community in other domains as well. Fractional
calculus has also found recent applications in computational
intelligence based control system design [6] with random time
delays e.g. in process control [7], nuclear reactor control [8]
etc. among others. Inspired by the successful applications in
these domains, the present paper investigates the applicability
of fractional order intelligent control for hybrid power system
or distributed energy generation. Other control approaches of
similar kind of system design include the standard PID
controller [9], robust H
∞
controller [10] etc.
The stochastic nature of the demand load and the renewable
generation terms i.e. solar and wind energies introduces
fluctuations in the system frequency [11-12]. The controller in
the hybrid power system tries to minimize the aberrations in
the system frequency so that the power quality is maintained.
This leads to the concept of AGC for grid frequency
oscillation damping in the context of distributed energy
generation [13]. This is done by sending an appropriate
control signal to the energy storage systems to absorb (release)
the surplus (deficit) power from (to) the grid. The controllers
are generally tuned in an output optimal fashion by
minimizing some error criterion [9]. However, these obtained
values of controller parameters may not be robust to the slight
variability during actual hardware implementation. This is
more relevant for fractional order controllers, since they are
realized in hardware using band limited approximations of
higher order transfer functions using different techniques like
Carlson, Matsuda, Continued Fraction Expansions (CFE) etc.
[5, 6, 14]. Under such circumstances, the optimal response as
obtained in the time domain simulations would vary
significantly. Therefore a robust optimization based controller
parameter tuning scheme is proposed in the present paper to
overcome this issue and facilitate practical implementation.
The controllers tuned with the robust algorithms show slight
variations in time domain performance, as opposed to the
drastic deterioration of performance with the optimally tuned
controllers, when the controller parameters are perturbed. In
the present work, different variants of PSO are used for robust
optimization. PSO has been used in other smart grid
applications like demand response and resource scheduling
[15], sizing of distributed generation and storage capacity
Fractional Order AGC for Distributed Energy
Resources Using Robust Optimization
Indranil Pan and Saptarshi Das
D
IEEE TRANSACTIONS ON SMART GRID
[16], wind power control [17], bi-directional energy trading
[18] etc. The present paper introduces a FOPID based
centralized AGC scheme for grid frequency oscillation
damping in a DER system. We use the PSO based robust
optimization technique for tuning the controller and report the
achievable parametric robustness of the hybrid power system.
II. DESCRIPTION OF THE DISTRIBUTED ENERGY SYSTEM
The schematic representation of the hybrid power system
using different energy generation/storage is illustrated in Fig.
1 with its various components described in Table I.
A. Models of Different Generation Subsystems
TABLE I
NOMINAL PARAMETERS OF THE COMPONENTS OF HYBRID POWER SYSTEM
Component Gain (K) Time constant (T)
Wind turbine generator (WTG) K
WTG
= 1 T
WTG
= 1.5
Aqua Electrolyzer (AE) K
AE
= 0.002 T
AE
= 0.5
Fuel Cell (FC) K
FC
= 0.01 T
FC
= 4
Flywheel energy storage system (FESS) K
FESS
= -0.01 T
FESS
= 0.1
Battery energy storage system (BESS) K
BESS
= -0.003 T
BESS
= 0.1
Diesel engine generator (DEG) K
DEG
= 0.003 T
DEG
= 2
Solar Photovoltaic (PV) K
PV
= 1 T
PV
= 1.8
For small signal analysis, the dynamics of the WTG, PV, FC
and DEG can be modeled by first order transfer functions (1)-
(4) with the associated gain and time constants given in Table
I [9, 19], where k represents the number of units. These
transfer functions represent the electrical power produced
(P
WTG
, P
PV
) from the renewable energy sources like wind
power (P
w
), solar irradiation (Φ) etc. In this study, a
centralized controller has been used for the hybrid energy
system (Fig. 1) as opposed to multiple decentralized
controllers [9] for each individual sub-systems like the battery,
flywheel and diesel. This helps in easier maintenance, reduced
wiring and also makes the design problem tractable by
reducing the number of controller parameters. However, there
would obviously be some deterioration in the performance, as
the same control signal is being used for all the sub-systems.
Nevertheless, here we show that the centralized controller
results in acceptable time domain performance. Different rate
limiters in each subsystem have been provided so that the
control signal is appropriately modified with respect to the
individual electromechanical characteristics of the
storage/generating devices.
1,1,2,3
k
WTG WTG WTG WTG W
GsK sT P Pk
(1)
1
PV PV PV PV
GsK Ts P
(2)
1,1,2
kk
FC FC FC FC AE
GsK sT P Pk
(3)
1
DEG DEG DEG DEG
GsK sT P u
(4)
B. Model of Aqua Electrolyzer
The aqua-electrolyzer produces hydrogen for the fuel cell
by using a part of the power generated from the renewable
sources like wind and/or solar. The dynamics of the AE can be
represented by the transfer function (5) [19] and it uses
1
n
K
fraction of the total power of WTG and PV to
produce hydrogen which is again used by two FCs to produce
power as an additional source to the grid.
11
AE AE AE AE WTG PV n
GsK sT P P P K
(5)
where,
,0.6
ntWTGPV n
KPP PK
.
Fig. 1 Schematic of the hybrid power system used in this study.
C. Models of Different Energy Storage Systems
In the hybrid energy system of Fig. 1, the FESS and the
BESS are connected in the feedback loop and are actuated by
the signal from the FOPID or PI
λ
D
μ
controller. These absorb
or release energy from or to the grid if there is a surplus or
deficit amount of power respectively. Their corresponding
dynamical models can be represented by (6)-(7) [19].
1
FESS FESS FESS FESS
GsK sT P u
(6)
1
BESS BESS BESS BESS
GsK sT P u
(7)
Here, the incremental control action of the FOPID controller is
represented as
CA
ut ut
, τ
CA
~ U(0.05, 0.15) [7]
that actuates the energy storage/generation devices. The
controller output
u
gets corrupted by the stochastic network
delay (
CA
) for the transmission of the signal from the
controller to actuator. Also, the models of the grid frequency
dependent energy storage/generating elements are considered
to have rate constraint nonlinearities as
0.9
FESS
P
,
0.2
BESS
P
,
0.01
DEG
P
respectively (Fig. 1),
such that all the three components operate in the nonlinear
zone for a wide range of controller values. The rate constraint
nonlinearities take care of the various electromechanical
constraints that these devices exhibit. However, the overall
system dynamics is dictated by the relative values of the gains
and time constants of the different components. For example,
the DEG has a larger time constant as compared to the FESS.
Therefore, as soon as a control signal is applied, the FESS will
respond more quickly and the DEG would take more time to
respond. Therefore, the overall system dynamics would be
governed by a combination of these fast and slow dynamics.
IEEE TRANSACTIONS ON SMART GRID
D. Power System Model Using Grid Frequency Deviation
The power system model can be represented as (8) which
describes the dynamics of power deficit/surplus (ΔP
e
) to the
grid frequency oscillation (Δf).
1
sys e L S
Gs fP fPP DMs
(8)
where,
M
and
D
are the equivalent inertia constant and
damping constant of the hybrid power/energy system [3] and
their typical values are considered as 0.4 and 0.03 respectively
for the present simulation study.
E. Wind Speed Model
The wind speed model should be able to capture the spatial
dependencies of wind flow like base fluctuation and small
stochastic components etc. To achieve this objective, the two
component wind model [19] is chosen as represented by (9).
WWBWN
VV V
(9)
where,
,
WB WN
VV
are the base wind component and noise wind
component respectively. The base wind component is the
constant component which is always present during the
operation of the wind turbine and has been considered as (10).
7.5 3 200 10.5 250
WB B
V K Ht Ht Ht
(10)
where,
B
K
is a constant (within a constant speed operation
regime) and
Ht
represents the Heaviside step function.
The noise component of the wind is expressed as (11).
2
1
2cos
N
WN V i i i
i
VS t
(11)
where,
12
i
i
, φ
i
~ U(0, 2π),
is the frequency
step to compute spectral density,
2
200
is variance of the
noise component and the spectral density function
Vi
S
is
given by (12).
43
2
22
21
Vi N i i
SKF F
(12)
Here,
0.004
N
K
is the surface drag coefficient,
2000F
is
the turbulence scale and
7.5
is the mean wind speed at
reference height. Here,
50N
and
0.5
rad/s are taken to
achieve an effective modeling accuracy.
F. Characteristics of Wind Turbine Model Output
The non-dimensional curves of the power coefficient
p
C
expressed as a function of tip speed ratio
(14) and blade
pitch angle
0.1745
is used to characterize the wind turbine
and is expressed as (13) [19].
3
0.44 0.0167 sin 0.0184 3
15 0.3
p
C
(13)
Here,
refers to the ratio of the speed at the blade tip of the
wind turbine to the wind speed and is expressed as (14).
blade blade W
RV
(14)
where,
23.5m
blade
R
is the radius of the wind turbine blades
and
3.14 rad/s
blade
is the rotational speed of the blades.
The output mechanical power of the wind turbine is given by
3
12
WrpW
P
AC V
(15)
where,
3
1.25 kg/m
refers to the air density and
2
1735 m
r
A
is the swept area of the blades.
G. Characteristics of PV Output Power and Demand Power
The power output of the photovoltaic system can be
represented by (16) as also done in [19].
10.005 25
PV a
PS T
(16)
where,
10%
is the conversion efficiency of the PV cells,
2
4084 mS
is the measured area of the PV array,
(17) in
kW/m
2
is the solar radiation on the surface of the PV cells and
o
25 C
a
T
is the ambient temperature.
0.5 0.3 25 0.3 75
0.3 150 , 0.1,0 1~.
nn
Ht Ht Ht
Ht t t U
(17)
Fig. 2 depicts one realization of the stochastic generation
components (
P
V
P
,
WTG
P
) along with the stochastic variations
in the load demand
L
P
(18) and the net power generated by the
renewable sources to the grid (
t
P
).
0.4 50 0.1 100 0.2 175
0.2 225 , 0.05,0.05~
L
LL
PHt Ht Ht Ht
Ht N N U
(18)
In all the three cases, there are sudden fluctuations in the
power levels with stochastic aberrations throughout in
WTG
P
,
P
V
P
and
L
P
which is representative of a realistic scenario.
Fig. 2 A single realization of the renewable generations and demand powers
which are independent of the controller structure.
H. Control Over Unreliable Communication Network
Since the different energy resources are located at different
places, they are assumed to communicate via a shared
communication medium [20, 21]. The use of a shared medium
introduces random delays in the control loop between the grid
frequency sensor to the controller (
SC
) and the controller to
the actuator (
CA
) [22]. A small amount of random delay can
IEEE TRANSACTIONS ON SMART GRID
induce instability while the control loop can be stable with a
larger amount of lumped static delay [6, 7]. Therefore these
stochastic delays must be considered in the optimization based
controller design procedure itself [8]. In the present
simulation,
SC
and
CA
are assumed to be randomly drawn
from a uniform distribution within a range of
0.05,0.15
. The
incremental control signal before and after the network are
represented as
u
and
u
. The whole hybrid power system
model in Fig. 2 thus can be viewed as a higher order stochastic
delay differential equation due to stochastic forcing and
stochastic delay terms and is numerically integrated with the
3
rd
order accurate Bogacki-Shampine formula over a total time
window of 300 sec with a fixed step-size of 0.01 sec.
III. BASICS OF FRACTIONAL ORDER CONTROLLER
Fractional calculus is an extension of the n
th
order
successive differ-integration of a function
f
t
having the
order as any real value
n
. For control system studies, the
non-integer order integro-differential operator is defined as
Caputo derivative (19) among three main definitions of
fractional calculus which under zero initial condition in
Laplace transform produce FO transfer functions [5].
0
1
0
1
,
,,1
m
t
t
m
Dft
Dft d
m
t
mZm m
(19)
The transfer function representation of a FOPID controller
is given in (20).
pi d
Cs K K s Ks
(20)
This typical controller structure has five independent
tuning knobs i.e. the three controller gains
,,
p
id
KKK
and
two FO integro-differential operators
,
.
For
1
and
1
the controller structure (20) reduces to the
classical PID controller in parallel structure [14, 6].
Various continuous and discrete time rational
approximation methods can be adopted to implement the FO
operators. In this paper, each guess value of the FO differ-
integrals
,
within the optimization process is continuously
rationalized with the Oustaloup’s 5
th
order rational
approximation (ORA) [14] within the chosen frequency range
of
22
10 ,10
rad/sec. Due to the fact that FO differ-
integrals represent infinite dimensional linear filters, their
band-limited realizations are necessary for implementation.
Here, each FO element has been rationalized using ORA given
by the equations (21) and (22). If it be assumed that the
expected fitting range or frequency range of controller
operation is
,
bh
, then the higher order filter which
approximates the FO element
s
can be written as (21).
()
N
k
f
kN
k
s
Gs s K
s
(21)
The poles, zeros and gain of the filter can be evaluated as (22).
0.5(1 ) 0.5(1 )
21 21
,,
kn kn
nn
kbhb kbhb h
K
(22)
In equations (21)-(22),
is the order of the differ-integration
and
21n
is the order of the realized analog filter. The
controller operates on the randomly delayed grid frequency
deviation signal to produce the control action (23).
, 0.05, 0.15~
pi d SCSC
ut K KD KD f t
(23)
IV. OPTIMIZATION ALGORITHMS AND CONTROL OBJECTIVES
A. The Concept of Robust Optimization
The difference between a robust and an optimal solution is
illustrated in Fig. 3. The optimum point has the lowest value
of the objective function. However, if the input design variable
has a certain variance, indicated by the first probability
distribution curve on the abscissa, then there is a
corresponding large deviation in the objective function value,
indicated by the probability distribution curve on the ordinate.
In case of the robust solution, the same variance in the input
variable produces a smaller variance in the objective function
value. Hence the latter solution is less sensitive to variation in
system parameters and is consequently a robust solution. The
robust solution has a higher value of objective function than
the optimal solution, but the worst case scenario for the robust
solution is much less severe than that of the corresponding
optimum solution [23].
There are various methods of assessing the robustness of
solutions (
x
) for the objective function
J
x
. The expected
fitness measure is used in this paper and is given by (24).
exp
.
J
x J x pdf d
(24)
where,
is the input variable fluctuations,
pdf
is the
probability distribution function of the occurrence of
over
the whole input variable space
,
N
, and
N
is the problem
dimension i.e. the number of decision variables.
Fig. 3 Schematic showing the difference between an optimum solution and a
robust solution.
B. Objective Function for Optimization Based Control
For effective functioning of the hybrid power system, the
controller gains and FO integro-differential orders need to be
tuned. For the controller design problem, the objective
function in (25) is considered. It consists of the integrals of
two weighted terms, which try to minimize the frequency
IEEE TRANSACTIONS ON SMART GRID
deviation in the hybrid power system (
f
), as well as the
incremental control signal (
u
) [7-8, 24].
max
22
0
1
T
n
J
wf wK u dt
(25)
where,
w
dictates the relative importance of the two objectives
(i.e. Integral of Squared Error – ISE and Integral of Squared
Deviation of Control Output – ISDCO) and
0.5w
is
considered to give equal weightage on both the parts of the
control objective.
5
10
n
K
is the normalizing constant to scale
ISE and ISDCO in an uniform scale and
max
300 sT
.
The objective function (25) is formulated in such a way that
along with the frequency oscillation (
f
), the incremental
control signals going to different actuators are also minimum.
This helps to limit the requirement of increased capacity for
the battery, reduces flywheel jerk and diesel consumption,
making the overall hybrid power system more cost-effective.
Also, due to the noisy nature of the frequency deviation
signal, the control signal gets amplified with the derivative
action of the PID/FOPID and even more due to the difference
operator (Δu) of the control signal. Therefore, there is a need
of a scale factor to bring the amplitude of the ISE and ISDCO
to a comparable platform and then only appropriate weights
could be assigned by a designer. An equal weightage is given
in the present case, as also reported in [7-8, 24]. In most of the
controller design problems, it is difficult to decide the weighs
a priori and a multi-objective optimization formalism should
be used to obtain the Pareto optimal trade-offs for different
weights as shown in [25]. The designer can then choose the
suitable weighting according to his requirements.
Nevertheless, the proposed methodology is still valid if the
weightings in the objective function is changed and our
simulations show one of such possible alternatives.
C.
Implementation Issues of Fractional Order Controller and
Need for Robust Optimization
Depending on the different methods of realization i.e.
analog (like Crone, Carlson, Matsuda, CFE with high and low
frequency approximations etc.) and digital (Tustin, Simpson,
backward difference, impulse response etc.) methods of a FO
element, the time/frequency domain characteristics of the filter
may be different [14]. Fig. 4 shows the phase response of the
band limited realization of
1
s
in the frequency range of
22
10 10
Hz. The corresponding time domain impulse
responses are shown in Fig. 5. It can be seen that there exists
significant differences in the time domain response among the
different realizations and also among different approximation
orders of the same realization. From the application point of
view, one can opt for any one of these realization techniques
and the filter order to implement a single FO operator without
paying much attention to the resulting frequency and time
domain discrepancies of the approximation. Therefore, the
design of the controller parameters themselves should be
robust enough to tolerate these imprecisions during the actual
hardware implementation, while still ensuring satisfactory
time domain performance. In addition, the components of the
hybrid energy system are generally modelled as low-pass first
order transfer functions, considering small signal stability
analysis [19], whereas in reality they may show more complex
nonlinear behavior which can be considered as a linear model
with uncertain parameters. The concept of robust optimization
in the present problem is introduced to handle both system
nonlinearities and FO controller implementation issues.
Fig. 4 Phase responses of different band-limited realizations of FO element
1
s
with different methods and order of approximation.
Fig. 5 Impulse responses of FO element
1
s
with different methods/orders.
D. Canonical PSO (CPSO) Optimizer and Its Variants
The CPSO algorithm tries to optimize an objective function
J
x
with respect to the design variable
n
x
as in (26).
n
x
minimize
J
x
(26)
where, the objective function
:
n
f
and the n-
dimensional search space
n
G
is pre-specified by the user.
The PSO algorithm consists of a swarm of particles
1, 2,...,
ip
x
in
with the maximum number of particles
p
n
specified by the user. The particles
i
x
search for an optimal
solution
n
x
of (26). The position of the
th
i
particle is
denoted by
,1 ,2 ,
:,,...,
T
n
iii in
xxx x
and the velocity is
denoted by
,1 , 2 ,
:,,...,
T
n
iii in
vvv v
, where
1, 2,...,
p
in
.
The position and velocity of the
th
i
particle,
n
i
x
is updated
in each iteration, based on equations (27)-(28) for
kZ
which indicates the iteration number.
11kkk
iii
x
xv
(27)