OPTIMAL REACTIVE POWER CONTROL OF GRID CONNECTED
PHOTOVOLTAIC RESOURCES
by
JOSHUA RYAN TRIMBLE
B.S., Georgia Institute of Technology, 2006
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Electrical Engineering
2013
This thesis for the Master of Science degree by
Joshua Ryan Trimble
has been approved for the
Electrical Engineering Program
by
Fernando MancillaDavid, Advisor
Titsa Papantoni
Dan Connors
December 12, 2013
11
Trimble, Joshua Ryan (M.S., Electrical Engineering)
Optimal Reactive Power Control of Grid Connected Photovoltaic Resources
Thesis directed by Assistant Professor Fernando MancillaDavid
ABSTRACT
As more photovoltaic distributed generation resources are installed on distribu
tion power systems, selective control of the inverters connecting the DC power sources
presents the opportunity to supply both real and reactive power at the point of com
mon coupling. This thesis presents a simulated distribution system with individually
controlled PV resources with the objective of minimizing total system losses while op
erating at the maximum power point and below the simulated rating of the associated
inverters. The control strategy assumes the characteristics of the distribution system
are known and solves for the optimal power flow operating point. The ability of each
PV source to provide real and reactive power varies instantaneously as irradiance
changes, so the operating point for each resource must be constantly recalculated and
adjusted. The assumption of known system paramaters can be justified in a Smart
Grid context, and a solution based on overall system power flow should be considered
as a benchmark for any other state estimation or local control approaches.
The form and content of this abstract are approved. I recommend its publication.
Approved: Fernando MancillaDavid
m
TABLE OF CONTENTS
Figures ..................................................................... vi
Tables..................................................................... viii
Chapter
1. Introduction.............................................................. 1
1.1 Reactive Supply Capabilities and Historical Limitations.................. 1
1.2 Economic Considerations and Operating Modes.............................. 2
2. Photovoltaic Array and Inverter Model..................................... 5
2.1 Model of the PV Cell and Array............................................ 5
2.2 Model for the Grid Connected PV System................................... 6
2.3 Voltage Source Inverter Model............................................ 7
3. Power Flow Algorithm..................................................... 10
3.1 Introduction............................................................ 10
3.2 Formulation of the Power Flow Equations ................................ 10
3.3 Modified Formulation for Reactive Power Control of PV Resources ... 15
3.3.1 Formulation of the Power flow equations............................... 16
3.3.2 Objective Function and System Constraints............................. 18
4. Power System Simulation.................................................. 20
4.1 Distribution Feeder Model............................................... 20
4.2 Control System Interface................................................ 22
4.3 Timed Breaker Logic..................................................... 23
5. Test Cases and Simulation Results........................................ 25
5.1 Test Case Constant Reactive Power Output ............................. 25
5.2 Test Case Variable Irradiance......................................... 26
5.3 Test Case Variable Load............................................... 28
5.4 Measurement of Objective Function....................................... 31
6. Conclusion............................................................... 33
iv
References . ,
35
Appendix
A. MATLAB Reactive Power Control Function.............................. 36
A.l Main Function QCTL............................................... 36
A.2 Sub Function readCF.............................................. 40
A.2.1 Sub Function findklk2.......................................... 49
A.3 Sub Function  Build Ybus .......................................... 49
A.4 Sub Function  ObjFun............................................... 50
A. 5 Sub Function  NonLinCon ........................................ 51
B. Fortran Interface Module to link PSCAD and MATLAB .................... 53
v
FIGURES
Figure
1.1 Operating Region of the Voltage Source Inverter.............................. 4
2.1 Equivalent circuit for (a) a single PV cell and (b) a PV array of an arbi
trary size.................................................................. 6
2.2 Circuit schematic of the power stage for the gridconnected PV system. . 7
2.3 Control sourcesbased dynamic equivalent circuit schematic................... 7
2.4 Control sourcesbased dynamic equivalent circuit schematic................... 9
4.1 Feeder Model Part 1......................................................... 21
4.2 Feeder Model Part 2......................................................... 21
4.3 PSCAD to MATLAB Interface Component ........................................ 24
5.1 Real and reactive power output of the Bus?, inverter under high irradiance
and zero reactive power command............................................ 25
5.2 Real power output of the Bus? inverter under low irradiance and zero
reactive power command..................................................... 26
5.3 Reactive power output of the Bus?, inverter under low irradiance and zero
reactive power command..................................................... 26
5.4 Real and reactive power output of the Bus? inverter with steadily increas
ing irradiance............................................................. 27
5.5 Real and reactive power output of the Bus? inverter with steadily decreas
ing irradiance............................................................. 28
5.6 Reactive power output at Bus? with steadily decreasing irradiance and
reactive power control signal.............................................. 28
5.7 Real and reactive loading at Bus4........................................... 29
5.8 Real and reactive loading at Bus?........................................... 29
vi
5.9 Reactive output of the inverter at Bus3 ............................. 30
5.10 Reactive output of the inverter at Buss ................................. 30
5.11 Reactive power supplied by the utility at Bus\ without inverter support 31
5.12 Reactive power supplied by the utility at Bus 1 with inverter support . . 32
vii
TABLES
Table
3.1 Summary of known and unknown quantities in the Power Flow problem. 11
3.2 Summary of known unknown, and monitored quantities in the Power Flow
problem modified for PV Control....................................... 18
4.1 Power Flow Bus Data................................................... 22
4.2 Power Flow Branch Data................................................ 22
4.3 Power Flow Generator Data............................................. 23
viii
1. Introduction
This thesis presents a reactive power control strategy for photovoltaic (PV) re
sources that are connected to a typical distribution feeder through a DCAC inverter.
The aim of the control strategy is to prioritize the transfer of real power, but utilize
any unused inverter capacity to supply reactive power for voltage support on a radial
distribution feeder. The amount of reactive power to be supplied is based on a power
flow solution that assumes the entire network topology is known at every moment in
time. While this was not previously practical, advancements in online monitoring
and utility System Control and Data Acquisition (SCADA) capabilities is narrowing
the theoretical gap. Furthermore, a solution based on a known system configuration
can serve as a benchmark for the same problem where state estimation techniques
are used in lieu of actual quantities.
This chapter discusses the potential operating modes of distributed generation
supplied power via an inverter and discusses the economic incentives of both the
utility and the power generator.
1.1 Reactive Supply Capabilities and Historical Limitations
In the case of photovoltaic sources, an inverter of some capacity will always be
associated with the connection to the utility grid to transform the DC power supplied
by the PV array to an AC output that can be connected to the utility system. The
complex power rating of the inverter will likely exceed the nameplate rating of the
PV system. Several maximum power point algorithms have been developed to adjust
the output of the inverter to supply as much real power as possible by considering the
VI curve of the solar array in question. During periods of reduced irradiance due to
cloud cover or dark periods of the day, a significant portion of the inverter capacity
goes unutilized.
1
The ability of inverter connected sources to supply reactive power has historically
been restricted by utility interconnection guidelines. This was done largely to ensure
that the interconnecting generators would not be adding harmonies to the utility grid.
Thus, most interconnecting parties have been required to maintain var neutrality at
the point of common coupling (PCC). While the work presented in this thesis largely
pertains to utility distribution systems, advancements in wind turbines are having
a similar effect at the transmission level. The capability of doubly fed induction
generators (DFIGs) to dynamically supply reactive power is giving utilities the option
of specifying a nonunity power factor for the injected power to help control system
voltages.
Furthermore, the optimal control of any reactive power source would require
online adjustment of power levels as the system operating point is constantly chang
ing. FACTS devices on the transmission level have been implemented with real time
SCADA supplied system data to optimize operation, but that level of metering and
communication has historically not been available at the distribution level. Switched
capacitors and inductors are most commonly controlled locally by monitoring vars at
the point of common coupling or simply switched when voltage levels exceed a preset
bandwidth. With the proliferation of SmartGrid driven initiatives to add SCADA to
medium voltage 3 phase feeder networks, smart control of distributed generation is
becoming much more commonplace.
1.2 Economic Considerations and Operating Modes
Utilizing the stranded capability of the PV resources to supply reactive power
will have a direct economic benefit for the utility. Any excess inverter capacity that
supplies reactive power as needed would offset the need for the utility to supply
switched or fixed reactive elements on the distribution system, and could potentially
save operations of voltage regulating devices such as load tap changers or voltage
2
regulators as they respond to low voltage conditions as power factor drifts away
from unity. The distributed generation provider is typically metered at the point
of common coupling strictly on a real power basis and required to maintain unity
power factor at all times. While the benefits to the utility presented above have a
real cost, it would be difficult if not impossible to attribute the reactive contribution
of a single energy provider to avoiding a direct cost if reactive power was being
supplied by multiple generators on a given feeder simultaneously. To account for this
fact, the control strategy presented in later chapters will give preference to supplying
the maximum amount of real power and supply reactive power only with whatever
inverter capacity remains. This restriction ensures that the generation provider is
compensated for all real power provided without losing any revenue for contributing
to the reactive compensation of the system.
Since real and reactive power add in quadrature, the maximum reactive power
that the inverter can supply Qinv is given by:
Qi
C*2
^inv
p2
mv
(l.i)
where the Sinv is the MVA rating of the inverter and Pinv is the active power being
supplied at any given moment.
Figure 1.1 represents the various operating modes of the inverter as a circle with
a radius equal to Sinv [7].
The control strategy presented in this thesis will allow Pinv to be determined by
irradiance, temperature, and a Maximum Power Point Tracking Algorithm (MPPT)
while the optimal Qintl will be first calculated by a power flow algorithm but limited
to the operating regions shown in figure 1.1. The inverter should operate in the left
hand plane where Pinv < 0 under normal operation when power is being supplied
to the grid, though some power will be absorbed during the initial energization to
charge the DC capacitor. Both positive and negative values of reactive power can be
commanded corresponding to a capacitive or inductive effect respectively as needed
3
Q
inv
Figure 1.1: Operating Region of the Voltage Source Inverter.
if sufficient inverter capacity is available.
4
2. Photovoltaic Array and Inverter Model
This chapter presents the modeling approaches for the PV array and inverter sys
tem. The final model was abstracted and used to represent independent photovoltaic
generators connected to a utility distribution system. The PV array model includes
inputs for temperature and irradiance to allow for real world conditions that will
change the output power of each individual distributed generator. The final model
was realized in PSCAD.
2.1 Model of the PV Cell and Array
Although the circuital model for a single PV cell and its generalization to a
number of cells in series is well established in terms of a current source, an anti
parallel diode, a series resistance and shunt resistance [5], the model used in this
thesis makes use of a modified circuit that replaces the antiparallel diode by an
external control current source as proposed in [3]. This model was explained in [6],
and it is illustrated in Fig. 2.1(a). A chief advantage of this model is that it allows
for including an arbitrary number of cells connected in series and/or parallel into a
single circuital representation including all details of each cell.
The output current of the PV cell of Fig. 2.1(a) may be expressed as,
Ipr lirr Iho Ip
(2.1)
where Ii is the photo current or irradiance current generated when the cell is exposed
to sunlight. Idio is the current flowing through the antiparallel diode and induces
the nonlinear characteristics of the PV cell. Ip is a shunt current due to the shunt
resistor Rp branch. Substituting relevant expressions for Idio and IP,
Vpv
lirr Ifi
exp
(j{vpv + ipvRs)
nkT
 1
Vpv + Rv^S
Rn
(2.2)
5
v
(b)
Figure 2.1: Equivalent circuit for (a) a single PV cell and (b) a PV array of an
arbitrary size
where q = 1.602 x 1019 C is the electrons electric charge, k = 1.3806503 x 1023 J/K
is the Boltzmann constant, T is the temperature of the cell, /0 is the diode saturation
current or cell reverse saturation current, n is the ideality factor or the ideal constant
of the diode, and Rs and Rp represent the series and shunt resistance, respectively
[5]. This model can be generalized to an arbitrary number of cells connected in series,
say Ns, and in parallel, say Np, to form an array of Ns x Np cells. The generalized
model is illustrated in Fig. 2.1(b).
2.2 Model for the Grid Connected PV System
The circuit schematic of Fig. 2.2 illustrates the power stage of the singlestage
gridconnected PV power plant used in simulation. It includes the cascade connection
of a dc capacitor that connects to the output terminals of an arbitrary size PV array;
a twolevel threephase voltage source inverter (VSI); a LChlter in series with a
coupling transformer that connects to the ac grid.
As suggested in Fig. 2.2, the overall control is realized through a twolayer
strategy. An outer layer tracks the PV arrays MPP utilizing a MPPT algorithm. An
6
Figure 2.2: Circuit schematic of the power stage for the gridconnected PV system.
inner loop regulates the VSI in current control mode in order to transfer the power
generated by the PV array into the ac grid at an arbitrary power factor. The ability
to operate at an arbitrary power factor enables the PV system to dynamically supply
reactive power to the grid as commanded by the control described in Chapter 3.
2.3 Voltage Source Inverter Model
The cascaded connection of the dc capacitor, VSI, LChlter and ac grid is modeled
in the d q reference frame using dynamic phasors. Fig. 2.3 illustrates the dynamic
equivalent circuit in terms of control sources [9, 8, 4],
Figure 2.3: Control sourcesbased dynamic equivalent circuit schematic.
The state space equations for the system shown in 2.3 are given by (2.3) through
(2.6).
d 1 / 3
dt^ = C ( ~ 2
(2.3)
7
(2.4)
(Jj J t  j j
i = ju>i + (mvpv ~(R + Rf) i Rfig vf)
jff = ~juvf + ^(? ta)
(2.5)
(2.6)
In (2.3) through (2.6), the vector quantities relate to their abc counterparts consider
ing a generic vector quantity x = x,i + jxq, with
modulating function.
In equations (2.4) through (2.6), ui and $ are the synchronizing frequency and
phase angle, respectively. The synchronizing quantities are measured at the LC
hlters output and synthesized through a phaselockedloop (PLL). The PLL com
putes the grid phase angle by sensing the grid voltage and projecting the correspond
ing space vector on the d and (/axis of a d q rotating frame. This d q frame
is then rotated in such a way which ensures that the daxis of the d q frame is
aligned with the grid voltage vector, that is, uq = 0. In steadystate, the d q frame
rotational speed equals the grid angular frequency and the extracted angle equals to
the grid voltage angle. This angle is used to synchronize the d q reference frame for
the current control.
(2.7)
where T^c is the Park transformation. In (2.3), the operator denotes the dot
product between two vectors, yielding a dc value. In (2.3) and (2.4), m is the VSPs
The overall control strategy for the VSI is shown in figure 2.4. As observed
in the figure, the control is broken down into four blocks: (i) a PLL utilized for
synchronization with the ac grid; (ii) a dc voltage control to track the PV arrays
MPP; (Hi) a powertocurrent transformation; and (iv) a conventional current control
to transfer the PV arrays power into the ac grid at an arbitrary power factor.
Figure 2.4: Control sourcesbased dynamic equivalent circuit schematic.
The value of Qref is specified to compensate for the reactive power consumed by
the LC filter and stepup transformer if strictly real power output is desired. Any
reactive power command to be injected by the inverter is added to these compensating
levels.
9
3. Power Flow Algorithm
3.1 Introduction
The basic aim of any power/load flow calculations is to determine the voltage at
every bus within a utility system and then analyze the flow of power through the grid.
The results can show the required ratings of conductors connecting various generator
and can be used to model future system load growth and/or system expansion effects.
Load flow studies were originally performed using smaller scale analog models of
utility systems called calculator boards where RLC elements were connected in the
same topology as the utility grid. Advancements in computing made it possible to
digitally simulate and solve the system of equations that represents the actual utility
system.
The load flow problem was extended to the optimal power flow (OPF) problem by
Carpentier in 1962 [1], which applies constraints to certain parameters (bus voltages
or angles for stability purposes, conductor ratings) while minimizing or maximizing
a desired objective function such as system losses or generation model costs.
3.2 Formulation of the Power Flow Equations
In order to solve for the voltage magnitude and angle at every bus, it is assumed
that the network topology consisting of the points below are known.
1. The total number of busses in the system N representing the connection points
for generators, loads, and switching stations, etc.
2. The transmission topology that connects the various busses of the system in
cluding conductor RLC properties.
3. The number of generators and their location within the system.
10
4. The number and size of loads served by the system and their location within
the system.
The busses of the power system are further classified as follows
1. If at least one generator connected to a bus, the bus is designated as a PV bus
where the real power and bus voltage magnitude are considered quantities that
can be specified by the generator control. The total number of PV busses in the
system is designated as NPy, and the total number of generator in the system
is given by NG
2. One of the generator busses is arbitrarily identified as a slack bus. The bus
voltage is assumed (both magnitude and angle) are assumed to be known. The
slack bus voltage angle is considered to be 0 and defines the reference for all
other bus voltage angles in the system designated as 0^. The power provided
by the slack bus is modeled to supply system losses only.
3. If no generators are connected to a bus it is identified as a PQ bus and is
typically considered a system load. The real and reactive power leaving the
system are assumed to be known based on load projection data. The total
number of PQ busses in the system is designated as NPq
A summary of the known and unknown quantities in the power system is given in
Table 3.1
Table 3.1: Summary of known and unknown quantities in the Power Flow problem.
Real Power Reactive Power Voltage Magnitude Voltage Angle
PQ bus known known unknown unknown
PV bus controlled unknown controlled unknown
slack bus unknown unknown controlled known
11
The unknown and controlled quantities of interest are represented by the variable
vector X.
X
\vN\
PNG
V j
(3.1)
The dimensions of X is equal to 2 (IV + No) 1. In order to find a unique solution,
we will need to write an equal number of equations that include the elements of X and
other known quantities. The basis for the system of equations will be the complex
power injected at every bus, Sbus, which includes both a real and imaginary part for
each element. Sbus can be defines as:
Sbus Vbus I I
bus
(3.2)
or rewritten as:
Sbus Vbus(YbusVbus)
(3.3)
where Ybus is an admittance matrix between all the busses of the system and takes
the form:
12
/
\
vn y 12 * yin
1/21 y22 1/2n
(3.4)
yynl yn2 ynn
Ybus contains the topology of the system and the electrical properties of each trans
mission line or branch. Ybus is of size NxN and has the following properties.
Ybug is symmetrical.
The diagonal elements, Yu, are equal to the sum of the admittances of all
elements connected to the ith node.
The offdiagonal elements, Yy, are equal to the negative of the admittance
connecting nodes i and j.
The ith element of S^w, can be written as
n
(3.5)
or
n
(3.6)
and finally as
n
(3.7)
where:
13
Yij Gij + jBij
0ij = Bl9j
hi = 1? 2,n
(3.8)
Equation (3.7) gives an expression for the complex power injected at every bus
within the system and can be separated into real and imaginary parts by applying
Eulers formula and grouping like terms as show in equations (3.9) through (3.12)
Sbusi = Y [^ll^l(cs% + jsm0ij)(Gij jB^)] i,j = 1,2(3.9)
j=i
$busi ^ ^ 11 11 by  (Gij cos dij I Bjj sin 0,j I jGij sin 0ij j B,j cos Bjj) \ (3.10)
3=1
Pbusi = '^2[\Vi\\Vj\(Gij cosOij + Bi:j sin9^ i, j = 1, 2,..., n (3.11)
3 = 1
QbuSi = [l^ll^il(Go sin% Bv cos^')] i.j = 1, 2,..., n (3.12)
3 = 1
where:
Pbusi ' The active power injection in bus i.
Qbusi The reactive power injection in bus i.
The final system of equations is written by observing that all complex power
into the system must be equal to the complex power out of the system plus losses to
preserve energy balance. Since the slack bus accounts for losses, equations for real
and reactive power at all other busses can be written as a function of the unknown
vector X as follows:
Pbusi = Pgi ~ PlA = f (V, B) Qbusi = Qgi ~ QlA = / (V, 0) (3.13)
14
where:
Pgi : The active power of the generator in the bus i.
Qgi :The reactive power of the generator in the bus i.
PLi ' The active power of the load in the bus i.
Qu The reactive power of the load in the bus i.
With the formulation of the power flow equations complete, the system can be
solved iteratively to determine the load flows in the system. Applying constraints on
the variable vector X pertaining to system stability and capacity will limit the results
to a practical operating condition. Typical system constraints include:
1. Limitations on bus voltage magnitude, \ Vmin\ < Tj < WjOX, to ensure safe and
effective power delivery
2. Limitations on voltage angle, 0min < 0i < Qmaxi to ensure system stability
3. Limitations on branch power flow to less than the ratings of the circuit,
Sij < SijBating, to prevent overloading.
The final step of the formulation is to impose an optional objective function to
minimize or maximize in order to achieve an optimal power flow from the solution
set. In the case of this simulation, the objective will be to minimize system losses. In
larger transmission networks, the objective function is usually representative of the
cost of dispatching the generators on the system, which, when minimized, results in
the dispatching of lowest cost that supplies the network.
3.3 Modified Formulation for Reactive Power Control of PV Resources
The formulation of the power flow equations for the work presented in this thesis
is a special case of the concepts presented in the previous section, but the fundamental
concepts remain the same. This section discusses the differences in the formulation of
the power flow equations. The system constrains imposed on the solution and overall
15
objective function for optimal control are also presented. The final set of equations
were programmed into MATLAB, and the complete code can be found in Appendix
A.
3.3.1 Formulation of the Power flow equations
As stated previously, the electrical system of interested for the simulation pre
sented in this thesis is a heavily loaded radial distribution feeder with multiple PV
distributed generators in operation. The primary objectives of the control algorithm
are as follows.
1. Supply real power from the PV generators whenever possible based on instan
taneous irradiance and temperature conditions.
2. Supply reactive power only when inverter capacity is available.
3. The amount of reactive power to be supplied will be determined by the optimal
power flow solution, but not to exceed the capacity of the inverter.
With these objectives in mind, it is clear that the system bus where the PV
generators connect cannot be treated as a traditonal PV bus in the formulation of
the power flow equations since the real power injected, Pgi, depends on temperature
and irradiance and cannot be externally commanded. In the special case of a radial
distribution feeder, the only true PV bus is the substation bus that feeds the entire
system. The vector of unknown quantities defined in equation (3.1) can be modified as
shown in (3.14). The dimensions of X0 is consequently reduced to 2(iV + 1) + Nq 1.
The active power injection of every PV generator becomes an input to the power flow
algorithm as opposed to a quantity to be solved for.
16
Xn
( dNl\
\vN\
Psubstatian
V Qng )
(3.14)
It is worth noting that this approach will require the active power output of the PV
generator to be monitored constantly and relayed via a SCADA system to enable
online computation of the reactive power command. The instantaneous real power
output of the PV generator is likely to be available locally as part of the inverter
control and monitoring system, but a SCADA link for distributed generators is usu
ally only required for larger installations to address islanding and synchronization
concerns.
Another key distinction between the classic power flow formulation and the spe
cial case presented is this thesis is the need to respond quickly to changes in load.
While load profiles at a transmission level may be very predictable based on historical
load profiles, weather, and time of year, distribution feeder loads can vary quickly as
commercial and residential customers change electric consumption. In order to calcu
late the optimal amount of reactive power to command from the PV generators, the
real and reactive power consumption at every load bus will need to be monitored and
supplied to the control algorithm online. The final list of quantities that will have
to be monitored online and supplied to the reactive power controller as inputs is as
follows.
1. The real power generated by every PV distributed generator. Assuming there
are no power sources other than PV distributed generatorand the utility, the
total number of inputs will equal Npy 1
2. The real power consumed at every PQ load bus, which implies NPq signals.
17
3. The reactive power consumed at every PQ load bus, which implies Npq signals.
As a result of the need to monitor the above mentioned quantities online, table 3.1
can be modified as shown in table 3.2 to redesignate some of the quantities previously
considered as known.
Table 3.2: Summary of known unknown, and monitored quantities in the Power
Flow problem modified for PV Control.
Real Power Reactive Power Voltage Magnitude Voltage Angle
PQ bus monitored monitored unknown unknown
PVSubstation bllS controlled unknown controlled unknown
PVlnverier bllS monitored unknown controlled unknown
slack bus unknown unknown controlled known
Online feeder monitoring is a key aim of SmartGrid topics and is becoming much
more prevalent. One of the challenges of real world implementation is the distributed
nature of feeder loads as opposed to discrete load busses. Several state estimation
techniques for lumped equivalents have been developed, but the aim of this thesis is
to provide an optimal benchmark assuming all of the electrical data is known.
3.3.2 Objective Function and System Constraints
The objective function chosen for the simulation presented in this thesis is mini
mization of the system losses. Given that the simulated system is a radial distribution
feeder, all complex power is delivered to the system from the utility at the substation
feeder exit plus the contribution of the distributed generators. The active power con
tribution of the distributed generators take priority over the reactive output and is
controlled by a MPPT algorithm. Since the utility has to supply all load not served
by the distributed generators, minimization of the transmission losses will protide the
most cost effective operating point.
18
The system constrains imposed on the solution to ensure a realistic solution are
as follows:
1. Limitations on bus voltage magnitude in per unit, 0.95T4cromaj < V) < 11.05^0^^1
2. Limitations on voltage angle, 50
3. Limitations on the complex power supplied by the inverter, Sinv = hQQkVA
As previously discussed, the inverters complex power rating Sinv is the quadrature
addition of the real and reactive power supplied by the inverter. The real power
component is controlled strictly by the maximum power point tracking algorithm as
a function of temperature and irradiance. The maximum reactive power component
is computed according to equation (1.1) where Sinv is limited to 500 kVAR and Pinv
is read as an input to the controller.
19
4. Power System Simulation
In order to test the reactive power control algorithms ability to respond to online
changes in system conditions and commands, a model of a distribution feeder was
developed in PSCAD. The PV array and voltage source inverter model presented
in Chapter 2 were abstracted as PSCAD components and applied to the simulated
distribution feeder as distributed generators. A PSCAD component was developed
to pass the online monitored signals to the MATLAB power flow algorithm and also
to return the reactive power commands. Finally, some timed breaker logic was added
to demonstrate the controllers ability to deal with changes in system loading
4.1 Distribution Feeder Model
As previously stated, the scope of the simulation presented in this thesis is a
radial medium voltage distribution feeder. The model for the feeder is show in 4.1
and 4.2. Note that the model has been split into two parts for clarity. Bus3 is shown
in both figures.
The three phase system is assumed to be balanced and does not include lateral ties
to other distribution feeders. The substation feeder exit is shown connected to Bus 1,
the PV distributed generators are shown connected at Bus3 and Bus\. Lumped loads
are shown connected between all of the power sources at Bus2, Bus4, and at the end
of the feeder at Busq. A branch impedance between all busses of 0.1 + j0.05 ohms is
included as well.
In order to simulate typical utility voltage regulation, the voltage source connected
at Bus\ that simulates the substation exit is modeled as a controlled voltage source.
The magnitude of the starting feeder voltage is controlled by the power flow algorithm
within a tolerance of +/ 5%, which is the ANSI limit for utility voltage tolerance.
20
Bus 1
Bus 2
Bus 3
5.0 [MW] 1.7 [MVAR]
Figure 4.1: Feeder Model Part 1
Bus 3
Bus 4
Bus 5
Bus 6
Qcmd2
Figure 4.2: Feeder Model Part 2
Distribution voltage regulation is typically achieved by the mechanical operation of a
load tap changer or external voltage regulator on the secondary side of the substation
21
power transformer,
In order to supply the reactive power control algorithm with the static known
system parameters, a text hie was developed that organizes the system parameters
in the IEEE standard format for power how data [2]. The data of interest for the
simulated feeder is shown in tables 4.1through 4.3
Table 4.1: Power Flow Bus Data.
Name Type Voltage (V) MW Gen Mvar Gen MW Load Mvar Load
Bus\ Slack 12470 0 0 0 0
BuS2 PQ 12470 0 0 5 1.7
BUS3 PV 12470 0 0.4 0 0
Bus A PQ 12470 0 0 3 1.0
Bu$5 PV 12470 0 0.4 0 0
Bus 3 PQ 12470 0 0 3 1.7
Table 4.2: Power Flow Branch Data.
From To Resistance (ohms) Reactance (ohms)
Busi Bus2 0.1 0.05
Bus2 Bus3 0.1 0.05
Bus3 Bus 4 0.1 0.05
Busa Buss 0.1 0.05
Bus$ Buse 0.1 0.05
4.2 Control System Interface
While the PV Array, voltage source inverter, and the feeder itself are all simu
lated in PSCAD, the power how equations and iterative solution are performed in
MATLAB. A custom component in PSCAD was developed to handle the exchange
22
Table 4.3: Power Flow Generator Data.
Location Im, (MW) PMin (MW) Qrnam (MVar) PMin (MVar)
Bus\ 20 20 20 20
Bu.s3 0.5 0.5 0.5 0.5
Bus5 0.5 0.5 0.5 0.5
of data between the two programs. As outlined in the previous chapter, the following
data has to be monitored online and fed to the controller as an input.
1. Supply real power from the PV generators whenever possible based on instan
taneous irradiance and temperature conditions.
2. Supply reactive power only when inverter capacity is available.
3. The amount of reactive power to be supplied will be determined by the optimal
power flow solution, but not to exceed the capacity of the inverter.
The output of the MATLAB power flow algorithm is the reactive power commands to
each voltage source inverter, which is also passed through the interface component.
The graphical representation of the interface as applied to the distribution feeder
model is show in 4.3.
The FORTRAN code that governs the components data transfer is shown in
Appendix B. Note that a triggering input is included in the component that governs
how frequently MATLAB will execute the control algorithm. For the purposes of the
test cases presented in this thesis, a period of 0.25 seconds was selected.
4.3 Timed Breaker Logic
In order to demonstrate the controllers ability to respond appropriately to
changes in feeder load, timed breaker logic was added to the PQ busses, Bus4 and
Bus6 as shown in figure 4.2.
23
Qcmd
Trig
Matlab Trigger pg
Run every 0.5 sec
G
P3
P5
4
o
P6
P4
K V
c > 4
Qcmdl Qcmd2
MATLAB
Pi
o
Q6
Q4
Vs lack
Vs lack
ql'
Figure 4.3: PSCAD to MATLAB Interface Component
The breaker logic is intended to change the load state quickly in between execu
tions of the reactive power control algorithm so that the subsequent reactive power
command can be evaluated. The details of each breaker operation will be detailed for
each test case in the next chapter.
24
5. Test Cases and Simulation Results
The models for the PV array, voltage source inverter, and distribution feeder were
modeled in PSCAD and verified by simulation. The behavior of the MATLAB power
flow algorithm can also be observed by considering the distribution feeder model,
the input parameters that are monitored online, the reactive power commands issued
from the controller, and the final system response. Several test cases are presented
in this chapter in order of least to most complex.
5.1 Test Case Constant Reactive Power Output
In order to first verify the simulation of the PV array and voltage source inverter,
the first case presented is simply a constant reactive power command. A case of
particular interest is operation at unity power factor. A constant reactive power
command of zero is supplied to the voltage source inverters connected at Bus3 and
Bus5 for the system shown in figures 4.1 and 4.2. Additionally, the PV array at Bu.s3
will be supplied with a large irradiance value to simulate full sun exposure, while the
PV array at Bus5 will be supplied with a minimal irradiance value.
Figure 5.1: Real and reactive power output of the Bus3 inverter under high irradi
ance and zero reactive power command
25
Figure 5.2: Real power output of the Bus^ inverter under low irradiance and zero
reactive power command
Main: Graphs
Figure 5.3: Reactive power output of the Bus$ inverter under low irradiance and
zero reactive power command
We can observe that the real power output of the PV Generator at Bus^ that was
subjected to a high irradiance value is nearly equal to the maximum inverter rating
Sinv of 500 kVA, while the real power output of the the PV Generator at Bus5 was
essentially zero. For both cases, the reactive power output was essentially zero. Both
results confirm that the invert models real power output is dependant on irradiance,
and a desired reactive power value can be successfully commanded from the MATLAB
algorithm.
5.2 Test Case Variable Irradiance
Another test case of interest is the real and reactive power response of the system
to steadily increasing or decreasing irradiance while the feeder is heavily loaded. The
goal of this simulation is to demonstrate the following:
26
1. Real power supplied by the PV array is directly proportional to irradiance.
2. The magnitude of reactive power supplied by the PV array is indirectly propor
tional to irradiance.
3. The magnitude of the reactive power supplied is limited by the inverter rating
Sinv as expressed in equation (1.1)
4. The system is able to respond to the quantities monitored online and subsequent
changing reactive power commands.
The simulations results are presented in figures 5.4 through 5.6
Figure 5.4: Real and reactive power output of the Bus3 inverter with steadily
increasing irradiance
We can see in figure 5.4 the expected response to steadily increasing irradiance,
and the quadrature relationship between Pinv and Qinv. The same effect is even more
pronounced in figure 5.5 where real power quickly drops off as the applied irradiance
decreases.
The online tracking of input parameters and response to changing reactive power
commands is illustrated in 5.6, which plots the reactive power command issued from
the MATLAB control algorithm on top of the measured reactive power output.
We can clearly see that the MATLAB interface component is executing the re
active power control algorithm and providing an updated control signal ever 0.25
seconds as expected and the subsequent response of the inverter control.
27
Figure 5.5: Real and reactive power output of the Bus5 inverter with steadily
decreasing irradiance
Figure 5.6: Reactive power output at Bus5 with steadily decreasing irradiance and
reactive power control signal
5.3 Test Case Variable Load
The next test case presented is intended to highlight the controllers ability to
respond to changes in load. Timed breaker logic is applied to the loads connected
at BuSi and Bus6 to simulate capacitor switching and abrupt changes in load that
can occur on a distribution system due to customer interaction or protective equip
ment operation. The timed breaker operations are designed to impose the following
conditions.
The initial load at Bus4 is 3 MW and 1.7 Mvar inductive.
The initial load at Busq is 3 MW and 1.7 Mvar inductive.
At t = 0.6sec, the load at Buse will quickly change to 3 MW and 1.2 Mvar
capacitive to simulate a capacitor being switched on.
At t = l.lsec, the load at Bus4 will quickly change to 1 MW and 0 MVar to
simulate a fuse operation.
The loading at Bus4 and Buse over the simulation period are shown in figures 5.7
and 5.8
Figure 5.7: Real and reactive loading at Bus.
Figure 5.8: Real and reactive loading at Busq.
The irradiance applied to both inverters will be held at practically zero to enable
the inverter to supply reactive power almost completely up to the inverter rating
Sinv = 500kVA .
The system response is shown in figures 5.9 and5.10
The following observations can be made at each 0.25 second interval when the
controller recalculates the reactive power commands:
29
Figure 5.9: Reactive output of the inverter at Bus3
Figure 5.10: Reactive output of the inverter at Bus5
At t = 0.25sec and t = 0.5sec, both of the inverters at Bus3 and Bus5 are
supplying maximum reactive power due to the heavy loading at nearby busses.
At t = 0.75sec, the control algorithm is called for the first time since the load
at Buse has changed significantly from the previous state. Less reactive power
is required to fully compensate the load and line loss in the rest of the feeder
from Buse
There is no significant change of state at t = 1.0sec.
At t = 1.25sec, the control algorithm is called for the first time since the load
at Busi has changed significantly from the previous state. The load at Bus^
30
has dropped to the point that minimal reactive power is required at B11S3, and
an excess of reactive power at the of the feeder has caused the inverter at Buss
to enter the inductive operating range and start absorbing vars.
5.4 Measurement of Objective Function
In order to demonstrate that the objective function is being probably realized,
a final test cases is presented. First, if the reactive power commands issued to the
inverter are designed to minimize system losses, any load demand for reactive power
should be passed directly to the inverter and supplied up to the inverters capacity.
To illustrate this, consider the same circuit shown in 4.1 and 4.2, without any timed
breaker operations and a low value of irradiance applied to every PV system. The
intent is to model a heavily loaded figure with the distributed generators able to supply
reactive power only. As a baseline, the reactive power supplied by the modeled utility
system a Bus\ is shown in figure 5.11.
Figure 5.11: Reactive power supplied by the utility at Bus\ without inverter sup
port
When the reactive power control algorithm is enabled, the reactive power supplied
by the system is shown in figure 5.12.
As expected, we can observe roughly a 1 Mvar reduction in the amount of reactive
power that the utility system has to supply.
31
Figure 5.12: Reactive power supplied by the utility at Bus\ with inverter support
32
6. Conclusion
This thesis has presented a comprehensive model of a distribution feeder with
multiple independent PV distributed generators, and a control strategy for the asso
ciated voltage source inverters that is capable of providing reactive compensation to
the utility system without compromising the amount of real power generated by the
PV systems. In other words, the control strategy takes advantage of any available
inverter capacity to help regulate the system without impacting the privately owned
distributed generator.
The challenges in realizing this system in practice is chiefly the required online
monitoring and adjustment of real power injection from the various generators and
load conditions. Prior to the implementation of cellular communications and the pro
liferation of fiber optic communication circuits, it would not have been possible to
achieve this level of control. Fiber optics are being increasingly installed on utility
circuits either as part of the neutral conductor or separately strung on the same over
head poles as the primary conductors as part of the SmartGrid initiative. Even with
the advances in communication capabilities on utility systems, the ability to monitor
load conditions comprehensively and instantaneously, especially at the distribution
level, will likely remain challenging.
The results presented in this thesis could be used as a benchmark for future work
on similar systems that assume less than perfect knowledge of the utility system
at every moment in time and apply state estimation techniques to extrapolate the
unknown quantities from particular data points. Additional functionality could also
be added to the reactive power control algorithm to deal with interruptions to the
flow of online input parameters when loss of communication occurs. Since the voltage
source inverter is capable of generating waveforms above the power frequency, control
33
could be added to provide active filtering of harmonics introduced by nonlinear loads
as well. This is readily achievable since the control of the inverter is already performed
in the dq reference frame. Finally, more sophisticated modeling of the distribution
feeder itself could lead to more interesting test cases.
34
REFERENCES
[1] J. Carpentier. Contribution aletude du dispatching economique. Belletin de la
Societe Frangaise des Electricens, 3(8):431447, January 1962.
[2] M. Cerroni, F. MancillaDavid, A. Arancibia, F. RigantiFulginei, and E. Muladi.
Common format for exchange of solved load flow data. Power Appuratus and
Systems, IEEE Transactions on, PAS92(6):19161925, November 1973.
[3] M. Cerroni, F. MancillaDavid, A. Arancibia, F. RigantiFulginei, and E. Mu
ladi. A maximum power point tracker variable dclink threephase inverter con
trol scheme for gridconnected pv panels. In Innovative Smart Grid Technologies
(ISGT Europe), 2012 3rd IEEE PES International Conference and Exhibition on,
pages 17, 2012.
[4] R. Kadri, J. P. Gaubert, and G. Champenois. An improved maximum power
point tracking for photovoltaic gridconnected inverter based on voltageoriented
control. IEEE Transactions on Industrial Electronics, 58(l):6675, January 2011.
[5] G. M. Masters. Renewable and Efficient Electric Power Systems. John Wiley &
Sons, New Jersey, 2004.
[6] H. Tian, F. MancillaDavid, K. Ellis, E. Muljadi, and P. Jenkins. A cell
tomoduletoarray detailed model for photovoltaic panels. Solar Energy,
86(9):26952706, September 2012.
[7] K. Turitsyn, P. Sulk, S. Backhaus, and M. Chertkox. Local control of reactive
power by distributed photovoltaic generators. In Distributed Control of Reactive
Power Flow in a Radially Distribution Circuit with High Photovoltaic Penetration,
pages 17, 2012.
[8] Y. Yang, M. Kazerani, and V. H. Quintana. Modeling, control and implementa
tion of threephase pwm converters. IEEE Transactions on Power Electronics,
18(3):857 864, May 2003.
[9] A. Yazdani and R. Iravani. VoltageSourced Converters in Power Systems. Wiley
IEEE Press, Reading, MA, 2010.
35
Appendix A. MATLAB Reactive Power Control Function
The C language code for the reactive power flow algorithm developed in MATLAB
is presented in this appendix. The code is organized as the main function called by
PSCAD during the simulation and also includes the referenced sub functions.
A.l Main Function QCTL
/Output QCMD reactive power commands.
/Input PG real power output of each inverter,
/Input PL real power consumed at each PQ Bus
/Input QL reactive power consumed at each PQ Bus
function [Qcmd, Vs] =Qctl(PG, PL, QL)
/Initial settings
j = sqrt(l);
S.base.S=1000000;
S.base.Vg=12470;
S.base.Zg=S.base.Vg~2/S.base.S;
S.base.Ig=S.base.S/S.base.Vg/sqrt(3);
S.base.Vi=270;
S.base.Zi=S.base.Vi~2/S.base.S;
S.base.Ii=S.base.S/S.base.Vi/sqrt(3);
36
/. %This initializes S
S.fname=' ';
% %This calls the fucntion read_cf
S=read_cf(S);
save Sfile
/(Update input values monitored online.
S.Bus,MWGen(S.Bus.PV) = PG;
S.Bus,MWLoad(S.Bus,PQ(2:length(S.Bus.PQ))) = PL;
S.Bus.MVARLoad(S.Bus.PQ(2:length(S.Bus.PQ))) = QL;
% Build the Ybus matrix
S=build_ybus(S);
/(Variables:
ID Angles at every bus( excluding slack )
/02) Voltages at every bus
#/04) Active Power at substation
/05) Reactive Power at every generator (PV + substation)
na=S.Bus.n1;
nv=S.Bus.n;
np=l;
nq=length(S.Gen.Bus);
37
tetamin=50;
tetamax=50;
Vmin=.95;
Vmax=l.05;
Slower bound
LB= [tetamin*ones(na,l)*pi/180; Vmin*ones(nv,1); inf; inf; ones(nq1,1)];
%upper bound
UB=[tetamax*ones(na,l)*pi/180; Vmax*ones(nv,l); inf; inf; ones(nq1,1)];
/(Initial condition
/ox0= [zeros(na, 1); ones(nv,l); 0.5*(S.Gen.Pmin+S.Gen.Pmax)/100; 0.5*(S.Gen.Qmin+S.G
x0=[zeros(na,1); ones(nv,l); 14000000/S.base.S; 1/S.base.S*[3000000; 500000; 50000
A=[] ;
b=[] ;
Aeq= [] ;
beq= [] ;
/(Defines optimset
OptOpf = optimset('Display,iter,Diagnostics',on,DerivativeCheck,on,...
'LargeScale, off GradConstr,off', GradObj, off, Hessian, off ,
'Maxlter, 5000, 'MaxFunEvals, 100000);
[x,f val,exitflag,output,lambda] = fmincon(@0bjFun,xO,A,b,Aeq,beq,LB,UB,ONonLinCon,
"/compute number of each variable
na=S.Bus.n1;
nv=S.Bus.n;
np=l;
nq=length(S.Gen.Bus);
"/.Updates Voltages
Vbus=x (na+1: na+nv) ; "/added
i_av=f ind(S.Bus.type~=3);
"/Vbus (i_av)=x(na+l:na+nv) ;
Vbus(S.Bus.Swing)=Vbus(S.Bus.Swing)*exp(j*0);
Vbus(i_av)=abs(Vbus(i_av)),*exp(j*x(l:na));
Sinj =Vbus.*conj(S.Ybus*Vbus);
Qcmd = x(na+nv+np+2:length(x))J;
Vs = x(na+l);
end
A.2 Sub Function readCF
39
The intent of this sub function is to build a MATLAB structure with all of the
static bus and branch data of the power system needed to build the power flow
equations.
function S=read_cf(S)
#/gets the name of the file
fname=JRTFeeder.cf;
S.fname=fname;
70opens the file for reading
fcf=fopen(fname,'r');
%gets the name of the case
s=fgetl(fcf);
while strcmp(s,), s=fgetl(fcf) ; end;
S.CaseName=s;
breads bus data
while strcmp(s(1:min(3,length(s))),BUS)~=1, % Find the start of bus data
s=fgetl(fcf);
end
s=fgetl(fcf);
while strcmp(s,), s=fgetl(fcf); end;
40
S.Bus.number=[];
S.Bus.name=[] ;
S.Bus.area=[] ;
S.Bus.zone=[];
S.Bus.type=[] ;
S.Bus.V=[] ;
S.Bus.angle= [] ;
S.Bus.MWLoad=[];
S. Bus. MVARLoad= [] ;
S Bus MWGen= [] ;
S.Bus.MVARGen=[];
%Extra stuff added
S.Bus.BasekV=[] ;
S.Bus.DesiredVolt=[] ;
S.Bus.MVarMax=[];
S.Bus.MVarMin=[];
S.Bus.G=[] ;
S.Bus.B=[] ;
while s(l)~=''
k2=l;
[kl k2]=find_klk2(s,k2);
41
n_aux=str2num(s(kl:k2));
S.Bus.number(n_aux,1)=n_aux;
[kl k2]=find_klk2(s,k2);
S.Bus.name=strvcat(S.Bus.name,s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Bus.area(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Bus.zone(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Bus.type(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Bus.V(n_aux,l)=str2num(s(kl:k2))/S.base.Vg;
[kl k2]=find_klk2(s,k2);
S.Bus.angle(n_aux,l)=str2num(s(kl:k2))*pi/180;
[kl k2]=find_klk2(s,k2);
42
S.Bus.MWLoad(n_aux,l)=str2num(s(kl:k2))/S.base.S;
[kl k2]=find_klk2(s,k2);
S.Bus.MVARLoad(n_aux,l)=str2num(s(kl:k2))/S.base.S;
[kl k2]=find_klk2(s,k2);
S.Bus.MWGen(n_aux,l)=str2num(s(kl:k2))/S.base.S;
[kl k2]=find_klk2(s,k2);
S.Bus .MVARGen(n_aux,l)=str2num(s(kl:k2))/S.base.S;
[kl k2]=find_klk2(s,k2);
S.Bus.BasekV(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Bus.DesiredVolt(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Bus.MVarMax(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Bus.MVarMin(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Bus.G(n_aux,l)=str2num(s(kl:k2));
43
[kl k2]=find_klk2(s,k2);
S. Bus B (n_aux, 1) =str2num(s (kl: k2) ) ;
s=fgetl(fcf);
end
%Find indices for various tyoe of buses
S.Bus.PQ=f ind(S.Bus.type==0) ;
S.Bus.PV=find(S.Bus.type==2);
S.Bus.Swing=f ind(S.Bus.type==3);
S.Bus.n=max(S.Bus.number);
%reads branch data
while strcmp(s(1:min(6,length(s))),BRANCH')~=1, % Find the start of bus data
s=fgetl(fcf);
end
s=fgetl(fcf);
while strcmp(s,J), s=fgetl(fcf); end;
S.Branch.From=[];
S .Branch.To= [] ;
S. Branch. area= [] ;
S Branch. zone= [] ;
44
S.Branch.circuit=[] ;
S.Branch.type=[];
S.Branch.r=[] ;
S.Branch.x=[] ;
S.Branch.b=[] ;
S.Branch.MVA1=[];
n_aux=0;
while s(l)~=,>
n_aux=n_aux+l;
k2=l;
[kl k2]=find_klk2(s,k2);
S.Branch.From(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Branch.To(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Branch.area(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Branch.zone(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Branch.circuit(n_aux,1)=str2num(s(kl:k2));
45
[kl k2]=find_klk2(s,k2);
S.Branch.type(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Branch.r(n_aux,l)=str2num(s(kl:k2))/S.base.Zg;
[kl k2]=find_klk2(s,k2);
S.Branch.x(n_aux,l)=str2num(s(kl:k2))/S.base.Zg;
[kl k2]=find_klk2(s,k2);
S.Branch.b(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Branch.MVA1(n_aux,l)=str2num(s(kl:k2))/S.base.S;
s=fgetl(fcf);
end
S.Branch.n=n_aux;
while strcmp(s(1:min(9,length(s))),GENERATOR)~=1, % Find the start of generator
s=fgetl(fcf) ;
end
s=fgetl(fcf);
46
while strcmp(s,'), s=fgetl(fcf) ; end;
if strcmp(s(1),('), s=fgetl(fcf); end;
while strcmp(s,''), s=fgetl(fcf); end;
S.Gen.Bus=[] ;
S.Gen.Pmax=[];
S.Gen.Pmin=[];
S.Gen.Qmax=[];
S.Gen.Qmin=[];
S.Gen.S=[] ;
S.Gen.a=[] ;
S.Gen.b=[] ;
S .Gen. c= [] ;
n_aux=0;
while s(l)~=,:
n_aux=n_aux+l;
k2=l;
[kl k2]=find_klk2(s,k2);
S.Gen.Bus(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Gen.Pmax(n_aux,l)=str2num(s(kl:k2));
47
[kl k2]=find_klk2(s,k2);
S.Gen.Pmin(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Gen.Qmax(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Gen.Qmin(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Gen.S(n_aux,l)=str2num(s(kl:k2))/S.base.S;
[kl k2]=find_klk2(s,k2);
S.Gen.a(n_aux,1)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Gen.b(n_aux,l)=str2num(s(kl:k2));
[kl k2]=find_klk2(s,k2);
S.Gen.c(n_aux,l)=str2num(s(kl:k2));
s=fgetl(fcf) ;
end
S.Gen.n=n_aux;
48
fclose(fcf);
A.2.1 Sub Function findklk2
The intent of this helper function is to find the start and stop coordinates of
numerical values in a text hie in order to avoid problems with excessive white space.
function [kl k2]=find_klk2(s,k2);
kl=k2+l;
while s(kl)==' kl=kl+l; end
k2=kl;
while s(k2)~=
if k2==length(s)
break
end
k2=k2+l;
end
A.3 Sub Function Build Ybus
This subfunction builds the sparse Ybus Matrix.
function S=build_ybus(S)
"/.Build ybus
"/.fills the diagonal elements from shunt compensators
S.Ybus=sparse(diag(S.Bus.G+j *S.Bus.B));
"/.fills diagonal elements from transmission lines and transformers
S.Ybus=S.Ybus+sparse(S.Branch.From,S.Branch.From,...
1./(S.Branch.r+j *S.Branch.x)+j *S.Branch.b/2,...
49
S.Bus.n,S.Bus.n);
S.Ybus=S.Ybus+sparse(S.Branch.To,S.Branch.To,...
1./(S.Branch.r+j*S.Branch.x)+j*S.Branch.b/2,...
S.Bus.n,S.Bus.n);
/.fills off diagonal elements
S.Ybus=S.Ybus+sparse(S.Branch.From,S.Branch.To,...
1./(S.Branch.r+j *S.Branch.x),...
S.Bus.n,S.Bus.n);
S.Ybus=S.Ybus+sparse(S.Branch.To,S.Branch.From,...
1./(S.Branch.r+j *S.Branch.x),...
S.Bus.n,S.Bus.n);
A.4 Sub Function ObjFun
This subfunction defines the objective function to be minimized or maximized by
MATLABs built in Fmincon function.
function f=ObjFun(x,S)
na=S.Bus.n1;
nv=S.Bus. n1;
np=l;
nq=length(S.Gen.Bus);
50
f=x(na+nv+l)+sum(S.Bus,MWGen(S.Bus.PV))sum(S.Bus.MWLoad(S.Bus.PQ));
A.5 Sub Function NonLinCon
This subfunction defines the nonlinear equality and inequality constraints for use
with MATLABs built in Fmincon function.
function [C, Ceq]=NonLinCon(x,S)
'/.compute number of each variable
na=S.Bus.n1;
nv=S.Bus.n; '/.took out 1
np=l;
nq=length(S.Gen.Bus);
'/.Updates Voltages
"/Vbus=S.Bus.V;
Vbus=x (na+1: na+nv) ; '/.added
Vbus(S.Bus.Swing)=Vbus(S.Bus.Swing)*exp(j*0);
i_av=find(S.Bus.type"=3);
'/.Vbus (i_av)=x(na+1:na+nv) ;
Vbus(i_av)=abs(Vbus(i_av)).*exp(j*x(l:na));
'/.Lines' Flow
Sfrom=Vbus(S.Branch.From).*conj((Vbus(S.Branch.From)Vbus(S.Branch.To))./(S.Branch
Vbus(S.Branch.From).*(j *S.Branch.b/2));
Sto=Vbus(S.Branch.To).*conj( (Vbus(S.Branch.To)Vbus(S.Branch.From))./(S.Branch.r+
Vbus(S.Branch.To).*(j *S.Branch.b/2) );
51
*/.Sinj=Vbus *conj (S Ybus*Vbus) ;
"/.update injected power
Pinj=S.Bus.MWLoad;
Qinj=S.Bus.MVARLoad;
Pinj(S.Bus.PV)=S.Bus.MWGen(S.Bus.PV);
Pinj(S.Bus.Swing)=x(na+nv+l);
Qinj(S.Gen.Bus)=x(na+nv+np+1:na+nv+np+nq);
Sinj=Pinj+j*Qinj;
Rate=S.Branch.MVA1;
"/.Inequality Constraints
C=[abs(Sfrom)Rate; abs(Sto)Rate; sqrt( real(Sinj(S.Gen.Bus)).~2+imag(Sinj(S.Gen.
"/.Power flow equations
Sv=Vbus.*conj(S.Ybus*Vbus);
"/.Equality constraint
Ceq=[Pinjreal(Sv); Qinjimag(Sv)];
52
Appendix B. Fortran Interface Module to link PSCAD and MATLAB
The following code was developed for the PSCAD component that links the sim
ulation to the MATLAB reactive power control algorithm. The code handles the
memory assignments of the various inputs and outputs.
#ST0RAGE REAL:8
#L0CAL INTEGER I_CNT
IF ($TRIG==1) THEN
DO I_CNT=1,2,1
STORF(NST0RF+I_CNT1) = $pg(I_CNT)
ENDDO
DO I_CNT=1,2,1
STORF(NST0RF+2+I_CNTl) = $pl(I_CNT)
ENDDO
DO I_CNT=1,2,1
STORF(NST0RF+4+I_CNTl) = $ql(I_CNT)
ENDDO
CALL MLAB_INT("$Path", "$Name","R(2) R(2) R(2)","R(2)")
DO I_CNT=1,2,1
$Qcmd(I_CNT) = STORF(NST0RF+6+I_CNTl)
ENDDO
END IF
53
NSTORF = NSTORF + 8
54

PAGE 1
OPTIMALREACTIVEPOWERCONTROLOFGRIDCONNECTEDPHOTOVOLTAICRESOURCESbyJOSHUARYANTRIMBLEB.S.,GeorgiaInstituteofTechnology,2006AthesissubmittedtotheFacultyoftheGraduateSchooloftheUniversityofColoradoDenverinpartialfulllmentoftherequirementsforthedegreeofMasterofScienceElectricalEngineering2013
PAGE 2
ThisthesisfortheMasterofSciencedegreebyJoshuaRyanTrimblehasbeenapprovedfortheElectricalEngineeringProgrambyFernandoMancillaDavid,AdvisorTitsaPapantoniDanConnorsDecember12,2013ii
PAGE 3
Trimble,JoshuaRyan(M.S.,ElectricalEngineering)OptimalReactivePowerControlofGridConnectedPhotovoltaicResourcesThesisdirectedbyAssistantProfessorFernandoMancillaDavidABSTRACT Asmorephotovoltaicdistributedgenerationresourcesareinstalledondistributionpowersystems,selectivecontroloftheinvertersconnectingtheDCpowersourcespresentstheopportunitytosupplybothrealandreactivepoweratthepointofcommoncoupling.ThisthesispresentsasimulateddistributionsystemwithindividuallycontrolledPVresourceswiththeobjectiveofminimizingtotalsystemlosseswhileoperatingatthemaximumpowerpointandbelowthesimulatedratingoftheassociatedinverters.Thecontrolstrategyassumesthecharacteristicsofthedistributionsystemareknownandsolvesfortheoptimalpowerowoperatingpoint.TheabilityofeachPVsourcetoproviderealandreactivepowervariesinstantaneouslyasirradiancechanges,sotheoperatingpointforeachresourcemustbeconstantlyrecalculatedandadjusted.TheassumptionofknownsystemparamaterscanbejustiedinaSmartGridcontext,andasolutionbasedonoverallsystempowerowshouldbeconsideredasabenchmarkforanyotherstateestimationorlocalcontrolapproaches. Theformandcontentofthisabstractareapproved.Irecommenditspublication.Approved:FernandoMancillaDavidiii
PAGE 4
TABLEOFCONTENTS Figures.......................................vi Tables........................................viii Chapter 1.Introduction...................................1 1.1ReactiveSupplyCapabilitiesandHistoricalLimitations..........1 1.2EconomicConsiderationsandOperatingModes..............2 2.PhotovoltaicArrayandInverterModel....................5 2.1ModelofthePVCellandArray.......................5 2.2ModelfortheGridConnectedPVSystem.................6 2.3VoltageSourceInverterModel........................7 3.PowerFlowAlgorithm.............................10 3.1Introduction..................................10 3.2FormulationofthePowerFlowEquations.................10 3.3ModiedFormulationforReactivePowerControlofPVResources...15 3.3.1FormulationofthePowerowequations.................16 3.3.2ObjectiveFunctionandSystemConstraints...............18 4.PowerSystemSimulation............................20 4.1DistributionFeederModel..........................20 4.2ControlSystemInterface...........................22 4.3TimedBreakerLogic.............................23 5.TestCasesandSimulationResults.......................25 5.1TestCaseConstantReactivePowerOutput...............25 5.2TestCaseVariableIrradiance.......................26 5.3TestCaseVariableLoad..........................28 5.4MeasurementofObjectiveFunction.....................31 6.Conclusion....................................33iv
PAGE 5
References ......................................35Appendix A.MATLABReactivePowerControlFunction.................36 A.1MainFunctionQCTL............................36 A.2SubFunctionreadCF............................40 A.2.1SubFunctionndk1k2..........................49 A.3SubFunctionBuildYbus..........................49 A.4SubFunctionObjFun............................50 A.5SubFunctionNonLinCon..........................51 B.FortranInterfaceModuletolinkPSCADandMATLAB..........53v
PAGE 6
FIGURES Figure 1.1OperatingRegionoftheVoltageSourceInverter..............4 2.1Equivalentcircuitfor(a)asinglePVcelland(b)aPVarrayofanarbitrarysize...................................6 2.2Circuitschematicofthepowerstageforthegrid{connectedPVsystem..7 2.3Controlsources{baseddynamicequivalentcircuitschematic........7 2.4Controlsources{baseddynamicequivalentcircuitschematic........9 4.1FeederModelPart1.............................21 4.2FeederModelPart2.............................21 4.3PSCADtoMATLABInterfaceComponent................24 5.1RealandreactivepoweroutputoftheBus3inverterunderhighirradianceandzeroreactivepowercommand......................25 5.2RealpoweroutputoftheBus5inverterunderlowirradianceandzeroreactivepowercommand...........................26 5.3ReactivepoweroutputoftheBus5inverterunderlowirradianceandzeroreactivepowercommand...........................26 5.4RealandreactivepoweroutputoftheBus3inverterwithsteadilyincreasingirradiance.................................27 5.5RealandreactivepoweroutputoftheBus5inverterwithsteadilydecreasingirradiance.................................28 5.6ReactivepoweroutputatBus5withsteadilydecreasingirradianceandreactivepowercontrolsignal.........................28 5.7RealandreactiveloadingatBus4......................29 5.8RealandreactiveloadingatBus6......................29vi
PAGE 7
5.9ReactiveoutputoftheinverteratBus3..................30 5.10ReactiveoutputoftheinverteratBus5..................30 5.11ReactivepowersuppliedbytheutilityatBus1withoutinvertersupport31 5.12ReactivepowersuppliedbytheutilityatBus1withinvertersupport..32vii
PAGE 8
TABLES Table 3.1SummaryofknownandunknownquantitiesinthePowerFlowproblem.11 3.2Summaryofknownunknown,andmonitoredquantitiesinthePowerFlowproblemmodiedforPVControl.......................18 4.1PowerFlowBusData.............................22 4.2PowerFlowBranchData...........................22 4.3PowerFlowGeneratorData.........................23viii
PAGE 9
1.Introduction Thisthesispresentsareactivepowercontrolstrategyforphotovoltaic(PV)resourcesthatareconnectedtoatypicaldistributionfeederthroughaDC{ACinverter.Theaimofthecontrolstrategyistoprioritizethetransferofrealpower,bututilizeanyunusedinvertercapacitytosupplyreactivepowerforvoltagesupportonaradialdistributionfeeder.Theamountofreactivepowertobesuppliedisbasedonapowerowsolutionthatassumestheentirenetworktopologyisknownateverymomentintime.Whilethiswasnotpreviouslypractical,advancementsinonlinemonitoringandutilitySystemControlandDataAcquisition(SCADA)capabilitiesisnarrowingthetheoreticalgap.Furthermore,asolutionbasedonaknownsystemcongurationcanserveasabenchmarkforthesameproblemwherestateestimationtechniquesareusedinlieuofactualquantities. Thischapterdiscussesthepotentialoperatingmodesofdistributedgenerationsuppliedpowerviaaninverteranddiscussestheeconomicincentivesofboththeutilityandthepowergenerator. 1.1ReactiveSupplyCapabilitiesandHistoricalLimitations Inthecaseofphotovoltaicsources,aninverterofsomecapacitywillalwaysbeassociatedwiththeconnectiontotheutilitygridtotransformtheDCpowersuppliedbythePVarraytoanACoutputthatcanbeconnectedtotheutilitysystem.ThecomplexpowerratingoftheinverterwilllikelyexceedthenameplateratingofthePVsystem.SeveralmaximumpowerpointalgorithmshavebeendevelopedtoadjusttheoutputoftheinvertertosupplyasmuchrealpoweraspossiblebyconsideringtheVIcurveofthesolararrayinquestion.Duringperiodsofreducedirradianceduetocloudcoverordarkperiodsoftheday,asignicantportionoftheinvertercapacitygoesunutilized.1
PAGE 10
Theabilityofinverterconnectedsourcestosupplyreactivepowerhashistoricallybeenrestrictedbyutilityinterconnectionguidelines.Thiswasdonelargelytoensurethattheinterconnectinggeneratorswouldnotbeaddingharmonicstotheutilitygrid.Thus,mostinterconnectingpartieshavebeenrequiredtomaintainvarneutralityatthepointofcommoncoupling(PCC).Whiletheworkpresentedinthisthesislargelypertainstoutilitydistributionsystems,advancementsinwindturbinesarehavingasimilareectatthetransmissionlevel.Thecapabilityofdoublyfedinductiongenerators(DFIG's)todynamicallysupplyreactivepowerisgivingutilitiestheoptionofspecifyinganonunitypowerfactorfortheinjectedpowertohelpcontrolsystemvoltages. Furthermore,theoptimalcontrolofanyreactivepowersourcewouldrequireonlineadjustmentofpowerlevelsasthesystemoperatingpointisconstantlychanging.FACTSdevicesonthetransmissionlevelhavebeenimplementedwithrealtimeSCADAsuppliedsystemdatatooptimizeoperation,butthatlevelofmeteringandcommunicationhashistoricallynotbeenavailableatthedistributionlevel.Switchedcapacitorsandinductorsaremostcommonlycontrolledlocallybymonitoringvarsatthepointofcommoncouplingorsimplyswitchedwhenvoltagelevelsexceedapresetbandwidth.WiththeproliferationofSmartGriddriveninitiativestoaddSCADAtomediumvoltage3phasefeedernetworks,smartcontrolofdistributedgenerationisbecomingmuchmorecommonplace. 1.2EconomicConsiderationsandOperatingModes UtilizingthestrandedcapabilityofthePVresourcestosupplyreactivepowerwillhaveadirecteconomicbenetfortheutility.Anyexcessinvertercapacitythatsuppliesreactivepowerasneededwouldosettheneedfortheutilitytosupplyswitchedorxedreactiveelementsonthedistributionsystem,andcouldpotentiallysaveoperationsofvoltageregulatingdevicessuchasloadtapchangersorvoltage2
PAGE 11
regulatorsastheyrespondtolowvoltageconditionsaspowerfactordriftsawayfromunity.Thedistributedgenerationprovideristypicallymeteredatthepointofcommoncouplingstrictlyonarealpowerbasisandrequiredtomaintainunitypowerfactoratalltimes.Whilethebenetstotheutilitypresentedabovehavearealcost,itwouldbedicultifnotimpossibletoattributethereactivecontributionofasingleenergyprovidertoavoidingadirectcostifreactivepowerwasbeingsuppliedbymultiplegeneratorsonagivenfeedersimultaneously.Toaccountforthisfact,thecontrolstrategypresentedinlaterchapterswillgivepreferencetosupplyingthemaximumamountofrealpowerandsupplyreactivepoweronlywithwhateverinvertercapacityremains.Thisrestrictionensuresthatthegenerationprovideriscompensatedforallrealpowerprovidedwithoutlosinganyrevenueforcontributingtothereactivecompensationofthesystem. Sincerealandreactivepoweraddinquadrature,themaximumreactivepowerthattheinvertercansupplyQinvisgivenby:Qinv=q S2inv)]TJ /F3 11.955 Tf 11.95 0 Td[(P2inv(1.1) wheretheSinvistheMVAratingoftheinverterandPinvistheactivepowerbeingsuppliedatanygivenmoment. Figure1.1representsthevariousoperatingmodesoftheinverterasacirclewitharadiusequaltoSinv[7]. ThecontrolstrategypresentedinthisthesiswillallowPinvtobedeterminedbyirradiance,temperature,andaMaximumPowerPointTrackingAlgorithm(MPPT)whiletheoptimalQinvwillberstcalculatedbyapowerowalgorithmbutlimitedtotheoperatingregionsshowningure1.1.TheinvertershouldoperateinthelefthandplanewherePinv<0undernormaloperationwhenpowerisbeingsuppliedtothegrid,thoughsomepowerwillbeabsorbedduringtheinitialenergizationtochargetheDCcapacitor.Bothpositiveandnegativevaluesofreactivepowercanbecommandedcorrespondingtoacapacitiveorinductiveeectrespectivelyasneeded3
PAGE 12
Figure1.1:OperatingRegionoftheVoltageSourceInverter.ifsucientinvertercapacityisavailable.4
PAGE 13
2.PhotovoltaicArrayandInverterModel ThischapterpresentsthemodelingapproachesforthePVarrayandinvertersystem.Thenalmodelwasabstractedandusedtorepresentindependentphotovoltaicgeneratorsconnectedtoautilitydistributionsystem.ThePVarraymodelincludesinputsfortemperatureandirradiancetoallowforrealworldconditionsthatwillchangetheoutputpowerofeachindividualdistributedgenerator.ThenalmodelwasrealizedinPSCAD. 2.1ModelofthePVCellandArray AlthoughthecircuitalmodelforasinglePVcellanditsgeneralizationtoanumberofcellsinseriesiswellestablishedintermsofacurrentsource,ananti{paralleldiode,aseriesresistanceandshuntresistance[5],themodelusedinthisthesismakesuseofamodiedcircuitthatreplacestheanti{paralleldiodebyanexternalcontrolcurrentsourceasproposedin[3].Thismodelwasexplainedin[6],anditisillustratedinFig.2.1(a).Achiefadvantageofthismodelisthatitallowsforincludinganarbitrarynumberofcellsconnectedinseriesand/orparallelintoasinglecircuitalrepresentationincludingalldetailsofeachcell. TheoutputcurrentofthePVcellofFig.2.1(a)maybeexpressedas,ipv=Iirr)]TJ /F3 11.955 Tf 11.95 0 Td[(Idio)]TJ /F3 11.955 Tf 11.95 0 Td[(Ip(2.1) whereIirristhephotocurrentorirradiancecurrentgeneratedwhenthecellisexposedtosunlight.Idioisthecurrentowingthroughtheanti{paralleldiodeandinducesthenonlinearcharacteristicsofthePVcell.IpisashuntcurrentduetotheshuntresistorRPbranch.SubstitutingrelevantexpressionsforIdioandIp,ipv=Iirr)]TJ /F3 11.955 Tf 11.96 0 Td[(I0expq(vpv+ipvRS) nkT)]TJ /F1 11.955 Tf 11.95 0 Td[(1)]TJ /F3 11.955 Tf 13.15 8.08 Td[(vpv+ipvRS Rp;(2.2)5
PAGE 14
Figure2.1:Equivalentcircuitfor(a)asinglePVcelland(b)aPVarrayofanarbitrarysize whereq=1:60210)]TJ /F4 7.97 Tf 6.59 0 Td[(19Cistheelectron'selectriccharge,k=1:380650310)]TJ /F4 7.97 Tf 6.58 0 Td[(23J/KistheBoltzmannconstant,Tisthetemperatureofthecell,I0isthediodesaturationcurrentorcellreversesaturationcurrent,nistheidealityfactorortheidealconstantofthediode,andRSandRPrepresenttheseriesandshuntresistance,respectively[5].Thismodelcanbegeneralizedtoanarbitrarynumberofcellsconnectedinseries,sayNS,andinparallel,sayNP,toformanarrayofNSNPcells.ThegeneralizedmodelisillustratedinFig.2.1(b). 2.2ModelfortheGridConnectedPVSystem ThecircuitschematicofFig.2.2illustratesthepowerstageofthesingle{stagegridconnectedPVpowerplantusedinsimulation.ItincludesthecascadeconnectionofadccapacitorthatconnectstotheoutputterminalsofanarbitrarysizePVarray;atwo{levelthree{phasevoltagesourceinverter(VSI);aLC{lterinserieswithacouplingtransformerthatconnectstotheacgrid. AssuggestedinFig.2.2,theoverallcontrolisrealizedthroughatwo{layerstrategy.AnouterlayertracksthePVarray'sMPPutilizingaMPPTalgorithm.An6
PAGE 15
Figure2.2:Circuitschematicofthepowerstageforthegrid{connectedPVsystem.innerloopregulatestheVSIincurrentcontrolmodeinordertotransferthepowergeneratedbythePVarrayintotheacgridatanarbitrarypowerfactor.TheabilitytooperateatanarbitrarypowerfactorenablesthePVsystemtodynamicallysupplyreactivepowertothegridascommandedbythecontroldescribedinChapter3. 2.3VoltageSourceInverterModel Thecascadedconnectionofthedccapacitor,VSI,LC{lterandacgridismodeledinthed)]TJ /F3 11.955 Tf 11.92 0 Td[(qreferenceframeusingdynamicphasors.Fig.2.3illustratesthedynamicequivalentcircuitintermsofcontrolsources[9,8,4]. Figure2.3:Controlsources{baseddynamicequivalentcircuitschematic. Thestatespaceequationsforthesystemshownin2.3aregivenby(2.3)through(2.6). d dtvpv=1 Cipv)]TJ /F1 11.955 Tf 13.15 8.09 Td[(3 2~m~i(2.3)7
PAGE 16
d dt~i=)]TJ /F3 11.955 Tf 9.3 0 Td[(j!~i+1 L(~mvpv)]TJ /F1 11.955 Tf 11.96 0 Td[((R+Rf)~i)]TJ /F3 11.955 Tf 11.95 0 Td[(Rf~ig)]TJ /F3 11.955 Tf 11.48 0 Td[(~vf)(2.4) d dt~vf=)]TJ /F3 11.955 Tf 9.3 0 Td[(j!~vf+1 Cf(~i)]TJ /F3 11.955 Tf 10.11 2.69 Td[(~ig)(2.5) d dt~ig=)]TJ /F3 11.955 Tf 9.3 0 Td[(j!~ig+1 Lg(~vf)]TJ /F1 11.955 Tf 11.96 0 Td[((Rg+Rf)~i+Rf~ig)]TJ /F3 11.955 Tf 11.48 0 Td[(~vg)(2.6) In(2.3)through(2.6),thevectorquantitiesrelatetotheirabccounterpartsconsideringagenericvectorquantity~x=xd+jxq,with xdxqT=TdqabcxaxbxcT;(2.7) whereTdqabcistheParktransformation.In(2.3),theoperator\"denotesthedotproductbetweentwovectors,yieldingadcvalue.In(2.3)and(2.4),~mistheVSI'smodulatingfunction. Inequations(2.4)through(2.6),!andarethesynchronizingfrequencyandphaseangle,respectively.ThesynchronizingquantitiesaremeasuredattheLC{lter'soutputandsynthesizedthroughaphase{locked{loop(PLL).ThePLLcomputesthegridphaseanglebysensingthegridvoltageandprojectingthecorrespondingspacevectoronthed{andq{axisofad)]TJ /F3 11.955 Tf 12.45 0 Td[(qrotatingframe.Thisd)]TJ /F3 11.955 Tf 12.45 0 Td[(qframeisthenrotatedinsuchawaywhichensuresthatthed{axisofthed)]TJ /F3 11.955 Tf 12.7 0 Td[(qframeisalignedwiththegridvoltagevector,thatis,vq=0.Insteady{state,thed)]TJ /F3 11.955 Tf 11.49 0 Td[(qframerotationalspeedequalsthegridangularfrequencyandtheextractedangleequalstothegridvoltageangle.Thisangleisusedtosynchronizethed)]TJ /F3 11.955 Tf 11.11 0 Td[(qreferenceframeforthecurrentcontrol. TheoverallcontrolstrategyfortheVSIisshowningure2.4.Asobservedinthegure,thecontrolisbrokendownintofourblocks:(i)aPLLutilizedfor8
PAGE 17
synchronizationwiththeacgrid;(ii)adcvoltagecontroltotrackthePVarray'sMPP;(iii)apower{to{currenttransformation;and(iv)aconventionalcurrentcontroltotransferthePVarray'spowerintotheacgridatanarbitrarypowerfactor. Figure2.4:Controlsources{baseddynamicequivalentcircuitschematic. ThevalueofQrefisspeciedtocompensateforthereactivepowerconsumedbytheLClterandstepuptransformerifstrictlyrealpoweroutputisdesired.Anyreactivepowercommandtobeinjectedbytheinverterisaddedtothesecompensatinglevels.9
PAGE 18
3.PowerFlowAlgorithm 3.1Introduction Thebasicaimofanypower/loadowcalculationsistodeterminethevoltageateverybuswithinautilitysystemandthenanalyzetheowofpowerthroughthegrid.Theresultscanshowtherequiredratingsofconductorsconnectingvariousgeneratorandcanbeusedtomodelfuturesystemloadgrowthand/orsystemexpansioneects.LoadowstudieswereoriginallyperformedusingsmallerscaleanalogmodelsofutilitysystemscalledcalculatorboardswhereRLCelementswereconnectedinthesametopologyastheutilitygrid.Advancementsincomputingmadeitpossibletodigitallysimulateandsolvethesystemofequationsthatrepresentstheactualutilitysystem. Theloadowproblemwasextendedtotheoptimalpowerow(OPF)problembyCarpentierin1962[1],whichappliesconstraintstocertainparameters(busvoltagesoranglesforstabilitypurposes,conductorratings)whileminimizingormaximizingadesiredobjectivefunctionsuchassystemlossesorgenerationmodelcosts. 3.2FormulationofthePowerFlowEquations Inordertosolveforthevoltagemagnitudeandangleateverybus,itisassumedthatthenetworktopologyconsistingofthepointsbelowareknown.1. ThetotalnumberofbussesinthesystemNrepresentingtheconnectionpointsforgenerators,loads,andswitchingstations,etc.2. ThetransmissiontopologythatconnectsthevariousbussesofthesystemincludingconductorRLCproperties.3. Thenumberofgeneratorsandtheirlocationwithinthesystem.10
PAGE 19
4. Thenumberandsizeofloadsservedbythesystemandtheirlocationwithinthesystem. Thebussesofthepowersystemarefurtherclassiedasfollows1. Ifatleastonegeneratorconnectedtoabus,thebusisdesignatedasaPVbuswheretherealpowerandbusvoltagemagnitudeareconsideredquantitiesthatcanbespeciedbythegeneratorcontrol.ThetotalnumberofPVbussesinthesystemisdesignatedasNPV,andthetotalnumberofgeneratorinthesystemisgivenbyNG.2. Oneofthegeneratorbussesisarbitrarilyidentiedasaslackbus.Thebusvoltageisassumed(bothmagnitudeandangle)areassumedtobeknown.Theslackbusvoltageangleisconsideredtobe0anddenesthereferenceforallotherbusvoltageanglesinthesystemdesignatedasN.Thepowerprovidedbytheslackbusismodeledtosupplysystemlossesonly.3. IfnogeneratorsareconnectedtoabusitisidentiedasaPQbusandistypicallyconsideredasystemload.Therealandreactivepowerleavingthesystemareassumedtobeknownbasedonloadprojectiondata.ThetotalnumberofPQbussesinthesystemisdesignatedasNPQ AsummaryoftheknownandunknownquantitiesinthepowersystemisgiveninTable3.1 Table3.1:SummaryofknownandunknownquantitiesinthePowerFlowproblem. RealPower ReactivePower VoltageMagnitude VoltageAngle PQbus known known unknown unknown PVbus controlled unknown controlled unknown slackbus unknown unknown controlled known 11
PAGE 20
TheunknownandcontrolledquantitiesofinterestarerepresentedbythevariablevectorX. X=0BBBBBBB@N)]TJ /F1 11.955 Tf 11.95 0 Td[(1jVNjPNGQNG1CCCCCCCA(3.1) ThedimensionsofXisequalto2(N+NG))]TJ /F1 11.955 Tf 11.96 0 Td[(1.Inordertondauniquesolution,wewillneedtowriteanequalnumberofequationsthatincludetheelementsofXandotherknownquantities.Thebasisforthesystemofequationswillbethecomplexpowerinjectedateverybus,Sbus,whichincludesbotharealandimaginarypartforeachelement.Sbuscanbedenesas: Sbus=VbusIbus(3.2) orrewrittenas: Sbus=Vbus(YbusVbus)(3.3) whereYbusisanadmittancematrixbetweenallthebussesofthesystemandtakestheform:12
PAGE 21
Ybus=0BBBBBBBBBBB@y11y12::y1ny21y22::y2n::::::::::yn1yn2::ynn1CCCCCCCCCCCA(3.4)Ybuscontainsthetopologyofthesystemandtheelectricalpropertiesofeachtransmissionlineorbranch.YbusisofsizeNxNandhasthefollowingproperties. Ybusissymmetrical. Thediagonalelements,Yii,areequaltothesumoftheadmittancesofallelementsconnectedtotheithnode. Theodiagonalelements,Yij,areequaltothenegativeoftheadmittanceconnectingnodesiandj. TheithelementofSbuscanbewrittenas Sbusi=Vbusi"nXj=1(yijVbusj)#i;j=1;2;:::;n(3.5) or Sbusi=nXj=1(VbusiVbusjyij)i;j=1;2;:::;n(3.6) andnallyas Sbusi=nXj=1jVijjVjjejij(Gij)]TJ /F3 11.955 Tf 11.95 0 Td[(jBij)i;j=1;2;:::;n(3.7) where:13
PAGE 22
Yij=Gij+jBijij=i)]TJ /F3 11.955 Tf 11.96 0 Td[(ji;j=1;2;:::;n(3.8) Equation(3.7)givesanexpressionforthecomplexpowerinjectedateverybuswithinthesystemandcanbeseparatedintorealandimaginarypartsbyapplyingEuler'sformulaandgroupingliketermsasshowinequations(3.9)through(3.12) Sbusi=nXj=1[jVijjVjj(cosij+jsinij)(Gij)]TJ /F3 11.955 Tf 11.96 0 Td[(jBij)]i;j=1;2;:::;n(3.9) Sbusi=nXj=1[jVijjVjj(Gijcosij+Bijsinij+jGijsinij)]TJ /F3 11.955 Tf 11.96 0 Td[(jBijcosij)](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. Thenalsystemofequationsiswrittenbyobservingthatallcomplexpowerintothesystemmustbeequaltothecomplexpoweroutofthesystempluslossestopreserveenergybalance.Sincetheslackbusaccountsforlosses,equationsforrealandreactivepoweratallotherbussescanbewrittenasafunctionoftheunknownvectorXasfollows: 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)14
PAGE 23
where:Pgi:Theactivepowerofthegeneratorinthebusi.Qgi:Thereactivepowerofthegeneratorinthebusi.PLi:Theactivepoweroftheloadinthebusi.QLi::Thereactivepoweroftheloadinthebusi. Withtheformulationofthepowerowequationscomplete,thesystemcanbesolvediterativelytodeterminetheloadowsinthesystem.ApplyingconstraintsonthevariablevectorXpertainingtosystemstabilityandcapacitywilllimittheresultstoapracticaloperatingcondition.Typicalsystemconstraintsinclude:1. Limitationsonbusvoltagemagnitude,jVminjjVijjVmaxj,toensuresafeandeectivepowerdelivery2. Limitationsonvoltageangle,minimax,toensuresystemstability3. Limitationsonbranchpowerowtolessthantheratingsofthecircuit,SijSijRating,topreventoverloading. Thenalstepoftheformulationistoimposeanoptionalobjectivefunctiontominimizeormaximizeinordertoachieveanoptimalpowerowfromthesolutionset.Inthecaseofthissimulation,theobjectivewillbetominimizesystemlosses.Inlargertransmissionnetworks,theobjectivefunctionisusuallyrepresentativeofthecostofdispatchingthegeneratorsonthesystem,which,whenminimized,resultsinthedispatchingoflowestcostthatsuppliesthenetwork. 3.3ModiedFormulationforReactivePowerControlofPVResources Theformulationofthepowerowequationsfortheworkpresentedinthisthesisisaspecialcaseoftheconceptspresentedintheprevioussection,butthefundamentalconceptsremainthesame.Thissectiondiscussesthedierencesintheformulationofthepowerowequations.Thesystemconstrainsimposedonthesolutionandoverall15
PAGE 24
objectivefunctionforoptimalcontrolarealsopresented.ThenalsetofequationswereprogrammedintoMATLAB,andthecompletecodecanbefoundinAppendixA. 3.3.1FormulationofthePowerowequations Asstatedpreviously,theelectricalsystemofinterestedforthesimulationpresentedinthisthesisisaheavilyloadedradialdistributionfeederwithmultiplePVdistributedgeneratorsinoperation.Theprimaryobjectivesofthecontrolalgorithmareasfollows.1. SupplyrealpowerfromthePVgeneratorswheneverpossiblebasedoninstantaneousirradianceandtemperatureconditions.2. Supplyreactivepoweronlywheninvertercapacityisavailable.3. Theamountofreactivepowertobesuppliedwillbedeterminedbytheoptimalpowerowsolution,butnottoexceedthecapacityoftheinverter. Withtheseobjectivesinmind,itisclearthatthesystembuswherethePVgeneratorsconnectcannotbetreatedasatraditonalPVbusintheformulationofthepowerowequationssincetherealpowerinjected,Pgi;dependsontemperatureandirradianceandcannotbeexternallycommanded.Inthespecialcaseofaradialdistributionfeeder,theonlytruePVbusisthesubstationbusthatfeedstheentiresystem.Thevectorofunknownquantitiesdenedinequation(3.1)canbemodiedasshownin(3.14).ThedimensionsofX0isconsequentlyreducedto2(N+1)+NG)]TJ /F1 11.955 Tf 11.96 0 Td[(1.TheactivepowerinjectionofeveryPVgeneratorbecomesaninputtothepowerowalgorithmasopposedtoaquantitytobesolvedfor.16
PAGE 25
X0=0BBBBBBB@N)]TJ /F1 11.955 Tf 11.96 0 Td[(1jVNjPSubstationQNG1CCCCCCCA(3.14) ItisworthnotingthatthisapproachwillrequiretheactivepoweroutputofthePVgeneratortobemonitoredconstantlyandrelayedviaaSCADAsystemtoenableonlinecomputationofthereactivepowercommand.TheinstantaneousrealpoweroutputofthePVgeneratorislikelytobeavailablelocallyaspartoftheinvertercontrolandmonitoringsystem,butaSCADAlinkfordistributedgeneratorsisusuallyonlyrequiredforlargerinstallationstoaddressislandingandsynchronizationconcerns. Anotherkeydistinctionbetweentheclassicpowerowformulationandthespecialcasepresentedisthisthesisistheneedtorespondquicklytochangesinload.Whileloadprolesatatransmissionlevelmaybeverypredictablebasedonhistoricalloadproles,weather,andtimeofyear,distributionfeederloadscanvaryquicklyascommercialandresidentialcustomerschangeelectricconsumption.InordertocalculatetheoptimalamountofreactivepowertocommandfromthePVgenerators,therealandreactivepowerconsumptionateveryloadbuswillneedtobemonitoredandsuppliedtothecontrolalgorithmonline.Thenallistofquantitiesthatwillhavetobemonitoredonlineandsuppliedtothereactivepowercontrollerasinputsisasfollows.1. TherealpowergeneratedbyeveryPVdistributedgenerator.AssumingtherearenopowersourcesotherthanPVdistributedgeneratorandtheutility,thetotalnumberofinputswillequalNPV)]TJ /F1 11.955 Tf 11.96 0 Td[(12. TherealpowerconsumedateveryPQloadbus,whichimpliesNPQsignals.17
PAGE 26
3. ThereactivepowerconsumedateveryPQloadbus,whichimpliesNPQsignals. Asaresultoftheneedtomonitortheabovementionedquantitiesonline,table3.1canbemodiedasshownintable3.2toredesignatesomeofthequantitiespreviouslyconsideredasknown. Table3.2:Summaryofknownunknown,andmonitoredquantitiesinthePowerFlowproblemmodiedforPVControl. RealPower ReactivePower VoltageMagnitude VoltageAngle PQbus monitored monitored unknown unknown PVSubstationbus controlled unknown controlled unknown PVInverterbus monitored unknown controlled unknown slackbus unknown unknown controlled known OnlinefeedermonitoringisakeyaimofSmartGridtopicsandisbecomingmuchmoreprevalent.Oneofthechallengesofrealworldimplementationisthedistributednatureoffeederloadsasopposedtodiscreteloadbusses.Severalstateestimationtechniquesforlumpedequivalentshavebeendeveloped,buttheaimofthisthesisistoprovideanoptimalbenchmarkassumingalloftheelectricaldataisknown. 3.3.2ObjectiveFunctionandSystemConstraints Theobjectivefunctionchosenforthesimulationpresentedinthisthesisisminimizationofthesystemlosses.Giventhatthesimulatedsystemisaradialdistributionfeeder,allcomplexpowerisdeliveredtothesystemfromtheutilityatthesubstationfeederexitplusthecontributionofthedistributedgenerators.TheactivepowercontributionofthedistributedgeneratorstakepriorityoverthereactiveoutputandiscontrolledbyaMPPTalgorithm.Sincetheutilityhastosupplyallloadnotservedbythedistributedgenerators,minimizationofthetransmissionlosseswillprotidethemostcosteectiveoperatingpoint.18
PAGE 27
Thesystemconstrainsimposedonthesolutiontoensurearealisticsolutionareasfollows:1. Limitationsonbusvoltagemagnitudeinperunit,0:95jVnominaljjVijj1:05Vnominalj.2. Limitationsonvoltageangle,)]TJ /F1 11.955 Tf 9.3 0 Td[(50i503. Limitationsonthecomplexpowersuppliedbytheinverter,Sinv=500kVA Aspreviouslydiscussed,theinverter'scomplexpowerratingSinvisthequadratureadditionoftherealandreactivepowersuppliedbytheinverter.Therealpowercomponentiscontrolledstrictlybythemaximumpowerpointtrackingalgorithmasafunctionoftemperatureandirradiance.Themaximumreactivepowercomponentiscomputedaccordingtoequation(1.1)whereSinvislimitedto500kVARandPinvisreadasaninputtothecontroller.19
PAGE 28
4.PowerSystemSimulation Inordertotestthereactivepowercontrolalgorithm'sabilitytorespondtoonlinechangesinsystemconditionsandcommands,amodelofadistributionfeederwasdevelopedinPSCAD.ThePVarrayandvoltagesourceinvertermodelpresentedinChapter2wereabstractedasPSCADcomponentsandappliedtothesimulateddistributionfeederasdistributedgenerators.APSCADcomponentwasdevelopedtopasstheonlinemonitoredsignalstotheMATLABpowerowalgorithmandalsotoreturnthereactivepowercommands.Finally,sometimedbreakerlogicwasaddedtodemonstratethecontroller'sabilitytodealwithchangesinsystemloading 4.1DistributionFeederModel Aspreviouslystated,thescopeofthesimulationpresentedinthisthesisisaradialmediumvoltagedistributionfeeder.Themodelforthefeederisshowin4.1and4.2.Notethatthemodelhasbeensplitintotwopartsforclarity.Bus3isshowninbothgures. Thethreephasesystemisassumedtobebalancedanddoesnotincludelateraltiestootherdistributionfeeders.ThesubstationfeederexitisshownconnectedtoBus1,thePVdistributedgeneratorsareshownconnectedatBus3andBus1.LumpedloadsareshownconnectedbetweenallofthepowersourcesatBus2,Bus4,andattheendofthefeederatBus6.Abranchimpedancebetweenallbussesof0:1+j0:05ohmsisincludedaswell. Inordertosimulatetypicalutilityvoltageregulation,thevoltagesourceconnectedatBus1thatsimulatesthesubstationexitismodeledasacontrolledvoltagesource.Themagnitudeofthestartingfeedervoltageiscontrolledbythepowerowalgorithmwithinatoleranceof+/5%,whichistheANSIlimitforutilityvoltagetolerance.20
PAGE 29
Figure4.1:FeederModelPart1 Figure4.2:FeederModelPart2Distributionvoltageregulationistypicallyachievedbythemechanicaloperationofaloadtapchangerorexternalvoltageregulatoronthesecondarysideofthesubstation21
PAGE 30
powertransformer. Inordertosupplythereactivepowercontrolalgorithmwiththestaticknownsystemparameters,atextlewasdevelopedthatorganizesthesystemparametersintheIEEEstandardformatforpowerowdata[2].Thedataofinterestforthesimulatedfeederisshownintables4.1through4.3 Table4.1:PowerFlowBusData. Name Type Voltage(V) MWGen MvarGen MWLoad MvarLoad Bus1 Slack 12470 0 0 0 0 Bus2 PQ 12470 0 0 5 1.7 Bus3 PV 12470 0 0.4 0 0 Bus4 PQ 12470 0 0 3 1.0 Bus5 PV 12470 0 0.4 0 0 Bus6 PQ 12470 0 0 3 1.7 Table4.2:PowerFlowBranchData. From To Resistance(ohms) Reactance(ohms) Bus1 Bus2 0.1 0.05 Bus2 Bus3 0.1 0.05 Bus3 Bus4 0.1 0.05 Bus4 Bus5 0.1 0.05 Bus5 Bus6 0.1 0.05 4.2ControlSystemInterface WhilethePVArray,voltagesourceinverter,andthefeederitselfareallsimulatedinPSCAD,thepowerowequationsanditerativesolutionareperformedinMATLAB.AcustomcomponentinPSCADwasdevelopedtohandletheexchange22
PAGE 31
Table4.3:PowerFlowGeneratorData. Location PMax(MW) PMin(MW) QMax(MVar) PMin(MVar) Bus1 20 20 20 20 Bus3 0.5 0.5 0.5 0.5 Bus5 0.5 0.5 0.5 0.5 ofdatabetweenthetwoprograms.Asoutlinedinthepreviouschapter,thefollowingdatahastobemonitoredonlineandfedtothecontrollerasaninput.1. SupplyrealpowerfromthePVgeneratorswheneverpossiblebasedoninstantaneousirradianceandtemperatureconditions.2. Supplyreactivepoweronlywheninvertercapacityisavailable.3. Theamountofreactivepowertobesuppliedwillbedeterminedbytheoptimalpowerowsolution,butnottoexceedthecapacityoftheinverter. TheoutputoftheMATLABpowerowalgorithmisthereactivepowercommandstoeachvoltagesourceinverter,whichisalsopassedthroughtheinterfacecomponent.Thegraphicalrepresentationoftheinterfaceasappliedtothedistributionfeedermodelisshowin4.3. TheFORTRANcodethatgovernsthecomponentsdatatransferisshowninAppendixB.NotethatatriggeringinputisincludedinthecomponentthatgovernshowfrequentlyMATLABwillexecutethecontrolalgorithm.Forthepurposesofthetestcasespresentedinthisthesis,aperiodof0.25secondswasselected. 4.3TimedBreakerLogic Inordertodemonstratethecontroller'sabilitytorespondappropriatelytochangesinfeederload,timedbreakerlogicwasaddedtothePQbusses,Bus4andBus6asshowningure4.2.23
PAGE 32
Figure4.3:PSCADtoMATLABInterfaceComponent Thebreakerlogicisintendedtochangetheloadstatequicklyinbetweenexecutionsofthereactivepowercontrolalgorithmsothatthesubsequentreactivepowercommandcanbeevaluated.Thedetailsofeachbreakeroperationwillbedetailedforeachtestcaseinthenextchapter.24
PAGE 33
5.TestCasesandSimulationResults ThemodelsforthePVarray,voltagesourceinverter,anddistributionfeederweremodeledinPSCADandveriedbysimulation.ThebehavioroftheMATLABpowerowalgorithmcanalsobeobservedbyconsideringthedistributionfeedermodel,theinputparametersthataremonitoredonline,thereactivepowercommandsissuedfromthecontroller,andthenalsystemresponse.Severaltestcasesarepresentedinthischapterinorderofleasttomostcomplex. 5.1TestCaseConstantReactivePowerOutput InordertorstverifythesimulationofthePVarrayandvoltagesourceinverter,therstcasepresentedissimplyaconstantreactivepowercommand.Acaseofparticularinterestisoperationatunitypowerfactor.AconstantreactivepowercommandofzeroissuppliedtothevoltagesourceinvertersconnectedatBus3andBus5forthesystemshowningures4.1and4.2.Additionally,thePVarrayatBus3willbesuppliedwithalargeirradiancevaluetosimulatefullsunexposure,whilethePVarrayatBus5willbesuppliedwithaminimalirradiancevalue. Figure5.1:RealandreactivepoweroutputoftheBus3inverterunderhighirradianceandzeroreactivepowercommand25
PAGE 34
Figure5.2:RealpoweroutputoftheBus5inverterunderlowirradianceandzeroreactivepowercommand Figure5.3:ReactivepoweroutputoftheBus5inverterunderlowirradianceandzeroreactivepowercommand WecanobservethattherealpoweroutputofthePVGeneratoratBus3thatwassubjectedtoahighirradiancevalueisnearlyequaltothemaximuminverterratingSinvof500kVA,whiletherealpoweroutputofthethePVGeneratoratBus5wasessentiallyzero.Forbothcases,thereactivepoweroutputwasessentiallyzero.Bothresultsconrmthattheinvertmodel'srealpoweroutputisdependantonirradiance,andadesiredreactivepowervaluecanbesuccessfullycommandedfromtheMATLABalgorithm. 5.2TestCaseVariableIrradiance Anothertestcaseofinterestistherealandreactivepowerresponseofthesystemtosteadilyincreasingordecreasingirradiancewhilethefeederisheavilyloaded.Thegoalofthissimulationistodemonstratethefollowing:26
PAGE 35
1. RealpowersuppliedbythePVarrayisdirectlyproportionaltoirradiance.2. ThemagnitudeofreactivepowersuppliedbythePVarrayisindirectlyproportionaltoirradiance.3. ThemagnitudeofthereactivepowersuppliedislimitedbytheinverterratingSinvasexpressedinequation(1.1)4. Thesystemisabletorespondtothequantitiesmonitoredonlineandsubsequentchangingreactivepowercommands. Thesimulationsresultsarepresentedingures5.4through5.6 Figure5.4:RealandreactivepoweroutputoftheBus3inverterwithsteadilyincreasingirradiance Wecanseeingure5.4theexpectedresponsetosteadilyincreasingirradiance,andthequadraturerelationshipbetweenPinvandQinv.Thesameeectisevenmorepronouncedingure5.5whererealpowerquicklydropsoastheappliedirradiancedecreases. Theonlinetrackingofinputparametersandresponsetochangingreactivepowercommandsisillustratedin5.6,whichplotsthereactivepowercommandissuedfromtheMATLABcontrolalgorithmontopofthemeasuredreactivepoweroutput. WecanclearlyseethattheMATLABinterfacecomponentisexecutingthereactivepowercontrolalgorithmandprovidinganupdatedcontrolsignalever0.25secondsasexpectedandthesubsequentresponseoftheinvertercontrol.27
PAGE 36
Figure5.5:RealandreactivepoweroutputoftheBus5inverterwithsteadilydecreasingirradiance Figure5.6:ReactivepoweroutputatBus5withsteadilydecreasingirradianceandreactivepowercontrolsignal 5.3TestCaseVariableLoad Thenexttestcasepresentedisintendedtohighlightthecontroller'sabilitytorespondtochangesinload.TimedbreakerlogicisappliedtotheloadsconnectedatBus4andBus6tosimulatecapacitorswitchingandabruptchangesinloadthatcanoccuronadistributionsystemduetocustomerinteractionorprotectiveequipmentoperation.Thetimedbreakeroperationsaredesignedtoimposethefollowingconditions. TheinitialloadatBus4is3MWand1.7Mvarinductive. TheinitialloadatBus6is3MWand1.7Mvarinductive. Att=0:6sec,theloadatBus6willquicklychangeto3MWand1.2Mvarcapacitivetosimulateacapacitorbeingswitchedon.28
PAGE 37
Att=1:1sec,theloadatBus4willquicklychangeto1MWand0MVartosimulateafuseoperation. TheloadingatBus4andBus6overthesimulationperiodareshowningures5.7and5.8 Figure5.7:RealandreactiveloadingatBus4. Figure5.8:RealandreactiveloadingatBus6. TheirradianceappliedtobothinverterswillbeheldatpracticallyzerotoenabletheinvertertosupplyreactivepoweralmostcompletelyuptotheinverterratingSinv=500kVA. Thesystemresponseisshowningures5.9and5.10 Thefollowingobservationscanbemadeateach0.25secondintervalwhenthecontrollerrecalculatesthereactivepowercommands:29
PAGE 38
Figure5.9:ReactiveoutputoftheinverteratBus3 Figure5.10:ReactiveoutputoftheinverteratBus5 Att=0:25secandt=0:5sec,bothoftheinvertersatBus3andBus5aresupplyingmaximumreactivepowerduetotheheavyloadingatnearbybusses. Att=0:75sec,thecontrolalgorithmiscalledforthersttimesincetheloadatBus6haschangedsignicantlyfromthepreviousstate.LessreactivepowerisrequiredtofullycompensatetheloadandlinelossintherestofthefeederfromBus6. Thereisnosignicantchangeofstateatt=1:0sec. Att=1:25sec,thecontrolalgorithmiscalledforthersttimesincetheloadatBus4haschangedsignicantlyfromthepreviousstate.TheloadatBus430
PAGE 39
hasdroppedtothepointthatminimalreactivepowerisrequiredatBus3,andanexcessofreactivepowerattheofthefeederhascausedtheinverteratBus5toentertheinductiveoperatingrangeandstartabsorbingvars. 5.4MeasurementofObjectiveFunction Inordertodemonstratethattheobjectivefunctionisbeingprobablyrealized,analtestcasesispresented.First,ifthereactivepowercommandsissuedtotheinverteraredesignedtominimizesystemlosses,anyloaddemandforreactivepowershouldbepasseddirectlytotheinverterandsupplieduptotheinverter'scapacity.Toillustratethis,considerthesamecircuitshownin4.1and4.2,withoutanytimedbreakeroperationsandalowvalueofirradianceappliedtoeveryPVsystem.Theintentistomodelaheavilyloadedgurewiththedistributedgeneratorsabletosupplyreactivepoweronly.Asabaseline,thereactivepowersuppliedbythemodeledutilitysystemaBus1isshowningure5.11. Figure5.11:ReactivepowersuppliedbytheutilityatBus1withoutinvertersupport Whenthereactivepowercontrolalgorithmisenabled,thereactivepowersuppliedbythesystemisshowningure5.12. Asexpected,wecanobserveroughlya1Mvarreductionintheamountofreactivepowerthattheutilitysystemhastosupply.31
PAGE 40
Figure5.12:ReactivepowersuppliedbytheutilityatBus1withinvertersupport32
PAGE 41
6.Conclusion ThisthesishaspresentedacomprehensivemodelofadistributionfeederwithmultipleindependentPVdistributedgenerators,andacontrolstrategyfortheassociatedvoltagesourceinvertersthatiscapableofprovidingreactivecompensationtotheutilitysystemwithoutcompromisingtheamountofrealpowergeneratedbythePVsystems.Inotherwords,thecontrolstrategytakesadvantageofanyavailableinvertercapacitytohelpregulatethesystemwithoutimpactingtheprivatelyowneddistributedgenerator. Thechallengesinrealizingthissysteminpracticeischieytherequiredonlinemonitoringandadjustmentofrealpowerinjectionfromthevariousgeneratorsandloadconditions.Priortotheimplementationofcellularcommunicationsandtheproliferationofberopticcommunicationcircuits,itwouldnothavebeenpossibletoachievethislevelofcontrol.FiberopticsarebeingincreasinglyinstalledonutilitycircuitseitheraspartoftheneutralconductororseparatelystrungonthesameoverheadpolesastheprimaryconductorsaspartoftheSmartGridinitiative.Evenwiththeadvancesincommunicationcapabilitiesonutilitysystems,theabilitytomonitorloadconditionscomprehensivelyandinstantaneously,especiallyatthedistributionlevel,willlikelyremainchallenging. Theresultspresentedinthisthesiscouldbeusedasabenchmarkforfutureworkonsimilarsystemsthatassumelessthanperfectknowledgeoftheutilitysystemateverymomentintimeandapplystateestimationtechniquestoextrapolatetheunknownquantitiesfromparticulardatapoints.Additionalfunctionalitycouldalsobeaddedtothereactivepowercontrolalgorithmtodealwithinterruptionstotheowofonlineinputparameterswhenlossofcommunicationoccurs.Sincethevoltagesourceinverteriscapableofgeneratingwaveformsabovethepowerfrequency,control33
PAGE 42
couldbeaddedtoprovideactivelteringofharmonicsintroducedbynonlinearloadsaswell.Thisisreadilyachievablesincethecontroloftheinverterisalreadyperformedinthed)61()]TJ /F3 11.955 Tf 19.34 0 Td[(qreferenceframe.Finally,moresophisticatedmodelingofthedistributionfeederitselfcouldleadtomoreinterestingtestcases.34
PAGE 43
REFERENCES[1] J.Carpentier.Contributional'etudedudispatchingeconomique.BelletindelaSocieteFrancaisedesElectricens,3(8):431{447,January1962.[2] M.Cerroni,F.MancillaDavid,A.Arancibia,F.RigantiFulginei,andE.Muladi.Commonformatforexchangeofsolvedloadowdata.PowerAppuratusandSystems,IEEETransactionson,PAS92(6):1916{1925,November1973.[3] M.Cerroni,F.MancillaDavid,A.Arancibia,F.RigantiFulginei,andE.Muladi.Amaximumpowerpointtrackervariabledclinkthreephaseinvertercontrolschemeforgridconnectedpvpanels.InInnovativeSmartGridTechnologies(ISGTEurope),20123rdIEEEPESInternationalConferenceandExhibitionon,pages1{7,2012.[4] R.Kadri,J.P.Gaubert,andG.Champenois.Animprovedmaximumpowerpointtrackingforphotovoltaicgridconnectedinverterbasedonvoltageorientedcontrol.IEEETransactionsonIndustrialElectronics,58(1):66{75,January2011.[5] G.M.Masters.RenewableandEcientElectricPowerSystems.JohnWiley&Sons,NewJersey,2004.[6] H.Tian,F.MancillaDavid,K.Ellis,E.Muljadi,andP.Jenkins.Acell{to{module{to{arraydetailedmodelforphotovoltaicpanels.SolarEnergy,86(9):2695{2706,September2012.[7] K.Turitsyn,P.Sulk,S.Backhaus,andM.Chertkox.Localcontrolofreactivepowerbydistributedphotovoltaicgenerators.InDistributedControlofReactivePowerFlowinaRadiallyDistributionCircuitwithHighPhotovoltaicPenetration,pages1{7,2012.[8] Y.Yang,M.Kazerani,andV.H.Quintana.Modeling,controlandimplementationofthree{phasepwmconverters.IEEETransactionsonPowerElectronics,18(3):857{864,May2003.[9] A.YazdaniandR.Iravani.VoltageSourcedConvertersinPowerSystems.Wiley{IEEEPress,Reading,MA,2010.35
PAGE 44
AppendixA.MATLABReactivePowerControlFunction TheClanguagecodeforthereactivepowerowalgorithmdevelopedinMATLABispresentedinthisappendix.ThecodeisorganizedasthemainfunctioncalledbyPSCADduringthesimulationandalsoincludesthereferencedsubfunctions. A.1MainFunctionQCTL %OutputQCMDreactivepowercommands. %InputPGrealpoweroutputofeachinverter, %InputPLrealpowerconsumedateachPQBus %InputQLreactivepowerconsumedateachPQBus function[Qcmd,Vs]=Qctl(PG,PL,QL) %Initialsettings j=sqrt(1); S.base.S=1000000; S.base.Vg=12470; S.base.Zg=S.base.Vg^2/S.base.S; S.base.Ig=S.base.S/S.base.Vg/sqrt(3); S.base.Vi=270; S.base.Zi=S.base.Vi^2/S.base.S; S.base.Ii=S.base.S/S.base.Vi/sqrt(3); %36
PAGE 45
%%ThisinitializesS S.fname=''; %%Thiscallsthefucntionread_cf S=read_cf(S); %saveSfile %Updateinputvaluesmonitoredonline. S.Bus.MWGen(S.Bus.PV)=PG; S.Bus.MWLoad(S.Bus.PQ(2:length(S.Bus.PQ)))=PL; S.Bus.MVARLoad(S.Bus.PQ(2:length(S.Bus.PQ)))=QL; %BuildtheYbusmatrix S=build_ybus(S); %Variables: %1)Anglesateverybus(excludingslack) %2)Voltagesateverybus %4)ActivePoweratsubstation %5)ReactivePowerateverygenerator(PV+substation) na=S.Bus.n1; nv=S.Bus.n; np=1; nq=length(S.Gen.Bus);37
PAGE 46
tetamin=50; tetamax=50; Vmin=.95; Vmax=1.05; %lowerbound LB=[tetamin*ones(na,1)*pi/180;Vmin*ones(nv,1);inf;inf;ones(nq1,1)]; %upperbound UB=[tetamax*ones(na,1)*pi/180;Vmax*ones(nv,1);inf;inf;ones(nq1,1)]; %Initialcondition %x0=[zeros(na,1);ones(nv,1);0.5*(S.Gen.Pmin+S.Gen.Pmax)/100;0.5*(S.Gen.Qmin+S.Gen.Qmax)/100]; x0=[zeros(na,1);ones(nv,1);14000000/S.base.S;1/S.base.S*[3000000;500000;500000]]; A=[]; b=[]; Aeq=[]; beq=[]; %Definesoptimset OptOpf=optimset('Display','iter','Diagnostics','on','DerivativeCheck','on',... 'LargeScale','off','GradConstr','off','GradObj','off','Hessian','off',...38
PAGE 47
'MaxIter',5000,'MaxFunEvals',100000); [x,fval,exitflag,output,lambda]=fmincon(@ObjFun,x0,A,b,Aeq,beq,LB,UB,@NonLinCon,OptOpf,S); %computenumberofeachvariable na=S.Bus.n1; nv=S.Bus.n; np=1; nq=length(S.Gen.Bus); %UpdatesVoltages Vbus=x(na+1:na+nv);%added i_av=find(S.Bus.type~=3); %Vbus(i_av)=x(na+1:na+nv); Vbus(S.Bus.Swing)=Vbus(S.Bus.Swing)*exp(j*0); Vbus(i_av)=abs(Vbus(i_av)).*exp(j*x(1:na)); Sinj=Vbus.*conj(S.Ybus*Vbus); Qcmd=x(na+nv+np+2:length(x))'; Vs=x(na+1); end A.2SubFunctionreadCF39
PAGE 48
TheintentofthissubfunctionistobuildaMATLABstructurewithallofthestaticbusandbranchdataofthepowersystemneededtobuildthepowerowequations. functionS=read_cf(S) %getsthenameofthefile fname='JRTFeeder.cf'; S.fname=fname; %opensthefileforreading fcf=fopen(fname,'r'); %getsthenameofthecase s=fgetl(fcf); whilestrcmp(s,''),s=fgetl(fcf);end; S.CaseName=s; %readsbusdata whilestrcmp(s(1:min(3,length(s))),'BUS')~=1,%Findthestartofbusdata s=fgetl(fcf); end s=fgetl(fcf); whilestrcmp(s,''),s=fgetl(fcf);end;40
PAGE 49
S.Bus.number=[]; S.Bus.name=[]; S.Bus.area=[]; S.Bus.zone=[]; S.Bus.type=[]; S.Bus.V=[]; S.Bus.angle=[]; S.Bus.MWLoad=[]; S.Bus.MVARLoad=[]; S.Bus.MWGen=[]; S.Bus.MVARGen=[]; %Extrastuffadded S.Bus.BasekV=[]; S.Bus.DesiredVolt=[]; S.Bus.MVarMax=[]; S.Bus.MVarMin=[]; S.Bus.G=[]; S.Bus.B=[]; whiles(1)~='' k2=1; [k1k2]=find_k1k2(s,k2);41
PAGE 50
n_aux=str2num(s(k1:k2)); S.Bus.number(n_aux,1)=n_aux; [k1k2]=find_k1k2(s,k2); S.Bus.name=strvcat(S.Bus.name,s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Bus.area(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Bus.zone(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Bus.type(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Bus.V(n_aux,1)=str2num(s(k1:k2))/S.base.Vg; [k1k2]=find_k1k2(s,k2); S.Bus.angle(n_aux,1)=str2num(s(k1:k2))*pi/180; [k1k2]=find_k1k2(s,k2);42
PAGE 51
S.Bus.MWLoad(n_aux,1)=str2num(s(k1:k2))/S.base.S; [k1k2]=find_k1k2(s,k2); S.Bus.MVARLoad(n_aux,1)=str2num(s(k1:k2))/S.base.S; [k1k2]=find_k1k2(s,k2); S.Bus.MWGen(n_aux,1)=str2num(s(k1:k2))/S.base.S; [k1k2]=find_k1k2(s,k2); S.Bus.MVARGen(n_aux,1)=str2num(s(k1:k2))/S.base.S; [k1k2]=find_k1k2(s,k2); S.Bus.BasekV(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Bus.DesiredVolt(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Bus.MVarMax(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Bus.MVarMin(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Bus.G(n_aux,1)=str2num(s(k1:k2)); 43
PAGE 52
[k1k2]=find_k1k2(s,k2); S.Bus.B(n_aux,1)=str2num(s(k1:k2)); s=fgetl(fcf); end %Findindicesforvarioustyoeofbuses S.Bus.PQ=find(S.Bus.type==0); S.Bus.PV=find(S.Bus.type==2); S.Bus.Swing=find(S.Bus.type==3); S.Bus.n=max(S.Bus.number); %readsbranchdata whilestrcmp(s(1:min(6,length(s))),'BRANCH')~=1,%Findthestartofbusdata s=fgetl(fcf); end s=fgetl(fcf); whilestrcmp(s,''),s=fgetl(fcf);end; S.Branch.From=[]; S.Branch.To=[]; S.Branch.area=[]; S.Branch.zone=[];44
PAGE 53
S.Branch.circuit=[]; S.Branch.type=[]; S.Branch.r=[]; S.Branch.x=[]; S.Branch.b=[]; S.Branch.MVA1=[]; n_aux=0; whiles(1)~='' n_aux=n_aux+1; k2=1; [k1k2]=find_k1k2(s,k2); S.Branch.From(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Branch.To(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Branch.area(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Branch.zone(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Branch.circuit(n_aux,1)=str2num(s(k1:k2));45
PAGE 54
[k1k2]=find_k1k2(s,k2); S.Branch.type(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Branch.r(n_aux,1)=str2num(s(k1:k2))/S.base.Zg; [k1k2]=find_k1k2(s,k2); S.Branch.x(n_aux,1)=str2num(s(k1:k2))/S.base.Zg; [k1k2]=find_k1k2(s,k2); S.Branch.b(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Branch.MVA1(n_aux,1)=str2num(s(k1:k2))/S.base.S; s=fgetl(fcf); end S.Branch.n=n_aux; whilestrcmp(s(1:min(9,length(s))),'GENERATOR')~=1,%Findthestartofgeneratordata s=fgetl(fcf); end s=fgetl(fcf);46
PAGE 55
whilestrcmp(s,''),s=fgetl(fcf);end; ifstrcmp(s(1),'('),s=fgetl(fcf);end; whilestrcmp(s,''),s=fgetl(fcf);end; S.Gen.Bus=[]; S.Gen.Pmax=[]; S.Gen.Pmin=[]; S.Gen.Qmax=[]; S.Gen.Qmin=[]; S.Gen.S=[]; S.Gen.a=[]; S.Gen.b=[]; S.Gen.c=[]; n_aux=0; whiles(1)~='' n_aux=n_aux+1; k2=1; [k1k2]=find_k1k2(s,k2); S.Gen.Bus(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Gen.Pmax(n_aux,1)=str2num(s(k1:k2)); 47
PAGE 56
[k1k2]=find_k1k2(s,k2); S.Gen.Pmin(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Gen.Qmax(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Gen.Qmin(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Gen.S(n_aux,1)=str2num(s(k1:k2))/S.base.S; [k1k2]=find_k1k2(s,k2); S.Gen.a(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Gen.b(n_aux,1)=str2num(s(k1:k2)); [k1k2]=find_k1k2(s,k2); S.Gen.c(n_aux,1)=str2num(s(k1:k2)); s=fgetl(fcf); end S.Gen.n=n_aux; 48
PAGE 57
fclose(fcf); A.2.1SubFunctionndk1k2 Theintentofthishelperfunctionistondthestartandstopcoordinatesofnumericalvaluesinatextleinordertoavoidproblemswithexcessivewhitespace. function[k1k2]=find_k1k2(s,k2); k1=k2+1; whiles(k1)=='',k1=k1+1;end k2=k1; whiles(k2)~='' ifk2==length(s) break end k2=k2+1; end A.3SubFunctionBuildYbus ThissubfunctionbuildsthesparseYbusMatrix. functionS=build_ybus(S) %Buildybus %fillsthediagonalelementsfromshuntcompensators S.Ybus=sparse(diag(S.Bus.G+j*S.Bus.B)); %fillsdiagonalelementsfromtransmissionlinesandtransformers S.Ybus=S.Ybus+sparse(S.Branch.From,S.Branch.From,... 1./(S.Branch.r+j*S.Branch.x)+j*S.Branch.b/2,...49
PAGE 58
S.Bus.n,S.Bus.n); S.Ybus=S.Ybus+sparse(S.Branch.To,S.Branch.To,... 1./(S.Branch.r+j*S.Branch.x)+j*S.Branch.b/2,... S.Bus.n,S.Bus.n); %fillsoffdiagonalelements S.Ybus=S.Ybus+sparse(S.Branch.From,S.Branch.To,... 1./(S.Branch.r+j*S.Branch.x),... S.Bus.n,S.Bus.n); S.Ybus=S.Ybus+sparse(S.Branch.To,S.Branch.From,... 1./(S.Branch.r+j*S.Branch.x),... S.Bus.n,S.Bus.n); A.4SubFunctionObjFun ThissubfunctiondenestheobjectivefunctiontobeminimizedormaximizedbyMATLAB'sbuiltinFminconfunction. functionf=ObjFun(x,S) na=S.Bus.n1; nv=S.Bus.n1; np=1; nq=length(S.Gen.Bus); 50
PAGE 59
f=x(na+nv+1)+sum(S.Bus.MWGen(S.Bus.PV))sum(S.Bus.MWLoad(S.Bus.PQ)); A.5SubFunctionNonLinCon ThissubfunctiondenesthenonlinearequalityandinequalityconstraintsforusewithMATLAB'sbuiltinFminconfunction. function[C,Ceq]=NonLinCon(x,S) %computenumberofeachvariable na=S.Bus.n1; nv=S.Bus.n;%tookout1 np=1; nq=length(S.Gen.Bus); %UpdatesVoltages %Vbus=S.Bus.V; Vbus=x(na+1:na+nv);%added Vbus(S.Bus.Swing)=Vbus(S.Bus.Swing)*exp(j*0); i_av=find(S.Bus.type~=3); %Vbus(i_av)=x(na+1:na+nv); Vbus(i_av)=abs(Vbus(i_av)).*exp(j*x(1:na)); %Lines'Flow Sfrom=Vbus(S.Branch.From).*conj((Vbus(S.Branch.From)Vbus(S.Branch.To))./(S.Branch.r+j*S.Branch.x)+... Vbus(S.Branch.From).*(j*S.Branch.b/2)); Sto=Vbus(S.Branch.To).*conj((Vbus(S.Branch.To)Vbus(S.Branch.From))./(S.Branch.r+j*S.Branch.x)... Vbus(S.Branch.To).*(j*S.Branch.b/2));51
PAGE 60
%Sinj=Vbus.*conj(S.Ybus*Vbus); %updateinjectedpower Pinj=S.Bus.MWLoad; Qinj=S.Bus.MVARLoad; Pinj(S.Bus.PV)=S.Bus.MWGen(S.Bus.PV); Pinj(S.Bus.Swing)=x(na+nv+1); Qinj(S.Gen.Bus)=x(na+nv+np+1:na+nv+np+nq); Sinj=Pinj+j*Qinj; Rate=S.Branch.MVA1; %InequalityConstraints C=[abs(Sfrom)Rate;abs(Sto)Rate;sqrt(real(Sinj(S.Gen.Bus)).^2+imag(Sinj(S.Gen.Bus)).^2)S.Gen.S]; %Powerflowequations Sv=Vbus.*conj(S.Ybus*Vbus); %Equalityconstraint Ceq=[Pinjreal(Sv);Qinjimag(Sv)];52
PAGE 61
AppendixB.FortranInterfaceModuletolinkPSCADandMATLAB ThefollowingcodewasdevelopedforthePSCADcomponentthatlinksthesimulationtotheMATLABreactivepowercontrolalgorithm.Thecodehandlesthememoryassignmentsofthevariousinputsandoutputs. #STORAGEREAL:8 #LOCALINTEGERI_CNT IF($TRIG==1)THEN DOI_CNT=1,2,1 STORF(NSTORF+I_CNT1)=$pg(I_CNT) ENDDO DOI_CNT=1,2,1 STORF(NSTORF+2+I_CNT1)=$pl(I_CNT) ENDDO DOI_CNT=1,2,1 STORF(NSTORF+4+I_CNT1)=$ql(I_CNT) ENDDO CALLMLAB_INT("$Path","$Name","R(2)R(2)R(2)","R(2)") DOI_CNT=1,2,1 $Qcmd(I_CNT)=STORF(NSTORF+6+I_CNT1) ENDDO ENDIF53
PAGE 62
NSTORF=NSTORF+854
