Citation
Optimal power flow in microgrid

Material Information

Title:
Optimal power flow in microgrid by Mirsadraddin Darakhshan Rokhsari
Creator:
Darakhshan Rokhsari, Mirsadraddin
Place of Publication:
Denver, CO
Publisher:
University of Colorado Denver
Publication Date:
Language:
English
Physical Description:
1 electronic file. : ;

Subjects

Subjects / Keywords:
Electric power production ( lcsh )
Distributed generation of electric power ( lcsh )
Distributed generation of electric power ( fast )
Electric power production ( fast )
Genre:
non-fiction ( marcgt )

Notes

Review:
This thesis is about the optimal power flow in microgrid with using multi-objective function. The thesis puposed new algorithm for optimal power flow plus new desicion about the determination of droop characteristic at optimized point.
Thesis:
Thesis (M.S.)--University of Colorado Denver. Electrical engineering
General Note:
Department of Electrical Engineering

Record Information

Source Institution:
|University of Colorado Denver
Holding Location:
|Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
857503788 ( OCLC )
ocn857503788

Downloads

This item has the following downloads:


Full Text
OPTIMAL POWER FLOW IN MICROGRID
by
Mirsadraddin Darakhshan Rokhsari
B.S.E.E University of Mashad, Iran, 1989
A thesis submitted to the
University of Colorado Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Applied Electrical
2012


This thesis for the Master of Science
degree by
Mirsadraddin Darakhshan Rokhsari
has been approved
by
Fernando Mancilla-David
Titsa Papantoni
Dan Connors
May fO, 20f2


Darakhshan Rokhsari, Mirsadraddin (Electrical Engineering)
Optimal Power Flow in Microgrid
Thesis directed by Professor Fernando Mancilla-David
ABSTRACT
A microgrid will play an essential role in the future, since it will reduce
the stress by the increasing electricity demands, as compared to the constraints
imposed by power plants and transmission lines. Hence, identifying the optimal
operating point of any distributed generators (DGs) in a microgrid leads them to
identifying the optimal operation of any generator in therm of cost minimization
efficiency maximization.
This thesis proposes a new algorithm for determination of the optimal oper-
ating point of distributed generators. In this algorithm, the new droop charac-
teristic slopes (frequency vs. active power and voltage vs. reactive power) of the
inverter will be detected of the optimized point. Conventionally, inverter-based
DGs droop characteristic slopes have a constant slope in their maximum and
minimum capacity in the different operation mode of microgrid (connect to or
disconnect from the grid). If any change is necessary, the slope will be chosen
arbitrary. To minimize the deviation of droop characteristic slopes and the other
objectives (defined by the decision maker), the multi-objective optimization al-
gorithm can be used. The multi-objective optimization is used in this paper,


can provide the reasonable decision to determine the new droop characteristic
slopes in optimized point with minimum deviation.


This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Approved: Fernando Mancilla-David


DEDICATION
This thesis is dedicated to my wife and son who have supported me all the
way since the beginning of my studies. They provide me with a great source
of motivation and inspiration. Also, this thesis is dedicated to my mother and
memories of my father for instilling in me the importance of hard work and
higher education. Finally, this thesis is dedicated to all those who believe in the
richness of learning.


ACKNOWLEDGMENT
I would like to express my deepest gratitude to my advisor, Professor Fernando
Mancilla-David, for his excellent guidance, caring, patience, and providing me
with an excellent atmosphere for doing research.
I greatly appreciate Professors Titsa Papantoni, and Professor Dan Connors
for forming part of my dissertation defense committee. Finally, I would like to
thank my wife and son; they were always there cheering me up and stood by me
through the good and bad times.
I would also like to thank my parents. They were always supporting me and
encouraging me with their best wishes.
Finally, I would like to thank my sister in law, Soraya for her supporting.


CONTENTS
Figures ................................................................ x
Tables.................................................................. xi
Chapter
1. Introduction......................................................... 1
2. Microgrid ........................................................... 3
2.1 Distributed Energy Resources unit (DERs) .......................... 3
2.2 Microgrid definition and characteristics........................... 3
2.2.1 Power management strategies for microgrid ........................ 6
3. Optimal Power Flow (OPF)............................................. 8
3.1 Introduction...................................................... 8
3.2 Formulation........................................................ 8
4. Droop Control........................................................ 13
4.1 Concept of droop control.......................................... 13
4.1.1 Droop in microgrid................................................ 15
4.2 Load sharing through P-f control.................................. 17
5. OPF with Droop Control............................................... 20
5.1 Mathematic model................................................. 20
5.2 Mathematical formula for optimization............................. 25
6. Case of Study........................................................ 30
6.1 Introduction...................................................... 30
viii


6.2 Microgrid connected to the grid.................................... 34
6.3 Microgrid disconnected from the grid (autonomous mode)............. 35
6.3.1 Sharing the power from droop characteristics .................. 35
6.3.2 Sharing the power with OPF .................................... 35
7. Conclusion........................................................... 40
Appendix
A. fmincon in MATLAB................................................. 41
B. The code of software in MATLAB...................................... 43
B.l Reading the hie.................................................... 43
B.2 Calculate Y-Bus ................................................... 51
B.3 Computational code ................................................ 52
B.4 Creating fun function.............................................. 58
B.5 creatin noncolon function....................................... 59
B.6 writing the hie after optimization................................. 62
References.............................................................. 66
IX


FIGURES
Figure
2.1 Microgrid structure ....................................................... 4
4.1 Power flowing through a line between two sources.......................... 13
4.2 Frequency droop characteristic ........................................... 16
4.3 Voltage droop characteristic.............................................. 17
4.4 Frequency droop characteristic for three DGs ......................... 18
5.1 Frequency droop characteristic for DG in bus i ....................... 21
5.2 Voltage droop characteristic for DG in bus i.............................. 22
5.3 Frequency droop characteristics for DG in bus i in different load . 25
5.4 Voltage droop characteristics for DG in bus i in different load ... 26
5.5 Direction of frequency droop characteristics changing for DG in bus
i for different load...................................................... 26
5.6 Direction of voltage droop characteristics changing for DG in bus i
for different load..................................................... 27
6.1 The case of study in microgrid ........................................... 30
6.2 Case of study with all buses.............................................. 31
6.3 Droop characteristic deviation for DG 1 38
6.4 Droop characteristic deviation for DG2 38
6.5 Droop characteristic deviation for DG3 ................................... 39
x


TABLES
Table
6.1 Amount of impedances between busi and busj in pu.................. 32
6.2 Data for types of bus ............................................... 33
6.3 Data for all of the generators.................................. 33
6.4 Data for all generators and buses in the connected mode with the grid 34
6.5 Data for microgrid in autonomous mode before optimization .... 36
6.6 The normalized of objective functions in autonomous mode after op-
timization
37


1. Introduction
The demand for electricity has increased steadily over the last few years due
to rapidly population growth. This trend will likely continue over the next few
decades [18]. Power generation in large centralized power plants such as nuclear,
fossil fuel (coal, gas), etc. has had an impact on the induced cost. However, the
customers frequently reside too far from the generating point of electricity. Elec-
tric power systems and the power grid are currently undergoing restructuring
due to a number of factors. These factors include: increasing demand, increasing
uncertainty caused by the integration of intermittent renewable energy sources
and further deregulation of the industry [5], [9].
Due to many economical, practical and technical reasons, the power grid can
not provide the necessary energy for the consumer. In providing a contingency
to the power shortage in many areas of demand, it has become necessary to
search for another way to generate the power required by the consumer. Any
new decision to produce energy should be researched in parallel with present
grid supply systems. In the past, all of the power plants have transferred energy
through transmission lines. However, restrictions in transmission lines force
the necessities of network design which allows generators to be placed close to
their loads. In conjunction with the advances in power electronics, the small-
scaled distributed power generation (non 50-60 HZ) is in the position to change
the conventional power system [11], Such change includes the deployment of
renewable energy sources (RES) such as wind, solar, fuel cells, etc.
1


The energy that is generated from natural resources such as wind, sunlight,
rain, tides and geothermal heat has been named renewable energy (RE), where
RE unit is an environmentally friendly way of generating power and a serious
contender to substituting of non renewable energies like coal and oil. Renewable
energy can be used for a wide range of environments as the resources of such
environments dictate. Such environments include villages as well as large rural
areas. Green power refers to a subset of renewable energy that characterized by
high environmental benefits.
A large number of research studies have been conducted in many countries
for substituting conventional energy with RES. For example, It has been re-
ported that from 1995 to 2000, the U.S. government had spent 1.456 billion
dollars on renewable energy research [19], while Japan is the rising star of green
power research. By adjusting their government policies and strong state finan-
cial support, Japan has made full use of renewable energy. It has also made
great achievements in wind, solar and ocean energy power generation [20]. Fur-
thermore, the utilization of green power has disseminated rapidly in Germany
over the past few years, where it is projected that it will gradually replace the
nuclear power generation in the next few years. [19].
2


2. Microgrid
The reduction of fossil fuel sources, pollution of the environment, constraints
in transmission lines, and the poor efficiency of energy drives us to a new type of
power generation. These kinds of generators should be located at the distribu-
tion voltage level and near residential districts. The types of generators in the
distribution network generate power by using non-conventional energy sources
(RES). This kind of generator is called a distributed generator (DG) and its
energy source is called distributed energy resources (DERs).
2.1 Distributed Energy Resources unit (DERs)
Distributed generators are small-range power generators that are connected
to the distribution network, or directly at the customer application [2], The
DERs refer to all technologies which can be used to provide energy close to the
consumer [8]. The high penetration of DGs in the network is able to alleviate the
pressure on the electric system by feeding into part of the local load. DER units
contain distributed generator (DG) and distributed storage (DS) with different
rates of flow [13].
2.2 Microgrid definition and characteristics
In 2001, Professor Robert H.Lasseter from the University of Wisconsin-
Madison proposed the definition of the microgrid [14] and created a microgrid
test bed. By this definition, the microgrid has a single unit role in a power
3


system. It is defined as the microgrid acting as a group of DER units in the
same way that distributed voltage can supply electricity and/or heat a group of
domestic load demand [18]. The proximity of the generator and the load (heating
and electrical) in the microgrid decrease the capital cost and reduce the losses in
the transmission lines. A traditional structure of microgrid is shown in Fig. 2.1.
The microgrid in Fig. 2.1 can operate in two modes:
Grid
Load
Distributed storage
Wind power
Load
Figure 2.1: Microgrid structure
1. In this mode, the microgrid is connected to the grid where the grid and
DER units share the generation of power sufficient for consumers in the
microgrid section.
2. In the second mode the microgrid is disconnected from the grid and acts as
4


an autonomous island in the microgrid. In this mode all of the DER units
that exist in the microgrid provide the energy for the consumer. There
are many reasons for going with this mode. The reasons can range from
dropped voltage in the grid, a lower price of energy in the microgrid, in
the event of fault in the grid, and so on. In this case the microgrid can
stay in stand-alone mode without requiring any power from the grid.
DER units in microgrid can be connected to the microgrid into two ways:
1. They are interfaced to the microgrid through rotating machines with re-
spect to slow dynamics and great inertia.
2. They are connected through a power electronic converter (PEC) with small
output impedance and fast dynamics in response to the changing of power.
So in the transition stage between the two different modes of microgrid
(stand-alone and connected to the grid), the two connection types can not have
the same dynamics response. In this situation, in order to increase response time
for a rotating machine, DS can be used for computational purposes, whereas in
a steady state condition the rotating machine can provide its power. Discussing
dynamics in detail is outside the scope of this project.
DER units and loads in the microgrid can be placed everywhere but not
necessarily located in the same physical location. All of the type of DER units
can contribute in the microgrid structure, such as diesel generators,fuel cells,
and all renewable energy sources, as solar cells (SC), or wind power generators,
and microturbines as well as other possible types (i.e., geothermal, biomass,
etc.), as depicted in Fig. 2.1.
5


Fig. 2.1 shows the existence of a DS that provides an uninterruptible supply
of electricity to the microgrid. Currently many types of DS, such as DC bat-
tery, Flywheel, Supercapacitor (SC), Superconduction Magnetic Energy Storage
(SMES), etc. are available for the microgrid.
The microgrid unit has an electrical connection point to the grid that is
referred to as the point of common coupling (PCC).
2.2.1 Power management strategies for microgrid
The objective of this project is to compute the amount of power generated by
any DER unit in the microgrid. To achieve this objective, a power management
strategy (PMS) is required. In the microgrid, every DER unit has a different
rate of flow with no dominant DG in the stand-alone mode of microgrid [12].
Therefore, the PMS must consider all types of DER units and create a control
methodology for energy flow in the microgrid in both modes (connected to and
disconnected from the grid). In this control, the PMS should supervise the
active and reactive power flow with high stability and reliability in the system.
The strategy for control should be to provide the peer-to-peer (P2P) and
plug-and-play (PnP) concepts in the microgrid.
In the context of the microgrid, P2P means there is no master control for any
DER unit, so if one DER unit is lost, the system can continue its operation with
an extra DER unit. PnP means, any DER units can be added in any location
of the microgrid without redesigning new control strategies for the system.
PMS has another responsibility which is to provide the decision in the dy-
namics in transition between the two modes of the microgrid, but this is outside
6


the scope of this project.
7


3. Optimal Power Flow (OPF)
3.1 Introduction
OPF was presented for the first time by Carpentier in 1962 [4]. OPF is an
algorithm that optimizes the objective function, such as the total cost of power
generation, magnitude and phase angle deviation of bus voltage, power losses,
and emission of pollutants, etc. The objective function can also be predefined
by the user. OPF finds the optimal set point of the electrical system which
satisfies the system power flow equations and all of the limit, constraints, and
contingencies, associated with the power network.
3.2 Formulation
In a power grid with N busses, if the voltage magnitude and the phase
angle of the busses are computed by utilizing the OPF formula then all of the
variables associated with the system, such as the active and reactive power
transferred between the buses and also the amount of the current in the lines,
will be satisfied. In the power system, there are three types of busses:
1. The PV bus or generator bus, with at least one generator connected to
the bus.
2. The PQ bus or load bus, where no generator is connected to the bus.
3. The slack bus or swing bus, this is defined by the reference for phase angle
of the voltage and can provide the necessary active and reactive power for
balancing the system at the optimized solution.
8


