2-STAGE APPROACH FOR THE OPTIMAL POWER FLOW OF
PHOTOVOLTAIC-PENETRATED DISTRIBUTION SYSTEMS
B.S., University of Colorado Denver, 2016
A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Electrical Engineering Program
This thesis for the Master of Science degree by Amin Bunaiyan has been approved for the Electrical Engineering Program by
Fernando Mancilla-David, Chair Jae Do Park Dan Connors
Date: May 13, 2017
Bunaiyan, Amin (M.S., Electrical Engineering)
2-Stage Approach for the Optimal Power Flow of Photovoltaic-Penetrated Distribution Systems
Thesis directed by Associate Professor Fernando Mancilla-David
The thesis presents a 2-stage model for achieving feasible optimal solution along with numerical validation for the exactness of the available second order cone relaxation in the distribution optimal power flow problem. To this end, the IEEE 4-bus, 13-bus, 34-bus and 123-bus distribution power systems are adapted by balancing the various feeder segments and by populating selected nodes with photovoltaic inverters. The different test feeders are studied under typical daily power profiles. Results generated by Matlab Software for Disciplined Convex Programming (CVX) and the nonlinear solver (fmincon) suggest that second order cone relaxation for the branch flow model is exact only for small radial systems but not necessary for large radial systems.
The form and content of this abstract are approved. I recommend its publication.
Approved: Fernando Mancilla-David
This thesis is dedicated to my family who have supported me since the beginning of studies. They provide me with a great source of motivation and inspiration. Finally, this thesis is dedicated to all those who believe that knowledge is power.
Greatest thanks to my parents who have supported me completely through good times and bad. Also, I would like to express my deepest gratitude to my advisor, Dr. Fernando Mancilla-David, for his guidance, valuable advice and support during this research.
TABLE OF CONTENTS
1. INTRODUCTION..................................................... 1
2. DISTRIBUTION OPTIMAL POWER FLOW BALANCED FEEDER . 5
2.1 Original Branch Flow Model.................................... 5
2.2 Convex Relaxation of the BFM ................................. 7
3. BRANCH ANALYSIS.................................................. 10
3.1 Unbalanced Series Impedance Matrices......................... 10
3.1.1 Overhead lines........................................ 10
3.1.2 Underground Lines..................................... 11
3.2 Unbalanced Shunt Admittance Matrices......................... 13
3.2.1 Overhead Lines........................................ 13
3.2.2 Underground Lines..................................... 14
3.3 Approximating The Line Matrix................................ 15
3.3.1 Approximating the Series Impedance Matrix............. 15
3.3.2 Approximating the Shunt Admittance Matrix ............ 17
4. CASE STUDIES....................................................... 21
4.1 Introduction................................................. 21
4.2 Error Computation............................................ 22
4.3 Preparing Case Studies....................................... 23
4.3.1 Approximating a Balanced Feeder....................... 23
4.3.2 Synthesizing Load Profiles............................ 24
4.3.3 Synthesizing PV Generation Profiles................... 26
4.4 Case Studies................................................. 27
4.4.1 Case Study IEEE 4-bus feeder.......................... 27
4.4.2 Case Study IEEE 13-bus feeder......................... 30
4.4.3 Case Study IEEE 34-bus feeder......................... 33
4.4.4 Case Study IEEE 123-bus feeder .................... 36
5. CONCLUSIONS.............................................. 39
A. MATHEMATICAL BACKGROUND......................................... 42
A.l Convexity................................................. 42
A. 1.1 Convex Set........................................ 42
A. 1.2 Convex Function .................................. 42
A. 1.3 Convexity conditions.............................. 43
A. 1.4 Relaxation........................................ 43
A. 2 Lagrangian Optimization................................... 44
A. 3 Duality................................................... 44
A.4 Karush-Kuhn-Tucker (KKT) conditions....................... 46
A.5 Semidefinite Programming SDP.............................. 47
A.6 Second Order Cone Programming SOCP......................... 47
LIST OF TABLES
3.1 Nomenclature and Definitions for Chapter 3................................. 20
4.1 IEEE 4-bus feeder: Line Configuration Data................................ 25
4.2 IEEE 123-bus feeder: Line Configuration Data.............................. 25
4.3 IEEE 4-bus feeder: tightness results...................................... 28
4.4 IEEE 4-bus feeder: duality gap and relative error of decision variables ... 29
4.5 IEEE 13-bus feeder: tightness results........................................ 31
4.6 IEEE 13-bus feeder: duality gap and relative error of decision variables ... 32
4.7 IEEE 34-bus feeder: tightness results........................................ 34
4.8 IEEE 34-bus feeder: duality gap and relative error of decision variables ... 35
4.9 IEEE 123-bus feeder: tightness results.................................... 37
4.10 IEEE 123-bus feeder: duality gap and relative error of decision variables . . 38
LIST OF FIGURES
2.1 Balanced model of distribution feeder line connecting two buses.................... 5
3.1 unbalanced line feeder............................................................ 10
4.1 2-Stage approach.................................................................. 21
4.2 The IEEE 4-bus test feeder network................................................ 27
4.3 IEEE 4-bus feeder power profiles: total demanded power and total available power
from PV generation............................................................. 27
4.4 The IEEE 13-bus test feeder network............................................... 30
4.5 IEEE 13-bus feeder power profiles: (a) total demanded power and total available
power from PV generation, (b) individual PV panel profiles..................... 30
4.6 The IEEE 34-bus test feeder network............................................... 33
4.7 IEEE 34-bus feeder power profiles: (a) total demanded power and total available
power from PV generation, (b) individual PV panel profiles..................... 33
4.8 The IEEE 123-bus test feeder network.............................................. 36
4.9 IEEE 123-bus feeder power profiles: (a) total demanded power and total available
power from PV generation, (b) individual PV panel profiles..................... 36
A.l Examples of convex and nonconvex sets ............................................ 42
A.2 Definition of convex function .................................................... 42
A.3 1st order condition of convex function............................................ 43
A.4 (a) relaxed set, (b) relaxed function............................................. 44
A.5 Epigraph geometric interpretation for an optimization problem..................... 46
The Optimal Power Flow (OPF) problem integrates the power flow (PF) problem and the economic dispatch (ED) problem that, known as a nonlinear problem , The ED finds the lowest generation cost that satisfies demand, while the PF counts for losses and the topology of the network. It can also identify the power transfer distribution function (PTDF). The OPF problem, which was introduced by J. Carpentier in French in 1962  is NP-hard because of its nonlinearity. The problem is NP-hard because of its nonlinearity. However, it is a highly desirable problem that can be solved in many ways for different applications. More specifically, there are many techniques for solving this problem by linearizing or relaxing the nonlinearity.
A decoupled (DC) power flow approximation was proposed in 1974 . The DC OPF is solved through linear programming (LP), which counts for the topology of the system but not for losses of the system. From an economic point of view, the DC OPF works fairly well for transmission networks but not for distribution networks. The reason is that DC approximation works better when the X/R ratio is being sufficiently high , Distribution systems were designed to be unidirectional  and their X/R ratio is not high enough to use DC approximation.
The problem has become more interesting on a distribution level as renewable resources have become more available to individual users. The AC OPF was solved using a nonlinear solver, but there is no guarantee its solution will converge to the global optimum; instead, it converges to a local minimum that is highly dependent on the initial points. Therefore, studies have taken the path of linearization and convex relaxation.
The semidefinite programming (SDP) relaxation technique is one way to obtain a solution to a convex problem. In , the AC OPF was solved using SDP in 2008. SDP relaxation for the OPF problem was proposed using primal-dual interior point
algorithm that is a convex relaxation. It states that SDP works out as linear programming but in its matrices form. This method takes advantage of the strong duality of the LP. It was proven that hrst-order KKT optimum conditions are the same for primal and dual OPF problems. The disadvantage is that, due to limitations in computer capacity, problems in extremely large systems cant be solved .
Further, ,  and  describe the use of SDP, whose globally optimal solution may be found with available numerical algorithms in polynomial time. In , SDP is used to obtain solutions to convex problems solvable in polynomial-time complexity. As a result, the method demonstrates the ability to reach the global optimal solution of the original nonconvex optimal power flow problem. In , SDP is used to obtain a computationally feasible convex reformulation. This is shown using real world photovoltaic (PV) generation and load profile data for an illustrative low-voltage residential distribution system. In , the OPF problem is solved by leveraging SDP through the use of a centralized computational device with reduced computational burden.
Second order cone programming (SOCP) is another way to solve convex problems. In , it was proven that SOCP is a special case of SDP. This means that SDP is a generalization of SOCP. However, SOCP is more tractable optimization class than SDP. Therefore, SOCP should not be solved as SDP.
Angle reparameterization along with second order cone (SOC) relaxation are additional techniques used to solve a convex OPF problem. In , the relaxation in the bus injection model (BIM) and the branch flow model (BFM) were presented along with a proven relationship between the two models. In , three types of sufficient conditions of exactness were investigated. The first condition involves power injections that require some of the active and reactive powers injected at the two ends of each line to not be bounded. In other words, upper and lower constraints of both ends of each line cannot be active at the two ends unless the objective function is strictly convex. The
second condition involves voltage magnitudes that require all the upper bounds of voltages to be inactive. The third condition refers to voltage angles, which requires small enough differences across each line.
One main advantage of solving the OPF problem using SOC relaxation finding a globally optimal solution to the problem , If the solution found using the relaxation is exact, then the optimal solution to the original nonconvex OPF problem can be computed. Through the study in , it has been found that the topology of the network is an important factor when solving the OPF problem. The angle relaxation is exact for radial systems but not for mesh systems. The conditions are insufficient for mesh systems; however, they could be sufficient if the phase shifters are controlled in a strategic way that makes the system behave like a radial system ,
Another investigation in  shows that an additional sufficient condition is necessary for ensuring the exactness of the SOC relaxation for distribution systems. The condition is related to the thermal limits and capacity of the branches. This condition states that there can only be active reversed power, only reactive reversed power, or no reverse power flows through a line. This condition is important due to the presence of renewable energy devices capable of injecting power in the lines.
The branch flow model, which involves relaxation and convexification, is one such design. The design is made to carry out optimization of mesh topologies in addition to radial networks. This phase includes reviewing the convex relationship, controlling load flow, checking on optimal power flow, and managing the entire power system. The power flow model focuses on different variables associated with the distributed power, which includes the voltages and injections of current together with power injections as to regulate power flow in the network .
The relaxed branch flow approach presented in . Solving the OPF problem is carried out using two stages of relaxation. The initial step involves the elimination of the current and voltage angles while the purpose of the second stage is to approximate
the outcoming problem via a conic program. The relaxation steps have been proven equal through the radial systems .
The main contributions of this thesis are 1.) to numerically validate the exactness of the SOC relaxation in the BFM of the optimal power flow problem at the distribution level (DOPF) for balanced feeders, 2.) to investigate the strength of the solution found by the relaxed model, and 3.) to seek a true solution for the DOPF for different test feeders. This thesis is organized into the following chapters: Chapter II describes the original BFM and the two-stage process of SOC relaxation for balanced networks. Chapter III discusses the process of approximating balanced feeders of the available IEEE test feeders that are unbalanced by nature. Chapter IV shows the numerical results of the tightness, duality gap, and decision variable comparison for different case studies. Chapter V provides an analysis and conclusion to the case studies.
DISTRIBUTION OPTIMAL POWER FLOW BALANCED FEEDER
The nature of the OPF problem is that it is nonlinear and nonconvex. Therefore, there are different approaches to solving this optimization problem. Also, there are different setups of the model; the model used in this research is the branch flow model.
In Fig. 2.1, Ijk represents the complex current flow through the line, Vj represents the complex voltages at bus j, the series line impedance is Zjk = Rjk + jXjk, and Sjk is the complex power flowing through the line.
The branch flow model is designed based on the computation of the power flow through the line. The optimal power flow optimization is used to minimize power losses.
Figure 2.1: Balanced model of distribution feeder line connecting two buses.
2.1 Original Branch Flow Model
The following model along with the relaxation are explained in :
min f(x) = Y
Vk Vj Zjkljk, (2.2)
Sjk = Vjl*k, (2.3)
sj = 4 Vi2l!f)
- Zyliul2 + j|Vj|2t|i), kiS\j (2.4)
|Vmin| < |Vj| < | V max 11 (2.5)
| Ijk | | Imax 11 (2.6)
Pf,min < Pf < Pjjmax (2.7)
dj.min *4j dj.max' (2.8)
In the above optimization problem, the objective function is described in (2.1), where it minimizes the produced active power that leads to minimizing the active power losses. Equation (2.2) represents Ohms law. Further, (2.3) represents the computation of apparent power flow through a line where denotes a conjugate.
In (2.4), Sj is the net power injected at bus j, which is the bus of interest, where the power flows from i to j and then from j to k. The net injected power is computed as follows:
>i U>j P-i ./b|j q-) (2.9)
where p^pj1 are the generated and demanded active powers, and q^qj1 are the generated and demanded reactive powers.
Furthermore, equations (2.5), (2.6), (2.7), (2.8) are the upper and lower bounds of the decision variables.
This model is a nonconvex and nonlinear model. Its decision variables are complex
quantities, which are Ijk, Vj, Sjk, and Sj. Equation (2.3), and the left inequality of (2.5) are the main sources of nonlinearity.
2.2 Convex Relaxation of the BFM
The first step of relaxation is performed by multiplying each side of (2.2) and (2.3) by their corresponding conjugates that attain (2.10) and (2.11), respectively.
Also, it is essential to modify (2.4) by separating the real and imaginary parts, shown in (2.9), that attain (2.12) and (2.13).
After applying the proposed step, the model of the optimization problem becomes:
min f(x) = J^pf (2.1)
|Vkp = |V,|2 (SjtZjt* + ZJt V) + |ZJt|2|IJt|2, (2.10)
Pjk + Q?k = |Vj|2|ijkl2, (2.11)
p-p?= Erv-IhPii-RiiN2)- (2.12)
4 d = E (Q* IvjI2w)
-E(ci'i-xlI'il2 + lvil2T)' i:i>j (2.13)
|Vmin|2 < |Vj|2 < |Vmax|2, (2.14)
IT 12 ^ IT 12 I pjk | Pmaxl i (2.15)
P? < p? < p? J- j^min l j J- j^max^ (2.7)
q? < q? < q? (2.8)
In the new formulation, the unknown variables become |Vj|2, |Ijk|2, Pjk, Qjk, Pj and q_j. The relaxation only considers the magnitudes of the complex quantities, called the angle-relaxation model. This is more of reparametrization than an actual mathematical relaxation. This relaxation is always exact for radial systems . This
model is implemented using fmincon solver, which is a nonlinear solver built in Matlab software.
However, the formulation remains nonconvex and nonlinear due to (2.11). Therefore, the second relaxation applied is the SOC, which relaxes the equality of (2.11) to an inequality (2.16). The model will become as follows:
min f(x) = J^pf (2.1)
|Vkp = |V,|2 (SjtZjt* + ZJt V) + |ZJt|2|IJt|2, (2.10)
Pjk + Qjl < |Vil2|iik|2. (2.16)
if if = E (Q* IvjI2w)
- Ew>i xiiN2 + lvjl2-|C i:i>j (2.13)
|Vmin|2 < |Vj|2 < |Vmax|2, (2.14)
IT 12 ^ IT 12 |-*-jk| Pmaxl i (2.15)
P? < p? < p? J- j^min l j i j^max? (2.7)
q? < q? < q? (2.8)
This model is called the SOC-relaxation model, which is a convex model, but its exactness is related to the explained sufficient condition. The model is implemented using CVX, which is a package for specifying and solving convex programs .
An important step is to recover the angle values after solving either the angle-relaxation model or the SOC-relaxation model. The distributed algorithm of
angle recovery was presented in , and it is as follows:
ZIjk = ZVj zsjk (2.17)
ZVj = Z( Vj 2 ZZjk*Sjk) (2.18)
Line segments are configured differently between overhead and underground lines. For example, the IEEE 123-bus test has eleven varying connection types for overhead lines and one configuration for the underground line. Impedance matrices are provided in  for each configuration and are unbalanced.
All test feeders in  are unbalanced feeders, which is numerically recognized by looking at the line segment data. Because the provided matrices are not ready to be approximated, a backward step is made by recalculating the matrices using Carsons equation and the Kron reduction.
3.1 Unbalanced Series Impedance Matrices
3.1.1 Overhead lines
Bus j Bus k
(sending end) (receiving end)
Figure 3.1: unbalanced line feeder .
Using the configuration data provided in  and applying it to Carsons equations (3.1), (3.2), the primitive impedance matrix is built  .
zu = P + re +j0.12134 x Zij = re + jO. 12134 x ln(^b)
The primitive impedance matrix is computed as follows:
^aa ^ab -^ac ^an
^ba ^bb %hc ^bn Zca Zcb Zcc Zcn
Zna Znb Znc Znn
From this matrix, the Kron reduction (3.3) is used to compute the phase impedance matrix [ZabC], which is an equivalent matrix that includes the effect of the neutral phase.
Zij = Zij ^ (3.3)
The phase impedance matrix is as follows:
Zaa Zab Zac
Zba Zbb Zbc
Zca Zcb Zcc
This matrix is unsymmetrical because its off-diagonal elements are not equivalent. The computed matrices do match all the provided matrices  for the overhead line configurations of all different test feeders.
3.1.2 Underground Lines
In some test feeders in , there are underground line configurations that are Concentric Neutral Cable, which has 3 phases/ 3 wires. Therefore, computing [Zabc] is not the same as overhead line configuration . Since each cable has its own conductors, [Zprim] is calculated using the modified Carsons equations (3.4) (3.11) as
[Zd [Z, [Z,] [Z
Zaa Zab Zac Zba ^bb Zbc Zca Zcb Zcc
zu = r + re + jO. 12134 x ln()
zy = re + jO. 12134 x ln(^)
Zaan Zabn ZaCn Zban Zbbn Zbcn Zcan Zcbn ZCcn
Ziin = r + re + jO. 12134 x ln()
zn = re +j0.12134 x ln(-^-)
Zana Zanb ZanC Zbna Zbnb Zbnc Zcna Zcnb Zcnc
zini = r + re +j0.12134 x ln()
re + jO. 12134 x ln(-!p-)
^t>n&n ^bnbn ^bnCn
Zinin = r + re +jO. 12134 x ln(^) Zinjn = re + jO. 12134 x ln(-j^-)
Similar approach is used for the overhead lines, from primitive matrix [Zprim], the Kron reduction (3.12) is used to compute the phase impedance matrix [Zabc] .
[Zbc] = [Z] [Zln][Zm,]-1[Zj]
3.2 Unbalanced Shunt Admittance Matrices
3.2.1 Overhead Lines
Using the same methodology of the computation of the series impedance matrices, a step back is needed to compute the shunt capacitive matrices before applying any approximation. According to , the primitive potential coefficient matrix [Ppi-im] is needed in the process of computing the line admittance matrix Yabc is computed.
[Pprim] has the size of the number of conductors by the number of conductors. For the overhead configurations, the lines are three phase-four wires. The primitive potential coefficient matrix is as follows:
Paa Pab Pac Pan Pba Pbb Pbc Pbn Pea Pcb Pcc Pen Pna Pnb Pnc Pnn
The elements of the primitive potential coefficient matrix are computed as follows:
* = 2^ X ^ to<|)
Using Kron reduction 3.15, the phase potential coefficient matrix [Pabc] is as follows:
Paa Pab Pac
Pba Pbb Pbc
Pea Pcb Pec
Now, it is possible to compute the phase capacitive matrix that is [Cabc] = [Pabc Then, the shunt admittance matrix [Yabc] is computed as follows:
[Yabc] = j 27t [Cabc]
Yaa Yab Yac
Yba Ybb Ybc
Yca Ycb Ycc
The results of this matrix do match the provided data of different feeders in . The matrix has unequal off-diagonal elements, which make the matrices asymmetrical.
3.2.2 Underground Lines
Again, some of the feeders have configurations that are underground. The computation of the Yabc matrix is different than the overhead lines. Since all phases of the concentric neutral strand are grounded, there is only the capacitance between each phase of ground Ypg to be computed.
j 27rf x
This implies that all phases have the same charge , which results in the following matrix:
y 1 pg 0 0
0 y 1 pg 0
0 0 Y 1 pg
Finally, this is the last sanity check of the computed matrices with the provided matrices in .
3.3 Approximating The Line Matrix
3.3.1 Approximating the Series Impedance Matrix
The goal is to approximate the phase impedance matrix with a balanced symmetric matrix that has all diagonal elements (self-impedance elements Zs) equal to each other and all off diagonal elements (mutual impedance elements Zs) equal to each other . This approach has been used to solve power flow and fault problems.
Therefore, recalculating an approximated balanced version of the primitive impedance matrix Zprim.b requires considering a geometric mean distance (GMD) for overhead and underground line configurations, which will make the Carsons equations (3.2), (3.5), (3.7), (3.9), and (3.11) change to:
zy = re + jO. 12134 x ln(^)
where GMD = vDabDbcDCa for phase-phase mutual impedances and GMD = v^DanDbnDcn for phase-neutral mutual impedances.
After applying the Kron reduction to Zpi.im.b, the balanced impedance matrix Zabc.b will be:
^aa Zab Zac 7 7 7
^ba Zfob Tjhc = 7 7 7 /-J'm
^ca ^ch Zcc 7 7 7
This matrix is symmetrical; all off-diagonal elements are equal to each other, and all diagonal elements are also equal to each other. This matrix can substitute Zabc in order to compute an approximated solution.
Out of Zabc.b) the sequence impedance matrix Z0i2 is calculated as the following:
Z0 0 0 0 Zi 0
0 0 z2
Zo Zs + 2Zm
Zi z2 zs ZE
From this matrix, three different models of approximated matrices could be built :
1. The Sequence Phase Impedance Model ZSeq
(2Zi + Z0) (Z0 + Zi) (Z0 + Zi)
(Zo + Zi) (2Zi + Z0) (Z0 + Zi)
(Zo + Zi) (Z0 + Zi) (2Zi + Z0)
This matrix happens to have exact entries of Zabcb, which is the matrix calculated using the GMD. Assuming a balanced power load and generation, this matrix could be used to solve one phase representation.
2. The Modified Sequence Phase Impedance Model ZMd
(2Zi + Zo) 0 0
0 (2Zi + Zo) 0 0 0 (2Zi + Zo)
The off-diagonal elements of this matrix, considering the mutual phase coupling between phases, can be eliminated. Also, it is one phase representation that has been used in the case study of the approximated model in the coming chapters.
3. The Positive Sequence Phase Impedance Model Zpos
Zj 0 0
[Zpos] = 0 Ta\ 0
This matrix only considers the positive sequence phase. Numerically, the values are slightly off. Therefore, it might not be the best choice for the case study.
3.3.2 Approximating the Shunt Admittance Matrix
According to , it is possible to approximate Pabc where Pabc.b has the following structure:
P P P
-1- s y m y m
P P P
-1- m -1- s -1- m
P P P
L m L m s
using the resultant matrix to compute the following:
[Yabc.b] = j 27T [Cabc] = j2yr[Pabc] 1
Y Y Y
-1- s ^ m ^ m
Y Y Y
-*-111 -1- S -*-111
Y Y Y
This is done using the same technique for Zabc by computing GMD and considering line transpose, which leads to the following matrix, called (sequence potential coefficient
matrix P012 ):
[P012] = [A]_ 1 [Pabc.b] [A]
1 1 1
1 a-2 a 1 a Q'2
Equation 3.20 could be computed following the same process of series impedances:
Po 0 0
[P012] = 0 Pi 0
Po = Ps + 2Pm Pi = ?2 = Ps Pm
using the resultant P012 to compute C012 and Y0i2 as the following:
Co 0 0
[C012] = [P012] 1 = 0 Ci 0
Y0 0 0
0 Yi 0
0 0 Y2
For the underground line, it is not necessary to approximate nor to have the sequence matrix since the matrix computed originally is diagonal with the equivalent diagonal element. Approximating the shunt admittance matrix for the overhead line will follow the same process used when approximating the series line impedance In order to have the same features at the approximated matrices.
1. The Sequence Phase Admittance Model Yseq
(2Y1+Y0) (Yo+Yi) (Yo+Yi)
(Yo+Yi) (2Y1+Y0) (Yo+Yi)
(Yo+Yi) (Yo+Yi) (2Y1+Y0)
Similar to the series impedance work, this matrix happens to have the same entries as Yabcb. Assuming a balanced power load and generation, this matrix could be used to solve one phase representation.
2. The Modified Sequence Phase Admittance Model Ymd
(2Y1+Y0) 0 0
0 (2Y1+Y0) 0
0 0 (2Zi + Yq)
Following the same assumption made for Zmd, this matrix considers the mutual phase coupling between phases where the off-diagonal elements can be eliminated. Ymd is used in the case study of the approximated model in future chapters.
3. The Positive Sequence Phase Admittance Model Ypos
Yi 0 0
[Ypos] = 0 Y[ 0
0 0 Yi
Similarly, this matrix considers only the positive sequence phase. Numerically, the values are slightly off. Therefore, it might not be the best choice for this case study.
Table 3.1: Nomenclature and Definitions for Chapter 3.
Zii (Q/mi) diagonal impedances
z (Q/mi) off-diagonal impedances
ri (Q/mi) conductor resistance
le = 1.588 10~3f(n/mi) resistance of Carsons equivalent of earth return conductor
De = 2160 (f )-* (ft) function of earth resistivity and frequency
D (ft) the distance between conductor i and j
GMD (ft) the geometric mean distance between phases
GMR (ft) conductor geometric mean radius
f = 60(Hz) frequency
P = 100(Q/m) average earth resistivity
[Zij] (Q/mi) phase phase conductor matrix
[Zin] (Q/mi) phase neutral of conductor matrix
[Znj] (Q/mi) neutral phase conductor matrix
[Znn] (Q/mi) neutral neutral of conductor matrix
k - number of concentric strands
dc (inch) diameter of the phase conductor
ds (inch) diameter of the phase strand
d-od (inch) outside diameter
R = ^(ft) distance between the center of the phase
conductor to its own strand
r (Q/mi) phase conductor resistance
rs (Q/mi) neutral strand resistance
^cn = ifT/mi) equivalent resistance of the concentric
GMRg (ft) strand geometric mean radius
GMRcn = dGMRskRk-!(ft) equivalent geometric mean radius of the
su (ft) distance between the phase and its image
s (ft) distance between phase i and image of phase j
RDi (inch) radius of the conductor
P (f) mutual potential coefficient
o = 1.424 10-2(si) permittivity of the air
Rb - dVd(ft) distance between the center of the
phase conductor to its own strand
RDC (inch) radius of the phase conductor
RDS (inch) radius of the neutral strand
er = 2.3074(^f) permittivity of the medium
CASE STUDIES 4.1 Introduction
Conditions under which SOC-based relaxations are exact for the OPF problem are currently being investigated in  and . In contrast, the first part of this case study investigates the numerical results of the relaxed model. The goal is to test the exactness of this SOC relaxation by checking if the inequality of (2.16) is tight at different power levels. The second part of the case study investigates the distance between the results and the true results by using the 2-stage model showing in Fig 4.1. This approach leads to finding a feasible point using the nonlinear solver because the relaxed solution may be infeasible.
Figure 4.1: 2-Stage approach.
Recalling Model 1: SOC-relaxation model and Model 2: angle-relaxation model, which are explained in detail in chapter 2. Model 1 is an exact model but nonlinear. This model is implemented using the fmincon solver in Matlab. Model 2 is a relaxed and convex model. This model is implemented using CVX in Matlab. The decision variables solution of Model 2 will be fed to the nonlinear solver as initial guess values for solving Model 2.
The CVX algorithm was configured with the following settings: semidefinite programming SDPT3 solver and precision controlled to a set tolerance called best . For the fmincon solver, the nonlinear tolerance was assigned as xlO-13 to insure tightness.
4.2 Error Computation
For each case study, three different types of errors were computed at different times.
Recall the following equations from chapter 2.
pjk + Qjk= |Vj|2|ijk|2, (2.11)
Pjk + Qjk<|Vj|2|ijk|2, (2.16)
The above equations are the only difference between the two models. The tightness compares the right-hand side with the left-hand side of the inequality in (2.16), where the left-hand side is LHS = P?k + Q?k, and the right-hand side is RHS = |Vj|2|Ijk|2-
The tightness relative error is computed as follows:
| LHS RHS |
| RHS |
This relative error is computed per line for each system at each time. And to attain exactness, the error should be small enough to imply that LHS = RHS.
2. Duality Gap:
The duality gap compares the objective values computed in the two models.
There will be one error calculation per time.
The solution computed by the relaxed model is the lower bound (L), where the solution computed by the exact model is the upper bound (U). The relative error is computed as follows:
DG = x 100. (4.2)
3. Decision Variables:
The decision variables x of both models are |Vj|2, |Ijk|2, Pjk, Qjk, Pj and qj. The computed error will compare both the solution found by the relaxed model (x0) and the solution found by the exact model (x). The relative error is computed as follows:
DX = X~X\ x fOO (4.3)
4.3 Preparing Case Studies
In order to conduct the case study, there are some assumptions and approximations that have been made. This is necessary because the test feeders provided by  are unbalanced distribution systems and the algorithm that is used only considers balanced systems. Therefore, preparing the benchmark is necessary for the relaxed BFM.
4.3.1 Approximating a Balanced Feeder
Data in  shows that some of the positions are consistently closed. An assumption is made that, for all closed switches, the impedance is negligible, creating a short circuit .
One important approximation is that of the balance branch impedances and admittance matrices. Approximation done for the branch impedance matrices was obtained by following  and . Following the same methodology, the approximation of the branch admittance matrix follows . These approximations were explained previously in chapter 3.
The first step of the proposed approximation is considering the geometric mean distance to Carsons equation as represented in  to obtain an approximated primitive matrix. After applying the Kron reduction, the resultant matrix is a balanced matrix as follows:
' 7 7
7 7 m
r 7 7
In this matrix, Zs is the self impedance element and Zm is the mutual impedance element. Matching the assumption made in , assuming a balanced power load and generation, this matrix could be used to solve one phase representation by not considering the mutual phase coupling of phases. This results in the elimination of the off-diagonal elements, which generates the matrix called the modified sequence phase impedance matrix [ZMd] :
Zs 0 0 0 Zs 0
0 0 zs
Consequently, the matrix called the modified sequence phase admittance matrix [YMd]  is:
Ys 0 0 0 Ys 0
0 0 Ys
For example, Tables 4.1, 4.2 show the approximated results for each configuration at the different test feeders where Zs = r + jx and Ys = jb.
4.3.2 Synthesizing Load Profiles
The test feeder data in  provides spot load data that contains different types of load connections. However, all loads are represented as kW or kVAR at each phase.
These case studies follow some approximations that are made in . One assumption is that the three phase loads are uniformly distributed. Due to this type of distribution,
Table 4.1: IEEE 4-bus feeder: Line Configuration Data
configuration r x b
number (Q / mi) (Q / mi) (S / mi)
1 0.4837 0.9932 5.763
Table 4.2: IEEE 123-bus feeder: Line Configuration Data
configuration number r (Q / mi) X (Q / mi) b (S / mi;
1 0.4837 0.9932 5.763
2 0.4837 0.9932 5.763
3 0.4837 0.9932 5.763
4 0.4837 0.9932 5.763
5 0.4837 0.9932 5.763
6 0.4837 0.9932 5.763
7 0.4737 1.025 4.8518
8 0.4737 1.025 4.8518
9 1.3292 1.3474 4.5193
10 1.3292 1.3474 4.5193
11 1.3292 1.3474 4.5193
12 1.5241 0.7406 67.2242
the three phases are identical, which allows them to be decoupled and presented as a single phase.
The IEEE 4-bus feeder is assumed to have the same load at all times, which is not realistic but serves the purpose of a test for such a small system. However, for all other test feeders, the load value provided is used as the peak value so that the lines wont be overloaded. Coinciding with the case study, the reactive power was computed with a
power factor of 0.9 as proposed in . The synthesized temporal distribution load profile, demonstrated in figures 4.3, 4.5, 4.7, and 4.9, for each individual case is reflected by a typical daily power demand similar to . The nodal load was distributed across the nodes proportionally to the load level described in .
4.3.3 Synthesizing PV Generation Profiles
In the original data, devices including voltage regulators and shunt capacitors have been used in order to maintain voltage levels. In these case studies, the voltage regulators and shunt capacitors will not be considered. However, adding the PV panels will substitute for the functionality of all the systems voltage regulators and shunt capacitors.
The location and active power rating of the PV components were chosen at random, and their individual profiles will be shown for every single case. The computation of the reactive power is as follows :
Srated 1 1 V P a. a /\ i. max (4.4)
Qmax . /c2 p2 V 'rated max (4.5)
4.4 Case Studies
4.4.1 Case Study IEEE 4bus feeder
12 3 4
Figure 4.2: The IEEE 4-bus test feeder network.
0 1000 Oh
0 5 10 15 20
Figure 4.3: IEEE 4-bus feeder power profiles: total demanded power and total available power
from PV generation
The IEEE 4-bus feeder is the smallest feeder provided in . The system is created to
be a test feeder with the possibility of choosing different specifications. In this case
study, a balanced close connection load has been chosen. The load is designed to be 1.8
MW at a power factor of 0.9 lagging. The voltage level is 12.47 kV at the substation
bus that is located at bus 1 as shown in Fig.4.2. A step-down transformer has been
chosen, where the secondary voltage is 4.16 kV. A 4-wire line configuration has been
chosen, and the approximated line impedances are shown in Table4.1.
As pointed out, a balanced close connection load has been chosen. The load is designed
to be 1.8 MW at a power factor of 0.9 lagging. This load has been kept as a constant at
all times for simplicity. Specifically, the feeder has been mapped close to Denver
International Airport (DIA), a commercial location that is open most of the time.
Table 4.3: IEEE 4-bus feeder: tightness results
Demand Total PV CVX fmincon
Time P (kW) P (kW) max(Â£rei) max(erei)
05 a. m. 1800 0.00 4.24 x io-10 3.05 x I013
06 a. m. 1800 0.00 4.24 x io-10 3.05 x IO-13
07 a.ni. 1800 0.00 4.24 x io-10 3.05 x I013
08 a. m. 1800 12.70 1.01 x io-10 1.07 x IO"12
09 a. m. 1800 113.60 1.86 x 10-9 6.91 x I013
10 a. m. 1800 370.33 2.78 x 10-9 5.92 x I014
If a.m. 1800 967.31 1.37 x io-11 6.05 x I014
12 p.ni. 1800 1259.59 1.45 x io-11 1.45 x IO11
01 p.ni. 1800 1315.37 4.55 x io-12 4.55 x IO"12
02 p.ni. 1800 1217.85 1.16 x io-11 1.16 x IO11
03 p.ni. 1800 1054.13 2.55 x io-11 2.02 x I014
04 p.ni. 1800 715.89 9.35 x 10-9 2.51 x IO"12
05 p.ni. 1800 323.39 2.32 x 10-9 2.02 x I014
06 p.ni. 1800 187.93 1.36 x 10-9 2.73 x IO-12
07 p.ni. 1800 43.91 1.34 x 10-9 1.46 x I013
According to the website of the U.S. Department of Energy, the PV farm next to DIA has a rating of 1.9 MW. The latitude and longitude of these PV generations were provided to the renewables ninja website (https://www.renewables.mnja/) to show data on available power on April 10, 2014.
The blue line in Fig.4.3 displays the load profile of active power during the day, while the red line represents the available active power from the PV generator. Table 4.3 shows the optimization results of fifteen case studies at different times throughout the day including the active power demanded, available active power from the PV panels,
Table 4.4: IEEE 4-bus feeder: duality gap and relative error of decision variables
Time Demand P (kW) Total PV P (kW) CVX L (kW) fmincon U (kW) DG (%) rnax(DX) (%)
05 a.m. 1800.00 0.00 1864.95 1864.95 0.000 0.024
06 a.m. 1800.00 0.00 1864.95 1864.95 0.000 0.024
07 a.m. 1800.00 0.00 1864.95 1864.95 0.000 0.024
08 a.m. 1800.00 12.70 1864.89 1864.89 0.000 0.005
09 a.m. 1800.00 113.60 1864.49 1864.49 0.000 0.005
10 a.m. 1800.00 370.33 1863.56 1863.56 0.000 0.000
If a.m. 1800.00 967.31 1861.95 1861.95 0.000 0.032
12 p.m. 1800.00 1259.59 1861.65 1861.65 0.000 0.000
01 p.m. 1800.00 1315.37 1861.70 1861.70 0.000 0.000
02 p.m. 1800.00 1217.85 1861.65 1861.65 0.000 0.000
03 p.m. 1800.00 1054.13 1861.80 1861.80 0.000 0.025
04 p.m. 1800.00 715.89 1862.53 1862.53 0.000 0.005
05 p.m. 1800.00 323.39 1863.72 1863.72 0.000 0.001
06 p.m. 1800.00 187.93 1864.20 1864.20 0.000 0.000
07 p.m. 1800.00 43.91 1864.77 1864.77 0.000 0.001
and the maximum relative error of the tightness for both models. Table 4.4 shows the the optimal values of the relaxed model and the exact model, the duality gap of those optimal values, and the relative error of the worst decision variables.
4.4.2 Case Study IEEE 13bus feeder
Figure 4.4: The IEEE 13-bus test, feeder network.
0 5 10 15 20
Figure 4.5: IEEE 13-bus feeder power profiles: (a) total demanded power and total available power from PV generation, (b) individual PV panel profiles.
The IEEE 13-bus test feeder presented in Fig. 4.4 is designed to operate at a nominal voltage of 4.16kV , The load profile used in this case study is obtained from ,
Table 4.5: IEEE 13-bus feeder: tightness results
Demand Total PV CVX fmincon
Time P (kW) P (kW) max (Del) max (Del)
05 a. m. 356.00 0.00 2.06 X 10- -08 2.78 X 10- -n
06 a. m. 302.60 90.78 9.40 X 10- -06 4.11 X 10- -n
07 a. m. 356.00 181.56 7.99 X 10- -06 1.69 X 10- -10
08 a. m. 445.00 242.08 1.51 X io- -05 9.79 X io- -13
09 a. m. 534.00 302.60 1.71 X io- -06 5.49 X io- -12
10 a. m. 551.80 347.99 5.10 X IO" -08 3.37 X io- -13
11 a. m. 516.20 363.12 5.83 X IO" -07 4.53 X io- -12
12 p.m. 498.40 363.12 5.70 X IO" -08 2.22 X io- -13
01 p.m. 551.80 363.12 7.21 X IO" -08 2.19 X io- -13
02 p.m. 569.60 332.86 6.98 X IO" -08 2.28 X io- -11
03 p.m. 587.40 272.34 4.33 X IO" -07 6.19 X io- -12
04 p.m. 587.40 211.82 5.96 X IO" -08 1.20 X io- -11
05 p.m. 747.60 151.30 2.74 X IO" -07 1.59 X io- -11
06 p.m. 1068.00 30.26 1.11 X IO" -08 6.41 X 10 -13
07 p.m. 640.80 0.00 3.15 X IO" -07 3.79 X io- -12
which shows the power flow result as 1.1 MW. This value is used as the peak value so that the lines will not be overloaded. The blue line in Fig. 4.5(a) displays the load profile of active power during the day. The synthesized temporal distribution load profile, demonstrated in Fig. 4.5(a), is reflected by a typical daily power demand similar to . The nodal load was distributed across the nodes proportionally to the load level described in . In Fig. 4.5(a), the red line represents the aggregated available active power from the PV generators. The individual profiles are shown in Fig. 4.5(b). Table 4.5 shows the optimization results of fifteen case studies at different times throughout
Table 4.6: IEEE 13-bus feeder: duality gap and relative error of decision variables
Time Demand P (kW) Total PV P (kW) cvx L (kW) fmincon U (kW) DG (%) max(DX) (%)
05 a.m. 356.00 0.00 357.59 357.59 0.000 0.000
06 a.m. 302.60 90.78 303.27 303.27 0.000 0.056
07 a.m. 356.00 181.56 356.59 356.59 0.000 0.706
08 a.m. 445.00 242.08 445.86 445.86 0.000 7.686
09 a.m. 534.00 302.60 535.20 535.20 0.000 0.007
10 a.m. 551.80 347.99 552.98 552.98 0.000 0.000
11 a.m. 516.20 363.12 517.12 517.12 0.000 0.208
12 p.m. 498.40 363.12 499.20 499.20 0.000 0.002
01 p.m. 551.80 363.12 552.97 552.97 0.000 0.002
02 p.m. 569.60 332.86 570.96 570.96 0.000 0.279
03 p.m. 587.40 272.34 589.21 589.21 0.000 0.016
04 p.m. 587.40 211.82 589.64 589.64 0.000 12.369
05 p.m. 747.60 151.30 752.61 752.61 0.000 0.004
06 p.m. 1068.00 30.26 1082.34 1082.34 0.000 0.000
07 p.m. 640.80 0.00 646.01 646.01 0.000 0.000
the day including the active power demanded, available active power from the PV panels, and the maximum relative error of the tightness for both models. Table 4.6 shows the optimal values of the relaxed model and the exact model, the duality gap of those optimal values, and the relative error of the worst decision variables.
4.4.3 Case Study IEEE 34bus feeder
Figure 4.6: The IEEE 34-bus test, feeder network.
Figure 4.7: IEEE 34-bus feeder power profiles: (a) total demanded power and total available power from PV generation, (b) individual PV panel profiles.
The IEEE 34-bus test feeder presented in Fig. 4.6 is designed to operate at a nominal voltage of 24.9kV . The load profile used in this case study was obtained from , which shows the power flow result as 240 kW. This value is used as the peak value so that the lines will not be overloaded. The blue line in Fig. 4.7(a) displays the load profile of active power during the day. The synthesized temporal distribution load
Table 4.7: IEEE 34-bus feeder: tightness results
Demand Total PV CVX fmincon
Time P (kW) P (kW) max (Del) max (Del)
05 a. m. 80.00 0.00 1.15 X 10- -03 6.32 X 10- -10
06 a. m. 68.00 27.60 5.16 X 10- -04 1.55 X 10- -10
07 a. m. 80.00 55.20 1.70 X 10- -04 2.40 X 10- -10
08 a. m. 100.00 73.60 1.27 X io- -04 9.60 X io- -10
09 a. m. 120.00 92.00 1.27 X io- -04 5.24 X io- -10
10 a. m. 124.00 105.80 1.33 X IO" -04 8.30 X io- -10
11 a. m. 116.00 110.40 1.50 X IO" -04 3.93 X io- -09
12 p.m. 112.00 110.40 1.43 X IO" -04 2.01 X io- -10
01 p.m. 124.00 110.40 1.39 X IO" -04 3.27 X io- -10
02 p.m. 128.00 101.20 5.01 X IO" -05 3.05 X io- -10
03 p.m. 132.00 82.80 9.19 X IO" -05 1.27 X io- -10
04 p.m. 132.00 64.40 1.02 X IO" -04 1.56 X io- -10
05 p.m. 168.00 46.00 1.43 X IO" -04 3.72 X io- -10
06 p.m. 240.00 9.20 1.95 X IO" -04 2.53 X 10 -10
07 p.m. 144.00 0.00 7.01 X IO" -03 1.88 X io- -10
profile, demonstrated in Fig. 4.7(a), is reflected by a typical daily power demand similar to . The nodal load was distributed across the nodes proportionally to the load level described in . In Fig. 4.7(a), the red line represents the aggregated available active power from the PV generators. The individual profiles are shown in Fig. 4.7(b). Table 4.7 shows the optimization results of fifteen case studies at different times throughout the day including the active power demand, available active power from the PV panels, and the maximum relative error of the tightness for both models. Table 4.8 shows the optimal values of the relaxed model and the exact model, the duality gap of those
Table 4.8: IEEE 34-bus feeder: duality gap and relative error of decision variables
Time Demand P (kW) Total PV P (kW) cvx L (kW) fmincon U (kW) DG (%) rnax(DX) (%)
05 a.m. 80.00 0.00 80.88 80.88 0.000 0.000
06 a.m. 68.00 27.60 68.61 68.61 0.000 0.471
07 a.m. 80.00 55.20 80.72 80.72 0.000 0.000
08 a.m. 100.00 73.60 101.03 101.03 0.000 0.000
09 a.m. 120.00 92.00 121.43 121.43 0.000 0.000
10 a.m. 124.00 105.80 125.52 125.52 0.000 0.238
If a.m. 116.00 110.40 117.35 117.35 0.000 7.594
12 p.m. 112.00 110.40 113.27 113.27 0.000 1.770
01 p.m. 124.00 110.40 125.52 125.52 0.000 0.000
02 p.m. 128.00 101.20 129.62 129.62 0.000 0.000
03 p.m. 132.00 82.80 133.76 133.76 0.000 0.001
04 p.m. 132.00 64.40 133.84 133.84 0.000 0.544
05 p.m. 168.00 46.00 171.21 171.21 0.000 10.051
06 p.m. 240.00 9.20 247.39 247.39 0.000 0.036
07 p.m. 144.00 0.00 146.63 146.63 0.000 0.000
optimal values, and the relative error of the worst decision variables.
Active Power (kW)
4.4.4 Case Study IEEE 123bus feeder
Figure 4.8: The IEEE 123-bus test, feeder network.
Figure 4.9: IEEE 123-bus feeder power profiles: (a) total demanded power and total available power from PV generation, (b) individual PV panel profiles.
The IEEE 123-bus test feeder presented in Fig. 4.8 is designed to operate at a nominal voltage of 4.16kV , The approximated line configurations are shown in Table 4.2.
Table 4.9: IEEE 123-bus feeder: tightness results
Demand Total PV CVX fmincon
Time P (kW) P (kW) max (Del) max (Del)
05 a. m. 400 0.00 1.34 X 10- -05 1.37 X 10- -09
06 a. m. 340 130.80 2.84 X 10- -03 1.30 X 10- -10
07 a. m. 400 261.60 5.31 X 10- -03 7.57 X 10- -10
08 a. m. 500 348.80 8.54 X io- -03 2.01 X io- -10
09 a. m. 600 436.00 4.38 X io- -04 2.65 X io- -10
10 a. m. 620 501.40 4.28 X IO" -03 7.29 X io- -10
11 a. m. 580 523.20 8.30 X IO" -03 2.21 X io- -09
12 p.m. 560 523.20 6.67 X IO" -03 4.02 X io- -10
01 p.m. 620 523.20 2.63 X IO" -03 1.54 X io- -10
02 p.m. 640 479.60 3.83 X IO" -05 5.49 X io- -10
03 p.m. 660 392.40 1.91 X IO" -03 7.23 X io- -10
04 p.m. 660 305.20 2.31 X IO" -04 3.88 X io- -10
05 p.m. 840 218.00 4.06 X IO" -05 7.48 X io- -10
06 p.m. 1200 43.60 1.21 X IO" -05 8.03 X 10 -11
07 p.m. 720 0.00 6.79 X IO" -07 2.49 X io- -10
The load profile used in this case study is obtained from , which shows the power flow result as 1.2 MW. An additional load was added to bus 610 because an old version of  shows that this bus is connected to an induction machine. This value is used as the peak value so that the lines will not be overloaded. The blue line in Fig. 4.9(a) displays the load profile of active power during the day. The synthesized temporal distribution load profile, demonstrated in Fig. 4.9(a), is reflected by a typical daily power demand similar to . The nodal load was distributed across the nodes proportionally to the load level described in . In Fig. 4.9(a), the red line represents
Table 4.10: IEEE 123-bus feeder: duality gap and relative error of decision variables
Time Demand P (kW) Total PV P (kW) cvx L (kW) fmincon U (kW) DG (%) max(DX) (%)
05 a.m. 400 0.00 401.52 401.52 0.000 0.000
06 a.m. 340 130.80 340.42 341.58 0.340 198.732
07 a.m. 400 261.60 400.22 401.50 0.317 196.209
08 a.m. 500 348.80 500.29 501.48 0.236 189.284
09 a.m. 600 436.00 600.38 601.46 0.179 192.776
10 a.m. 620 501.40 620.32 620.95 0.100 197.101
11 a.m. 580 523.20 580.25 580.46 0.037 199.160
12 p.ni. 560 523.20 560.22 561.11 0.159 197.287
01 p.ni. 620 523.20 620.32 622.99 0.429 199.266
02 p.ni. 640 479.60 640.41 641.42 0.159 197.316
03 p.ni. 660 392.40 660.76 661.54 0.118 199.715
04 p.ni. 660 305.20 661.23 661.88 0.098 196.002
05 p.ni. 840 218.00 843.67 843.67 0.000 7.981
06 p.ni. 1200 43.60 1212.92 1212.92 0.000 0.321
07 p.ni. 720 0.00 724.95 724.95 0.000 0.000
the aggregated available active power from the PV generators. The individual profiles are shown in Fig. 4.9(b). Table 4.9 shows the optimization results of fifteen case studies at different times throughout the day including the active power demand, available active power from the PV panels, and the maximum relative error of the tightness for both models. Table 4.10 shows the optimal values of the relax model and the exact model, the duality gap of those optimal values, and the relative error of the worst decision variables.
Based on the results of the approximated balanced IEEE 4-bus, 13-bus, 34-bus, and 123-bus test feeders presented in chapter 4, the SOC relaxation provides a strong solution when compared to the tightness results. This indicates that the numerical results of the tightness for all test feeders at different load and generation levels are close to 0% for the radial modified systems. However, results show that the smaller the test feeder, the better is the tightness.
The duality gap results computed for the different feeders indicate that the objective values were exact for all cases except for the IEEE 123-bus test feeder that has the highest gap to be 0.43%. This indicates that the relaxed model finds an objective function that is close enough to a feasible local optimum which is most likely to be the global optimum.
However, the percent error of decision variables demonstrates that feasible solutions are not necessarily the same as those computed by the relaxed model. With the exception of the IEEE 123-bus test feeder, the error values for all tests were below 15%, which relatively small. However, the error computed for the IEEE 123-bus test feeder were large enough that both the total dispatch and the power flow were different.
With the help of the 2-stage model, results indicate that, for small systems, the relaxation is strong enough and the second stage is able to determine the global optimum. However, the 2-stage model was able to determine a feasible local minimum that may or may not be the global optimum for the IEEE 123-bus test feeder even if the objective values are similar in the two models.
It worth pointing out that SOC relaxation is a powerful tool that facilitated the development of the 2-stage model by providing the initial values to the nonlinear solver. Otherwise, the solver would start from a flat point with no guarantee of finding global optimality. Therefore, further investigation is needed to approach global solution.
 A. J. Wood, B. F. Wollenberg, and G. B. Shebl, Power Generation, Operation, and Control, 3rd Edition. Hoboken, New Jersey, USA: Wiley, 2014.
 J. Carpentier, Contribution to the economic dispatch problem, Bulletin de la Societe Francoise des Electriciens, vol. 3, no. 8, pp. 431-447, 1962.
 B. Stott and O. Alsac, Fast decoupled load flow, IEEE Transactions on Power Apparatus and Systems, vol. PAS-93, no. 3, pp. 859-869, May 1974.
 S. Huang, Q. Wu, J. Wang, and H. Zhao, A sufficient condition on convex relaxation of ac optimal power flow in distribution networks, IEEE Transactions on Power Systems, vol. 32, no. 2, pp. 1359-1368, Mar. 2017.
 E. DalFAnese, S. V. Dhople, and G. B. Giannakis, Optimal dispatch of photovoltaic inverters in residential distribution systems, IEEE Transactions on Sustainable Energy, vol. 5, no. 2, pp. 487-497, Apr. 2014.
 X. Bai, H. Wei, K. Fujisawa, and Y. Wang, Semidehnite programming for optimal power flow problems, Inti J. of Electrical Power & Energy Systems, 208.
 E. DallAnese, H. Zhu, and G. B. Giannakis, Distributed optimal power flow for smart microgrids, IEEE Transactions on Smart Grid, vol. 4, no. 3, pp. 1464-1475, Sep. 2013.
 E. DalFAnese, S. V. Dhople, B. B. Johnson, and G. B. Giannakis, Decentralized optimal dispatch of photovoltaic inverters in residential distribution systems, IEEE Transactions on Energy Conversion, vol. 29, no. 4, pp. 957-967, Dec. 2014.
 J. A. Taylor, Convex Optimization of Power Systems. Padstow, Cornwall, UK: Cambridge University Press, 2015.
 S. H. Low, Convex relaxation of optimal power flow; Part I: Formulations and equivalence, IEEE Transactions on Control of Network Systems, vol. 1, no. 1, pp. 15-27, Mar. 2014.
 S. H. Low, Convex relaxation of optimal power flow; Part II: Exactness, IEEE Transactions on Control of Network Systems, vol. 1, no. 2, pp. 177-189, Jun. 2014.
 M. Farivar and S. H. Low, Branch flow model: Relaxations and convexihcation (parts I, II), IEEE Transactions on Power Systems, vol. 28, no. 3, pp. 2554-2572, Aug. 2013.
 M. Grant and S. Boyd, CVX: Matlab software for disciplined convex programming, version 2.1, http://cvxr.com/cvx, 2014.
 W. H. Kersting, Radial distribution test feeders, in IEEE Power Engineering Society Winter Meeting, vol. 2, Jan. 2001, pp. 908-912.
 W. H. Kersting and W. H. Phillips, Distribution feeder line models, IEEE Transactions on Industry Applications, vol. 31, no. 4, pp. 715-720, Jul. 1995.
 W. H. Kersting, Distribution System Modeling and Analysis. Boca Raton, FL, USA: CRC Press, 2002.
 T. Gonen, Modern Power System Analysis. Boca Raton, FL, USA: CRC Press, 2013.
 L. Gan, N. Li, U. Topcu, and S. H. Low, Exact convex relaxation of optimal power flow in radial networks, IEEE Transactions on Automatic Control, vol. 60, no. 1, pp. 72-87, Jan. 2015.
 A. OConnell and A. Keane, Multi-period three-phase unbalanced optimal power flow, in IEEE PES Innovative Smart Grid Technologies, Europe, Oct. 2014, pp. 1-6.
 V. Kekatos, G. Wang, A. J. Conejo, and G. B. Giannakis, Stochastic reactive power management in microgrids with renewables, IEEE Transactions on Power Systems, vol. 30, no. 6, pp. 3386-3395, Nov. 2015.
 L. Gan and S. H. Low, Convex relaxations and linear approximation for optimal power flow in multiphase radial networks, in Power Systems Computation Conference, Aug. 2014, pp. 1-9.
 S. Boyd and L. Vandenberghe, Convex Optimization. UK: Cambridge University Press, 2009.
This appendix highlights the mathematical background needed for the study performed in this thesis. For a comprehensive understanding of the materials in this appendix, reading  is highly recommended.
A. 1.1 Convex Set
S C E is convex if
x,y E S ==>- Xx + jjy G S for all A,fi > 0 with A + fi = 1
Figure A.l: Examples of convex and nonconvex sets
A.1.2 Convex Function
/ : E ^ E is convex if the domain is a convex set and f(ex + (i-e)y)
Figure A.2: Definition of convex function
A.1.3 Convexity conditions
1st order condition:
If f(x) is differentiable, and its gradient is \/f(x), then f(x) is convex if: f(y) > f(x) + V/(^)T(y x) for ah x, y G dom/.
2nd order condition:
If f(x) is twice differentiable, and its Hessian is \/2f(x), then f(x) is convex if: \/2f(x) A 0 Positive Semidehnite (PSD)
Figure A.3: 1st order condition of convex function.
To sum up, an equality set has to be an affine (linear) function to be convex, which is is known as hyperplane. For a convex inequality set, this can be an affine called halfspace or any function that satisfies one of the two conditions. For example, a quadratic function f(x) = xTPx + 2qTx + r is convex only if P >z 0.
U = min fo(x)
x G X
L = min f0(x)
x G Xr
Relaxed if X C Xr then L < U.
Figure A.4: (a) relaxed set, (b) relaxed function
A.2 Lagrangian Optimization
Consider the following optimization function:
min f0(x) (A.5)
fi(x) < 0 (A.6)
hj(x) = 0 (A.7)
The lagrangian function of the problem is as follows:
L(x, A,p.) = f0(x) + ^(/'ACC) (A.8)
the scalar A* and fij are called lagrangian multipliers. The vector A and fi are called dual variables.
In order to explain duality, we will be using Linear Programming (LP): Consider a Primal function in the standard of LP
p = min C x (A.9)
Ax = b (A.10)
x < 0 (A.11)
The lagrangian function of the Primal problem is as follows:
L(x, A, n) = hTp + (ATp A + C)T;r
The dual function in the standard of LP is as follows:
d* = max bT fi
Ax A + C = 0
A > 0
The lagrangian function of the dual problem is as follows:
g(A, n) = hTp + inf((ATp, A + C)T;r)
Lower bound property states that if A > 0
then g(A, //.) < p*
where A and p are the dual variables,
and p* is the optimal value of the primal problem.
When solving the concave dual problem, we are finding the max A (A*).
1. Dual is infeasible
2. Primal is infeasible
3. Both dual and primal are infeasible
4. Both are feasible (typically)
Solution of the dual problem is d*.
Weak duality d* < p*. This hold for convex and nonconvex problems.
Strong duality d* = p*. This does not generally hold but usually holds for convex
Using Slaters Constraints Qualification, which are exempt from the linear inequality constraints but must satisfy nonlinear inequality constraints, strong duality for SDP is solvable.
The epigraph geometric interpretation in fig. A.5.
Figure A.5: Epigraph geometric interpretation for an optimization problem.
A.4 Karush-Kuhn-Tucker (KKT) conditions
Suppose fi(x) and hj(x) are differentiable and x*, A* and p* are a local optimum solution:
1. fi(x*) < 0 hj(x*) = 0 (Primal Feasibility)
2. A* > 0 (Dual Feasibility)
3. :>L ')A|A 41 ^ = 0 (Lagrangian Optimality)
4. A*fi(x*) = 0 (Complementary Slackness)
A.5 Semidefinite Programming SDP
The SDP can be similar to linear programming, but it is in its matrix form. The following example has A0 and Aj as symmetric constant matrices. The trace function sums up the diagonal elements of a matrix.
min tr(Aox) (A.17)
tr(Ajx) = bj, i = 1,..., m (A.18)
The Dual SDP:
max bTy (A.20)
A0 Ai) h 0 (A.21)
A.6 Second Order Cone Programming SOCP
The SOCP is proven to be a special type of SDP. However, the SOCP is more tractable than SDP, which is an advantage. The following example shows a typical way to express the SOCP.
fT min j x (A.22)
AiX + hi 2 < cfx + di i = 1,..., m, (A.23)