If we assume the power system includes n busses and ng generators, the
dimension of the admittance matrix Ybus will be a n x n" matrix which can be
denoted as:
( \
y 11 yi2 yin
2/21 ?/22 2/2n
Yh
bus
(3-1)
y2/ral 2/ra2 ynn y
where the elements of i/ij, i,j = 1,2,n of Ybus are constructed as defined in
[3] as follows:
The Ybus is symmetric.
The diagonal elements Yu of matrix are equal to the sum of the admittances
of all the components connected to the ith node.
The off-diagonal elements Yij are equal to the negative of primitive admit-
tances of all components connected between nodes i and j.
For this system the complex power injected matrix Sbus is defined as:
Si,, = (3.2)
where the Vbus is the vector of the bus voltage and Ibus is a vector and denotes
the bus current injected in the bus. Manipulating the equations, the results will
be as follows:
with Ibus (3.3)
9


(3.4)
Shi s Vbus{YbusVtyas)
Sbusi Vbusi /* '(yij^busj) -3 = 1 hi = 1,2 ,...,n (3.5)
n 7=1 i,j = 1,2 ,...,n (3.6)
C = he* Vj = \ Vj\ePj i,j = l,2,...,n (3.7)
Yij Gij T jBij b)-ij 9i-9j i,j = 1,2,...,n (3.8)
n 3 = 1 ~jBij)] h j = l,2,...,n (3.9)
where:
Sbug. : The complex power injection in bus i.
n
Sbusi = ^2 [|e|lCl(cos^i + jsin9ij)(Gij jB^)} i,j = 1,2 (3.10)
3 = 1
n
Pbusi = ^[{ViWVjKGij cosdij + BijSinflij)] %,j = 1,2,...,n (3.11)
i=i
n
Qbusi = ^2 [lCHCl(Gl7sinC? W'Cos^)] h3 = r 2, ...,ra (3.12)
3 = 1
where:
P&uSi : The active power injection in bus i.
Qbusi ' The reactive power injection in bus i.
In the power flow problem, the load is a known variable. In setting up the
formula, the active and reactive power injection to every bus should consider
the amount of the load and generation as follows:
PbuSi = P9r ~ Pu = f(v, 9) QbuSt = Qgi Qu = f(V, 9) (3.13)
10


where:
Pgi : The active power of the generator in the bus i.
Qgi :The reactive power of the generator in the bus i.
Pu : The active power of the load in the bus i.
Qia : : The reactive power of the load in the bus i.
According to the optional objective function and its constraints, alternative
formulas can be presented. The OPF will solve quadratic problem include mini-
mizing a quadratic function subject to linear and nonlinear equality, inequality,
that is of interest in this paper.
The formula in which the objective function minimizes the cost (C'j(P),
i = 1, 2,..., rig) or/and emissions (Totemmissioni, i = 1, 2,..., ng) or/and total losses
between the lines (Totiosses) in the system or/and any other objective function
can be set up as follows:
ng
Minimize: ^^Ci(P) i=l,2,...,ng
i 1
or/and
ng
Minimize: E Totemmissioni ^ 1; 2, ..., Tig
i= 1
or/and
Minimize: TotiOSses (3.16)
or/and . .
subject to:
inequality constraints:
(3.14)
(3.15)
11


$min ^ ^ @max i 1,2,..., n 1 (3.17)
Vmin | | V'max *=1,2,..., n (3.18)
P < P < P 1 mini 1 gi 1 maxi * = 1,2,... ,n9 (3.19)
Qmini Qgi Qmaxi * = 1,2,. ..,ng (3.20)
\SijFrom\ < RatC hj = !,2 ,..., n, i ^ j (3.21)
^'To Rat& hj = !,2 (3.22)
where:
\SijFrom \ and \SijTo \ are the magnitude of complex power which transferred from
bus(i) to bus(j) at the sending and receiving point,
equality constraints:
n
Pbusi = Pgi PLl ^ [\Vi\\Vj\(Gij cosOij + Bij sin%)] (3.23)
j=i
n
Qbusi = Qgi -Qu = ^2 [\Vi\\Vj\{Gij sindij Bij cos^)] (3.24)
3=1
From equations (3.14) to (3.24), the user can find the variables for the objective
function.
12


4. Droop Control
In the power system, active power and frequency are dependent on each
other. In a rotating machine, increasing the load will result in an increase
in torque. The torque is related to the rotational speed with relation to the
frequency of the system. This means that the lower the frequency, the higher
the load will be. Therefore PMS can try to determine a way to control the
system with these characteristics when the DER unit in microgrid is interfaced
to the system with an inverter.
4.1 Concept of droop control
To understand the concept of droop control, we return to the power system
for delivering complex power through a transmission line between two buses.
XL
Z = R+ JXL
Figure 4.1: Power flowing through a line between two sources
Fig. 4.1 shows the model of the transmission line between bus i and bus j with
13


impedance Z where:
Z = R + jXL
and the power flow from the bus i to bus j is defined as:
Pij I jQij
Sij =
s, = vXXX
LJ 1 y*
by using this definition:
'13 Vi Uj
(4.1)
(4.2)
(4.3)
(4.4)
(4.5)
with Eulers formula:
Sa

z*
(4.6)
e0ij = cos 9ij + j sin 6b,- (4.7)
we can separate the active and reactive power to this form:
Pit = xJffjpW Ml + |l)|VLsiny) (4.8)
Qv = A.yg(AV|V)| |V)|cosy |V'|flEin%) (4.9)
If in a transmission line the resistance of the line is much smaller than the
inductance, then equation (4.8) and equation (4.9) will be changed to this form:
= (4.10)
14


and
(4.11)
Furthermore, if the difference in angle % is very small, in computing the
radian, cos 66? = 1 and sin % = 6b/. Therefore:
From equation (4.12) the difference of the phase angle between two busses (6b,)
is dependent on the active power delivered between the buses, and from equation
(4.13) the magnitude voltages (|Ib| and \ Vi\) in two buses is influenced by reactive
power. This implies that by changing 6b/, the active power, and changing the
| lb | and | Vi |, the reactive power can be controlled.
4.1.1 Droop in microgrid
In relation to droop control in the microgrid, the frequency is used instead of
the phase angle for controlling the active power. The microgrid includes two or
more DGs which use droop characteristics to share the active and reactive power,
These DGs are interfaced with the inverter. When the load in the microgrid is
changed, PMS must have specific rules to determine the set point of every DG
for the sharing of active and reactive power at the optimized point.
Fig. 4.2 shows the frequency droop characteristic. In this figure the /0 is
the nominal frequency and P0 is the momentary set point for active power of
(4.12)
and
mm\ w) ~ xLQtj
(4.13)
15


Figure 4.2: Frequency droop characteristic
the inverter [7]. From the Fig. 4.2, this equation can be writhen as follows:
./Wo + ^Pi-Po) (4.14)
Where Kp is the slope of the characteristic and /i is the the frequency of
the inverter for generating the power Pi from frequency droop characteristics.
In Fig.4.3, the voltage droop characteristic is shown.
From Fig.4.3, the equation can be written as follows:
I Vi | = |Vo| + Kq{Qi Qq) (4.15)
Where Vo is the nominal voltage and Qo is the reactive power of the inverter.
From equations (4.14) and (4.15) by adjusting P and Q independently, the
frequency and magnitude of the voltage can be determined [7].
4.2 Load sharing through P-f control
16


Voltage magnitude
Figure 4.3: Voltage droop characteristic
In a microgrid with two or more DGs, changing the load in both modes
(connectd to the grid or disconnected from grid), PMS uses the frequency droop
characteristics to change the set-up of the operating point to access the local
power balance in the new situation of loading. As an example, a microgrid which
has three DGs with droop characteristics is shown in Fig.4.4.
Every DG has a different slope in frequency droop characteristics (KPl, KP2
and Kp3).
When the microgrid is connected to the grid, the main grid is large and ro-
bust enough to compensate any variation in the microgrid [11]. So the frequency
in the microgrid will be the nominal frequency of the grid (/o). From Fig. 4.4
and their droop characteristics, the three DGs provide the load Pl0 as follows:
Pl = Pi + PS + PS (4.16)
17


Frequency
Figure 4.4: Frequency droop characteristic for three DGs
The Pl0 equals:
PL0=PL-Pgrid (4-17)
The Pgrid denotes the amount of load in the microgrid shared by the grid,
and the Pl is the total loads in the microgrid. Now if the microgrid for any
reason is disconnected from the grid, with the same load Pl, the frequency of
the system will be changed to /i and the power share for every DG will be P\,
Pf and Pf, where Pl equals:
PLl = Pl + Pf + Pf (4.18)
In this case:
PLl = PL (4.19)
If the equation (4.14) is developed for every DG, this equation can be written
as follows:
h=h + K,H(Pi-n) (4.20)
18


with:
3 3
Pl, =J2P! fi = £Po A Pl = PLi-Pl (4.21)
i 1 i= 1
Using equation (4.16), equation (4.18), and equation (4.21) the amount of
power for every DG will be:
pi
A PL^~
____APi
3 1
Ef
1 ^Pi
1=1
pi
r0
(4.22)
19


5. OPF with Droop Control
An experimental microgrid based on distributed generators (DG) consists
of several modules: the utility power grid, inverters and loads. This chapter
will address inverter management and control. Every inverter has autonomous
droop control characteristics for sharing the active and reactive power amongst
DCs to meet the load.
5.1 Mathematic model
Chapter 4 explained the load sharing process through the p-f controller and
equation (4.22) shows the amount of power for every DG. In this calculation, the
total of the power is equal to the loads, but the reality the system is nonlinear
and the losses should be considered as the equation will be:
ng
Y P9i = Pioad + Losses (5.1)
i= 1
The Losses are not constant and are dependent on the voltage of the buses.
In the autonomous microgrid mode the Losses vary with frequency because the
impedance is dependent on the frequency.
To connect many inverters parallel in the microgrid, the frequency of the
whole system should be by making them the same. The frequency vs. active
power characteristics for the DG in bus i plot is shown in the Fig. 5.1.
According to the Fig. 5.1 the following equation for Frequency-droop char-
acteristic can be written:
20


Figure 5.1: Frequency droop characteristic for DG in bus i
f = KpiPgi + bPi (5-2)
In this equation, the amount of power is dependent on the frequency, so with the
constant slope in the lower frequency the output power is larger. In stand alone
mode where the grid is disconnected, every DG should produce more power in
the same lower frequency.
The plotting of the voltage vs. reactive power is shown in Fig. 5.2. And
equation (5.3) shows the relationship of magnitude of the voltage with respect
to the reactive power.
Ilflil = KqiQgi + bqi (5-3)
Where / and \ Vgi\ are the frequency and amplitude of the voltage of the DG
in bus i respectively. Kpi and Kqi are defined by the corresponding slopes of
droop.
21


Voltage magnitude
li
\
Figure 5.2: Voltage droop characteristic for DG in bus i
When the microgrid is connected to the grid, the frequency of the system will
be the same as the utility grid. To share the active and reactive power between
the DGs at the optimized point, the slopes Kpi and Kqi can be adjusted. Hence,
the impedance between two buses i and j can be written as :
Zij = Rij + 2n f'Lij where: / = 60 (5.4)
In autonomous microgrid mode the load will be different, and in order to
get excellent power sharing, the frequency and magnitude of the DG voltage will
be changed. So, the change of impedance between two buses i and j caused by
changing the frequency can be expressed as follows:
X = 0L and vu) = ^+:)vfa no
For a network with n independent busses, and ng generators, the dimension
of admittance matrix will be a n x n and this matrix will be a function of /
22


as follows (see Chapter 3 and equation 5.5):
^yuU) yi2(f) yin(f)^
V21 (/) 1/22(/) y2n{f)
Ybus(f) =
(5.6)
y Uni (/) Vn2 (/) Vnn{f) y
for this network with n buses the following equations can be written:
= vbU8rbus
where Ibus is the bus current injection vector.
Sbus = Vbus[Ybus(f)Vbus]* (5.8)
&bu.Si Vbusi -i=i i,j = 1,2,...,n (5.9)
n = £ j=i [lU^I/y(/)] bj = 1,2,...,7i (5.10)
15 = 15 eJ< l-j = \Vj\e bj = 1,2,...,7i (5.11)
Yiiif) = Gn(f) +jBiAf) % = 0i~ h3 = 1, 2, ...,n (5.12)
n
su, = £ [m\V,\e-{Gi,U) - h) = 1,2,..., (5.13)
i= i
where:
Stas* : The complex power injection in bus i.
n
Sbus, = Y [|15||15l(cos^ +isin%)(Gq(/) jBij(f))} hi = 1,2,...,71
3 = 1
(5.14)
23


Pbusi = ^[iViWVjKGijif) cos dij + Bijif) sindij)] i,j = 1,2 (5.15)
j=i
n
Qbusi = WVi\\Vj\(Gij(f) sinij Bijif) cos dij)] i,j = 1,2 (5.16)
3 = 1
In the power flow problem, the load demands are known variables, so we
can define the following bus power injections as:
Pbasi = Pgi PlA = fiV, e, /) = Qgl QLl = f(V, 9, /) (5.17)
where:
Pgi : The active power of the DG in the bus i.
Qgi : The reactive power of the DG in the bus i.
Pn : The active power of the load in the bus i.
Qn : The reactive power of the load in the bus i.
From equation (5.2) and equation (5.3) for every DG, the output power (ac-
tive and reactive) can be computed corresponding to their droop characteristics
P,ji = BP Qsi = IMp** (5.18)
In optimal power flow the main objective function is the total generation
cost which is defined as follow:
n
Total^ = ^2 Ci(P) (5-19)
i= 1
where:
Ci(P) is the cost function for the DG in bus i When the grid is connected to
24


the microgrid, all of the DGs provide their energy froman optimal point at the
frequency of grid in Kpio and Kqio. Those are defined by the corresponding slopes
of their droops. During stand-alone microgrid mode the grid is disconnected
from the microgrid and the power delivery from each DG will be changed. In
this variation the amount of power and reactive power sharing will be changed
according to equation (5.2) and equation (5.3) in the new frequency /. The
change is shown in Fig. 5.3 and Fig. 5.4.
Reactive Power
Figure 5.3: Frequency droop characteristics for DG in bus i in different load
However, these points are not optimized points. After calculating the opti-
mized point, the slope and Y-intersect of the droop will be changed to the new
value Kpi, Kqi, bPi and bqi. Fig. 5.5 and Fig. 5.6, depict the direction of changes
for both droops.
5.2 Mathematical formula for optimization
To set up the equations in an optimal power flow, other objectives can be
added to the optimization problem to minimize the deviation of droop charac-
teristics. This means the minimizing of AKpi, AKqi,Abpi and Abqi.
25


Voltage magnitude
Figure 5.4: Voltage droop characteristics for DG in bus i in different load
Figure 5.5: Direction of frequency droop characteristics changing for DG in
bus i for different load
The equation can be written as follows:
n
Minimize: Totalcost = E Ci{P) i = 1,2, ...,ng (5.20)
i= 1
and
26


Voltage magnitude

Figure 5.6: Direction of voltage droop characteristics changing for DG in bus
i for different load
Minimize: (AK^, Ab^, AK^, Ab^fj i = l,2,...,ng (5-21)
Such that:
$min @max i = 1,2,... ,n (5.22)
Vmin A VI Kma: i = = 1,2,. ..,n (5.23)
P 1 mini Pgi < P 1 maxi i = = 1,2,. ,ng (5.24)
Qmini < Qgi Qmaxi i = = 1,2, ...,ng (5.25)
/ min j Pi fmax (5.26)
~^Pmin 27


Kqmin Fqi T Kqrnax i = 1, 2,..., ng (5.28)
bpmin. bpi T bPrnax ^ = 1,2, , j ng (5.29)
bqmm bqi < bqmax 3 = 1,2, ...,ng (5.30)
1 ^^3 From < Rate i,j = 1,2,... ,n,i j (5.31)
\SHTo\ < Rate i,j - = 1,2,..., n, ^ j (5.32)
n ElNIMGaCOcosfly + Bij(f) sin 9ij)\ * = 1,2,..., n
3 = 1
(5.33)
WVi\\Vj\ (Go(/)sin% - BM) COS $ij)\ * = 1,2, ...,**
3 = 1
(5.34)
p f bPi ^ 91 Kpi = 1,2,... ,ng (5.35)
Qgi 1 Fgi | bqi Kqi *=1,2,. -,ng (5.36)
AKpi x K lvpl 1^pi 0 * = 1,2 , ...,ng (5.37)
o & 1 * = 1,2 , ...,ng (5.38)
Abpi bpi bpio i = 1,2,... .,ng (5.39)
3^-bqi hqi bqi0 *=1,2,. ..,ng (5.40)
The objective functions can be evaluated by using the weighted sum method
that explained in [10],[15],[17],[16],[6]. The weighted objective function for the
system can be expressed as follows:
F = W\Fi + W2F2 + ... + WiFi = 1,2, ...,4ng + 1
(5.41)
28


where: Fi = E Ci(P),
i= 1
F2 = ^Klv Fs =
F* = cm ^ <1 f5 = CM <1
Fa-2 = F l'\i = ^4*+l CM 5? <4 II
and,
W\ + W2 + W3 + ... + Wa+i 1
(5.42)
29


6. Case of Study
6.1 Introduction
The case that we picked up for this study is shown in the Fig.6.1 [12]. Fig.
Grid Utility source
1000 MWA
X/R=22.2
Vbase=13.8 KV
Sbase= 100MVA
, 69/13.8
' 15MVA
DG2
0 732+jO 095/
13 8/0
1 5MVA
6 48+J38 3

L5 L4
3 8/0 4!
25 MV,
6+J48 0%

L3
DG1
13 8/2
2 5 MV A
3 29+J2 3%
4 vXj:
\ /'pi
T ?
L2
3 8/0 48
0MVA
21+J57 5%
L1
Figure 6.1: The case of study in rnicrogrid
6.1 is a single-diagram of a 13.8-KV distribution system used to investigate a
possible microgrid PMS. This system has three DG units, with different rating,
ie.,DGl (1.8-MVA), DG2 (2.5-MVA) and DG3 (1.5-MVA). All of the DG units
30


can be dispatched [12]. To study this case, all of the buses for the loads and DGs
are shown in Fig. 6.2. As seen as in Fig. 6.2, the system has fourteen buses,
Figure 6.2: Case of study with all buses
and three DGs: DG1, DG2 and DG3 connect respectively to Bus6, Bus 13, and
Bus 12. The values of the impedance between the lines are stated in Table 6.1.
This microgrid is connected to the grid through Bus\\. In Table. 6.2 the
type of the buses with the amounts of the loads (active and reactive power) are
listed. In the column Type of bus of this table, the No.O means that there
is no generator connected to the bus, the No.2 means the DG is connected to
the bus and the No.3 means that the bus is connected to the grid.
31


Table 6.1: Amount of impedances between busi and busj in pu.
Busi Busj Rij XLl3
1 7 0.5314 5.7635
1 5 0.4022 0.2395
1 6 0.21 1.094
1 14 0.3976 0.5127
2 14 0.6141 0.3066
3 14 0.6065 1.015
3 8 0.816 4.8332
3 9 0.6903 0.000
14 11 0.0818 0.5826
12 4 0.822 0.48332
13 3 0.822 0.48332
4 14 0.3564 0.2661
10 4 0.822 4.8332
The data for all of the generators are listed in the Table 6.3.
The generator cost curves are presented by quadratic functions [1]. The param-
eters a, b, and c in Table. 6.3, are the coefficient of the quadratic function and
are expressed in the following form:
Q = OiPl + btPm + bt i = 1,2, ...,ng (6.1)
Some renewable energy sources such as solar and wind power do not have
any fuel costs because the resources are free, so in the OPF equation, this type
of distributed generator resource unit should not be considered for dispatch. For
32


Ta ole 6.2: Data for types of bus
Bus No. Type of bus Pl(MW) Ql{MVAR)
Bust 0 0 0
Bus2 0 0.6 0.42
Bus3 0 0 0
Bus4 0 0 0
Bus5 0 0.61 0.42
Bus6 2 0 0
Bus7 0 0.600 0.390
Bus8 0 0.7 0.45
Bus9 0 0.8 0.5
BuslO 0 0.9 0.61
Busll 3 0 0
Busl2 2 0 0
Busl3 2 0 0
Busl4 2 0 0
Table 6.3: Data for all of the generators
DGs.Name BusNo. Type Pmax(Mws) Pmin(MW) Qmax(MV AR) Qmin(MVAR) a c Kp K, K
DG1 6 1.80 0.000 1.80 -1.80 0.00482 7.970 78.00 -0.0282 60.0137 0.0056 1.0306
DG2 13 2.50 0.000 2.50 -2.50 0.00156 7.920 561.00 -0.0536 60.0453 -0.1488 1.0183
DG3 12 1.50 0.000 1.50 -1.50 0.00194 7.850 310.00.0 -0.2537 60.2417 -0.2157 1.0255
Grid 11 1000 0.000 1000 -1000 0.00139 7.060 500.00 -
coding in the Table 6.3, if the generator is able to dispatch the type then it is
(1), otherwise it is (0). Droop characteristic parameters are stated in the last
four columns of Table 6.3 (Kp, bp, Kq and bq). The loads L\ in Bust, L2 in
Bust,, L3 in Bus2, L4 in Bu.sg, L5 in Busg and Lq in Busio are stated and the
33


value of the loads are expressed in Table 6.2.
6.2 Microgrid connected to the grid
In this case, the grid is connected to the microgrid, and the power necessary
for the loads is shared by DGs and the grid. Table 6.4 shows the values of the
magnitude, the phase angle of the voltage of each bus, and the amount of the
power and reactive power shared by any generators, as well as parameters of
droop characteristics for every DG at the optimized point.
Table 6.4: Data for all generators and buses in the connected mode with the
grid
DGs.Name BusNo. vu ^-^deg Pl(MW) QUMVAR) Pg(MW) Q,(MVAR) Kp bp \ 1
Busl 1.062 -0.804 0 0 0 0 - - -
Bus2 1.058 -0.537 0.600 0.420 0 0 - -
Bus3 1.058 -0.887 0 0 0 0 -
Bus4 1.063 -0.585 0 0 0 0 -
Bus5 1.059 -0.792 0.61 0.420 0 0 - -
Bus6 DG1 1.075 -0.610 0 0 0.5753 1.1608 -0.0293 60.0169 0.0656 0.9987
Bus? 1.037 -2.496 0.600 0.390 0 0 -
Bus8 1.031 -2.472 0.700 0.450 0 0 -
Bus9 1.033 -2.319 0.800 0.500 0 0 - - -
BuslO 1.026 -2.607 0.900 0.610 0 0 - - -
Busll Grid 1.066 0.000 0 0 1.9840 0.2543 -
Busl2 DG3 1.073 -0.660 0 0 0.8708 0.6927 -1.0648 60.9272 0.3767 0.8116
Busl3 DG2 1.068 -1.050 0 0 0.8375 0.8847 -0.1078 60.0903 0.6452 0.4973
Busl4 1.063 -0.575 0 0 0 0 - - -
34


6.3 Microgrid disconnected from the grid (autonomous mode)
If the grid is disconnected from the microgrid, it will not contribute to
sharing any generated power with the DGs. In this situation, all of the DGs in
the microgrid should provide the necessary power required for the load without
the grid.
6.3.1 Sharing the power from droop characteristics
Initially, when the grid is disconnected from the microgrid, the droop char-
acteristics of all inverters will not be changed. New amounts of power generated
for every DG are computed by equation 4.22 in the same frequency. In Table.6.5,
the amounts of active and reactive power are generated by any DG without any
change in the parameters of droop characteristic.
6.3.2 Sharing the power with OPF
The amount of power for every DG in previous section will not be at op-
timized point. According to the algorithm expressed in chapter 5, the multi-
objective function optimization with weighted parameter can be used as follows:
13
(6.2)
i= 1
where:
n
fwJTvp)
i= 1
35


Table 6.5: Data for microgrid in autonomous mode before optimization
DGs.Name BusNo. ivw Pl(MW) Ql(MVAR) P,(MW) Q,(MVAS) Kp bp K, K
Busl - 1.073 0.125 0 0 0 0
Bus2 1.063 0.037 0.600 0.420 0 0 -
Bus3 1.065 -0.102 0 0 0 0 -
Bus4 1.067 -0.007 0 0 0 0 -
Bus5 1.070 0.137 0.61 0.420 0 0
Bus6 DG1 1.091 1.110 0 0 2.1114 1.4093 -0.0293 60.0169 0.0656 0.9987
Bus? 1.048 -1.530 0.600 0.390 0 0 - - -
BusS 1.038 -1.665 0.700 0.450 0 0 -
Bus9 1.041 -1.514 0.800 0.500 0 0
BuslO 1.031 -2.009 0.900 0.610 0 0
Busll Grid 1.067 0.000 0 0 0.000 0.000 - -
Busl2 DG3 1.078 -0.077 0 0 0.9130 0.7061 -1.0648 60.9272 0.3767 0.8116
Busl3 DG2 1.078 -0.169 0 0 1.2548 0.9004 -0.1078 60.0903 0.6452 0.4973
Busl4 1.067 0.000 0 0 0 0 -
f2 = AK. F3 =
f4 = f5 = ^K2q2
F6 = F7 =
Al 00 J A62 Pi f9 = ^2qi
F\o - -- Ab2 V2 Fn = Ab2 Q2
to II -- Ab2 V3 Fl3 = Ab2 qs
and
W1 + W2 + W3+W4 + W5 + W6+W7 + W8 + W9 + Ww + Wn + W12 + W13 = 1 (6.3)
For the case study, the problem solved with multi-objective optimization by
fmincon" command in MAT LAB (see appendix A). The values of objectives
36


F1} F2, F3, F4j F5, F6, F7, Fg, F9, F10, Fn, F12 and Fi3 are normalized to be
consistent with the weights assigned to the objective functions. The functions
are normalized using the form:
Fi-Fi .
q ^ 0 vmin 1
- ~p. _ /'. -
'Lmax 'I'min
(6.4)
where the Fj is the ith objective function; id- is minimum value of vh objective
function; and Fimax is the maximum value of ith objective function in the different
weighted, that is assigned to the objective function. Results are shown in Table
6.6.
Table 6.6: The normalized of objective functions in autonomous mode after
optimization
Ft F2 Fi fa F't Ft Fr Fs F* Fl0 Eli Fi2 Fi
0,085106383 0.0143888 34 0,088864526 0,697111083 1.07868E-18 3.87819E-05 1,000012322 0128438349 0,124272989 0.135580785 1.12556E-17 1,000139336
0,9 0,00833 0,687043262 0.032172178 1.245 37E-18 0,841746085 0,27036738 0.096737545 0,119596111 0,032 13545 1.000000015 0.999949561 0,658812061 0.788466167
O.S 0.0167 0,709219858 0.022172178 0.416655789 0.142153052 0,315088821 0.128846533 0,581670718 0011908363 0,03949486 2.4568 IE-17 0,781715352 0,902958985
0.7 0,025 0,680851064 0.099674757 0.10919729 0,480503545 0.000129537 3.87819E-05 0.035064372 0,418399402 1.000002212 0,970045689 0.078618498 0,027269824 0,090394067
0.6 0.0333 0,695035461 3.02 4 38E-19 0,164773485 0,547370822 1,000010934 0.215431531 0.040584768 -4.87464E-06 2.70814E-05 0,049037947 0.644746978 0,028824834 0.210738986
O.S 0,04166 0,758865248 0.424917685 1.00003586 0,196919165 0,44534645 0,146438493 0.013043907 0,550473308 -1.50154E-17 0,013757267 0,291117184 0,97142 7611 0,166015378
0,4 0,05 0,680851064 1.000028609 0,567440913 0,3605 10891 0.051668461 4.28715E-18 0.015638986 0.105922775 0.009402929 0,033540092 0.008405018 1,000210091 -1,1165 IE-17
0.3 0.0583 0,737588652 0.688 3845 24 0,354241764 0,218181902 1.67047001 Q, 12638824 0.011437077 0,267534592 0.010167426 0,003169549 1.047614681 0,982890181 0,000332506
0.2 0.0667 0,723404255 0.926835785 0,516317125 0,118091182 0018364031 0,009272961 0.003304351 0.27557615 0.008642708 0,01149258 0.048573027 0,039976134 0,018240316
0,1 0,075 0,872340426 0,27919063 0.10761944 0,015645135 0100306981 0,291417765 1.97157E-06 0,040762902 0.000224224 0,001106104 Q,021120473 0,014092911 0,055126289
0.0833 0.057076946 0,018571209 8,50237E-08 0.047681697 Q,00042456 0.005651012 0,007778609 0.000230869 6.62201E-Q7 Q,027824141 0,008042287 0,062017075
The droop characteristics of all DGs are depicted in Fig. 6.3, Fig. 6.4 and Fig.
6.5 for two conditions: no optimization and optimization with W\ = 0.7, because
this value of W\ = 0.7 is the most suitable weighted coefficient for optimization.
37


Frequency-droop characteristic for DG1 Voltage-droop characteristic for DG1
Active power(kw) Reactive power(kVar)
Figure 6.3: Droop characteristic deviation for DG1
Frequency-droop characteristic for DG2
Voltage-droop characteristic for DG2
Figure 6.4: Droop characteristic deviation for DG'2
From this algorithm, the
operator can define the new slope of droop char-
acteristics at the optimized point.
38


Frequency-droop characteristic for DG3
Voltage-droop characteristic for DG3
Figure 6.5: Droop characteristic deviation for DG3
39


7. Conclusion
This study has presented a methodology for determining optimal power
flow in a microgrid in a steady state, for either one of two modes of operation,
i.e., connected to the grid and disconnected from the grid, acting as a stand-
alone isolated microgrid. This study develops an algorithm which defines the
droop characteristic parameters of an inverter. The goal is to find the output in
determining the optimal power to the user. This algorithm assists the user to
dispatch power in the microgrid.
If the load changes in the microgrid, this results in changes in the droop
characteristics in both the slope ( Kp and Kq) and the Y-intersect ( bp and
bq) deviating at the optimized point. In addition, when the microgrid is dis-
connected from the grid, the frequency is different from the frequency of the
grid. From the frequency-droop characteristics the amount of power of DG will
be changed. The droop characteristic parameters of the inverter-based DGs has
computed to reach proper load sharing at optimized point rather than to be used
a fixed value. The power necessary for the system is provided by the optimal
power flow that minimizes the cost of fuel.
Transition modes from the connected operating mode to the isolated oper-
ating mode in a microgrid were not covered in this paper, but only the steady
state condition in each mode is addressed in this paper. The dynamic transition
between two modes can be left to a future discussion as an extension of this
study scope.
40


APPENDIX A. fmincon in MATLAB
fmincon function in MATLAB is proposed to find the minimum of a con-
strained non-linear multi-variable function. Otherwise fmincon finds a con-
strained minimum of a function of several variables. The equation can be solved
by this function expressed as:
minf(x)
suchthat :
c{x) < 0
ceq(x) = 0
Ax < b
Aeqx = beq
lb < x x, b, beq, lb and ub are vectors, A and Aeq are matrices, c(x) and ceq(x) are
function that return vectors (@nonlincon), and f(x) is a function that returns
a scalar (@fun). f(x), c(x), and ceq(x) can be nonlinear function. The function
of fmincon in MATLAB is expressed as:
[x, fval, exit flag, output, lambda] = fmincon(@fun, xO, A, b, Aeq, beq, lb, ub, @nonlcon)
(A.l)
41


where
objective Objective function
xO Initial point for x
A Matrix for linear inequality constraints
b Vector for linear inequality constraints
Aeq Matrix for linear equality constraints
beq Vector for linear equality constraints lb Vector
nonicon Nonlinear constraint function
lb Lower bound
up Upper bound
X optimized value
fval the value of the objective function at the soluti
lambda Lagrangean multipliers
exitflag specification at algorithm terminatio
solver 'fmincon'
42


APPENDIX B. The code of software in MATLAB
B.l Reading the file
function S=Read_filecf(S)
filename=uigetfile('*.cf') ;
S.f ilename=f ilename;
filecf=fopen(S.filename,'r') ;
s=fgetl(filecf) ;
while strcmp(s,'')
s=fgetl(filecf);
end;
S.CaseName=s;
while strcmp(s(1:min(3,length(s))),'BUS')~=1
s=fgetl(filecf);
end
s=fgetl(filecf);
while strcmpCs,'')
s=fgetl(filecf) ;
end;
S.Bus.number= [];
S. Bus .name= [] ;
/ S.Bus.zone=[] ;
43


S. Bus type= [] ;
S. Bus V= [] ;
S. Bus ang= [] ;
S. Bus MWL= [] ;
S. Bus MVARL= [] ;
S. Bus MWG= [] ;
S. Bus MVARG= [] ;
S. Bus Bkv= [] ;
S. Bus DV= [] ;
S. Bus ,max= [] ;
S. Bus ,min= [] ;
S. Bus G= [] ;
S. Bus B= [] ;
while 3(1:1)='-'
k2=0;
[kl k2]=find_k(s,k2);
ite=str2num(s(kl:k2));
S.Bus.number(ite,l)=ite;
[kl k2] =f ind_k(s,k2) ;
S.Bus.name=strvcat(S.Bus.name,s(kl:k2));
[kl k2] =f ind_k(s,k2) ;
S.Bus.type(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
44


S.Bus.V(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.Bus.ang(ite,l)=str2num(s(kl:k2))*pi/180;
[kl k2]=find_k(s,k2);
S.Bus,MWL(ite,l)=str2num(s(kl:k2))/100;
[kl k2]=find_k(s,k2);
S.Bus.MVARL(ite,l)=str2num(s(kl:k2))/100;
[kl k2]=find_k(s,k2);
S.Bus.MWG(ite,l)=str2num(s(kl:k2))/I00;
[kl k2]=find_k(s,k2);
S.Bus.MVARG(ite,l)=str2num(s(kl:k2))/100;
[kl k2]=find_k(s,k2);
S.Bus.Bkv(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.Bus.DV(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.Bus.max(ite,l)=str2num(s(kl:k2));
[kl k2] =f ind_k(s,k2) ;
S.Bus,min(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.Bus.G(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.Bus.B(ite,l)=str2num(s(kl:k2));
45


s=fgetl(filecf);
end
while strcmp(s(1:min(6,length(s))),'BRANCH')
s=fgetl(filecf);
end
s=fgetl(filecf);
while strcmp(s,'')
s=fgetl(filecf) ;
end;
S.Branch.fr=[] ;
S.Branch.to=[] ;
S.Branch.r=[] ;
S.Branch.x=[] ;
S.Branch.b=[] ;
S. Branch. Mvar= [] ;
ite=l;
while s(l:l)~='-'
k2=0;
[kl k2]=find_k(s,k2);
S.Branch.fr(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.Branch.to(ite,1)=str2num(s(kl:k2));
46


[kl k2]=find_k(s,k2);
S.Branch.r(ite,1)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.Branch.x(ite,1)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.Branch.b(ite,1)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.Branch.Mvar(ite,l)=str2num(s(kl:k2))/100;
s=fgetl(filecf) ;
ite=ite+l;
end
while strcmp(s(1:min(9,length(s))),'GENERATOR')
s=fgetl(filecf);
end
s=fgetl(filecf) ;
while strcmp(s, ')
s=fgetl(filecf) ;
end;
if strcmp(s(l),'(');
s=fgetl(filecf);
end;
while strcmp(s,'');
s=fgetl(filecf);
end
47


S.G.Bus=[] ;
S.G.Type=[] ;
S.G.Kind=[] ;
S.G.Pmax=[] ;
S.G.Pmin=[] ;
S. G. Qmax= [] ;
S. G. Qmin= [] ;
S. G. a= [] ;
S.G.b=[] ;
S. G. c= [] ;
% S.G.f=[] ;
S.G.Kp=[] ;
S. G. Yp= [] ;
S.G.Kq=[] ;
S. G. Yq= [] ;
S.G.P_ava_UnDis=[];
ite=l;
while 3(1:1)='-'
k2=0;
[kl k2]=find_k(s,k2);
S.G.Bus(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.G.Type(ite,l)=str2num(s(kl:k2)) ;
[kl k2]=find_k(s,k2);
48


S.G.Kind(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.G.Pmax(ite,l)=str2num(s(kl:k2))/100;
[kl k2]=find_k(s,k2);
S.G.Pmin(ite,l)=str2num(s(kl:k2))/100;
[kl k2]=find_k(s,k2);
S.G.Qmax(ite,l)=str2num(s(kl:k2))/100;
[kl k2]=find_k(s,k2);
S.G.Qmin(ite,l)=str2num(s(kl:k2))/100;
[kl k2]=find_k(s,k2);
S.G.a(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.G.b(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.G.c(ite,l)=str2num(s(kl:k2));
*/. [kl k2]=find_k(s,k2);
% S.G.f(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.G.Kp(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.G.Yp(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.G.Kq(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
49


S.G.Yq(ite,l)=str2num(s(kl:k2));
[kl k2]=find_k(s,k2);
S.G.P_ava_UnDis(ite,l)=str2num(s(kl:k2))/I00;
s=fgetl(filecf) ;
ite=ite+l;
end
while strcmp(s(1:min(9,length(s))),'FREQUENCY')
s=fgetl(filecf);
end
s=fgetl(filecf);
while strcmpCs,'')
s=fgetl(filecf) ;
end;
S.f=[] ;
ite=l;
k2=0;
[kl k2]=find_k(s,k2);
S.f=str2num(s(kl:k2)) ;
S.Bus.pq=find(S.Bus.type==0);
S.Bus.pv=find(S.Bus.type==2);
S.Bus.s=f ind(S.Bus.type==3);
S.G.Dis=find(S.G.Kind==l);
S.G.UnDis=f ind(S.G.Kind==0);
50


/ S. G. Inv=f ind (S. G. Type==0) ;
N=length(S.G.Bus);
Nl=N-length(S.G.UnDis) ;
a=max(S.Bus.number);
S.G.Pmax_Dis=S.G.Pmax(f ind(S.G.Kind==l));
S.G.Pmin_Dis=S.G.Pmin(find(S.G.Kind==l));
fclose(filecf);
B.2 Calculate Y-Bus
function S=calculate_YBus(S)
% S .filename^
% S=Read_filecf(S);
S.YB=sparse(diag(S.Bus.G+j*S.Bus.B));
S.YB=S.YB+sparse(S.Branch.fr,S.Branch.fr,...
1./(S.Branch.r+j *S.Branch.x)+j *S.Branch.b/2,
max(S.Bus.number),max(S.Bus.number));
S.YB=S.YB+sparse(S.Branch.to,S.Branch.to,...
1./(S.Branch.r+j *S.Branch.x)+j *S.Branch.b/2,
max(S.Bus.number),max(S.Bus.number));
S.YB=S.YB+sparse(S.Branch.fr,S.Branch.to,...
51


-1./(S.Branch.r+j*S.Branch.x),max(S.Bus.number),max(S.Bus.number));
S.YB=S.YB+sparse(S.Branch.to,S.Branch.fr,...
-1./(S.Branch.r+j*S.Branch.x),max(S.Bus.numb)
B.3 Computational code
clear all
close all
clc
/(starting the reading the data.
S.filename=' ';
S=Read_filecf(S);
%For putting the bus that is connected to the grid in the lasr raw:
S=0rg(S);
%NN is the number of generators without the grid.
NN=length(S.G.Bus)-l;
xl = input('IS THE SYSTEM CONNECTED TO THE GRID ? ','s')
if (xl=='n')
S.G.Busl=S.G.Bus(1:NN);
S.G.Kindl=S.G. Kind (CNN) ;
S.G.Typel=S.G. Type (CNN) ;
S.G.Pmaxl=S.G.Pmax(1:NN);
S.G.Pminl=S.G.Pmin(l:NN);
S.G.Qmaxl=S.G.Qmax(1:NN);
52


S.G.Qminl=S.G.Qmin(l:NN);
S.G.al=S.G.a(l:NN);
S.G.bl=S.G.b(l:NN);
S.G.cl=S.G.c(l:NN);
/ S.G.f=[] ;
S.G.Kpl=S.G.Kp(1:NN);
S.G.Ypl=S.G.Yp(1:NN);
S.G.Kql=S.G.Kq(l:NN);
S.G.Yql=S.G.Yq(l:NN);
S.G.Dis=S.G.Bus(f ind(S.G.Kindl==l));
S.G.UnDis=S.G.Bus(f ind(S.G.Kindl==0));
S.G.Invl=f ind(S.G.Typel==l);
Nl=length(S.G.Invl);
%Ng is the number of generator,
%Nd is the number of generator that can be dispatchable
%and Nb is the number of busses.
Ng=length(S.G.Bus1);
Nd=length(S.G.Dis);
Nb=max(S.Bus.number);
/S.G.Pmax_Dis is the matrix of maximum reactive power of dispatching
%generators and S.G.Pmin_Dis is the minimum reactive power of dispatching
%generators.
53


S. G.Pmax_Dis=S.G.Pmaxl(find(S.G.Kindl==l));
y.
S.G.Pmin_Dis=S.G.Pminl(find(S.G.Kindl==l));
S.G.Kp2=S.G.Kpl(S.G.Invl);
S.G.Yp2=S.G.Ypl(S.G.Invl);
S.G.Kq2=S.G.Kql(S.G.Invl);
S.G.Yq2=S.G.Yql(S.G.Invl);
*/. Initials value
xO=[zeros(Nb,1);ones(Nb,1);(S.G.Pmin_Dis+S.G.Pmax_Dis);(S.G.Qminl+S.G.Qmax
A= [] ;
b= [] ;
Aeq=[] ;
beq= [] ;
/0upper bond of variable.
ub=[pi*ones(Nb,1);1.l*ones(Nb,1);S.G.Pmax_Dis;S.G.Qmaxl;60.8;Inf*ones(16,1
%lower bond of variable.
lb=[-pi*ones(Nb,1);0.9*ones(Nb,1);S.G.Pmin_Dis;S.G.Qminl;59.2;-Inf*ones(16
options=optimset('DerivativeCheck','on','Diagnostics','on','Display',...
54


'iter','FunValCheck','off','GradObj','off' ,'MaxFunEvals',30000,...
'Maxiter,10000,'Hessian','off ');
[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(@fun,x0,A,b,Aeq,beq,It
Nd=length(S.G.Bus1)-length(S.G.UnDis);
Ng=Nd+length(S.G.UnDis);
Nb=max(S.Bus.number);
S.Bus.V=x(Nb+1:2*Nb);
S.Bus.ang=x(1:Nb);
S. Bus MWG= [] ;
S. Bus MVARG= [] ;
S.Bus.MWG(S.G.Dis)=x(2*Nb+l:2*Nb+Nd)';
S.Bus.MWG(S.G.UnDis)=S.G.P_ava_UnDis(find(S.G.Kindl==0))';
S.Bus.MVARG(S.G.Busl)=x(2*Nb+Nd+l:2*Nb+Nd+Ng)';
f=x(2*Nb+Ng+Nd+l);
S.G.Kp(S.G.Inv1)=x(2*Nb+Ng+Nd+2:2*Nb+Nd+Ng+Nl+1);
S.G.Yp(S.G.Invl)=x(2*Nb+Ng+Nd+Nl+2:2*(Nb+Nl)+Nd+Ng+1);
S.G.Kq(S.G.Invl)=x(2*(Nb+Nl)+Nd+Ng+2:2*(Nb+Nl)+Nd+Ng+Nl+1);
S.G.Yq(S.G.Invl)=x(2*(Nb+Nl)+Nd+Ng+Nl+2:2*(Nb+Nl)+Nd+Ng+Nl+Nl+1);
S.f=x(2*Nb+Ng+Nd+l);
55


S.Bus.MWL(f ind(S.Bus.type==2))=0;
S.Bus.MVARL(find(S.Bus.type==2))=0;
else
S.G.Inv=find(S.G.Type==l);
Ng=length(S.G.Bus);
Nd=length(S.G.Dis);
Nb=max(S.Bus.number);
S=calculate_YBus(S);
S.G.Dis=S.G.Bus(f ind(S.G.Kind==l));
S.G.UnDis=S.G.Bus(f ind(S.G.Kind==0));
NN=length(S.G.Bus)-1;
S.G.Kpl=S.G.Kp(S.G.Inv);
S.G.Ypl=S.G.Yp(S.G.Inv);
S.G.Kql=S.G.Kq(S.G.Inv);
S.G.Yql=S.G.Yq(S.G.Inv);
x0=[zeros(Nb,1);ones(Nb,1);ones(Nd,1);(S.G.Qmin+S.G.Qmax)/2;S.G.Kpl;S.G.Yp
A= [] ;
b= [] ;
Aeq= [] ;
beq= [] ;
ub=[pi*ones(Nb,1);1.l*ones(Nb,1);S.G.Pmax_Dis;S.G.Qmax;[]];
lb= [-pi*ones(Nb,1);0.9*ones(Nb,1);S.G.Pmin_Dis;S.G.Qmin;[]];
56


options=optimset('DerivativeCheck','on','Diagnostics','on','Display',...
'iter','FunValCheck','off','GradObj','off','MaxFunEvals',5000,...
'Maxiter,10000,'Hessian','off ');
[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(@funi,x0,A,b,Aeq,beq,1
Nl=length(S.G.Inv);
S.Bus.V=x(Nb+1:2*Nb);
S.Bus.ang=x(1:Nb);
S.G.Kp(S.G.Inv)=x(2*Nb+Ng+Nd+l:2*Nb+Nd+Ng+N1);
S.G.Yp(S.G.Inv)=x(2*Nb+Ng+Nd+Nl+l:2*Nb+Ng+Nd+Nl+Nl);
S.G.Kq(S.G.Inv)=x(2*(Nb+Nl)+Nd+Ng+1:2*(Nb+Nl)+Nd+Ng+N1);
S.G.Yq(S.G.Inv)=x(2*(Nb+Nl)+Nd+Ng+Nl+1:2*(Nb+Nl)+Nd+Ng+Nl+Nl);
S. Bus MWG= [] ;
S. Bus MVARG= [] ;
S.Bus.MWG(S.G.Dis)=x(2*Nb+l:2*Nb+Nd)';
S.Bus.MWG(S.G.UnDis)=S.G.P_ava_UnDis(find(S.G.Kind==0))';
S.Bus.MVARG(S.G.Bus)=x(2*Nb+Nd+l:2*Nb+Nd+Ng)';
S.f=60;
57


S.Bus.MWL(f ind(S.Bus.type==2))=0;
S.Bus.MVARL(find(S.Bus.type==2))=0;
end
S=Write_filecf(S);
B.4 Creating fun function
function opt=fun2(x,S)
/ opt=0;
% x=ones(100,1);
Ng=length(S.G.Bus1) ;
Nd=Ng-length(S.G.UnDis);
Nb=max(S.Bus.number);
Nc=Nb-l;
Nl=length(S.G.Invl);
p=100*x(Nc+Nb+l:Nc+Nb+Nd);
kp=x(Nb+Nc+Ng+Nd+2:Nb+Nc+Nd+Ng+Nl+1);
yp=x(N c+Nb+Ng+Nd+Nl+2:N c+Nb+2*Nl+Nd+Ng+l);
kq=x(Nb+Nc+Nd+Ng+Nl+Nl+2:Nb+Nc+Nd+Ng+Nl+Nl+l+Nl);
yq=x(Nb+Nc+2*Nl+Nd+Ng+Nl+2:Nb+Nc+2*N1+Nd+Ng+Nl+Nl+1);
sl=S.G.a(l:Nd);
s2=S.G.b(l:Nd);
58


s3=S.G.c(1:Nd);
c=sl.*(p.~2)+s2.*p+s3;
Fl=sum(c
w=l:2*Nl+l;
w(l)=0.6;
for i=2:5*N1+1;
w(i)=(l-w(l))/(4*N1);
end
/ w(2)=l-w(l);
op=0;
for ii=l:Nl;
op=op+w(4*11-2)*((kp(ii)-S.G.Kp(ii))~2)+w(4*ii-l)*((kq(ii)-S.G.Kq(ii))~2
end
opt=op+Fl;
B.5 creatin noncolon function
function [C, Ceq]=nonlcon(x,S)
Nd=length(S.G.Bus1)-length(S.G.UnDis);
Ng=Nd+length(S.G.UnDis);
Nb=max(S.Bus.number);
Vb=x(Nb+l:2*Nb);
59


ii=l:Nb;
Vb(ii)=abs(Vb(ii)).*exp(j*x(l:Nb));
f=x(2*Nb+Ng+Nd+l);
Nl=length(S.G.Invl);
Kp=x(2*Nb+Ng+Nd+2:2*Nb+Nd+Ng+Nl+l) ;
Yp=x(2*Nb+Ng+Nd+Nl+2:2*Nb+Ng+Nd+Nl+Nl+l);
Kq=x(2*(Nb+Nl)+Nd+Ng+2:2*(Nb+Nl)+Nd+Ng+Nl+1);
Yq=x(2*(Nb+Nl)+Nd+Ng+Nl+2:2*(Nb+Nl)+Nd+Ng+Nl+Nl+1);
/ f =x (2* (a+3*N)+1: 2* (a+3*N)+N) ;
f0=S.f*ones(N1,1);
fl=f*ones(N1,1);
R=S.Branch.Mvar;
SF=Vb(S.Branch.fr).*conj(((Vb(S.Branch.fr)-Vb(S.Branch.to))./...
(S.Branch.r+j *S.Branch.x*f/60))+Vb(S.Branch.fr).*(j *(S.Branch.b/2)*f/6C
ST=Vb(S.Branch.to).*conj(((Vb(S.Branch.to)-Vb(S.Branch.fr))./...
(S.Branch.r+j *S.Branch.x*f/60))-Vb(S.Branch.to).*(j *(S.Branch.b/2)*f/6C
C=[abs(SF)-R;abs(ST)-R] ;
C=full(C);
S.YB=sparse(diag(S.Bus.G+j*S.Bus.B*f/60)) ;
S.YB=S.YB+sparse(S.Branch.fr,S.Branch.fr,...
60


1./(S.Branch.r+j*S.Branch.x*f/60)+j*S.Branch.b/2*f/60,...
max(S.Bus.number),max(S.Bus.number));
S.YB=S.YB+sparse(S.Branch.to,S.Branch.to,...
1./(S.Branch.r+j*S.Branch.x*f/60)+j*S.Branch.b/2*f/60,...
max(S.Bus.number),max(S.Bus.number));
S.YB=S.YB+sparse(S.Branch.fr,S.Branch.to,...
-1./(S.Branch.r+j*S.Branch.x*f/60),max(S.Bus.number),max(S.Bus.number));
S.YB=S.YB+sparse(S.Branch.to,S.Branch.fr,...
-1./(S.Branch.r+j*S.Branch.x*f/60),max(S.Bus.number),max(S.Bus.number));
S.YB=full(S.YB);
Ib=S.YB*Vb;
Si=Vb.*conj(lb);
Pi=-S.Bus.MWL;
Qi=-S.Bus.MVARL;
Pi(S.G.Dis)=x(2*Nb+l:2*Nb+Nd);
Pi(S.G.UnDis)=S.G.P_ava_UnDis(f ind(S.G.Kindl==0));
Qi(S.G.Busl)=x(2*Nb+Nd+l:2*Nb+Nd+Ng);
f2=Kp.*Pi(S.G.Busl(find(S.G.Typel==l)))+Yp;
% Pr=(fl-Yp)./Kp;
*/. Qr= (abs (Vb (S Bus pv) ) -Yq) /Kq- (S Bus V (S. Bus pv) -S G. Yq) /S G. Kq+S Bus MV
/. Qr=(abs(Vb(S.G.Bus))-Yq) ./Kq;
V2=Kq. *Qi(S.G.Bus1(find(S.G.Typel==l)))+Yq;
61


Ceq=[Pi-real(Si);Qi-imag(Si);f2-f1;V2-abs(Vb(S.G.Busl(find(S.G.Typel==l))))
/ ;Pi(S.G. Bus) -Pr (S G. Bus) ;Qi(S.G. Bus) -Qr (S G. Bus)
Ceq=full(Ceq);
B.6 writing the file after optimization
function S=Write_filecf(S)
S.filename=uiputfile( *.cf ) ;
filecf=fopen(S.filename,'w');
fprintf(filecf,'%s\n',S.CaseName) ;
fprintf(filecf,'BUS DATA F0LL0WS\n');
for ii=l:length(S.Bus.number)
fprintf(filecf,'%4.Of ',S.Bus.number(ii));
fprintf(filecf,'%s',S.Bus.name(ii,:));
y.
/ fprintf (f ilecf ,1. Of ,S. Bus area(ii));
/ fprintf (f ilecf ,'% 1. Of ,S.Bus .zone(ii));
fprintf(filecf,1.Of',S.Bus.type(ii));
fprintf(filecf,5.3f ',S.Bus.V(ii));
fprintf(filecf,'% 6.3f ',S.Bus.ang(ii)*180/pi);
fprintf(filecf,7.3f',(S.Bus.MWL(ii))*100);
62


fprintf(filecf,7.3f',(S.Bus.MVARL(ii))*100);
fprintf(filecf,7,4f,(S.Bus.MWG(ii))*100);
fprintf(filecf,8.4f,(S.Bus.MVARG(ii))*100);
fprintf(filecf,'%2.Of,S.Bus.Bkv(ii));
fprintf(filecf,%2.Of',S.Bus.DV(ii));
fprintf (f ilecf ,0/02. Of ,S Bus .max(ii) ) ;
fprintf (f ilecf ,0/2. Of 1 ,S Bus .min(ii) ) ;
fprintf(filecf,2.Of',S.Bus.G(ii));
fprintf(filecf,2.If',S.Bus.B(ii));
fprintf(filecf,'\n') ;
end
fprintf(filecf,'-999\n') ;
fprintf(filecf,'BRANCH DATA F0LL0WS\n');
for ii=l:length(S.Branch.r);
fprintf(filecf,3.Of ',S.Branch.fr(ii));
fprintf(filecf,3.Of ',S.Branch.to(ii));
% fprintf (f ilecf, ,0/0 3.Of ,S.Branch.ar(ii)) ;
% fprintf(filecf,3.Of ,S.Branch.zone(ii));
% fprintf(filecf,3.Of ',S.Branch.type(ii));
% fprintf(filecf,3.Of ,S.Branch.circuit(ii));
fprintf(filecf,5.3f ',S.Branch.r(ii));
fprintf(filecf,% 5.3f ',S.Branch.x(ii));
fprintf(filecf,5.3f ',S.Branch.b(ii));
fprintf (f ilecf ,0/0 5.3f ,S.Branch.Mvar (ii)*100) ;
63


fprintf(filecf,\n);
end
fprintf(filecf,,-999\nJ);
fprintf(filecf,'GENERATOR DATA FOLLOWSW)
fprintf(filecf,(Bus kind Pmax Pmin Qmax Qmin a b c Kp Yp Kq Yq )\nJ)
for ii=l:length(S.G.Bus);
fprintf(filecf,4.Of',S.G.Bus(ii));
fprintf (f ilecf ,0/ 4. Of' S. G .Type (ii));
fprintf(filecf,% 4.Of,S.G.Kind(ii));
fprintf(filecf,'Z 7.2f ,(S.G.Pmax(ii))*100);
fprintf(filecf,7.2f ,(S.G.Pmin(ii))*100);
fprintf (f ilecf ,0/ 7.2f ', (S .G. Qmax(ii))*100) ;
fprintf(filecf,7.2f ,(S.G.Qmin(ii))*100);
fprintf (f ilecf, ,0/ 7.5f ,S.G.a(ii))
fprintf(filecf,% 6.3f ,S.G.b(ii))
fprintf(filecf,6.2f ,S.G.c(ii))
fprintf(filecf,7.4f ',S.G.Kp(ii));
fprintf (f ilecf ,5/ 7.4f ,S.G.Yp(ii)) ;
fprintf (f ilecf,7.4f ,S.G.Kq(ii));
fprintf(filecf,% 7.4f ,S.G.Yq(ii));
fprintf (f ilecf, ,0/ 5.2f ,S.G.P_ava_UnDis(ii)*100);
fprintf(f ilecf,\n);
end
fprintf(filecf,-99\n) ;
64


fprintf(filecf,FREQUENCY IN WHOLE SYSTEM\n');
fprintf(filecf,7.3f ,S.f);
fprintf(filecf,J\n);
fprintf(filecf,'-99\n');
fclose(filecf) ;
65


REFERENCES
[1] M.A. Abido. Environmental/economic power dispatch using multiobjec-
tive evolutionary algorithms: a comparative study. In Power Engineering
Society General Meeting, 2003, IEEE, volume 1, page 4 vol. 2666, july 2003.
[2] T. Ackermann, G. Andersson, and L. Soder. Electricity market regulations
and their impact on distributed generation. In Electric Utility Deregulation
and Restructuring and Power Technologies, 2000. Proceedings. DRPT 2000.
International Conference on, pages 608 -613, 2000.
[3] Vijay Vital Artur R.Bergen. power system analysis,. Department of Elec-
trical and Computer Engineering, pages 306-334, 2009.
[4] J. Carpentier. Contribution e letude do dispatching economique. Bull.
Soc.Franc. Elect., pages 431-447, 1962.
[5] J.M. Carrasco, L.G. Franquelo, J.T. Bialasiewicz, E. Galvan, R.C.P.
Guisado, Ma.A.M. Prats, J.I. Leon, and N. Moreno-Alfonso. Power-
electronic systems for the grid integration of renewable energy sources: A
survey. Industrial Electronics, IEEE Transactions on, 53(4):1002 -1016,
june 2006.
[6] D. Cvetkovic and I.C. Parmee. Genetic algorithm-based multi-objective
optimisation and conceptual engineering design. In Evolutionary Compu-
tation, 1999. CEC 99. Proceedings of the 1999 Congress on, volume 1, pages
3 vol. (xxxvii+2348), 1999.
[7] K.. De Brabandere, B.. Bolsens, J.. Van den Keybus, A.. Woyte, J.. Driesen,
and R.. Belmans. A voltage and frequency droop control method for parallel
inverters. Power Electronics, IEEE Transactions on, 22(4):1107 1115, july
2007.
[8] K. De Brabandere, K. Vanthournout, J. Driesen, G. Deconinck, and R. Bel-
mans. Control of microgrids. In Power Engineering Society General Meet-
ing, 2007. IEEE, pages 1 -7, june 2007.
[9] D. Gayme and U. Topcu. Optimal power flow with distributed energy
storage dynamics. In American Control Conference (ACC), 2011, pages
1536 -1542, 29 2011-july 1 2011.
66


[10] Jong hyun Ryu, Sujin Kim, and Hong Wan. Pareto front approximation
with adaptive weighted sum method in multiobjective simulation optimiza-
tion. In Simulation Conference (WSC), Proceedings of the 2009 Winter,
pages 623 -633, dec. 2009.
[11] Hyun-Koo Kang, Seon-Ju Ahn, and Seung-Il Moon. A new method to de-
termine the droop of inverter-based dgs. In Power Energy Society General
Meeting, 2009. PES 09. IEEE, pages 1 -6, july 2009.
[12] F. Katiraei and M.R. Iravani. Power management strategies for a microgrid
with multiple distributed generation units. Power Systems, IEEE Trans-
actions on, 21(4): 1821 -1831, nov. 2006.
[13] F. Katiraei, R. Iravani, N. Hatziargyriou, and A. Dimeas. Microgrids man-
agement. Power and Energy Magazine, IEEE, 6(3):54 -65, may-june 2008.
[14] B. Lasseter. Microgrids [distributed power generation]. In Power Engineer-
ing Society Winter Meeting, 2001. IEEE, volume 1, pages 146 -149 vol.l,
jan-1 feb 2001.
[15] Guanghong Liu, Gang Wu, Tao Zheng, and Qing Ling. Integrating pref-
erence based weighted sum into evolutionary multi-objective optimization.
In Natural Computation (ICNC), 2011 Seventh International Conference
on, volume 3, pages 1251 -1255, july 2011.
[16] A. Singhee and P. Castalino. Pareto sampling: Choosing the right weights
by derivative pursuit. In Design Automation Conference (DAC), 2010 47th
ACM/IEEE, pages 913 -916, june 2010.
[17] Xiaobo Tang and Guoqing Tang. Multi-objective planning for distributed
generation in distribution network. In Electric Utility Deregulation and
Restructuring and Power Technologies, 2008. DRPT 2008. Third Interna-
tional Conference on, pages 2664 -2667, april 2008.
[18] Zhutian Wang, Xinhong Huang, and Jin Jiang. Design and implementation
of a control system for a microgrid involving a fuel cell power module. In
Electrical Power Conference, 2007. EPC 2007. IEEE Canada, pages 207
-212, oct. 2007.
[19] Wu Zhe, Wan Yiru, He Chuan, Yu Jianhui, and Zhou Hao. Development
status of chinas renewable energy power generation. In Sustainable Power
Generation and Supply, 2009. SUPERGEN 09. International Conference
on, pages 1 -5, april 2009.
67


[20] Cao Yijia Zhou Xiao-xin. The future of china to develop non-hydro renew-
able sources power generation[j]. In Journal of Electric Power Science And
Technology, pages 2-7, 23(1): 2008.
68


Full Text

PAGE 1

OPTIMALPOWERFLOWINMICROGRIDbyMirsadraddinDarakhshanRokhsariB.S.E.EUniversityofMashad,Iran,1989AthesissubmittedtotheUniversityofColoradoDenverinpartialfulllmentoftherequirementsforthedegreeofMasterofScienceAppliedElectrical2012

PAGE 2

ThisthesisfortheMasterofSciencedegreebyMirsadraddinDarakhshanRokhsarihasbeenapprovedbyFernandoMancilla-DavidTitsaPapantoniDanConnorsMay10,2012

PAGE 3

DarakhshanRokhsari,Mirsadraddin(ElectricalEngineering)OptimalPowerFlowinMicrogridThesisdirectedbyProfessorFernandoMancilla-DavidABSTRACT Amicrogridwillplayanessentialroleinthefuture,sinceitwillreducethestressbytheincreasingelectricitydemands,ascomparedtotheconstraintsimposedbypowerplantsandtransmissionlines.Hence,identifyingtheoptimaloperatingpointofanydistributedgenerators(DGs)inamicrogridleadsthemtoidentifyingtheoptimaloperationofanygeneratorinthermofcostminimizationeciencymaximization. Thisthesisproposesanewalgorithmfordeterminationoftheoptimaloper-atingpointofdistributedgenerators.Inthisalgorithm,thenewdroopcharac-teristicslopes(frequencyvs.activepowerandvoltagevs.reactivepower)oftheinverterwillbedetectedoftheoptimizedpoint.Conventionally,inverter-basedDGsdroopcharacteristicslopeshaveaconstantslopeintheirmaximumandminimumcapacityinthedierentoperationmodeofmicrogrid(connecttoordisconnectfromthegrid).Ifanychangeisnecessary,theslopewillbechosenarbitrary.Tominimizethedeviationofdroopcharacteristicslopesandtheotherobjectives(denedbythedecisionmaker),themulti-objectiveoptimizational-gorithmcanbeused.Themulti-objectiveoptimizationisusedinthispaper,

PAGE 4

canprovidethereasonabledecisiontodeterminethenewdroopcharacteristicslopesinoptimizedpointwithminimumdeviation.

PAGE 5

Thisabstractaccuratelyrepresentsthecontentofthecandidate'sthesis.Irecommenditspublication.Approved:FernandoMancilla-David

PAGE 6

DEDICATIONThisthesisisdedicatedtomywifeandsonwhohavesupportedmeallthewaysincethebeginningofmystudies.Theyprovidemewithagreatsourceofmotivationandinspiration.Also,thisthesisisdedicatedtomymotherandmemoriesofmyfatherforinstillinginmetheimportanceofhardworkandhighereducation.Finally,thisthesisisdedicatedtoallthosewhobelieveintherichnessoflearning.

PAGE 7

ACKNOWLEDGMENTIwouldliketoexpressmydeepestgratitudetomyadvisor,ProfessorFernandoMancilla-David,forhisexcellentguidance,caring,patience,andprovidingmewithanexcellentatmospherefordoingresearch. IgreatlyappreciateProfessorsTitsaPapantoni,andProfessorDanConnorsforformingpartofmydissertationdefensecommittee.Finally,Iwouldliketothankmywifeandson;theywerealwaystherecheeringmeupandstoodbymethroughthegoodandbadtimes. Iwouldalsoliketothankmyparents.Theywerealwayssupportingmeandencouragingmewiththeirbestwishes. Finally,Iwouldliketothankmysisterinlaw,Sorayaforhersupporting.

PAGE 8

CONTENTS Figures....................................x Tables.....................................xi Chapter 1.Introduction................................1 2.Microgrid.................................3 2.1DistributedEnergyResourcesunit(DERs).............3 2.2Microgriddenitionandcharacteristics................3 2.2.1Powermanagementstrategiesformicrogrid............6 3.OptimalPowerFlow(OPF).......................8 3.1Introduction...............................8 3.2Formulation...............................8 4.DroopControl...............................13 4.1Conceptofdroopcontrol........................13 4.1.1Droopinmicrogrid..........................15 4.2LoadsharingthroughP-fcontrol...................17 5.OPFwithDroopControl.........................20 5.1Mathematicmodel...........................20 5.2Mathematicalformulaforoptimization................25 6.CaseofStudy...............................30 6.1Introduction...............................30viii

PAGE 9

6.2Microgridconnectedtothegrid....................34 6.3Microgriddisconnectedfromthegrid(autonomousmode)......35 6.3.1Sharingthepowerfromdroopcharacteristics...........35 6.3.2SharingthepowerwithOPF....................35 7.Conclusion.................................40Appendix A.\fmincon"inMATLAB..........................41 B.ThecodeofsoftwareinMATLAB....................43 B.1Readingthele.............................43 B.2CalculateY-Bus............................51 B.3Computationalcode..........................52 B.4Creatingfunfunction..........................58 B.5creatin"noncolon"function......................59 B.6writingtheleafteroptimization...................62 References ...................................66ix

PAGE 10

FIGURES Figure 2.1Microgridstructure..........................4 4.1Powerowingthroughalinebetweentwosources..........13 4.2Frequencydroopcharacteristic....................16 4.3Voltagedroopcharacteristic......................17 4.4FrequencydroopcharacteristicforthreeDGs............18 5.1FrequencydroopcharacteristicforDGinbusi...........21 5.2VoltagedroopcharacteristicforDGinbusi.............22 5.3FrequencydroopcharacteristicsforDGinbusiindierentload.25 5.4VoltagedroopcharacteristicsforDGinbusiindierentload...26 5.5DirectionoffrequencydroopcharacteristicschangingforDGinbusifordierentload...........................26 5.6DirectionofvoltagedroopcharacteristicschangingforDGinbusifordierentload............................27 6.1Thecaseofstudyinmicrogrid....................30 6.2Caseofstudywithallbuses......................31 6.3DroopcharacteristicdeviationforDG1...............38 6.4DroopcharacteristicdeviationforDG2...............38 6.5DroopcharacteristicdeviationforDG3...............39x

PAGE 11

TABLES Table 6.1Amountofimpedancesbetweenbusiandbusjinpu..........32 6.2Datafortypesofbus.........................33 6.3Dataforallofthegenerators.....................33 6.4Dataforallgeneratorsandbusesintheconnectedmodewiththegrid34 6.5Dataformicrogridinautonomousmodebeforeoptimization....36 6.6Thenormalizedofobjectivefunctionsinautonomousmodeafterop-timization................................37xi

PAGE 12

1.Introduction Thedemandforelectricityhasincreasedsteadilyoverthelastfewyearsduetorapidlypopulationgrowth.Thistrendwilllikelycontinueoverthenextfewdecades[18].Powergenerationinlargecentralizedpowerplantssuchasnuclear,fossilfuel(coal,gas),etc.hashadanimpactontheinducedcost.However,thecustomersfrequentlyresidetoofarfromthegeneratingpointofelectricity.Elec-tricpowersystemsandthepowergridarecurrentlyundergoingrestructuringduetoanumberoffactors.Thesefactorsinclude:increasingdemand,increasinguncertaintycausedbytheintegrationofintermittentrenewableenergysourcesandfurtherderegulationoftheindustry[5],[9]. Duetomanyeconomical,practicalandtechnicalreasons,thepowergridcannotprovidethenecessaryenergyfortheconsumer.Inprovidingacontingencytothepowershortageinmanyareasofdemand,ithasbecomenecessarytosearchforanotherwaytogeneratethepowerrequiredbytheconsumer.Anynewdecisiontoproduceenergyshouldberesearchedinparallelwithpresentgridsupplysystems.Inthepast,allofthepowerplantshavetransferredenergythroughtransmissionlines.However,restrictionsintransmissionlinesforcethenecessitiesofnetworkdesignwhichallowsgeneratorstobeplacedclosetotheirloads.Inconjunctionwiththeadvancesinpowerelectronics,thesmall-scaleddistributedpowergeneration(non50-60HZ)isinthepositiontochangetheconventionalpowersystem[11].Suchchangeincludesthedeploymentofrenewableenergysources(RES)suchaswind,solar,fuelcells,etc.1

PAGE 13

Theenergythatisgeneratedfromnaturalresourcessuchaswind,sunlight,rain,tidesandgeothermalheathasbeennamedrenewableenergy(RE),whereREunitisanenvironmentallyfriendlywayofgeneratingpowerandaseriouscontendertosubstitutingofnonrenewableenergieslikecoalandoil.Renewableenergycanbeusedforawiderangeofenvironmentsastheresourcesofsuchenvironmentsdictate.Suchenvironmentsincludevillagesaswellaslargeruralareas.Greenpowerreferstoasubsetofrenewableenergythatcharacterizedbyhighenvironmentalbenets. AlargenumberofresearchstudieshavebeenconductedinmanycountriesforsubstitutingconventionalenergywithRES.Forexample,Ithasbeenre-portedthatfrom1995to2000,theU.S.governmenthadspent1.456billiondollarsonrenewableenergyresearch[19],whileJapanistherisingstarofgreenpowerresearch.Byadjustingtheirgovernmentpoliciesandstrongstatenan-cialsupport,Japanhasmadefulluseofrenewableenergy.Ithasalsomadegreatachievementsinwind,solarandoceanenergypowergeneration[20].Fur-thermore,theutilizationofgreenpowerhasdisseminatedrapidlyinGermanyoverthepastfewyears,whereitisprojectedthatitwillgraduallyreplacethenuclearpowergenerationinthenextfewyears.[19].2

PAGE 14

2.Microgrid Thereductionoffossilfuelsources,pollutionoftheenvironment,constraintsintransmissionlines,andthepooreciencyofenergydrivesustoanewtypeofpowergeneration.Thesekindsofgeneratorsshouldbelocatedatthedistribu-tionvoltagelevelandnearresidentialdistricts.Thetypesofgeneratorsinthedistributionnetworkgeneratepowerbyusingnon-conventionalenergysources(RES).Thiskindofgeneratoriscalledadistributedgenerator(DG)anditsenergysourceiscalleddistributedenergyresources(DERs). 2.1DistributedEnergyResourcesunit(DERs) Distributedgeneratorsaresmall-rangepowergeneratorsthatareconnectedtothedistributionnetwork,ordirectlyatthecustomerapplication[2].TheDERsrefertoalltechnologieswhichcanbeusedtoprovideenergyclosetotheconsumer[8].ThehighpenetrationofDGsinthenetworkisabletoalleviatethepressureontheelectricsystembyfeedingintopartofthelocalload.DERunitscontaindistributedgenerator(DG)anddistributedstorage(DS)withdierentratesofow[13]. 2.2Microgriddenitionandcharacteristics In2001,ProfessorRobertH.LasseterfromtheUniversityofWisconsin-Madisonproposedthedenitionofthemicrogrid[14]andcreatedamicrogridtestbed.Bythisdenition,themicrogridhasasingleunitroleinapower3

PAGE 15

system.ItisdenedasthemicrogridactingasagroupofDERunitsinthesamewaythatdistributedvoltagecansupplyelectricityand/orheatagroupofdomesticloaddemand[18].Theproximityofthegeneratorandtheload(heatingandelectrical)inthemicrogriddecreasethecapitalcostandreducethelossesinthetransmissionlines.AtraditionalstructureofmicrogridisshowninFig.2.1.ThemicrogridinFig.2.1canoperateintwomodes: Figure2.1:Microgridstructure1. Inthismode,themicrogridisconnectedtothegridwherethegridandDERunitssharethegenerationofpowersucientforconsumersinthemicrogridsection.2. Inthesecondmodethemicrogridisdisconnectedfromthegridandactsas4

PAGE 16

anautonomousislandinthemicrogrid.InthismodealloftheDERunitsthatexistinthemicrogridprovidetheenergyfortheconsumer.Therearemanyreasonsforgoingwiththismode.Thereasonscanrangefromdroppedvoltageinthegrid,alowerpriceofenergyinthemicrogrid,intheeventoffaultinthegrid,andsoon.Inthiscasethemicrogridcanstayinstand-alonemodewithoutrequiringanypowerfromthegrid. DERunitsinmicrogridcanbeconnectedtothemicrogridintotwoways:1. Theyareinterfacedtothemicrogridthroughrotatingmachineswithre-specttoslowdynamicsandgreatinertia.2. Theyareconnectedthroughapowerelectronicconverter(PEC)withsmalloutputimpedanceandfastdynamicsinresponsetothechangingofpower. Sointhetransitionstagebetweenthetwodierentmodesofmicrogrid(stand-aloneandconnectedtothegrid),thetwoconnectiontypescannothavethesamedynamicsresponse.Inthissituation,inordertoincreaseresponsetimeforarotatingmachine,DScanbeusedforcomputationalpurposes,whereasinasteadystateconditiontherotatingmachinecanprovideitspower.Discussingdynamicsindetailisoutsidethescopeofthisproject. DERunitsandloadsinthemicrogridcanbeplacedeverywherebutnotnecessarilylocatedinthesamephysicallocation.AllofthetypeofDERunitscancontributeinthemicrogridstructure,suchasdieselgenerators,fuelcells,andallrenewableenergysources,assolarcells(SC),orwindpowergenerators,andmicroturbinesaswellasotherpossibletypes(i.e.,geothermal,biomass,etc.),asdepictedinFig.2.1.5

PAGE 17

Fig.2.1showstheexistenceofaDSthatprovidesanuninterruptiblesupplyofelectricitytothemicrogrid.CurrentlymanytypesofDS,suchasDCbat-tery,Flywheel,Supercapacitor(SC),SuperconductionMagneticEnergyStorage(SMES),etc.areavailableforthemicrogrid. Themicrogridunithasanelectricalconnectionpointtothegridthatisreferredtoasthepointofcommoncoupling(PCC). 2.2.1Powermanagementstrategiesformicrogrid TheobjectiveofthisprojectistocomputetheamountofpowergeneratedbyanyDERunitinthemicrogrid.Toachievethisobjective,apowermanagementstrategy(PMS)isrequired.Inthemicrogrid,everyDERunithasadierentrateofowwithnodominantDGinthestand-alonemodeofmicrogrid[12].Therefore,thePMSmustconsideralltypesofDERunitsandcreateacontrolmethodologyforenergyowinthemicrogridinbothmodes(connectedtoanddisconnectedfromthegrid).Inthiscontrol,thePMSshouldsupervisetheactiveandreactivepowerowwithhighstabilityandreliabilityinthesystem. Thestrategyforcontrolshouldbetoprovidethepeer-to-peer(P2P)andplug-and-play(PnP)conceptsinthemicrogrid. Inthecontextofthemicrogrid,P2PmeansthereisnomastercontrolforanyDERunit,soifoneDERunitislost,thesystemcancontinueitsoperationwithanextraDERunit.PnPmeans,anyDERunitscanbeaddedinanylocationofthemicrogridwithoutredesigningnewcontrolstrategiesforthesystem. PMShasanotherresponsibilitywhichistoprovidethedecisioninthedy-namicsintransitionbetweenthetwomodesofthemicrogrid,butthisisoutside6

PAGE 18

thescopeofthisproject.7

PAGE 19

3.OptimalPowerFlow(OPF) 3.1Introduction OPFwaspresentedforthersttimebyCarpentierin1962[4].OPFisanalgorithmthatoptimizestheobjectivefunction,suchasthetotalcostofpowergeneration,magnitudeandphaseangledeviationofbusvoltage,powerlosses,andemissionofpollutants,etc.Theobjectivefunctioncanalsobepredenedbytheuser.OPFndstheoptimalsetpointoftheelectricalsystemwhichsatisesthesystempowerowequationsandallofthelimit,constraints,andcontingencies,associatedwiththepowernetwork. 3.2Formulation InapowergridwithNbusses,ifthevoltagemagnitudeandthephaseangleofthebussesarecomputedbyutilizingtheOPFformulathenallofthevariablesassociatedwiththesystem,suchastheactiveandreactivepowertransferredbetweenthebusesandalsotheamountofthecurrentinthelines,willbesatised.Inthepowersystem,therearethreetypesofbusses:1. ThePVbusorgeneratorbus,withatleastonegeneratorconnectedtothebus.2. ThePQbusorloadbus,wherenogeneratorisconnectedtothebus.3. Theslackbusorswingbus,thisisdenedbythereferenceforphaseangleofthevoltageandcanprovidethenecessaryactiveandreactivepowerforbalancingthesystemattheoptimizedsolution.8

PAGE 20

Ifweassumethepowersystemincludesnbussesandnggenerators,thedimensionoftheadmittancematrixYbuswillbea\nn"matrixwhichcanbedenotedas:Ybus=0BBBBBBBBBBB@y11y12::y1ny21y22::y2n::::::::::yn1yn2::ynn1CCCCCCCCCCCA(3.1)wheretheelementsofyi;j,i;j=1;2;:::;nofYbusareconstructedasdenedin[3]asfollows: TheYbusissymmetric. ThediagonalelementsYiiofmatrixareequaltothesumoftheadmittancesofallthecomponentsconnectedtotheithnode. Theo-diagonalelementsYijareequaltothenegativeofprimitiveadmit-tancesofallcomponentsconnectedbetweennodesiandj. ForthissystemthecomplexpowerinjectedmatrixSbusisdenedas:Sbus=VbusIbus(3.2) wheretheVbusisthevectorofthebusvoltageandIbusisavectoranddenotesthebuscurrentinjectedinthebus.Manipulatingtheequations,theresultswillbeasfollows:withIbus=YbusVbus(3.3)9

PAGE 21

Sbus=Vbus(YbusVbus)(3.4)Sbusi=Vbusi"nXj=1(yijVbusj)#i;j=1;2;:::;n(3.5)Sbusi=nXj=1(VbusiVbusjyij)i;j=1;2;:::;n(3.6)Vi=jVijejiVj=jVjjejji;j=1;2;:::;n(3.7)Yij=Gij+jBijij=i)]TJ /F3 11.955 Tf 11.96 0 Td[(ji;j=1;2;:::;n(3.8)Sbusi=nXj=1jVijjVjjejij(Gij)]TJ /F3 11.955 Tf 11.95 0 Td[(jBij)i;j=1;2;:::;n(3.9) where:Sbusi:Thecomplexpowerinjectioninbusi.Sbusi=nXj=1[jVijjVjj(cosij+jsinij)(Gij)]TJ /F3 11.955 Tf 11.95 0 Td[(jBij)]i;j=1;2;:::;n(3.10)Pbusi=nXj=1[jVijjVjj(Gijcosij+Bijsinij)]i;j=1;2;:::;n(3.11)Qbusi=nXj=1[jVijjVjj(Gijsinij)]TJ /F3 11.955 Tf 11.95 0 Td[(Bijcosij)]i;j=1;2;:::;n(3.12)where:Pbusi:Theactivepowerinjectioninbusi.Qbusi:Thereactivepowerinjectioninbusi. Inthepowerowproblem,theloadisaknownvariable.Insettinguptheformula,theactiveandreactivepowerinjectiontoeverybusshouldconsidertheamountoftheloadandgenerationasfollows:Pbusi=Pgi)]TJ /F3 11.955 Tf 11.96 0 Td[(PLi=f(V;)Qbusi=Qgi)]TJ /F3 11.955 Tf 11.95 0 Td[(QLi=f(V;)(3.13)10

PAGE 22

where:Pgi:Theactivepowerofthegeneratorinthebusi.Qgi:Thereactivepowerofthegeneratorinthebusi.PLi:Theactivepoweroftheloadinthebusi.QLi::Thereactivepoweroftheloadinthebusi. Accordingtotheoptionalobjectivefunctionanditsconstraints,alternativeformulascanbepresented.TheOPFwillsolvequadraticproblemincludemini-mizingaquadraticfunctionsubjecttolinearandnonlinearequality,inequality,thatisofinterestinthispaper. Theformulainwhichtheobjectivefunctionminimizesthecost(Ci(P),i=1;2;:::;ng)or/andemissions(Totemmissioni,i=1;2;:::;ng)or/andtotallossesbetweenthelines(Totlosses)inthesystemor/andanyotherobjectivefunctioncanbesetupasfollows:Minimize:ngXi=1Ci(P)i=1;2;:::;ng(3.14) or/andMinimize:ngXi=1Totemmissionii=1;2;:::;ng(3.15) or/andMinimize:Totlosses(3.16) or/and...subjectto:inequalityconstraints:11

PAGE 23

minimaxi=1;2;:::;n)]TJ /F1 11.955 Tf 11.95 0 Td[(1(3.17)VminjVijVmaxi=1;2;:::;n(3.18)PminiPgiPmaxii=1;2;:::;ng(3.19)QminiQgiQmaxii=1;2;:::;ng(3.20)jSijFromjRatei;j=1;2;:::;n;i6=j(3.21)jSijTojRatei;j=1;2;:::;n;i6=j(3.22) where:jSijFromjandjSijTojarethemagnitudeofcomplexpowerwhichtransferredfrombus(i)tobus(j)atthesendingandreceivingpoint.equalityconstraints:Pbusi=Pgi)]TJ /F3 11.955 Tf 11.96 0 Td[(PLi=nXj=1[jVijjVjj(Gijcosij+Bijsinij)](3.23)Qbusi=Qgi)]TJ /F3 11.955 Tf 11.96 0 Td[(QLi=nXj=1[jVijjVjj(Gijsinij)]TJ /F3 11.955 Tf 11.96 0 Td[(Bijcosij)](3.24)Fromequations(3.14)to(3.24),theusercanndthevariablesfortheobjectivefunction.12

PAGE 24

4.DroopControl Inthepowersystem,activepowerandfrequencyaredependentoneachother.Inarotatingmachine,increasingtheloadwillresultinanincreaseintorque.Thetorqueisrelatedtotherotationalspeedwithrelationtothefrequencyofthesystem.Thismeansthatthelowerthefrequency,thehighertheloadwillbe.ThereforePMScantrytodetermineawaytocontrolthesystemwiththesecharacteristicswhentheDERunitinmicrogridisinterfacedtothesystemwithaninverter. 4.1Conceptofdroopcontrol Tounderstandtheconceptofdroopcontrol,wereturntothepowersystemfordeliveringcomplexpowerthroughatransmissionlinebetweentwobuses. Figure4.1:PowerowingthroughalinebetweentwosourcesFig.4.1showsthemodelofthetransmissionlinebetweenbusiandbusjwith13

PAGE 25

impedanceZwhere:Z=R+jXL(4.1) andthepowerowfromthebusitobusjisdenedas:Sij=Pij+jQij(4.2)Sij=ViIij(4.3)Sij=Vi(Vi)]TJ /F3 11.955 Tf 11.96 0 Td[(Vj) Z(4.4) byusingthisdenition:ij=i)]TJ /F3 11.955 Tf 11.96 0 Td[(j(4.5) Sij=jVij2)-222(jVijjVjjeij Z(4.6) withEuler'sformula: eij=cosij+jsinij(4.7) wecanseparatetheactiveandreactivepowertothisform: Pij=jVij XL2+R2(RjVij)-223(jVjjcosij+jVjjXLsinij)(4.8)Qij=jVij XL2+R2(XLjVij)-222(jVjjcosij)-222(jVjjRsinij)(4.9) Ifinatransmissionlinetheresistanceofthelineismuchsmallerthantheinductance,thenequation(4.8)andequation(4.9)willbechangedtothisform:Pij=jVijjVjj XLsinij(4.10)14

PAGE 26

andQij=jVij2 XL)]TJ 13.16 8.09 Td[(jVijjVjj XLcosij(4.11) Furthermore,ifthedierenceinangleijisverysmall,incomputingtheradian,cosij=1andsinij=ij.Therefore:sinij=ij=XLPij jVijjVjj(4.12) andjVij(jVij)-222(jVjj)=XLQij(4.13) Fromequation(4.12)thedierenceofthephaseanglebetweentwobusses(ij)isdependentontheactivepowerdeliveredbetweenthebuses,andfromequation(4.13)themagnitudevoltages(jVijandjVij)intwobusesisinuencedbyreactivepower.Thisimpliesthatbychangingij,theactivepower,andchangingthejVijandjVij,thereactivepowercanbecontrolled. 4.1.1Droopinmicrogrid Inrelationtodroopcontrolinthemicrogrid,thefrequencyisusedinsteadofthephaseangleforcontrollingtheactivepower.ThemicrogridincludestwoormoreDGswhichusedroopcharacteristicstosharetheactiveandreactivepower,TheseDGsareinterfacedwiththeinverter.Whentheloadinthemicrogridischanged,PMSmusthavespecicrulestodeterminethesetpointofeveryDGforthesharingofactiveandreactivepowerattheoptimizedpoint. Fig.4.2showsthefrequencydroopcharacteristic.Inthisgurethef0isthenominalfrequencyandP0isthemomentarysetpointforactivepowerof15

PAGE 27

Figure4.2:Frequencydroopcharacteristictheinverter[7].FromtheFig.4.2,thisequationcanbewrithenasfollows:f1=f0+Kp(P1)]TJ /F3 11.955 Tf 11.95 0 Td[(P0)(4.14) WhereKpistheslopeofthecharacteristicandf1isthethefrequencyoftheinverterforgeneratingthepowerP1fromfrequencydroopcharacteristics. InFig.4.3,thevoltagedroopcharacteristicisshown. FromFig.4.3,theequationcanbewrittenasfollows:jV1j=jV0j+Kq(Q1)]TJ /F3 11.955 Tf 11.96 0 Td[(Q0)(4.15) WhereV0isthenominalvoltageandQ0isthereactivepoweroftheinverter.Fromequations(4.14)and(4.15)byadjustingPandQindependently,thefrequencyandmagnitudeofthevoltagecanbedetermined[7]. 4.2LoadsharingthroughP-fcontrol16

PAGE 28

Figure4.3:Voltagedroopcharacteristic InamicrogridwithtwoormoreDGs,changingtheloadinbothmodes(connectdtothegridordisconnectedfromgrid),PMSusesthefrequencydroopcharacteristicstochangetheset-upoftheoperatingpointtoaccessthelocalpowerbalanceinthenewsituationofloading.Asanexample,amicrogridwhichhasthreeDGswithdroopcharacteristicsisshowninFig.4.4. EveryDGhasadierentslopeinfrequencydroopcharacteristics(Kp1,Kp2andKp3). Whenthemicrogridisconnectedtothegrid,themaingridislargeandro-bustenoughtocompensateanyvariationinthemicrogrid[11].Sothefrequencyinthemicrogridwillbethenominalfrequencyofthegrid(f0).FromFig.4.4andtheirdroopcharacteristics,thethreeDGsprovidetheloadPL0asfollows:PL0=P10+P20+P30(4.16)17

PAGE 29

Figure4.4:FrequencydroopcharacteristicforthreeDGs ThePL0equals:PL0=PL)]TJ /F3 11.955 Tf 11.96 0 Td[(Pgrid(4.17) ThePgriddenotestheamountofloadinthemicrogridsharedbythegrid,andthePListhetotalloadsinthemicrogrid.Nowifthemicrogridforanyreasonisdisconnectedfromthegrid,withthesameloadPL,thefrequencyofthesystemwillbechangedtof1andthepowershareforeveryDGwillbeP11,P21andP31,wherePLequals: PL1=P11+P21+P31(4.18) Inthiscase:PL1=PL(4.19) Iftheequation(4.14)isdevelopedforeveryDG,thisequationcanbewrittenasfollows:f1=f0+Kpi(Pi1)]TJ /F3 11.955 Tf 11.95 0 Td[(Pi0)(4.20)18

PAGE 30

with:PL1=3Xi=1Pi1PL0=3Xi=1Pi0PL=PL1)]TJ /F3 11.955 Tf 11.96 0 Td[(PL0(4.21) Usingequation(4.16),equation(4.18),andequation(4.21)theamountofpowerforeveryDGwillbe:Pi1=PL1 Kpi 3Xi=11 kpi+Pi0(4.22)19

PAGE 31

5.OPFwithDroopControl Anexperimentalmicrogridbasedondistributedgenerators(DG)consistsofseveralmodules:theutilitypowergrid,invertersandloads.Thischapterwilladdressinvertermanagementandcontrol.EveryinverterhasautonomousdroopcontrolcharacteristicsforsharingtheactiveandreactivepoweramongstDGstomeettheload. 5.1Mathematicmodel Chapter4explainedtheloadsharingprocessthroughthep-fcontrollerandequation(4.22)showstheamountofpowerforeveryDG.Inthiscalculation,thetotalofthepowerisequaltotheloads,buttherealitythesystemisnonlinearandthelossesshouldbeconsideredastheequationwillbe:ngXi=1Pgi=PLoad+Losses(5.1)TheLossesarenotconstantandaredependentonthevoltageofthebuses.IntheautonomousmicrogridmodetheLossesvarywithfrequencybecausetheimpedanceisdependentonthefrequency. Toconnectmanyinvertersparallelinthemicrogrid,thefrequencyofthewholesystemshouldbebymakingthemthesame.Thefrequencyvs.activepowercharacteristicsfortheDGinbusiplotisshownintheFig.5.1. AccordingtotheFig.5.1thefollowingequationforFrequency-droopchar-acteristiccanbewritten:20

PAGE 32

Figure5.1:FrequencydroopcharacteristicforDGinbusif=KpiPgi+bpi(5.2)Inthisequation,theamountofpowerisdependentonthefrequency,sowiththeconstantslopeinthelowerfrequencytheoutputpowerislarger.Instandalonemodewherethegridisdisconnected,everyDGshouldproducemorepowerinthesamelowerfrequency. Theplottingofthevoltagevs.reactivepowerisshowninFig.5.2.Andequation(5.3)showstherelationshipofmagnitudeofthevoltagewithrespecttothereactivepower.jVgij=KqiQgi+bqi(5.3) WherefandjVgijarethefrequencyandamplitudeofthevoltageoftheDGinbusirespectively.KpiandKqiaredenedbythecorrespondingslopesofdroop.21

PAGE 33

Figure5.2:VoltagedroopcharacteristicforDGinbusi Whenthemicrogridisconnectedtothegrid,thefrequencyofthesystemwillbethesameastheutilitygrid.TosharetheactiveandreactivepowerbetweentheDGsattheoptimizedpoint,theslopesKpiandKqicanbeadjusted.Hence,theimpedancebetweentwobusesiandjcanbewrittenas:Zij=Rij+2fLijwhere:f=60(5.4) Inautonomousmicrogridmodetheloadwillbedierent,andinordertogetexcellentpowersharing,thefrequencyandmagnitudeoftheDGvoltagewillbechanged.So,thechangeofimpedancebetweentwobusesiandjcausedbychangingthefrequencycanbeexpressedasfollows:Xij=260LijandZij(f)=Rij+jXijf 60(5.5) Foranetworkwithnindependentbusses,andnggenerators,thedimensionofadmittancematrixwillbea\nn"andthismatrixwillbeafunctionoff22

PAGE 34

asfollows(seeChapter3andequation5.5):Ybus(f)=0BBBBBBBBBBB@y11(f)y12(f)::y1n(f)y21(f)y22(f)::y2n(f)::::::::::yn1(f)yn2(f)::ynn(f)1CCCCCCCCCCCA(5.6) forthisnetworkwithnbusesthefollowingequationscanbewritten: Sbus=VbusIbus(5.7) whereIbusisthebuscurrentinjectionvector.Sbus=Vbus[Ybus(f)Vbus](5.8)Sbusi=Vbusi"nXj=1[yij(f)Vbusj]#i;j=1;2;:::;n(5.9)Sbusi=nXj=1[VbusiVbusjyij(f)]i;j=1;2;:::;n(5.10)Vi=jVijejiVj=jVjjejji;j=1;2;:::;n(5.11)Yij(f)=Gij(f)+jBij(f)ij=i)]TJ /F3 11.955 Tf 11.95 0 Td[(ji;j=1;2;:::;n(5.12)Sbusi=nXj=1jVijjVjjejij[Gij(f))]TJ /F3 11.955 Tf 11.95 0 Td[(jBij(f)]i;j=1;2;:::;n(5.13) where:Sbusi:Thecomplexpowerinjectioninbusi.Sbusi=nXj=1[jVijjVjj(cosij+jsinij)(Gij(f))]TJ /F3 11.955 Tf 11.96 0 Td[(jBij(f))]i;j=1;2;:::;n(5.14)23

PAGE 35

Pbusi=nXj=1[jVijjVjj(Gij(f)cosij+Bij(f)sinij)]i;j=1;2;:::;n(5.15)Qbusi=nXj=1[jVijjVjj(Gij(f)sinij)]TJ /F3 11.955 Tf 11.95 0 Td[(Bij(f)cosij)]i;j=1;2;:::;n(5.16) Inthepowerowproblem,theloaddemandsareknownvariables,sowecandenethefollowingbuspowerinjectionsas: Pbusi=Pgi)]TJ /F3 11.955 Tf 11.96 0 Td[(PLi=f(V;;f)Qbusi=Qgi)]TJ /F3 11.955 Tf 11.95 0 Td[(QLi=f(V;;f)(5.17) where:Pgi:TheactivepoweroftheDGinthebusi.Qgi:ThereactivepoweroftheDGinthebusi.PLi:Theactivepoweroftheloadinthebusi.QLi:Thereactivepoweroftheloadinthebusi. Fromequation(5.2)andequation(5.3)foreveryDG,theoutputpower(ac-tiveandreactive)canbecomputedcorrespondingtotheirdroopcharacteristics:Pgi=f)]TJ /F3 11.955 Tf 11.96 0 Td[(bpi KpiQgi=jVgij)]TJ /F3 11.955 Tf 17.93 0 Td[(bqi Kqi(5.18) Inoptimalpowerowthemainobjectivefunctionisthetotalgenerationcostwhichisdenedasfollow: Totalcost=nXi=1Ci(P)(5.19) where:Ci(P)isthecostfunctionfortheDGinbusiWhenthegridisconnectedto24

PAGE 36

themicrogrid,alloftheDGsprovidetheirenergyfromanoptimalpointatthefrequencyofgridinKpi0andKqi0.Thosearedenedbythecorrespondingslopesoftheirdroops.Duringstand-alonemicrogridmodethegridisdisconnectedfromthemicrogridandthepowerdeliveryfromeachDGwillbechanged.Inthisvariationtheamountofpowerandreactivepowersharingwillbechangedaccordingtoequation(5.2)andequation(5.3)inthenewfrequencyf.ThechangeisshowninFig.5.3andFig.5.4. Figure5.3:FrequencydroopcharacteristicsforDGinbusiindierentload However,thesepointsarenotoptimizedpoints.Aftercalculatingtheopti-mizedpoint,theslopeandY-intersectofthedroopwillbechangedtothenewvalueKpi,Kqi,bpiandbqi.Fig.5.5andFig.5.6,depictthedirectionofchangesforbothdroops. 5.2Mathematicalformulaforoptimization Tosetuptheequationsinanoptimalpowerow,otherobjectivescanbeaddedtotheoptimizationproblemtominimizethedeviationofdroopcharac-teristics.ThismeanstheminimizingofKpi,Kqi,bpiandbqi.25

PAGE 37

Figure5.4:VoltagedroopcharacteristicsforDGinbusiindierentload Figure5.5:DirectionoffrequencydroopcharacteristicschangingforDGinbusifordierentload Theequationcanbewrittenasfollows:Minimize:Totalcost=nXi=1Ci(P)i=1;2;:::;ng(5.20) and26

PAGE 38

Figure5.6:DirectionofvoltagedroopcharacteristicschangingforDGinbusifordierentload Minimize:(K2pi;b2pi;K2qi;b2qi)i=1;2;:::;ng(5.21) Suchthat: minimaxi=1;2;:::;n(5.22)VminjVijVmaxi=1;2;:::;n(5.23)PminiPgiPmaxii=1;2;:::;ng(5.24)QminiQgiQmaxii=1;2;:::;ng(5.25)fminffmax(5.26)KpminKpiKpmaxi=1;2;:::;ng(5.27)27

PAGE 39

KqminKqiKqmaxi=1;2;:::;ng(5.28)bpminbpibpmaxi=1;2;:::;ng(5.29)bqminbqibqmaxi=1;2;:::;ng(5.30) jSijFromjRatei;j=1;2;:::;n;i6=j(5.31)jSijTojRatei;j=1;2;:::;n;i6=j(5.32)Pbusi=Pgi)]TJ /F3 11.955 Tf 11.96 0 Td[(PLi=nXj=1[jVijjVjj(Gij(f)cosij+Bij(f)sinij)]i=1;2;:::;n(5.33)Qbusi=Qgi)]TJ /F3 11.955 Tf 10.06 0 Td[(QLi=nXj=1[jVijjVjj(Gij(f)sinij)]TJ /F3 11.955 Tf 11.95 0 Td[(Bij(f)cosij)]i=1;2;:::;n(5.34)Pgi=f)]TJ /F3 11.955 Tf 11.96 0 Td[(bpi Kpii=1;2;:::;ng(5.35)Qgi=jVgij)]TJ /F3 11.955 Tf 17.93 0 Td[(bqi Kqii=1;2;:::;ng(5.36)Kpi=Kpi)]TJ /F3 11.955 Tf 11.96 0 Td[(Kpi0i=1;2;:::;ng(5.37)Kqi=Kqi)]TJ /F3 11.955 Tf 11.96 0 Td[(Kqi0i=1;2;:::;ng(5.38)bpi=bpi)]TJ /F3 11.955 Tf 11.95 0 Td[(bpi0i=1;2;:::;ng(5.39)bqi=bqi)]TJ /F3 11.955 Tf 11.96 0 Td[(bqi0i=1;2;:::;ng(5.40) Theobjectivefunctionscanbeevaluatedbyusingtheweightedsummethodthatexplainedin[10],[15],[17],[16],[6].Theweightedobjectivefunctionforthesystemcanbeexpressedasfollows: F=W1F1+W2F2+:::+WiFii=1;2;:::;4ng+1(5.41)28

PAGE 40

where:F1=ngXi=1Ci(P),F2=K2p1,F3=K2q1,F4=b2p1,F5=b2q1......F4i)]TJ /F8 7.97 Tf 6.58 0 Td[(2=K2pi,F4i)]TJ /F8 7.97 Tf 6.59 0 Td[(1=K2qiF4i=b2pi,F4i+1=b2qiand,W1+W2+W3+:::+W4i+1=1(5.42)29

PAGE 41

6.CaseofStudy 6.1Introduction ThecasethatwepickedupforthisstudyisshownintheFig.6.1[12].Fig. Figure6.1:Thecaseofstudyinmicrogrid6.1isasingle-diagramofa13.8-KVdistributionsystemusedtoinvestigateapossiblemicrogridPMS.ThissystemhasthreeDGunits,withdierentrating,ie.,DG1(1.8-MVA),DG2(2.5-MVA)andDG3(1.5-MVA).AlloftheDGunits30

PAGE 42

canbedispatched[12].Tostudythiscase,allofthebusesfortheloadsandDGsareshowninFig.6.2.AsseenasinFig.6.2,thesystemhasfourteenbuses, Figure6.2:CaseofstudywithallbusesandthreeDGs:DG1,DG2andDG3connectrespectivelytoBus6,Bus13,andBus12.ThevaluesoftheimpedancebetweenthelinesarestatedinTable6.1. ThismicrogridisconnectedtothegridthroughBus11.InTable.6.2thetypeofthebuseswiththeamountsoftheloads(activeandreactivepower)arelisted.Inthecolumn\Typeofbus"ofthistable,the\No.0"meansthatthereisnogeneratorconnectedtothebus,the\No.2"meanstheDGisconnectedtothebusandthe\No.3"meansthatthebusisconnectedtothegrid.31

PAGE 43

Table6.1:Amountofimpedancesbetweenbusiandbusjinpu. Busi Busj Rij XLij 1 7 0.5314 5.7635 1 5 0.4022 0.2395 1 6 0.21 1.094 1 14 0.3976 0.5127 2 14 0.6141 0.3066 3 14 0.6065 1.015 3 8 0.816 4.8332 3 9 0.6903 0.000 14 11 0.0818 0.5826 12 4 0.822 0.48332 13 3 0.822 0.48332 4 14 0.3564 0.2661 10 4 0.822 4.8332 ThedataforallofthegeneratorsarelistedintheTable6.3.Thegeneratorcostcurvesarepresentedbyquadraticfunctions[1].Theparam-etersa,b,andcinTable.6.3,arethecoecientofthequadraticfunctionandareexpressedinthefollowingform:Ci=aiP2gi+biPgi+bii=1;2;:::;ng(6.1) Somerenewableenergysourcessuchassolarandwindpowerdonothaveanyfuelcostsbecausetheresourcesarefree,sointheOPFequation,thistypeofdistributedgeneratorresourceunitshouldnotbeconsideredfordispatch.For32

PAGE 44

Table6.2:Datafortypesofbus BusNo. Typeofbus PL(MW) QL(MVAR) Bus1 0 0 0 Bus2 0 0.6 0.42 Bus3 0 0 0 Bus4 0 0 0 Bus5 0 0.61 0.42 Bus6 2 0 0 Bus7 0 0.600 0.390 Bus8 0 0.7 0.45 Bus9 0 0.8 0.5 Bus10 0 0.9 0.61 Bus11 3 0 0 Bus12 2 0 0 Bus13 2 0 0 Bus14 2 0 0 Table6.3:Dataforallofthegenerators DGs:Name BusNo: Type Pmax(Mws) Pmin(MW) Qmax(MVAR) Qmin(MVAR) a b c Kp bp Kq bq DG1 6 1 1.80 0.000 1.80 -1.80 0.00482 7.970 78.00 -0.0282 60.0137 0.0056 1.0306 DG2 13 1 2.50 0.000 2.50 -2.50 0.00156 7.920 561.00 -0.0536 60.0453 -0.1488 1.0183 DG3 12 1 1.50 0.000 1.50 -1.50 0.00194 7.850 310.00.0 -0.2537 60.2417 -0.2157 1.0255 Grid 11 1 1000 0.000 1000 -1000 0.00139 7.060 500.00 codingintheTable6.3,ifthegeneratorisabletodispatchthetypethenitis(1),otherwiseitis(0).DroopcharacteristicparametersarestatedinthelastfourcolumnsofTable6.3(Kp,bp,Kqandbq).TheloadsL1inBus7,L2inBus5,L3inBus2,L4inBus8,L5inBus9andL6inBus10arestatedandthe33

PAGE 45

valueoftheloadsareexpressedinTable6.2. 6.2Microgridconnectedtothegrid Inthiscase,thegridisconnectedtothemicrogrid,andthepowernecessaryfortheloadsissharedbyDGsandthegrid.Table6.4showsthevaluesofthemagnitude,thephaseangleofthevoltageofeachbus,andtheamountofthepowerandreactivepowersharedbyanygenerators,aswellasparametersofdroopcharacteristicsforeveryDGattheoptimizedpoint. Table6.4:Dataforallgeneratorsandbusesintheconnectedmodewiththegrid DGs:Name BusNo: jVjpu ]Vdeg PL(MW) QL(MVAR) Pg(MW) Qg(MVAR) Kp bp Kq bq Bus1 1.062 -0.804 0 0 0 0 Bus2 1.058 -0.537 0.600 0.420 0 0 Bus3 1.058 -0.887 0 0 0 0 Bus4 1.063 -0.585 0 0 0 0 Bus5 1.059 -0.792 0.61 0.420 0 0 Bus6 DG1 1.075 -0.610 0 0 0.5753 1.1608 -0.0293 60.0169 0.0656 0.9987 Bus7 1.037 -2.496 0.600 0.390 0 0 Bus8 1.031 -2.472 0.700 0.450 0 0 Bus9 1.033 -2.319 0.800 0.500 0 0 Bus10 1.026 -2.607 0.900 0.610 0 0 Bus11 Grid 1.066 0.000 0 0 1.9840 0.2543 Bus12 DG3 1.073 -0.660 0 0 0.8708 0.6927 -1.0648 60.9272 0.3767 0.8116 Bus13 DG2 1.068 -1.050 0 0 0.8375 0.8847 -0.1078 60.0903 0.6452 0.4973 Bus14 1.063 -0.575 0 0 0 0 34

PAGE 46

6.3Microgriddisconnectedfromthegrid(autonomousmode) Ifthegridisdisconnectedfromthemicrogrid,itwillnotcontributetosharinganygeneratedpowerwiththeDGs.Inthissituation,alloftheDGsinthemicrogridshouldprovidethenecessarypowerrequiredfortheloadwithoutthegrid. 6.3.1Sharingthepowerfromdroopcharacteristics Initially,whenthegridisdisconnectedfromthemicrogrid,thedroopchar-acteristicsofallinverterswillnotbechanged.NewamountsofpowergeneratedforeveryDGarecomputedbyequation4.22inthesamefrequency.InTable.6.5,theamountsofactiveandreactivepoweraregeneratedbyanyDGwithoutanychangeintheparametersofdroopcharacteristic. 6.3.2SharingthepowerwithOPF TheamountofpowerforeveryDGinprevioussectionwillnotbeatop-timizedpoint.Accordingtothealgorithmexpressedinchapter5,themulti-objectivefunctionoptimizationwithweightedparametercanbeusedasfollows: F=13Xi=1WiFi(6.2) where:F1=nXi=1Ci(P),35

PAGE 47

Table6.5:Dataformicrogridinautonomousmodebeforeoptimization DGs:Name BusNo: jVjpu ]Vdeg PL(MW) QL(MVAR) Pg(MW) Qg(MVAR) Kp bp Kq bq Bus1 1.073 0.125 0 0 0 0 Bus2 1.063 0.037 0.600 0.420 0 0 Bus3 1.065 -0.102 0 0 0 0 Bus4 1.067 -0.007 0 0 0 0 Bus5 1.070 0.137 0.61 0.420 0 0 Bus6 DG1 1.091 1.110 0 0 2.1114 1.4093 -0.0293 60.0169 0.0656 0.9987 Bus7 1.048 -1.530 0.600 0.390 0 0 Bus8 1.038 -1.665 0.700 0.450 0 0 Bus9 1.041 -1.514 0.800 0.500 0 0 Bus10 1.031 -2.009 0.900 0.610 0 0 Bus11 Grid 1.067 0.000 0 0 0.000 0.000 Bus12 DG3 1.078 -0.077 0 0 0.9130 0.7061 -1.0648 60.9272 0.3767 0.8116 Bus13 DG2 1.078 -0.169 0 0 1.2548 0.9004 -0.1078 60.0903 0.6452 0.4973 Bus14 1.067 0.000 0 0 0 0 F2=K2p1,F3=K2q1,F4=K2p2,F5=K2q2F6=K2p3F7=K2q3F8=b2p1F9=b2q1F10=b2p2F11=b2q2F12=b2p3F13=b2q3and W1+W2+W3+W4+W5+W6+W7+W8+W9+W10+W11+W12+W13=1(6.3) Forthecasestudy,theproblemsolvedwithmulti-objectiveoptimizationby\fmincon00commandinMATLAB(seeappendixA).Thevaluesofobjectives36

PAGE 48

F1,F2,F3,F4,F5,F6,F7,F8,F9,F10,F11,F12andF13arenormalizedtobeconsistentwiththeweightsassignedtotheobjectivefunctions.Thefunctionsarenormalizedusingtheform:0Fi)]TJ /F3 11.955 Tf 11.96 0 Td[(Fimin Fimax)]TJ /F3 11.955 Tf 11.96 0 Td[(Fimin1(6.4) wheretheFiistheithobjectivefunction;Fiminisminimumvalueofithobjectivefunction;andFimaxisthemaximumvalueofithobjectivefunctioninthedierentweighted,thatisassignedtotheobjectivefunction.ResultsareshowninTable6.6. Table6.6:Thenormalizedofobjectivefunctionsinautonomousmodeafteroptimization W1 W2:W13 F1 F2 F3 F4 F5 F6 F7 F8 F9 F10 F11 F12 F13 1 0 0.085106383 0.014388834 0.088864536 0.697111063 1.97866E-18 3.87819E-05 1.000012322 1 0.128438349 0.124272989 0.135580785 1.12556E-17 1.000139336 0.9 0.00833 0.687943262 0.032172178 1.24537E-18 0.841746065 0.27036738 1 0.096737545 0.119596111 0.03213545 1.000000015 0.999949561 0.658812061 0.788466167 0.8 0.0167 0.709219858 0.032172178 0.416655789 1 0.142153052 0.315088821 0.128846533 0.581670718 0.011908363 0.03949486 2.45661E-17 0.781715352 0.902958985 0.7 0.025 0.680851064 0.099674757 0.10919729 0.480503545 0.000139537 3.87819E-05 0.035064372 0.418399402 1.000002212 0.970045669 0.078616496 0.027269824 0.090394067 0.6 0.0333 0.695035461 3.92438E-19 0.164773485 0.547379823 1.000019934 0.215431531 0.040584768 -4.67464E-06 2.70814E-05 0.049037947 0.644746978 0.038824834 0.210738986 0.5 0.04166 0.758865248 0.424917685 1.00003586 0.196919165 0.44534645 0.146438493 0.013043907 0.550473308 -1.50154E-17 0.013757267 0.291117184 0.971427611 0.166015378 0.4 0.05 0.680851064 1.000038609 0.567440913 0.360510891 0.051668461 4.38715E-18 0.015638986 0.105922775 0.009492929 0.033540092 0.008405018 1.000210091 -1.11651E-17 0.3 0.0583 0.737588652 0.688384534 0.354241764 0.218181903 1.67047901 0.12638824 0.011437077 0.267534592 0.010167436 0.003169549 1.047614681 0.982890181 0.000332506 0.2 0.0667 0.723404255 0.926835785 0.516317125 0.118091183 0.018364031 0.009272961 0.003304351 0.27557615 0.008642708 0.01149358 0.048573027 0.039976134 0.018240316 0.1 0.075 0.872340426 0.37919063 0.10761944 0.015645135 0.100306981 0.391417765 1.97157E-06 0.040762902 0.000224234 0.001106104 0.021129473 0.014092911 0.055126289 0 0.0833 1 0.057976946 0.018571299 8.50237E-08 0.047681697 0.00042456 0.005651012 0.007778609 0.000230869 6.62301E-07 0.027824141 0.008042287 0.062017075 ThedroopcharacteristicsofallDGsaredepictedinFig.6.3,Fig.6.4andFig.6.5fortwoconditions:nooptimizationandoptimizationwithW1=0:7,becausethisvalueofW1=0:7isthemostsuitableweightedcoecientforoptimization.37

PAGE 49

Figure6.3:DroopcharacteristicdeviationforDG1 Figure6.4:DroopcharacteristicdeviationforDG2 Fromthisalgorithm,theoperatorcandenethenewslopeofdroopchar-acteristicsattheoptimizedpoint.38

PAGE 50

Figure6.5:DroopcharacteristicdeviationforDG339

PAGE 51

7.Conclusion Thisstudyhaspresentedamethodologyfordeterminingoptimalpowerowinamicrogridinasteadystate,foreitheroneoftwomodesofoperation,i.e.,connectedtothegridanddisconnectedfromthegrid,actingasastand-aloneisolatedmicrogrid.Thisstudydevelopsanalgorithmwhichdenesthedroopcharacteristicparametersofaninverter.Thegoalistondtheoutputindeterminingtheoptimalpowertotheuser.Thisalgorithmassiststheusertodispatchpowerinthemicrogrid. Iftheloadchangesinthemicrogrid,thisresultsinchangesinthedroopcharacteristicsinboththeslope(KpandKq)andtheY-intersect(bpandbq)deviatingattheoptimizedpoint.Inaddition,whenthemicrogridisdis-connectedfromthegrid,thefrequencyisdierentfromthefrequencyofthegrid.Fromthefrequency-droopcharacteristicstheamountofpowerofDGwillbechanged.Thedroopcharacteristicparametersoftheinverter-basedDGshascomputedtoreachproperloadsharingatoptimizedpointratherthantobeusedaxedvalue.Thepowernecessaryforthesystemisprovidedbytheoptimalpowerowthatminimizesthecostoffuel. Transitionmodesfromtheconnectedoperatingmodetotheisolatedoper-atingmodeinamicrogridwerenotcoveredinthispaper,butonlythesteadystateconditionineachmodeisaddressedinthispaper.Thedynamictransitionbetweentwomodescanbelefttoafuturediscussionasanextensionofthisstudyscope.40

PAGE 52

APPENDIXA.\fmincon"inMATLAB fminconfunctioninMATLABisproposedtondtheminimumofacon-strainednon-linearmulti-variablefunction.Otherwisefminconndsacon-strainedminimumofafunctionofseveralvariables.Theequationcanbesolvedbythisfunctionexpressedas:minf(x)suchthat:c(x)0ceq(x)=0AxbAeqx=beqlbxubx,b,beq,lbandubarevectors,AandAeqarematrices,c(x)andceq(x)arefunctionthatreturnvectors(@nonlincon),andf(x)isafunctionthatreturnsascalar(@fun).f(x),c(x),andceq(x)canbenonlinearfunction.Thefunctionof"fmincon"inMATLABisexpressedas:[x;fval;exitflag;output;lambda]=fmincon(@fun;x0;A;b;Aeq;beq;lb;ub;@nonlcon)(A.1)41

PAGE 53

where: objectiveObjectivefunction x0Initialpointforx AMatrixforlinearinequalityconstraints bVectorforlinearinequalityconstraints AeqMatrixforlinearequalityconstraints beqVectorforlinearequalityconstraintslbVectoroflowerboundsubVectorofupperbounds nonlconNonlinearconstraintfunction lbLowerbound upUpperbound xoptimizedvalue fvalthevalueoftheobjectivefunctionatthesolutionX. lambdaLagrangeanmultipliers exitflagspecificationatalgorithmterminatio solver'fmincon' 42

PAGE 54

APPENDIXB.ThecodeofsoftwareinMATLAB B.1Readingthele functionS=Read_filecf(S) filename=uigetfile('*.cf'); S.filename=filename; filecf=fopen(S.filename,'r'); s=fgetl(filecf); whilestrcmp(s,'') s=fgetl(filecf); end; S.CaseName=s; whilestrcmp(s(1:min(3,length(s))),'BUS')~=1 s=fgetl(filecf); end s=fgetl(filecf); whilestrcmp(s,'') s=fgetl(filecf); end; S.Bus.number=[]; S.Bus.name=[]; %S.Bus.zone=[];43

PAGE 55

S.Bus.type=[]; S.Bus.V=[]; S.Bus.ang=[]; S.Bus.MWL=[]; S.Bus.MVARL=[]; S.Bus.MWG=[]; S.Bus.MVARG=[]; S.Bus.Bkv=[]; S.Bus.DV=[]; S.Bus.max=[]; S.Bus.min=[]; S.Bus.G=[]; S.Bus.B=[]; whiles(1:1)~='-' k2=0; [k1k2]=find_k(s,k2); ite=str2num(s(k1:k2)); S.Bus.number(ite,1)=ite; [k1k2]=find_k(s,k2); S.Bus.name=strvcat(S.Bus.name,s(k1:k2)); [k1k2]=find_k(s,k2); S.Bus.type(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2);44

PAGE 56

S.Bus.V(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.Bus.ang(ite,1)=str2num(s(k1:k2))*pi/180; [k1k2]=find_k(s,k2); S.Bus.MWL(ite,1)=str2num(s(k1:k2))/100; [k1k2]=find_k(s,k2); S.Bus.MVARL(ite,1)=str2num(s(k1:k2))/100; [k1k2]=find_k(s,k2); S.Bus.MWG(ite,1)=str2num(s(k1:k2))/100; [k1k2]=find_k(s,k2); S.Bus.MVARG(ite,1)=str2num(s(k1:k2))/100; [k1k2]=find_k(s,k2); S.Bus.Bkv(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.Bus.DV(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.Bus.max(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.Bus.min(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.Bus.G(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.Bus.B(ite,1)=str2num(s(k1:k2)); 45

PAGE 57

s=fgetl(filecf); end whilestrcmp(s(1:min(6,length(s))),'BRANCH')~=1 s=fgetl(filecf); end s=fgetl(filecf); whilestrcmp(s,'') s=fgetl(filecf); end; S.Branch.fr=[]; S.Branch.to=[]; S.Branch.r=[]; S.Branch.x=[]; S.Branch.b=[]; S.Branch.Mvar=[]; ite=1; whiles(1:1)~='-' k2=0; [k1k2]=find_k(s,k2); S.Branch.fr(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.Branch.to(ite,1)=str2num(s(k1:k2));46

PAGE 58

[k1k2]=find_k(s,k2); S.Branch.r(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.Branch.x(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.Branch.b(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.Branch.Mvar(ite,1)=str2num(s(k1:k2))/100; s=fgetl(filecf); ite=ite+1; end whilestrcmp(s(1:min(9,length(s))),'GENERATOR')~=1 s=fgetl(filecf); end s=fgetl(filecf); whilestrcmp(s,'') s=fgetl(filecf); end; ifstrcmp(s(1),'('); s=fgetl(filecf); end; whilestrcmp(s,''); s=fgetl(filecf); end47

PAGE 59

S.G.Bus=[]; S.G.Type=[]; S.G.Kind=[]; S.G.Pmax=[]; S.G.Pmin=[]; S.G.Qmax=[]; S.G.Qmin=[]; S.G.a=[]; S.G.b=[]; S.G.c=[]; %S.G.f=[]; S.G.Kp=[]; S.G.Yp=[]; S.G.Kq=[]; S.G.Yq=[]; S.G.P_ava_UnDis=[]; ite=1; whiles(1:1)~='-' k2=0; [k1k2]=find_k(s,k2); S.G.Bus(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.G.Type(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2);48

PAGE 60

S.G.Kind(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.G.Pmax(ite,1)=str2num(s(k1:k2))/100; [k1k2]=find_k(s,k2); S.G.Pmin(ite,1)=str2num(s(k1:k2))/100; [k1k2]=find_k(s,k2); S.G.Qmax(ite,1)=str2num(s(k1:k2))/100; [k1k2]=find_k(s,k2); S.G.Qmin(ite,1)=str2num(s(k1:k2))/100; [k1k2]=find_k(s,k2); S.G.a(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.G.b(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.G.c(ite,1)=str2num(s(k1:k2)); %[k1k2]=find_k(s,k2); %S.G.f(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.G.Kp(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.G.Yp(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.G.Kq(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2);49

PAGE 61

S.G.Yq(ite,1)=str2num(s(k1:k2)); [k1k2]=find_k(s,k2); S.G.P_ava_UnDis(ite,1)=str2num(s(k1:k2))/100; s=fgetl(filecf); ite=ite+1; end whilestrcmp(s(1:min(9,length(s))),'FREQUENCY')~=1 s=fgetl(filecf); end s=fgetl(filecf); whilestrcmp(s,'') s=fgetl(filecf); end; S.f=[]; ite=1; k2=0; [k1k2]=find_k(s,k2); S.f=str2num(s(k1:k2)); S.Bus.pq=find(S.Bus.type==0); S.Bus.pv=find(S.Bus.type==2); S.Bus.s=find(S.Bus.type==3); S.G.Dis=find(S.G.Kind==1); S.G.UnDis=find(S.G.Kind==0);50

PAGE 62

%S.G.Inv=find(S.G.Type==0); N=length(S.G.Bus); N1=N-length(S.G.UnDis); a=max(S.Bus.number); S.G.Pmax_Dis=S.G.Pmax(find(S.G.Kind==1)); S.G.Pmin_Dis=S.G.Pmin(find(S.G.Kind==1)); fclose(filecf); B.2CalculateY-Bus functionS=calculate_YBus(S) %S.filename=''; %S=Read_filecf(S); S.YB=sparse(diag(S.Bus.G+j*S.Bus.B)); S.YB=S.YB+sparse(S.Branch.fr,S.Branch.fr,... 1./(S.Branch.r+j*S.Branch.x)+j*S.Branch.b/2,... max(S.Bus.number),max(S.Bus.number)); S.YB=S.YB+sparse(S.Branch.to,S.Branch.to,... 1./(S.Branch.r+j*S.Branch.x)+j*S.Branch.b/2,... max(S.Bus.number),max(S.Bus.number)); S.YB=S.YB+sparse(S.Branch.fr,S.Branch.to,...51

PAGE 63

-1./(S.Branch.r+j*S.Branch.x),max(S.Bus.number),max(S.Bus.number)); S.YB=S.YB+sparse(S.Branch.to,S.Branch.fr,... -1./(S.Branch.r+j*S.Branch.x),max(S.Bus.numb) B.3Computationalcode clearall closeall clc %startingthereadingthedata. S.filename=''; S=Read_filecf(S); %Forputtingthebusthatisconnectedtothegridinthelasrraw: S=Org(S); %NNisthenumberofgeneratorswithoutthegrid. NN=length(S.G.Bus)-1; x1=input('ISTHESYSTEMCONNECTEDTOTHEGRID?','s') if(x1=='n') S.G.Bus1=S.G.Bus(1:NN); S.G.Kind1=S.G.Kind(1:NN); S.G.Type1=S.G.Type(1:NN); S.G.Pmax1=S.G.Pmax(1:NN); S.G.Pmin1=S.G.Pmin(1:NN); S.G.Qmax1=S.G.Qmax(1:NN);52

PAGE 64

S.G.Qmin1=S.G.Qmin(1:NN); S.G.a1=S.G.a(1:NN); S.G.b1=S.G.b(1:NN); S.G.c1=S.G.c(1:NN); %S.G.f=[]; S.G.Kp1=S.G.Kp(1:NN); S.G.Yp1=S.G.Yp(1:NN); S.G.Kq1=S.G.Kq(1:NN); S.G.Yq1=S.G.Yq(1:NN); S.G.Dis=S.G.Bus(find(S.G.Kind1==1)); S.G.UnDis=S.G.Bus(find(S.G.Kind1==0)); S.G.Inv1=find(S.G.Type1==1); N1=length(S.G.Inv1); %Ngisthenumberofgenerator, %Ndisthenumberofgeneratorthatcanbedispatchable %andNbisthenumberofbusses. Ng=length(S.G.Bus1); Nd=length(S.G.Dis); Nb=max(S.Bus.number); %S.G.Pmax_Disisthematrixofmaximumreactivepowerofdispatching %generatorsandS.G.Pmin_Disistheminimumreactivepowerofdispatching %generators.53

PAGE 65

S.G.Pmax_Dis=S.G.Pmax1(find(S.G.Kind1==1)); % S.G.Pmin_Dis=S.G.Pmin1(find(S.G.Kind1==1)); S.G.Kp2=S.G.Kp1(S.G.Inv1); S.G.Yp2=S.G.Yp1(S.G.Inv1); S.G.Kq2=S.G.Kq1(S.G.Inv1); S.G.Yq2=S.G.Yq1(S.G.Inv1); %Initialsvalue x0=[zeros(Nb,1);ones(Nb,1);(S.G.Pmin_Dis+S.G.Pmax_Dis);(S.G.Qmin1+S.G.Qmax1);59.5;S.G.Kp2/1.2;S.G.Yp2/1.2;S.G.Kq2/1.2;S.G.Yq2/1.2]; A=[]; b=[]; Aeq=[]; beq=[]; %upperbondofvariable. ub=[pi*ones(Nb,1);1.1*ones(Nb,1);S.G.Pmax_Dis;S.G.Qmax1;60.8;Inf*ones(16,1)]; %lowerbondofvariable. lb=[-pi*ones(Nb,1);0.9*ones(Nb,1);S.G.Pmin_Dis;S.G.Qmin1;59.2;-Inf*ones(16,1)]; options=optimset('DerivativeCheck','on','Diagnostics','on','Display',...54

PAGE 66

'iter','FunValCheck','off','GradObj','off','MaxFunEvals',30000,... 'MaxIter',10000,'Hessian','off'); [x,fval,exitflag,output,lambda,grad,hessian]=fmincon(@fun,x0,A,b,Aeq,beq,lb,ub,@nonlcon,options,S); Nd=length(S.G.Bus1)-length(S.G.UnDis); Ng=Nd+length(S.G.UnDis); Nb=max(S.Bus.number); S.Bus.V=x(Nb+1:2*Nb); S.Bus.ang=x(1:Nb); S.Bus.MWG=[]; S.Bus.MVARG=[]; S.Bus.MWG(S.G.Dis)=x(2*Nb+1:2*Nb+Nd)'; S.Bus.MWG(S.G.UnDis)=S.G.P_ava_UnDis(find(S.G.Kind1==0))'; S.Bus.MVARG(S.G.Bus1)=x(2*Nb+Nd+1:2*Nb+Nd+Ng)'; f=x(2*Nb+Ng+Nd+1); S.G.Kp(S.G.Inv1)=x(2*Nb+Ng+Nd+2:2*Nb+Nd+Ng+N1+1); S.G.Yp(S.G.Inv1)=x(2*Nb+Ng+Nd+N1+2:2*(Nb+N1)+Nd+Ng+1); S.G.Kq(S.G.Inv1)=x(2*(Nb+N1)+Nd+Ng+2:2*(Nb+N1)+Nd+Ng+N1+1); S.G.Yq(S.G.Inv1)=x(2*(Nb+N1)+Nd+Ng+N1+2:2*(Nb+N1)+Nd+Ng+N1+N1+1); S.f=x(2*Nb+Ng+Nd+1);55

PAGE 67

S.Bus.MWL(find(S.Bus.type==2))=0; S.Bus.MVARL(find(S.Bus.type==2))=0; else S.G.Inv=find(S.G.Type==1); Ng=length(S.G.Bus); Nd=length(S.G.Dis); Nb=max(S.Bus.number); S=calculate_YBus(S); S.G.Dis=S.G.Bus(find(S.G.Kind==1)); S.G.UnDis=S.G.Bus(find(S.G.Kind==0)); NN=length(S.G.Bus)-1; S.G.Kp1=S.G.Kp(S.G.Inv); S.G.Yp1=S.G.Yp(S.G.Inv); S.G.Kq1=S.G.Kq(S.G.Inv); S.G.Yq1=S.G.Yq(S.G.Inv); x0=[zeros(Nb,1);ones(Nb,1);ones(Nd,1);(S.G.Qmin+S.G.Qmax)/2;S.G.Kp1;S.G.Yp1;S.G.Kq1;S.G.Yq1]; A=[]; b=[]; Aeq=[]; beq=[]; ub=[pi*ones(Nb,1);1.1*ones(Nb,1);S.G.Pmax_Dis;S.G.Qmax;[]]; lb=[-pi*ones(Nb,1);0.9*ones(Nb,1);S.G.Pmin_Dis;S.G.Qmin;[]];56

PAGE 68

options=optimset('DerivativeCheck','on','Diagnostics','on','Display',... 'iter','FunValCheck','off','GradObj','off','MaxFunEvals',5000,... 'MaxIter',10000,'Hessian','off'); [x,fval,exitflag,output,lambda,grad,hessian]=fmincon(@fun1,x0,A,b,Aeq,beq,lb,ub,@nonlcon1,options,S); N1=length(S.G.Inv); S.Bus.V=x(Nb+1:2*Nb); S.Bus.ang=x(1:Nb); S.G.Kp(S.G.Inv)=x(2*Nb+Ng+Nd+1:2*Nb+Nd+Ng+N1); S.G.Yp(S.G.Inv)=x(2*Nb+Ng+Nd+N1+1:2*Nb+Ng+Nd+N1+N1); S.G.Kq(S.G.Inv)=x(2*(Nb+N1)+Nd+Ng+1:2*(Nb+N1)+Nd+Ng+N1); S.G.Yq(S.G.Inv)=x(2*(Nb+N1)+Nd+Ng+N1+1:2*(Nb+N1)+Nd+Ng+N1+N1); S.Bus.MWG=[]; S.Bus.MVARG=[]; S.Bus.MWG(S.G.Dis)=x(2*Nb+1:2*Nb+Nd)'; S.Bus.MWG(S.G.UnDis)=S.G.P_ava_UnDis(find(S.G.Kind==0))'; S.Bus.MVARG(S.G.Bus)=x(2*Nb+Nd+1:2*Nb+Nd+Ng)'; S.f=60;57

PAGE 69

S.Bus.MWL(find(S.Bus.type==2))=0; S.Bus.MVARL(find(S.Bus.type==2))=0; end S=Write_filecf(S); B.4Creatingfunfunction functionopt=fun2(x,S) %opt=0; %x=ones(100,1); Ng=length(S.G.Bus1); Nd=Ng-length(S.G.UnDis); Nb=max(S.Bus.number); Nc=Nb-1; N1=length(S.G.Inv1); p=100*x(Nc+Nb+1:Nc+Nb+Nd); kp=x(Nb+Nc+Ng+Nd+2:Nb+Nc+Nd+Ng+N1+1); yp=x(Nc+Nb+Ng+Nd+N1+2:Nc+Nb+2*N1+Nd+Ng+1); kq=x(Nb+Nc+Nd+Ng+N1+N1+2:Nb+Nc+Nd+Ng+N1+N1+1+N1); yq=x(Nb+Nc+2*N1+Nd+Ng+N1+2:Nb+Nc+2*N1+Nd+Ng+N1+N1+1); s1=S.G.a(1:Nd); s2=S.G.b(1:Nd);58

PAGE 70

s3=S.G.c(1:Nd); c=s1.*(p.^2)+s2.*p+s3; F1=sum(c w=1:2*N1+1; w(1)=0.6; fori=2:5*N1+1; w(i)=(1-w(1))/(4*N1); end %w(2)=1-w(1); op=0; forii=1:N1; op=op+w(4*ii-2)*((kp(ii)-S.G.Kp(ii))^2)+w(4*ii-1)*((kq(ii)-S.G.Kq(ii))^2)+w(4*ii)*((yp(ii)-S.G.Yp(ii))^2)+w(4*ii+1)*((yq(ii)-S.G.Yq(ii))^2); end opt=op+F1; B.5creatin"noncolon"function function[C,Ceq]=nonlcon(x,S) Nd=length(S.G.Bus1)-length(S.G.UnDis); Ng=Nd+length(S.G.UnDis); Nb=max(S.Bus.number); Vb=x(Nb+1:2*Nb);59

PAGE 71

ii=1:Nb; Vb(ii)=abs(Vb(ii)).*exp(j*x(1:Nb)); f=x(2*Nb+Ng+Nd+1); N1=length(S.G.Inv1); Kp=x(2*Nb+Ng+Nd+2:2*Nb+Nd+Ng+N1+1); Yp=x(2*Nb+Ng+Nd+N1+2:2*Nb+Ng+Nd+N1+N1+1); Kq=x(2*(Nb+N1)+Nd+Ng+2:2*(Nb+N1)+Nd+Ng+N1+1); Yq=x(2*(Nb+N1)+Nd+Ng+N1+2:2*(Nb+N1)+Nd+Ng+N1+N1+1); %f=x(2*(a+3*N)+1:2*(a+3*N)+N); f0=S.f*ones(N1,1); f1=f*ones(N1,1); R=S.Branch.Mvar; SF=Vb(S.Branch.fr).*conj(((Vb(S.Branch.fr)-Vb(S.Branch.to))./... (S.Branch.r+j*S.Branch.x*f/60))+Vb(S.Branch.fr).*(j*(S.Branch.b/2)*f/60)); ST=Vb(S.Branch.to).*conj(((Vb(S.Branch.to)-Vb(S.Branch.fr))./... (S.Branch.r+j*S.Branch.x*f/60))-Vb(S.Branch.to).*(j*(S.Branch.b/2)*f/60)); C=[abs(SF)-R;abs(ST)-R]; C=full(C); S.YB=sparse(diag(S.Bus.G+j*S.Bus.B*f/60)); S.YB=S.YB+sparse(S.Branch.fr,S.Branch.fr,...60

PAGE 72

1./(S.Branch.r+j*S.Branch.x*f/60)+j*S.Branch.b/2*f/60,... max(S.Bus.number),max(S.Bus.number)); S.YB=S.YB+sparse(S.Branch.to,S.Branch.to,... 1./(S.Branch.r+j*S.Branch.x*f/60)+j*S.Branch.b/2*f/60,... max(S.Bus.number),max(S.Bus.number)); S.YB=S.YB+sparse(S.Branch.fr,S.Branch.to,... -1./(S.Branch.r+j*S.Branch.x*f/60),max(S.Bus.number),max(S.Bus.number)); S.YB=S.YB+sparse(S.Branch.to,S.Branch.fr,... -1./(S.Branch.r+j*S.Branch.x*f/60),max(S.Bus.number),max(S.Bus.number)); S.YB=full(S.YB); Ib=S.YB*Vb; Si=Vb.*conj(Ib); Pi=-S.Bus.MWL; Qi=-S.Bus.MVARL; Pi(S.G.Dis)=x(2*Nb+1:2*Nb+Nd); Pi(S.G.UnDis)=S.G.P_ava_UnDis(find(S.G.Kind1==0)); Qi(S.G.Bus1)=x(2*Nb+Nd+1:2*Nb+Nd+Ng); f2=Kp.*Pi(S.G.Bus1(find(S.G.Type1==1)))+Yp; %Pr=(f1-Yp)./Kp; %Qr=(abs(Vb(S.Bus.pv))-Yq)./Kq-(S.Bus.V(S.Bus.pv)-S.G.Yq)./S.G.Kq+S.Bus.MVARG(S.Bus.pv); %Qr=(abs(Vb(S.G.Bus))-Yq)./Kq; V2=Kq.*Qi(S.G.Bus1(find(S.G.Type1==1)))+Yq;61

PAGE 73

Ceq=[Pi-real(Si);Qi-imag(Si);f2-f1;V2-abs(Vb(S.G.Bus1(find(S.G.Type1==1))))]; %;Pi(S.G.Bus)-Pr(S.G.Bus);Qi(S.G.Bus)-Qr(S.G.Bus) Ceq=full(Ceq); B.6writingtheleafteroptimization functionS=Write_filecf(S) S.filename=uiputfile('*.cf'); filecf=fopen(S.filename,'w'); fprintf(filecf,'%s\n',S.CaseName); fprintf(filecf,'BUSDATAFOLLOWS\n'); forii=1:length(S.Bus.number) fprintf(filecf,'%4.0f',S.Bus.number(ii)); fprintf(filecf,'%s',S.Bus.name(ii,:)); % %fprintf(filecf,'%1.0f',S.Bus.area(ii)); %fprintf(filecf,'%1.0f',S.Bus.zone(ii)); fprintf(filecf,'%1.0f',S.Bus.type(ii)); fprintf(filecf,'%5.3f',S.Bus.V(ii)); fprintf(filecf,'%6.3f',S.Bus.ang(ii)*180/pi); fprintf(filecf,'%7.3f',(S.Bus.MWL(ii))*100);62

PAGE 74

fprintf(filecf,'%7.3f',(S.Bus.MVARL(ii))*100); fprintf(filecf,'%7.4f',(S.Bus.MWG(ii))*100); fprintf(filecf,'%8.4f',(S.Bus.MVARG(ii))*100); fprintf(filecf,'%2.0f',S.Bus.Bkv(ii)); fprintf(filecf,'%2.0f',S.Bus.DV(ii)); fprintf(filecf,'%2.0f',S.Bus.max(ii)); fprintf(filecf,'%2.0f',S.Bus.min(ii)); fprintf(filecf,'%2.0f',S.Bus.G(ii)); fprintf(filecf,'%2.1f',S.Bus.B(ii)); fprintf(filecf,'\n'); end fprintf(filecf,'-999\n'); fprintf(filecf,'BRANCHDATAFOLLOWS\n'); forii=1:length(S.Branch.r); fprintf(filecf,'%3.0f',S.Branch.fr(ii)); fprintf(filecf,'%3.0f',S.Branch.to(ii)); %fprintf(filecf,'%3.0f',S.Branch.ar(ii)); %fprintf(filecf,'%3.0f',S.Branch.zone(ii)); %fprintf(filecf,'%3.0f',S.Branch.type(ii)); %fprintf(filecf,'%3.0f',S.Branch.circuit(ii)); fprintf(filecf,'%5.3f',S.Branch.r(ii)); fprintf(filecf,'%5.3f',S.Branch.x(ii)); fprintf(filecf,'%5.3f',S.Branch.b(ii)); fprintf(filecf,'%5.3f',S.Branch.Mvar(ii)*100);63

PAGE 75

fprintf(filecf,'\n'); end fprintf(filecf,'-999\n'); fprintf(filecf,'GENERATORDATAFOLLOWS\n') fprintf(filecf,'(BuskindPmaxPminQmaxQminabcKpYpKqYq)\n') forii=1:length(S.G.Bus); fprintf(filecf,'%4.0f',S.G.Bus(ii)); fprintf(filecf,'%4.0f',S.G.Type(ii)); fprintf(filecf,'%4.0f',S.G.Kind(ii)); fprintf(filecf,'%7.2f',(S.G.Pmax(ii))*100); fprintf(filecf,'%7.2f',(S.G.Pmin(ii))*100); fprintf(filecf,'%7.2f',(S.G.Qmax(ii))*100); fprintf(filecf,'%7.2f',(S.G.Qmin(ii))*100); fprintf(filecf,'%7.5f',S.G.a(ii)); fprintf(filecf,'%6.3f',S.G.b(ii)); fprintf(filecf,'%6.2f',S.G.c(ii)); fprintf(filecf,'%7.4f',S.G.Kp(ii)); fprintf(filecf,'%7.4f',S.G.Yp(ii)); fprintf(filecf,'%7.4f',S.G.Kq(ii)); fprintf(filecf,'%7.4f',S.G.Yq(ii)); fprintf(filecf,'%5.2f',S.G.P_ava_UnDis(ii)*100); fprintf(filecf,'\n'); end fprintf(filecf,'-99\n');64

PAGE 76

fprintf(filecf,'FREQUENCYINWHOLESYSTEM\n'); fprintf(filecf,'%7.3f',S.f); fprintf(filecf,'\n'); fprintf(filecf,'-99\n'); fclose(filecf);65

PAGE 77

REFERENCES[1] M.A.Abido.Environmental/economicpowerdispatchusingmultiobjec-tiveevolutionaryalgorithms:acomparativestudy.InPowerEngineeringSocietyGeneralMeeting,2003,IEEE,volume1,page4vol.2666,july2003.[2] T.Ackermann,G.Andersson,andL.Soder.Electricitymarketregulationsandtheirimpactondistributedgeneration.InElectricUtilityDeregulationandRestructuringandPowerTechnologies,2000.Proceedings.DRPT2000.InternationalConferenceon,pages608{613,2000.[3] VijayVitalArturR.Bergen.powersystemanalysis,.DepartmentofElec-tricalandComputerEngineering,pages306{334,2009.[4] J.Carpentier.Contributioneletudedodispatchingeconomique.Bull.Soc.Franc.Elect.,pages431{447,1962.[5] J.M.Carrasco,L.G.Franquelo,J.T.Bialasiewicz,E.Galvan,R.C.P.Guisado,Ma.A.M.Prats,J.I.Leon,andN.Moreno-Alfonso.Power-electronicsystemsforthegridintegrationofrenewableenergysources:Asurvey.IndustrialElectronics,IEEETransactionson,53(4):1002{1016,june2006.[6] D.CvetkovicandI.C.Parmee.Geneticalgorithm-basedmulti-objectiveoptimisationandconceptualengineeringdesign.InEvolutionaryCompu-tation,1999.CEC99.Proceedingsofthe1999Congresson,volume1,pages3vol.(xxxvii+2348),1999.[7] K..DeBrabandere,B..Bolsens,J..VandenKeybus,A..Woyte,J..Driesen,andR..Belmans.Avoltageandfrequencydroopcontrolmethodforparallelinverters.PowerElectronics,IEEETransactionson,22(4):1107{1115,july2007.[8] K.DeBrabandere,K.Vanthournout,J.Driesen,G.Deconinck,andR.Bel-mans.Controlofmicrogrids.InPowerEngineeringSocietyGeneralMeet-ing,2007.IEEE,pages1{7,june2007.[9] D.GaymeandU.Topcu.Optimalpowerowwithdistributedenergystoragedynamics.InAmericanControlConference(ACC),2011,pages1536{1542,292011-july12011.66

PAGE 78

[10] JonghyunRyu,SujinKim,andHongWan.Paretofrontapproximationwithadaptiveweightedsummethodinmultiobjectivesimulationoptimiza-tion.InSimulationConference(WSC),Proceedingsofthe2009Winter,pages623{633,dec.2009.[11] Hyun-KooKang,Seon-JuAhn,andSeung-IlMoon.Anewmethodtode-terminethedroopofinverter-baseddgs.InPowerEnergySocietyGeneralMeeting,2009.PES'09.IEEE,pages1{6,july2009.[12] F.KatiraeiandM.R.Iravani.Powermanagementstrategiesforamicrogridwithmultipledistributedgenerationunits.PowerSystems,IEEETrans-actionson,21(4):1821{1831,nov.2006.[13] F.Katiraei,R.Iravani,N.Hatziargyriou,andA.Dimeas.Microgridsman-agement.PowerandEnergyMagazine,IEEE,6(3):54{65,may-june2008.[14] B.Lasseter.Microgrids[distributedpowergeneration].InPowerEngineer-ingSocietyWinterMeeting,2001.IEEE,volume1,pages146{149vol.1,jan-1feb2001.[15] GuanghongLiu,GangWu,TaoZheng,andQingLing.Integratingpref-erencebasedweightedsumintoevolutionarymulti-objectiveoptimization.InNaturalComputation(ICNC),2011SeventhInternationalConferenceon,volume3,pages1251{1255,july2011.[16] A.SingheeandP.Castalino.Paretosampling:Choosingtherightweightsbyderivativepursuit.InDesignAutomationConference(DAC),201047thACM/IEEE,pages913{916,june2010.[17] XiaoboTangandGuoqingTang.Multi-objectiveplanningfordistributedgenerationindistributionnetwork.InElectricUtilityDeregulationandRestructuringandPowerTechnologies,2008.DRPT2008.ThirdInterna-tionalConferenceon,pages2664{2667,april2008.[18] ZhutianWang,XinhongHuang,andJinJiang.Designandimplementationofacontrolsystemforamicrogridinvolvingafuelcellpowermodule.InElectricalPowerConference,2007.EPC2007.IEEECanada,pages207{212,oct.2007.[19] WuZhe,WanYiru,HeChuan,YuJianhui,andZhouHao.Developmentstatusofchina'srenewableenergypowergeneration.InSustainablePowerGenerationandSupply,2009.SUPERGEN'09.InternationalConferenceon,pages1{5,april2009.67

PAGE 79

[20] CaoYijiaZhouXiao-xin.Thefutureofchinatodevelopnon-hydrorenew-ablesourcespowergeneration[j].InJournalofElectricPowerScienceAndTechnology,pages2{7,23(1):2008.68