DISTRIBUTION SYSTEM RELIABILITY ENHANCEMENT

A Thesis
Presented to
The Academic Faculty

by

Xuebei Yu

In Partial Fulfillment
of the Requirements for the Degree
Master of Science in the
School of Electrical and Computer Engineering

Georgia Institute of Technology

August 2011

COPYRIGHT 2011 BY XUEBEI YU

DISTRIBUTION SYSTEM RELIABILITY ENHANCEMENT

Approved by:

Dr. A.P. Meliopoulos, Advisor

School of Electrical and Computer Engineering
Georgia Institute of Technology

Dr. Deepakraj M Divan

School of Electrical and Computer Engineering
Georgia Institute of Technology

Dr. Carlos S Grijalva

School of Electrical and Computer Engineering
Georgia Institute of Technology

Date Approved: May 16th 2011

 ACKNOWLEDGEMENTS

I would like to thank my advisor Dr. A.P. Meliopoulos for an opportunity to work

on this project and his continued guidance. I would also like to thank the committee, Dr.

Deepakraj M Divan and Dr. Carlos S Grijalva for their time.

During the course of my Masters, I have had the opportunity to exchange ideas

with many of my fellow students, and this has helped me a great deal in my research. For

this, I would like to thank Vangelis Farantatos, Ye Tao, Renke Huang, Peter Suh, Ravi

Shankar, Joosung Kang and Anupama Keeli. I would particularly like to thank Renke

Huang and Vangelis Farantatos for taking the time to help me during the course of my

research, and also for providing many valuable pointers.

I would like to thank my parents for their love and support, this has continued to

be a source of encouragement.

iv
 TABLE OF CONTENTS

Page

ACKNOWLEDGEMENTS iv

LIST OF TABLES vii

LIST OF FIGURES ix

LIST OF SYMBOLS AND ABBREVIATIONS xi

SUMMARY xiii

CHAPTER

1 Introduction 1

1.2 Research Objectives Statement 1

1.2 Research Objectives 3

2 Background Information 5

2.1 System Drives 5

2.2 Power Outage Definition 5

2.3 System Indices 6

3 Initial Response 9

3.1 Introduction 9

3.2 Relays 10

3.2.1 Adaptive Relay Scheme – Directional Overcurrent Relay 11

3.2.2 Adaptive Relay Scheme – Impedance Relay 24

3.2.3 Coordination of Overcurrent Relays and Distance Relays 26

3.3 Fuse 34

3.3.1 Characteristic Parameters 35

3.3.2 Fuse Operation 37

v
 3.4 Coordination of Relays and Fuses 39

3.5 Test System and Calculation 41

3.5.1 System Description 41

3.5.2 Problem Formulation 43

3.5.3 Solution 48

4 Automatic System Response 50

4.1 Introduction 50

4.2 Switch Placement Optimization 51

4.2.1 Problem Formulation 53

4.2.2 Particle Swarm Optimization 55

4.2.3 Genetic Algorithm 72

4.2.4 Results and Conclusions 78

4.3 Automatic Network Reconfiguration 79

4.3.1 Simple Rule-based Reconfiguration 79

4.3.2 Simulation Results of the Optimal Reconfiguration 84

5 Conclusion and Future Work 85

APPENDIX A: Basic Reliability Calculations 87

APPENDIX B: Description of the Test System 89

REFERENCES 99

vi
 LIST OF TABLES

Page

Table 3.1: Constants and exponents for standard characteristics of overcurrent relays 15

Table 3.2: Constraints of the pickup current of each relay in example system 19

Table 3.3: Fault current at each relay to different fault locations 20

Table 3.4: Results for the MINLP and LP formulation 24

Table 3.5: Fault current at each relay to 8 different fault locations 30

Table 3.6: Results for the coordination of overcurrent and distance relay setting 33

Table 3.7: Constraints of the pickup current in IEEE 123 node test system 43

Table 3.8: Fault current in IEEE 123 node test system 44

Table 3.9: Results for the optimal relay setting in IEEE 123 node system 49

Table 4.1: Number of function evaluation for continuous case 61

Table 4.2: Objective function for continuous case 62

Table 4.3: Number of function evaluation for integer case 62

Table 4.4: Objective function for integer case 62

Table 4.5: Comparison of 1-switch and 2-switch simulations 64

Table 4.6: Fval, SAIDI Reduction, & Cost of 1-Switch Verification Case 65

Table 4.7: Number of Function Evaluations of 1-Switch PSO Optimization 67

Table 4.8: Objective Function of 1-Switch PSO Optimization 67

Table 4.9: Number of Function Evaluations of 2-Switch PSO Optimization 68

Table 4.10: Objective Function of 2-Switch PSO Optimization 68

Table 4.11: Optimization Results of PSO with 10 Fault Profiles 70

Table 4.12: Objective Function, SAIDI Reduction, and Cost of 120-Switch Case 71

Table 4.13: Comparison of PSO and GA 78

vii
Table 4.14: Optimal GA Switch Configuration Results 78

Table 4.15: Effects of Types of Switches on Reliability 84

Table B.1: Line segment data 90

Table B.2: Three Phase Switches 91

Table B.3: Overhead Line Configuratioions 91

Table B.4: Underground Line Configuration 91

Table B.5: Shunt Capacitors data 91

Table B.6: Overhead Line Spacing 92

Table B.7: Underground Line Spacing 92

Table B.8: Load data 93

Table B.9: Switch type and location in the test system 95

Table B.10: Load data in the test system 96

viii
 LIST OF FIGURES

Page

Figure 3.1: Example distribution system 18

Figure 3.2: Fault simulation in the example distribution system 20

Figure 3.3: Distance relay protection zone for a radial system 25

Figure 3.4: Coordination of overcurrent relays 28

Figure 3.5: Coordination of overcurrent relays and distance relays 28

Figure 3.6: Test system for coordination of overcurrent relays and distance relays 29

Figure 3.7: Minimum Melting Time and Total Clearing Time for a Specific Fuse 36

Figure 3.8: Fault time current curve 38

Figure 3.9: Time current curve for fuse 40K 39

Figure 3.10: Fuse blow scheme 40

Figure 3.11: Fuse saving scheme 40

Figure 3.12: IEEE 123-Node Test Feeder 42

Figure 3.13: IEEE 123-Node Test Feeder in WinIGS - F 43

Figure 3.14: IEEE 123-Node Test Feeder with fault simulation 44

Figure 4.1: Switch Placement Optimization 52

Figure 4.2: Contour plot of sample objective function 61

Figure 4.3: Validation of PSO Algorithm 63

Figure 4.4: Verification of PSO Algorithm with 1 Switch 69

Figure 4.5: Objective Function of 1-Switch PSO Optimization 71

Figure 4.6: SAIDI Reduction of 1-Switch PSO Optimization 72

Figure 4.7: Selection - Elitest vs. Roulette Wheel 74

Figure 4.8: Two-Point Crossover of GA 75

ix
Figure 4.9: Mutation of GA 76

Figure 4.10: Optimal switch locations for 1 year payback period 79

Figure 4.11: Data Structure for Node 13 80

Figure 4.12: Paths finding and switch commands determination 81

Figure 4.13: Optimal switch locations for automatic reconfiguration 82

Figure 4.14: Network Plot @ t=12s with a Fault at Node 19 84

Figure B.1: IEEE 123-Node Test Feeder 89

Figure B.2: Overhead Line Spacing 92

Figure B.3: Underground Line Spacing 92

Figure B.4: IEEE 123-Node Test Feeder in WinIGS - F 94

x
 LIST OF SYMBOLS AND ABBREVIATIONS

SAIDI System Average Interruption Duration Index

CAIDI Customer Average Interruption Duration Index

SAIFI System Average Interruption Frequency Index

MAIFI Momentary Average Interruption Frequency Index

FDIR Fault detection, isolation and restoration

DOE Department of Energy

TDS Time dial setting

Ip Pickup current

LP Linear Programming

NLP Nonlinear Programming

MINLP Mixed Integer Nonlinear Programming

IEEE Institute of Electrical and Electronics Engineers

WinIGS Windows Integrated Grounding System

LINGO Linear Interactive and General Optimizer

TCC Time Current Curve

PSO Particle Swarm Optimization

GA Genetic Algorithm

FVAL Function Value

i Nearest relay to the fault i=F(j)

Ti , j Operation time of relay i for a fault at location j

∆T Time interval between primary relay and backup relay

I ij Fault current at relay i for a fault at place j,

xi
∆T ' Time interval between overcurrent relay and distance relay

Tzi Time setting for zone i of distance relay

IN Rated current of fuse

VN Rated Voltage of fuse

NT Total number of customers served in an area

xki +1 Position of the particle i at the iteration k+1

vki +1 Velocity of the particle i at the iteration k+1

pi Best position found by particle i

r Independent random numbers between 0 and 1

w Weight

xii
 SUMMARY

Practically all everyday life tasks from economic transactions to entertainment

depend on the availability of electricity. Some customers have come to expect a higher

level of power quality and availability from their electric utility. Federal and state

standards are now mandated for power service quality and utilities may be penalized if

the number of interruptions exceeds the mandated standards. In order to meet the

requirement for safety, reliability and quality of supply in distribution system, adaptive

relaying and optimal network reconfiguration are proposed. By optimizing the system to

be better prepared to handle a fault, the end result will be that in the event of a fault, the

minimum number of customers will be affected. Thus reliability will increase.

The main function of power system protection is to detect and remove the faulted

parts as fast and as selectively as possible. The problem of coordinating protective relays

in electric power systems consists of selecting suitable settings such that their

fundamental protective function is met under the requirements of sensitivity, selectivity,

reliability, and speed. In the proposed adaptive relaying approach, weather data will be

incorporated as follows. By using real-time weather information, the potential area that

might be affected by the severe weather will be determined. An algorithm is proposed for

adaptive optimal relay setting (relays will optimally react to a potential fault). Different

types of relays (and relay functions) and fuses will be considered in this optimization

problem as well as their coordination with others. The proposed optimization method is

based on mixed integer programming that will provide the optimal relay settings

including pickup current, time dial setting, and different relay functions and so on.

xiii
 The main function of optimal network reconfiguration is to maximize the power

supply using existing breakers and switches in the system. The ability to quickly and

flexibly reconfigure the power system of an interconnected network of feeders is a key

component of Smart Grid. New technologies are being injected into the distribution

systems such as advanced metering, distribution automation, distribution generation and

distributed storage. With these new technologies, the optimal network reconfiguration

becomes more complicated. The proposed algorithms will be implemented and

demonstrated on a realistic test system. The end result will be improved reliability. The

improvements will be quantified with reliability indexes such as SAIDI.

xiv
 CHAPTER 1

INTRODUCTION

1.1 Problem Statement

The reliability of electric distribution systems in the US is one which can affect

any business or other entity dependent on electricity for any reason. Customers deal with

problems as small as a light flickers to extended outages which may take hours or even

days to fix. Some businesses that lose power for even a few seconds may completely lose

a product in development costing them hundreds of thousands of dollars. Outages have

cost the U.S. economy an average of $1.5 billion each week – $80 billion each year, and

customers and utilities alike are looking for a cost effective solution.

In order to improve reliability of the power system, many smart grid technologies

are being and will be implemented in the distribution system. Presently, most distribution

systems are designed based on a main trunk three phase feeder with single-phase laterals.

The main trunk carries most power away from the substation through the center of the

feeder service territory. Single phase laterals are used to connect the main trunk to

customer locations. Actual distribution systems have branching, normally-open loops,

and other complexities, but the overarching philosophy remains the same [1]. Traditional

distribution systems use time-current coordination for protection devices. These devices

assume that faster devices are topologically further from the substation.

Distribution systems will not resemble the distribution systems of today in several

years. New technologies are being injected into the distribution systems such as advanced

1
metering, distribution automation, distribution generation and distributed storage. An

example test system in built in Appendix B which includes some of the new smart grid

technologies. Through the integrated use of these technologies, smart grid will be able to

self heal, provide high reliability and power quality, be resistant to cyber attacks, operate

with multi-directional power flow, increase equipment utilization.

A Smart Grid does not just try to connect substations to customers for the lowest

cost. Instead, a Smart Grid is an enabling system that can be quickly and flexibly

reconfigured. Therefore, future distribution systems will be an integrated Grid of

distribution lines, with the grid being connected to multiple power sources. In a Smart

Grid, topology is flexible and the assumption that faster devices are topologically further

from the substation is problematic. System topology and system protection will have to

be coordinated properly to ensure proper protection coordination for a variety of

configurations. In order to take advantages of these changes, reengineering and

optimizing the existing distribution systems is necessary. For example, adaptive relay

programming is proposed to combat sustained power outages resulting from wind and

lightning related weather events, and it is claimed that the SAIDI reliability metric can be

improved by as much as 30% annually.

The ability to quickly and flexibly reconfigure an interconnected network of

feeders is a key component of Smart Grid, and the most important aspects are in the areas

of protection and switching (often integrated into the same device like breakers and

reclosers). Switches (including breakers, load break switches, tie switches) play

important roles in Fault Detection, Isolation and Restoration (FDIR). At a permanent

2
fault occurrence, the following operations are taken. 1) An automatic feeder

breaker/recloser will trip in real time operation to cut off the fault from the feeder, 2) the

feeders are broken up into sections isolated by switches or breakers 3) power is then

restored to non faulted sections.

1.2 Research Objectives

The present project addresses the above issues. We propose three specific focus

areas. Two related to real time operations and the third addressing the planning issue.

These are defined below.

Focus Area 1: Initial response

By using real-time weather information, the potential area that might be affected

by the severe weather will be determined. An algorithm will be built for adaptive relay

setting for the specific configuration of the system so that the relays can optimally react

to a potential fault. Different types of relays (and relay functions) and fuses will be

considered in this optimization problem.

Focus Area 2: Sectionalization and power restoration

An algorithm will be developed to detect the section of the feeder with fault

occurred, quickly isolate that feeder section by operating the isolating switches or

breakers and restoring power to the non faulted sections. The objective is to maximum

power supply to customers considering the existed smart grid technologies such as

distributed energy resources and distribution automation devices.

3
 Focus Area 3: Switch Placement Optimization

It is obvious that adding more switches in the power system can improve power

system reliability. There exists a trade-off between power reliability improvements and

switch costs. Therefore, several algorithms are built to decide the optimal number of

switches added to the system and the best locations to add those switches.

4
 CHAPTER 2

BACKGROUND AND LITERATURE REVIEW

2.1 Reliability in Power System

2.1.1 System Drives

With today's high electricity costs, customers have come to expect a higher level

of power quality from their electric utility. Higher power quality requirements have

become even more critical with the proliferation of electronic devices such as VCRs,

computers and clock radios, which are intolerant of even the smallest interruption of

power. These issues have become so important to utility customers that federal and state

standards are now mandated for power quality and utilities are now penalized if the

number of interruptions exceeds these mandated standards.

The primary reason for many outages is the aging of the electric grid, as well as

the lack of redundancy, due to the generally centralized nature of the system. In fact, over

80% of transformers, as well as 80% of breakers are 30 years old or older (DOE).

Uncontrollable events such as animal contact or human interaction on the line may also

lead to power outages. Weather can be especially hazardous, sometimes causing

cascading faults leading to millions of dollars of expense. Tree branches may touch lines

during a storm potentially causing an outage, leading to tree trimming programs.

2.1.2 Power Outage Definition

5
 Two types of outages may occur and can be categorized as momentary or

sustained outages. A momentary outage is defined as being any outage lasting from 1

second to 1 minute. Afterwards, it is considered to be a sustained outage. Anything less

than 1 second is a power quality issue. Other definitions for momentary outages consider

it to be any outage lasting less than 5 minutes. It is debatable, but the point is that an

outage, however short, still affects customers and is not desired.

2.1.3 System Indices

The objective of many distribution utilities today is to improve overall customer

satisfaction through quality of service, which is to maximize the reliability of power

distribution system. Quality of service in part is monitored through system performance

indices such as Momentary Average Interruption Frequency Index (MAIFI), Sustained

Average Interruption Duration Index (SAIDI), Customer Average Interruption Duration

Index (CAIDI) and Sustained Average Interruption Frequency Index (SAIFI). Definition

of these indices can be found in Appendix A. Utility companies are measured by these

indices on a yearly basis, and in some states may receive bonuses or penalties depending

on how far they are away from target values or how many standard deviations they are

away from the average of the past few years. Improving these metrics may be a highly

profitable project for utilities in many states to undertake.

2.2 Optimization Algorithms Used in Reliability Enhancement

The focus of this thesis is optimal relay setting and optimal network

reconfiguration in power distribution system. The main function of protection on the

6
power system is to detect and remove the faulted parts as fast and as selectively as

possible. The problem of coordinating protective relays in electric power systems consists

of selecting their suitable settings such that their fundamental protective function is met

under the requirements of sensitivity, selectivity, reliability, and speed. These

requirements must be met for a variety of system conditions and configurations, and can

be translated into conditions such as: (i) a variety of fault conditions must be detected by

the appropriate relays, (ii) the relays located closer to the fault should have priority of

operation, (iii) if a primary relay fails, a backup relay should operate, and (iv) the

operation of the relay should be as fast as possible to prevent equipment damage, and

must occur only in the presence of abnormal operating conditions which jeopardize the

system integrity [2]. The devices mostly used for distribution system protection are:

breakers and reclosers with different relay functions and fuses.

Directional overcurrent relays are widely used in distribution system. Directional

overcurrent relays have two types of settings: time dial setting (TDS) and pickup current

setting ( Ip ).In general, the protective relay coordination problem was formulated in

previous work either as a linear, nonlinear, or a mixed integer nonlinear programming

problem depending on the type of variables in the problem. And the pickup current is the

variable that determines the type of the problem according to the relay characteristic

functions. When Ip is fixed, the coordination problem becomes a Linear Programming

(LP) problem [4-5]. For continuous Ip values, the problem becomes a Nonlinear

Programming (NLP) problem, and when discrete values of Ip are taken into account, the

problem turns to be a Mixed Integer Nonlinear Programming (MINLP) problem.

7
 In the previous work, due to the complexity of nonlinear optimal programming

techniques, the coordination problem was formulated as an LP problem and was solved

using LP techniques. These techniques include Simplex, Dual Simplex and Two Phase

Simplex. The main disadvantage of this formulation is that the pickup currents and relay

characteristic curves are assumed predetermined. Since each relay can allow for various

Ip values and several relay characteristic curves, there could be a better pickup current

setting for each relay that would lead to a better objective function.

In this project, I am going to discuss about the digital relays, since both the time

dial setting and pickup current setting are continuous, the coordination problem is a

Nonlinear Programming (NLP) problem. And for each relay, there will be several options

for curve characteristics selection. Hence, the problem becomes a Mixed Integer

Nonlinear Programming (MINLP) problem. In addition, coordination of the overcurrent

relays and distance relays, and fuse saving scheme will be also considered in the project

so that more constraints are put into the optimization problem to make the formulation

more complicated and complete. Final optimal relay setting solution is obtained using

global solver in GAMS. And in the second half of the thesis, cost-benefit analysis of

adding more switches into the system are analyzed.

8
 CHAPTER 3

INITIAL RESPONSE

3.1 Introduction

Associated with the distribution networks themselves are a number of ancillary

systems to assist in meeting the requirement for safety, reliability and quality of supply.

The most important of these are the protection systems which are installed to clear faults

and limit any damage to distribution equipment. Among the principle causes of faults are

lightning discharges, the deterioration of insulation, vandalism, and tree branches and

animal contacting the electricity circuits. The majority of faults is of a transient nature

and can often be cleared with no loss of supply, or just the shortest of interruptions,

whereas permanent faults can result in longer outages. In order to avoid damage, suitable

and reliable protection should be installed on all circuits and electrical equipment.

Protective devices initiate the isolation of faulted sections of the network in order to

maintain supplies elsewhere on the system. This then leads to an improved electricity

service with better continuity and quality of supply.

A properly coordinated protection system is vital to ensure that an electricity

distribution network can operate within preset requirement for safety for individual items

of equipment, staff and public, and the network overall. Automatic operation is necessary

to isolate faults on the networks as quickly as possible in order to minimize damage. The

economic costs and the benefits of a protection device should also be considered to arrive

at a suitable balance between the requirements of the scheme and the available financial

resources.

9
 The main function of protection on the power system is to detect and remove the

faulted parts as fast and as selectively as possible. The problem of coordinating protective

relays in electric power systems consists of selecting their suitable settings such that their

fundamental protective function is met under the requirements of sensitivity, selectivity,

reliability, and speed. These requirements must be met for a variety of system conditions

and configurations, and can be translated into conditions such as: (i) a variety of fault

conditions must be detected by the appropriate relays, (ii) the relays located closer to the

fault should have priority of operation, (iii) if a primary relay fails, a backup relay should

operate, and (iv) the operation of the relay should be as fast as possible to prevent

equipment damage, and must occur only in the presence of abnormal operating

conditions which jeopardize the system integrity[2]. The devices most used for

distribution system protection are: over current relays, reclosers and fuses. Their

characteristics and coordination between each other will be discussed in the following

sections.

3.2 Relays

The evolution of protection relays started with the electromechanical relays and

the design of the protection relays has changed significantly over the past years with the

advancement in microprocessor and signal processing technology. The first commercial

microprocessor based relay was introduced in 1984. Subsequently in the 1990s

microprocessor technology, along with the improvements in mathematical algorithms,

spurred the development of the so-called numerical relays which are extremely popular

for their multifunctional capabilities, low prices and reliability.

10
 Traditional protective relays are usually set offline and kept constant during

operation. Settings are normally calculated based on the maximum generation mode in

the power system and may not be at optimal performance. In the worst case, when fault

occurs under minimum generation, the performance of protection relays will decrease

because of the poor sensitivity of the protection relays or even failure of operation may

take place. With enlargement of power system, network structure of power supply and

distribution line becomes complex and traditional protection of it cannot satisfy the

requirement of power system. Appearance of adaptive protection offers an effective

method to solve this problem. It can adapt to the changes of different topology and fault

types of power system and then change the setting of the relays to an optimal value.

Adaptive relay programming are proposed to combat sustained power outages resulting

from wind and lightning related weather events, and it is claimed that the SAIDI

reliability metric can be improved by as much as 30% annually.

Directional over current relay and distance relay are the most commonly used

protective relays in power distribution systems. These two types of relays will be

discussed separately in the following sections. And optimization algorithms will be built

to find the optimal settings.

3.2.1 Adaptive Relay Scheme - Directional Over current relay

3.2.1.1 Introduction

The following work is based on the assumption that the potential fault area is

determined by the weather information, and relay settings are going to be optimized so

that the protection devices in the system can react optimally to the fault.

11
 Directional relays operate only when the fault current flows in the specific

tripping direction, they avoid compromising line protection and provide selectivity.

Directional over current relays have two types of settings: time dial setting (TDS) and

pickup current setting ( Ip ). The settings should be chosen to minimize the overall time of

operation of relays while maintaining selectivity and reliability. Thus, the directional over

current relay coordination problem involves optimization, where the solution is the

optimal settings of each relay.

Generally speaking, the over current relay coordination problem was formulated

in previous work as a linear, nonlinear or a mixed integer nonlinear programming

problem. And the pickup current is the variable that determines the type of the problem

according to the relay characteristic functions. When Ip is fixed, the coordination

problem becomes a Linear Programming (LP) problem [4-5]. For continuous Ip values,

the problem becomes a Nonlinear Programming (NLP) problem, and when discrete

values of Ip are taken into account, the problem turns to be a Mixed Integer Nonlinear

Programming (MINLP) problem. In this project, we are going to discuss about the digital

relays, since both the time dial setting and pickup current setting are continuous, the

coordination problem is a Nonlinear Programming (NLP) problem. In addition, for each

relay, there will be several options for curve characteristics selection. Hence, the problem

becomes a Mixed Integer Nonlinear Programming (MINLP) problem.

3.2.1.2 MINLP Formulation for Optimal Over Current Relay Setting

The relay coordination algorithm is to calculate the time dial setting (TDS) and

pickup current setting ( Ip ) with the objective of minimizing the primary relay operation

12
 time with backup relay coordinated at all possible fault locations. The objective function

can be stated as follows:

=z min ∑ Ti , j ∀j ∈ L
j

Where

Ti , j refers to the operation time of relay i for a fault at location j

L is a set of possible fault locations

i indicates the nearest relay to the fault, i=F(j)

The objective function is to be achieved under the following constraints:

A. Coordination criteria

Ti −1, j − Ti , j ≥ ∆T

Where Ti −1, j is the operation time of the upstream relay of relay Ri for a given fault at

place j. ∆T is the coordination time interval between primary relay and backup relay

which includes the breaker operating time, relay overrun time after the fault has been

cleared and a safety margin to compensate for possible deviations. In numerical relays

there is no overrun, and therefore the margin could be chosen as 0.2s [6].

B. Bounds on relay settings

TDSi min ≤ TDSi ≤ TDSi max ∀i

Ipi min ≤ Ipi ≤ Ipi max ∀i

Where TDS is the time dial setting and Ip is the pickup current. In most modern relays

the time dial settings can start from values as low as 0.1s, and goes to as high as 20s. In

distribution systems where it is possible to increase the loading on feeders under

13
 emergency conditions, the minimum pickup current setting is usually 2 times the nominal

circuit current. And in order to make sure that the relay is able to trip under fault

condition, the maximum pick up current is usually chosen to be half of the smallest single

phase to ground fault current. Then the constraints of the relay setting can be rewritten as

0.1 ≤ TDSi ≤ 20 ∀i
1
2 I norm ,i ≤ Ipi ≤ I ph _ g ,i ∀i
2

C. Relay characteristics

All relays are assumed identical and with characteristic function approximated by:

Ak
=Ti , j ( + Bk ) × TDSi
( I ij / Ipi )Ck − 1

Where I ij is the fault current passing through the relay i for a fault at place j, and A, B

and C are constant values to provide selected curve characteristics. It can be seen from

the equation that the non-linearity comes from the pickup current term. By adjusting the

physical design parameters different time current characteristic can be implemented. The

time current characteristics are classified according to the steepness of the curves and the

overall time delay as:

 Short time

 Long time

 Definite Minimum Time

 Moderately Inverse

 Inverse

14
  Very Inverse, and

 Extremely Inverse

According to IEEE Std C37.112-1996, the IEEE Standard Inverse-Time

Characteristic Equations for over current relays, the constant and exponents in Table 3.1

define the shape of the standard Moderately Inverse, Very Inverse, and Extremely Inverse

trip characteristics.

Table 3.1- Constants and exponents for standard characteristics of overcurrent relays

Characteristic A B C

Moderately Inverse 0.0515 0.1140 0.0200

Very Inverse 19.61 0.491 2.0000

Extremely Inverse 28.2 0.1217 2.0000

As a result, an additional parameter should be added to represent the choice of the

time current characteristics curve, whether it is moderately inverse (1), very inverse (2) or

extremely inverse (3).

=K 1, 2or 3 ∀i

In this research area, only these three shapes will be considered. And different k

indicates a different group of A, B and C. In order to translate this constraint into a

mathematical way, the time current characteristics curve selection can be reformulated as

follows

A =0.0515 + K1 (19.61 − 0.0515) + K 2 (28.2 − 0.0515)

B =0.114 + K1 (0.491 − 0.114) + K 2 (0.1217 − 0.114)
C = 0.02 + K1 (2 − 0.02) + K 2 (2 − 0.02)

15
 Where

K1 + K 2 ≤ 1
K1 ,K 2 =0 or 1

So that each time only one of the time current characteristics curves is chosen. For

example, when K1=0 and K2=0, then A= 0.0515, B=0.114 and C=0.02, so moderately

inverse is picked. When K1=1 and K2=0, then A=19.61, B=0.491 and C= 2, so very

inverse is chosen. When K1=0 and K2=1, then A=28.2, B=0.1217 and C=2, which means

extremely inverse is chosen. If there are more time current characteristics curves, they

can all be written in this way to be formulated into the optimization problem.

In summary, the optimization problem, in term of the relay operation times has

the following basic form:

=z min ∑ Ti , j ∀j ∈ L
j

Subject to:

Ti −1, j − Ti , j ≥ ∆T

TDSi min ≤ TDSi ≤ TDSi max ∀i

Ipi min ≤ Ipi ≤ Ipi max ∀i

Ak
=Ti , j ( + Bk ) × TDSi
( I ij / Ipi )Ck − 1

Where
Ti , j refers to the operation time of relay i for a fault at location j

16
L is a set of possible fault locations

i indicates the nearest relay to the fault, i=F(j)

Ti −1, j is the operation time of the upstream relay of relay i for a given fault at place j

∆T is the coordination time interval between primary relay and backup relay

TDS is the time dial setting and Ip is the pickup current

3.2.1.3 Case Study

A simple distribution system is built in WinIGS in Fig. 3.1. There is one

substation with two feeders. Feeder 3.1 is protected by a breaker at the substation, two

reclosers at the indicated location and fuse at the indicated lateral. Feeder 3.2 is

connected to feeder 1 with a normally open switch and has one breaker and one recloser

along the transmission line. Assume that both the breaker and the recloser are equipped

with numerical relays which have over current protection with one of the following

selections:

0.0515
=
Moderately Inverse: T ( + 0.114) × TDSi
( I ij / Ipi )0.02 − 1

19.61
=
Very Inverse: T ( + 0.491) × TDSi
( I ij / Ipi ) 2 − 1

28.2
=
Extremely Inverse: T ( + 0.1217) × TDSi
( I ij / Ipi ) 2 − 1

17
 Figure 3.1 Example distribution system

The locations of the breakers and reclosers are shown in Fig. 3.1. According to

the requirements and rules for relay settings which are discussed in 1.2.1, the

coordination problem can be stated below.

z = min ∑T
j
i, j

s.t.

(Bounds on relay settings)

0.1 ≤ TDSi ≤ 20 for i =

1, 2, 3, 4, 5
1
2 I norm ,i ≤ Ipi ≤ I ph _ g ,i for i =
1, 2, 3, 4, 5
2

(Coordination criteria)

T1, j − T2, j ≥ ∆T

T2, j − T3, j ≥ ∆T

18
 T4, j − T5, j ≥ ∆T

(Relay characteristics)

Ai
=Ti , j ( + Bi ) × TDS
= for i 1, 2,3, 4,5
/ Ipi )Ci − 1
i
( Ii, j

Where

Ai =0.0515 + K i ,1 (19.61 − 0.0515) + K i ,2 (28.2 − 0.0515)

Bi =0.114 + K i ,1 (0.491 − 0.114) + K i ,2 (0.1217 − 0.114)
Ci =0.02 + K i ,1 (2 − 0.02) + K i ,2 (2 − 0.02)

K i ,1 + K i ,2 ≤ 1
K i ,1 , K i ,2 = 0 or 1

The algorithm is implemented to solve the relay optimization problem in Fig. 3.1.

Assume the coordination time between the upstream and downstream relay is 0.2s as is

discussed in section 3.1.2.1. The lower bond and upper bond of the pickup current of

each relay can be read from the simulation result of WinIGS and is listed in Table 3.2.

After simulating faults at all the possible locations 1, 2, 3, 4 and 5, as is shown in Fig3.2,

fault current at each relay are recorded in Table 3.3.

Table 3.2 - Constraints of the pickup current of each relay

Relay Lower bound Upper bound

1 340 375

2 180 355

3 60 350

4 320 395

5 160 380

19
 Table 3.3 – Fault current at each relay to different fault locations

Location 1 Location 2 Location 3 Location 4 Location 5

Fault Current at R1 1875 870 752

Fault Current at R2 837 714

Fault Current at R3 685

Fault Current at R4 1444 790

Fault Current at R5 760

Figure 3.2 Fault simulation in the example distribution system

The detailed equations of the problem are as follows:

20
 min T1,1 + T2,2 + T 3,3 +T4,4 + T5,5
s.t.

(Bounds on relay settings)

0.1 ≤ TDSi ≤ 20 i=1, 2, 3, 4 and 5

340 ≤ Ip1 ≤ 375

180 ≤ Ip2 ≤ 355
60 ≤ Ip3 ≤ 350
320 ≤ Ip4 ≤ 395
160 ≤ Ip5 ≤ 380

(Case 1: When fault occurs at place 1)

Ak1
=T1,1 ( + Bk1 ) × TDS1
−1
Ck1
( I1,1 / Ip1 )

(Case2: When fault occurs at place 2)

Ak1
=T1,2 ( + Bk1 ) × TDS1
−1
Ck1
( I1,2 / Ip1 )
Ak2
=T2,2 ( + Bk2 ) × TDS 2
−1
Ck 2
( I 2,2 / Ip2 )
T1,2 − T2,2 ≥ 0.2

(Case 3: When fault occurs at place 3)

21
 Ak1
=T1,3 ( + Bk1 ) × TDS1
−1
Ck1
( I1,3 / Ip1 )
Ak2
=T2,3 ( + Bk2 ) × TDS 2
−1
Ck2
( I 2,3 / Ip2 )
Ak3
=T3,3 ( + Bk3 ) × TDS3
−1
Ck3
( I 3,3 / Ip2 )
T1,3 − T2,3 ≥ 0.2
T2,3 − T3,3 ≥ 0.2

(Case 4: When fault occurs at place 4)

Ak4
=T4,4 ( + Bk4 ) × TDS 4
−1
Ck4
( I 4,4 / Ip4 )

Case 5: When fault occurs at place 5

Ak4
=T4,5 ( + Bk4 ) × TDS 4
−1
Ck4
( I 4,5 / Ip4 )
Ak5
=T5,5 ( + Bk5 ) × TDS5
−1
Ck5
( I 5,5 / Ip5 )
T4,5 − T5,5 ≥ 0.2

(Relay time current characteristic curve selection constraints)

A=i 0.0515 + Ki ,1 × (19.61 - 0.0515) + Ki ,2 × (28.2 - 0.0515) for =

i 1, 2, 3, 4 and 5
Bi= 0.114 + Ki ,1 × (0.491 - 0.114) + Ki ,2 × (0.1217 - 0.114) for i = 1, 2, 3, 4 and 5
Ci = 0.02 + Ki ,1 × (2.0 - 0.02) + Ki ,2 × (2.0 - 0.02) for i =
1, 2, 3, 4 and 5
Ki ,1 + Ki ,2 ≤ 1 for i =
1, 2, 3, 4 and 5
Ki ,1 ,Ki ,2 =0 or 1

22
 The optimization problem is built in LINGO. Using MINLP function, the result is

shown in Table 3.4. The codes for the programming and the results from LINGO are in

Appendix C. In order to demonstrate the effectiveness of the proposed algorithm, the

widely used LP is also tested in the same distribution system. In the LP formulation, the

objective function and the coordination criteria are the same as the ones in MINLP

formulation. The differences come from the constraints of pickup current and the

selection of relay functions. In LP formulation, the pickup up current and time current

characteristic curve of each relay are determined by experience, and only the TDS is

optimized. The formulation is listed below. And by using LP function in LINGO, the

result is also shown in Table 3.4.

z = min ∑T j
i, j

s.t.

(Bounds on relay settings)

0.1 ≤ TDSi ≤ 20 for i =

1, 2,3, 4,5

(Coordination criteria)

T1, j − T2, j ≥ ∆T

T2, j − T3, j ≥ ∆T

T4, j − T5, j ≥ ∆T

(Relay characteristics)

Ai
=Ti , j ( + Bi ) × TDS
= for i 1, 2,3, 4,5
/ Ipi )Ci − 1
i
( Ii, j

23
Where

Ai , Bi , Ci and Ipi are determined

Table 3.4 – Results for the MINLP and LP formulation

Relay TDS Pickup Current TC curve Objective Value

Relay 1 0.112 375 Extremely Inverse

Relay 2 0.1 201 Extremely Inverse

MINLP Relay 3 0.1 60 Very Inverse 0.4248

Relay 4 0.1 331.1 Extremely Inverse

Relay 5 0.1 160 Extremely Inverse

Relay 1 0.178 360 Extremely Inverse

Relay 2 0.1 300 Extremely Inverse

LP Relay 3 0.1 200 Extremely Inverse 0.9052

Relay 4 0.126 370 Extremely Inverse

Relay 5 0.1 250 Extremely Inverse

Compare the results from MINLP and LP formulation of the optimal relay setting,

it can be seen clearly from the table that MINLP can realize relay operation much faster

than LP by optimizing the relay settings.

3.2.2 Adaptive Relay Scheme - Distance Relay

Distance relay are set on the basis of the positive sequence impedance from the

relay location up to the point on the line to be protected. Line impedances are

proportional to the line lengths and it is this property that is used to determine the

position of the fault, starting from the location of the relay.

24
 Normally, three protection zones in the direction of the fault are used in order to

cover a section of line and to provide back-up protection to remote sections. In the

majority of cases the setting of the reach of the three main protection zones is made in

accordance with the following criteria [6] (as is shown in Figure 3.3):

Zone 1: this is set to cover between 80 and 85 per cent of the length of the protected

line;

Zone 2: this is set to cover all the protected line plus 50 per cent of the shortest next

line;

Zone 3: this is set to cover all the protected line plus 100 per cent of the second

longest line, plus 25 per cent of the shortest next line.

Figure 3.3 Distance relay protection zone for a radial system

In addition to the unit for setting the reach, each zone unit has a timer unit. The

operating time for zone 1, t1, is normally set by the manufacturer to trip instantaneously

since any fault on the protected line detected by the zone 1 unit should be cleared

immediately without the need to wait for any other device to operate. This operation is

fast with just a small delay (two to three cycles) to avoid tripping on transients. The

25
operating time for zone 2 is usually of the order of 0.25 to 0.4 s, and that of zone 3 is in

the range of 0.6 to 1.0 s.

Since the tripping produced by zone 1 is instantaneous, it should not reach as far as

the bus at the end of the first line so it is set to cover only 80-85% of the protected line.

The remaining 20-15% provides a factor of safety in order to mitigate against errors

introduced by the measurement transformers and line impedance calculations. The 20-

15% to the end of the line is protected by zone 2, which operates in t2 s. Zone 3 provides

the back-up and operates with a delay of t3 s. Since the reach and therefore the operating

time of the distance relays are fixed, their co-ordination is much easier than that for

overcurrent relays.

3.2.3 Coordination of over current relays and distance relays

It has been shown that when the line protection schemes are composed of distance

relays and overcurrent relays, the setting of the relays must be computed considering both

relays. Separate relay computation would lead to loss of selectivity. Thus, it is useful to

include distance relays parameters in the process of computing the time dial settings of

over current relays in distribution system which have a mixed scheme with overcurrent

relays and distance relays.

3.2.3.1 Formulation for coordination of overcurrent and distance relay

The problem of finding the time dial setting of directional overcurrent relays has

been stated and solved in section 3.2.1 using Mixed Integer Linear Programming

technique. The technique is based on the mathematical statement of the sensitivity, speed

26
 and selectivity conditions associated with the relay coordination problem. The basic

optimization problem, in term of the relay operation times has the following basic form:

=z min ∑ Ti , j ∀j ∈ L
j

Subject to:

Ti −1, j − Ti , j ≥ ∆T

TDSi min ≤ TDSi ≤ TDSi max ∀i

Ipi min ≤ Ipi ≤ Ipi max ∀i

Ak
=Ti , j ( + Bk ) × TDSi
( I ij / Ipi )Ck − 1

Where
Ti , j refers to the operation time of relay i for a fault at location j

L is a set of possible fault locations

i indicates the nearest relay to the fault, i=F(j)

Ti −1, j is the operation time of the upstream relay of relay i for a given fault at place j

∆T is the coordination time interval between primary relay and backup relay

TDS is the time dial setting and Ip is the pickup current

Figure 3.4 shows the case of two overcurrent relays with similar time current

characteristic, where the relevant faults are determined by points F1 and F2[7].

27
 Fig. 3.4 Coordination of overcurrent relays

The derivation of the coordination constraints for system with overcurrent relays

and distance relays is explained in Figure 3.5. Second zone of distance relay associated

with circuit breaker b must be slower than the overcurrent relays associated with main

circuit breaker m, which can be stated as follows:

Tzb , F 4 − Tm , F 4 ≥ ∆T '

Fig. 3.5 Coordination of overcurrent relays and distance relays

Where ∆T ' is a time coordination interval used in the selectivity constraints

between overcurrent relays and distance relays. For this constraint, the operation time of

28
main overcurrent relay is evaluated at point F4, which corresponds to the ohmic reach of

the second zone of the distance relay associated with the backup circuit breaker.

Overcurrent relay associated with circuit breaker b must be slower that the second zone

of distance relay associated with main circuit breaker m:

Tb , F 3 − Tzm , F 3 ≥ ∆T '

3.2.3.2 Case study for coordination of overcurrent and distance relay

All the relays are assumed to have multi functions including the overcurrent relay

function and distance relay function. The same example distribution system is used for

the optimal relay setting for overcurrent relays and distance relays. Eight faults, as is

shown in Figure 3.6, are simulated in the test system. Short circuit currents for the faults

are calculated and recorded in Table 3.5.

Fig. 3.6 Test system for coordination of overcurrent relay and distance relay

29
 Table 3.5 – Fault current at each relay to 8 different fault locations
Fault Fault Fault Fault Fault
Current(A) Current(A) Current(A) Current(A) Current(A)
/impedance(Ω) /impedance(Ω) /impedance(Ω) /impedance(Ω) /impedanc(Ω)e
at R1 at R2 at R3 at R4 at R5
Location 1 1875
Location 2 870/9.16 837/8.19
Location 3 752/10.56 714/7.69 685/3.86
Location 4 1444
Location 5 790/9.64 760/6.99
Location 6 1798
Location 7 799
Location 8 1380

The detailed equations of the problem are as follows:

min T1,1 +T2,2 +T3,3 +T4,4 +T5,5 +T2,6 +T3,7 +T5,8

Subject to

(Bounds on relay settings)

0.1 ≤ TDSi ≤ 20 i=1, 2, 3, 4 and 5

0.2 ≤ Tzi ≤ 0.6 i=1, 2, 3, 4 and 5
340 ≤ Ip1 ≤ 375
180 ≤ Ip2 ≤ 355
60 ≤ Ip3 ≤ 350
320 ≤ Ip4 ≤ 395
160 ≤ Ip5 ≤ 380

(Case 1: When fault occurs at place 1)

Ak1
=T1,1 ( + Bk1 ) × TDS1
−1
Ck1
( I1,1 / Ip1 )

30
 (Case2: When fault occurs at place 2)

Ak1
=T1,2 ( + Bk1 ) × TDS1
−1
Ck1
( I1,2 / Ip1 )
Ak2
=T2,2 ( + Bk2 ) × TDS 2
−1
Ck 2
( I 2,2 / Ip2 )
T1,2 − T2,2 ≥ 0.2
T1,2 − Tz2,2 ≥ 0.3

(Case 3: When fault occurs at place 3)

Ak1
=T1,3 ( + Bk1 ) × TDS1
−1
Ck1
( I1,3 / Ip1 )
Ak2
=T2,3 ( + Bk2 ) × TDS 2
−1
Ck2
( I 2,3 / Ip2 )
Ak3
=T3,3 ( + Bk3 ) × TDS3
−1
Ck3
( I 3,3 / Ip2 )
T1,3 − T2,3 ≥ 0.2
T2,3 − T3,3 ≥ 0.2
T2,3 − Tz3,3 ≥ 0.3

(Case 4: When fault occurs at place 4)

Ak4
=T4,4 ( + Bk4 ) × TDS4
−1
Ck4
( I 4,4 / Ip4 )

31
 (Case 5: When fault occurs at place 5)

Ak4
=T4,5 ( + Bk4 ) × TDS 4
−1
Ck4
( I 4,5 / Ip4 )
Ak5
=T5,5 ( + Bk5 ) × TDS5
−1
Ck5
( I 5,5 / Ip5 )
T4,5 − T5,5 ≥ 0.2
T4,5 − Tz5,5 ≥ 0.3

(Case 6: When fault occurs at place 6)

Ak2
=T2,6 ( + Bk2 ) × TDS 2
−1
Ck2
( I 2,6 / Ip2 )
Tz1,6 − T2,6 ≥ 0.3

(Case 7: When fault occurs at place 7)

Ak3
=T3,7 ( + Bk3 ) × TDS3
−1
Ck3
( I 3,7 / Ip3 )
Tz2,7 − T3,7 ≥ 0.3

(Case 8: When fault occurs at place 8)

Ak5
=T5,8 ( + Bk5 ) × TDS5
−1
Ck5
( I 4,8 / Ip4 )
Tz4,8 − T5,8 ≥ 0.3

32
(Relay time current characteristic curve selection constraints)

A=i 0.0515 + Ki ,1 × (19.61 - 0.0515) + Ki ,2 × (28.2 - 0.0515) for =

i 1, 2, 3, 4 and 5
Bi= 0.114 + Ki ,1 × (0.491 - 0.114) + Ki ,2 × (0.1217 - 0.114) for i = 1, 2, 3, 4 and 5
Ci = 0.02 + Ki ,1 × (2.0 - 0.02) + Ki ,2 × (2.0 - 0.02) for i =
1, 2, 3, 4 and 5
Ki ,1 + Ki ,2 ≤ 1 for i =
1, 2, 3, 4 and 5
Ki ,1 ,Ki ,2 =0 or 1

Results

The optimization problem is built in LINGO. Using MINLP function, the result is

shown in Table 3.6. In paper [7], the researchers used the same second zone time setting

for all the distance relays in the system. In this project, second zone time setting of each

distance relay is optimized separately. And after comparing the results of case 1 and case

2, it shows that with different second zone timing for each relay, the result is better than

the one using one second zone time setting for all the relays. The operation time of relays

can save up to 20.73%.

Table 3.6 – Results for the coordination of overcurrent and distance relay setting

Relay TDS Pickup TC curve 2nd zone Objective

Current timing Value
Relay 0.251 375 Extremely 0.388
1 Inverse
Case1 Relay 0.1 291.1 Extremely 0.323 0.9024
2 Inverse
Relay 0.1 60 Very Inverse 0.2
3
Relay 0.126 395 Extremely 0.351
4 Inverse
Relay 0.1 160 Extremely 0.2
5 Inverse

33
 Table 3.6 continued

Relay 0.289 375 Extremely 0.42

1 Inverse
Case 2 Relay 0.1 338.8 Extremely 0.42 1.1384
2 Inverse
Relay 0.1 60 Very Inverse 0.42
3
Relay 0.181 395 Extremely 0.42
4 Inverse
Relay 0.1 160 Extremely 0.42
5 Inverse

3.3 Fuse

Fuses act as both a protective and a disconnecting device, and are the earliest and

simplest overcurrent protection device in the power system. They basically consist of

metallic elements which melt in a time depending on the level of electric current when

the current exceeds a certain value. The size and construction of the element is

determined so that the heat produced for a normal current does not cause the element to

operate. During normal load conditions, the fuse must carry the normal operating current

without nuisance openings. However, when an overcurrent occurs the fuse must interrupt

the overcurrent and withstand the voltage across the fuse after internal arcing.

Fuse possesses one important advantage over mechanical interrupting devices

such as circuit breakers. It operates so quickly that it limits the possible damage to cable

and other equipment. The larger the current, the quicker the element melts. More

precisely, fuse has the ability to interrupt very large currents in a much shorter time – so

short that the current will be ‘cut off’ before it reaches its peak value, which is less than 5

34
ms in a 50Hz system. On the other hand, replacing a fuse is inconvenient, because it takes

longer to replace a fuse than to reclose a circuit breaker.

3.3.1 Characteristic Parameters

Rated Current I N

I N is the maximum current that the fuse can continuously conduct without interrupting

the circuit. The current rating of a fuse identifies its current carrying capacity based on a

controlled set of test conditions.

Rated Voltage VN

Voltage rating of the fuse must be greater than or equal to what would become the open

circuit voltage. Rated voltage should be larger than the maximum voltage source it would

have to disconnect.

Temperature

All electrical characteristics of a fuse are rated and validated at an ambient temperature of

25°C. Both higher and lower ambient temperatures will affect the fuse’s opening and

current carrying characteristics.

Melting Integral

∫ i dt is
2
called the Joule integral and is usually abbreviated to I 2t [1]. It is a most

convenient way of estimating the heating effect on the fuse due to the overcurrent.

Time Current Curve (TCC)

35
The time current characteristic curves of a protective device describe how fast the device

responds to the overcurrent condition. Time current characteristics are used in

coordinating multiple fuses and other protective equipment installed on the same

distribution feeder. Figure 3.7 shows the minimum melting time of a certain class of fuses

versus the fault current through the fuse. Similar curves provide the maximum melting

time of a fuse versus fault current through the fuse.

Figure 3.7 Minimum Melting Time and Total Clearing Time for a Specific Fuse [3]

Melting time

From the figure, we see the time required by a fuse to melt during a fault has a

distribution. The difference between the minimum and maximum times can be substantial.

Figure 3.3 shows a typical distribution between minimum and maximum melting times. It

is easy to understand the distribution in melting times if one considers the following two

facts: (a) the initial temperature of the fuse and its housing will vary and it will affect the

time at which the filament will reach the melting temperature, and (b) the size (cross

36
section) of the filament is not constant because of manufacturing imperfections and

therefore the thinner parts will reach melting temperature faster than the other parts. Once

the filament of the fuse melts, an electric arc is generated in the location of the filament.

The housing of the fuse element is designed to quench this electric arc. Various designs

are available for this purpose. Once the electric arc has been extinguished, then the fault

has been cleared. The time between the initiation of the fault and the extinguishing of the

electric arc is the total clearing time of the fuse [3]. The variability of the total clearing

time must be considered when fuses are coordinated with other protection devices.

3.3.2 Fuse Operation

In the project, the magnitude of the fault current will be provided for each

particular fuse at different locations. Given the fuse size and manufacture data, the time

current curve of the fuse is known. Our Objective is to determine whether and when the

fuse will melt.

A trial and error method is developed to achieve this objective using the following

procedure.

1. Given the fault time-current curve, initiate fuse starting melting time t1 , as is shown

in figure 3.8

2. Use time current curve for the particular fuse, t1 = f min ( I1 ) , calculate the value of I1 ,

as is shown in figure 3.9

t1

3. Compare the joule energy of ∫i

o
fault dt and I12t1

If the former is bigger than the latter, decrease t1 , and then go to step 2

37
 If the former is smaller than the latter, increase t1 , and then go to step 2

t1

4. When ∫i
o
fault dt equals I12t1 , then t1 is the starting melting time of the fuse

Figure 3.8 Fault time current curve

Total clearing time calculation follows the same step as the calculation of starting

melting time.

1. Initiate total clearing time t2 , and t2 > t1 , as is shown in figure 3.8

2. Use time current curve for the particular fuse, t2 = f max ( I 2 ) , calculate the value of I 2 ,

as is shown in figure 3.9

t 2

3. Compare the joule energy of ∫ i fault dt and I 22t2

o

3.1 If the former is bigger than the latter, decrease t2 , and then go to step 2

3.2 If the former is smaller than the latter, increase t2 , and then go to step 2
t 2

4. When ∫i
o
fault dt equals I 22t2 , then t2 is the starting melting time of the fuse, and

∆t = t2 − t1 is the fuse melting time

38
 Figure 3.9 Time current curve for fuse 40K

3.4 Coordination of Relays and Fuses

Fuse blow and fuse saving scheme are widely used in distribution system. For

fuse blow scheme, breaker is slow to trip so that fuse is allowed to blow for most faults,

as is shown in Figure 3.10. It is used primarily to minimize momentary interruptions

(reduce MAIFI) and is useful in high short circuit current areas and more suitable for

industrial type customers having high sensitive loads. For fuse saving scheme, breaker is

tripping on fast or instantaneous pickup to clear temporary fault without blowing fuse, as

is shown in Figure 3.11. It is used to minimize customer interruption time, and it can

reduce SAIDI but at the same time increase MAIFI. Fuse saving scheme works well in

most areas especially for residential and small commercial customers, however, not

suitable for certain industrial customers that cannot tolerate immediate reclosing.

39
 Many utilities use both schemes for a variety of reasons. They use fuse saving on

overhead and fuse blow on underground taps; they use fuse saving on rural and fuse blow

in urban; they use fuse save on stormy days and fuse blow on nice days; they use fuse

saving on some circuits and fuse blow on others depending on customer desires.

Figure 3.10 Fuse blow scheme

Figure 3.11 Fuse saving scheme

40
 In this project, the research is emphasized on the improvement of SAIDI in bad

weather conditions. So fuse saving scheme is recommended as most of the faults

occurred in severe weather conditions are momentary which are caused by lightning,

wind, tree branches and so on.

Take the test system in Figure 3.6 as an example, recloser 2 should operate faster

than the fuse in lateral to prevent the fuse from blowing. After putting the constraints into

the optimal relay setting algorithm, the new settings of the relays in the system will be

generated. A more detailed example is listed in the next section.

3.5 Test System

This section provides an overview of the adaptive relay programming approach

used for optimal relay setting for an example test system of IEEE 123 Node Test Feeder.

3.5.1 System Description

Figure 3.12 shows the original IEEE 123 node test system. The IEEE 123 node

test feeder operates at a nominal voltage of 13.8 kV. There are enough switches in the

feeder so that optimal configuration procedures can be tested. This is the most

comprehensive feeder and is characterized by:

1. Overhead and underground line segments with various phasing.

2. Unbalanced loading with all combinations of load types (PQ, constant I,

constant Z)

3. All loads are “spot loads” located at a node

4. Four step-type voltage regulators

41
 5. Shunt capacitor banks

6. Switching to provide alternate paths of power-flow

This feeder is well behaved and does not have a convergence problem. It provides

a test of the modeling of the phasing of the lines. The four voltage regulators provide a

good test to assure that the changing of individual regulator taps is coordinated with the

other regulators. A complete description of the ratings of all the equipment in the system

is given in Appendix B.

Figure 3.12 IEEE 123-Node Test Feeder

Figure 3.13 is the same test system built in WinIGS – F for simulation and

computation. Fuses and relays are added to the original system for the purpose of optimal

relay programming.

42
 Figure 3.13 IEEE 123-Node Test Feeder in WinIGS - F

3.5.2 Problem Formulation

The lower bond and upper bond of the pickup current of each relay can be read

from the simulation result of WinIGS – F, and is listed in Table 3.7. After simulating

faults at all the possible locations, as is shown in Figure 3.14, fault current at each relay is

recorded in Table3.8.

Table 3.7 Constraints of the pickup current in IEEE 123 node test system

Relay Lower Bond Upper Bond

1 300 500

2 100 400

3 150 300

43
 Figure 3.14 IEEE 123-Node Test Feeder with fault simulation

Table 3.8 Fault current in IEEE 123 node test system

Location Fault Current at Fault Current at Fault Current at Fault Current

R1 R2 R3 at Fuse

1 1002 945

2 1016 980

3 685 630

4 726 700

5 583 522

6 662 626 618

7 422 381 362

8 454 400 380

9 432 372 352

44
 The detailed equations of the problem are as follows:

min T1,1 +T1,2 +T1,3 +T1,4 +T2,5 +T2,6 +T3,7 +T3,8 +T3,9

Subject to

(Bounds on relay settings)

0.1 ≤ TDSi ≤ 20 i=1, 2 and 3

300 ≤ Ip1 ≤ 400

60 ≤ Ip2 ≤ 300
120 ≤ Ip3 ≤ 200

(Case 1: When fault occurs at place 1)

Ak1
=T1,1 ( + Bk1 ) × TDS1
−1
Ck1
( I1,1 / Ip1 )

(Case2: When fault occurs at place 2)

Ak1
=T1,2 ( + Bk1 ) × TDS1
−1
Ck1
( I1,2 / Ip1 )

(Case 3: When fault occurs at place 3)

Ak1
=T1,3 ( + Bk1 ) × TDS1
−1
Ck1
( I1,3 / Ip1 )

(Case 4: When fault occurs at place 4)

45
 Ak4
=T1,4 ( + Bk1 ) × TDS1
−1
Ck1
( I1,4 / Ip4 )

(Case 5: When fault occurs at place 5)

Ak2
=T2,5 ( + Bk2 ) × TDS 2
−1
Ck2
( I 2,5 / Ip2 )
Ak1
=T1,5 ( + Bk1 ) × TDS1
−1
Ck1
( I1,5 / Ip1 )
T1,5 − T2,5 ≥ 0.2

(Case 6: When fault occurs at place 6)

Ak2
=T2,6 ( + Bk2 ) × TDS 2
−1
Ck2
( I 2,6 / Ip2 )
Ak1
=T1,6 ( + Bk1 ) × TDS1
−1
Ck1
( I1,6 / Ip1 )
T1,6 − T2,6 ≥ 0.2

(Case 7: When fault occurs at place 7)

46
 Ak3
=T3,7 ( + Bk3 ) × TDS3
−1
Ck3
( I 2,7 / Ip2 )
Ak1
=T1,7 ( + Bk1 ) × TDS1
−1
Ck1
( I1,7 / Ip1 )
T1,7 − T3,7 ≥ 0.2

(Case 8: When fault occurs at place 8)

Ak3
=T3,8 ( + Bk3 ) × TDS3
−1
Ck3
( I 3,8 / Ip3 )
Ak1
=T1,8 ( + Bk1 ) × TDS1
−1
Ck1
( I1,8 / Ip1 )
T1,8 − T3,8 ≥ 0.2

(Case 9: When fault occurs at place 9)

Ak3
=T3,9 ( + Bk3 ) × TDS3
−1
Ck3
( I 2,9 / Ip3 )
Ak1
=T1,9 ( + Bk1 ) × TDS1
−1
Ck1
( I1,9 / Ip1 )
T1,9 − T3,9 ≥ 0.2

(Relay time current characteristic curve selection constraints)

47
A=i 0.0515 + Ki,1 × (19.61 - 0.0515) + Ki,2 × (28.2 - 0.0515) for =i 1, 2 and 3
Bi= 0.114 + Ki,1 × (0.491 - 0.114) + Ki,2 × (0.1217 - 0.114) for i = 1, 2 and 3
Ci = 0.02 + Ki,1 × (2.0 - 0.02) + Ki,2 × (2.0 - 0.02) for i = 1, 2and 3
Ki,1 + Ki,2 ≤ 1 for i =
1, 2 and 3
Ki,1,Ki,2 =0 or 1

(Fuse saving scheme)

T3,7 ≤ 0.06
T3,8 ≤ 0.05
T3,9 ≤ 0.06

3.5.3 Solution

The optimization problem is built in LINGO. Using MINLP formulation, the

result is shown in Table 3.9. Considering the fuse saving scheme, recloser 2 can react in

0.02s which is much faster than the fuse blow process. And recloser 3 is capable to

operate within 0.1s to avoid fuse blowing. With fuse saving process, reclosers and

breakers are more sensitive to fault current so that long term power outages are greatly

reduced at the sacrifice of increasing the momentary power outages.

48
Table 3.9 – Results for the optimal relay setting in IEEE 123 node system
Relay TDS Pickup TC curve Objective
Current Value
Relay 1 0.1 173.8 Extremely
Inverse
Relay 2 0.1 35 Extremely
0.9604
Inverse
Relay 3 0.1 70 Extremely
Inverse

49
 CHAPTER 4

AUTOMATIC SYSTEM RESPONSE

4.1 Introduction

Distribution system is the network which transfers the power from the

transmission lines to the load centers over long distance. The distribution network

consisting of feeders, transformers, loads, lines and protective devices is generally built

as interconnected mesh network. While in operation, the network is arranged in the form

of radial line structures, indicating that the systems are divided into a number of

subsystems of radial feeders. Each feeder is divided into several load sections with closed

sectionalizing switches and has connections to other feeders via several open tie switches.

The main purpose of the sectionalizing switches is the isolation or restoration of loads

during an outage or maintenance. As an outage occurs on the system, distribution

networks must be reconfigured as quickly as possible to restore as much out-of-service

loads as possible. Service restoration after an outage usually refers to an emergency

situation.

Implementation of smart grid technologies will result in (a) more switching

devices deployed in distribution systems and (b) a communication infrastructure that will

provide information on status of switches / breakers and will also enable remote control

of switches. The remote-controlled switches have become economically viable due to the

large amount of automation suppliers and the arrival of new communication technologies.

This arrangement will be ideal to support optimization and reconfiguration for the

purpose of restoring power to customers.

50
 Thus the objective of network reconfiguration is to maximize the amount of

power restored and reduce the amount of time to restore the energy supply while the

following constraints are satisfied,

(i) voltage constraints,

(ii) radiality constraints,

(iii) line losses

(iv) loading constraints and

Thus service restoration is a complicated combinatorial optimization problem

because there are a great number of switches in the distribution system. It may take a

long time to reach a feasible restoration plan which satisfies all the requirements.

Therefore, the dispatchers at many utilities tend to use their experience to narrow down

and reach a proper restoration plan in a short period.

4.2 Switch Placement Optimization

It is obvious that smart switches are required to improve power system reliability,

i.e., to reduce reliability indices. However, adding smart switches to the system is directly

related to the increase in costs. Therefore, there exists a trade-off between power

reliability improvement and costs. Then, the question is: What the optimal number of

smart switches that leads to the minimal cost? Here, optimization techniques come in to

play.

Figure 4.1 shows how the optimization of switch placement works. Initially, an

optimization algorithm generates a set of random switch locations and feed it to the

51
power model. Then, pre-defined, random fault profiles are applied to the power system,

which will eventually return back to normal, un-faulted conditions through the control

model’s switching algorithm. Once the power system returns to its normal conditions,

SAIDI savings, which is the reduction in SAIDI scores by adding new switches, and

related switch costs are calculated, both of which consist of the objective function. Based

upon the value of the objective function, the optimization algorithm updates switch

locations again and feeds them back to the power model until it reaches an optimum

placement or until it is terminated by some criteria. We use two heuristic optimization

techniques: Particle Swarm Optimization (PSO) and Genetic Algorithm (GA)

Optimization.

Figure 4.1 Switch Placement Optimization

Most general-purpose optimizations make use of gradient-based algorithms,

mainly due to their computational efficiency. However, non-gradient-based, probabilistic

optimization algorithms have attracted much attention from the research community.

52
Although these probabilistic optimization algorithms generally require many more

function evaluations to find an optimum solution, as compared to gradient-based

algorithms, they do provide several advantages. These algorithms are generally easy to

program, can efficiently make use of large numbers of processors, do not require

continuity in the problem definition, and are generally better suited for finding a global,

or near global, solution. In particular, these algorithms are ideally suited for solving

discrete and/or combinatorial-type optimization problems.

There are several modern heuristic optimization techniques such as evolutionary

computation, simulated annealing, tabu search, particle swarm, etc. Recently, genetic

algorithm (GA) and particle swarm optimization (PSO) techniques appeared as promising

algorithms for handling the optimization problems because of their versatility and ability

to optimize in complex multimodal search spaces applied to non-differentiable objective

functions. Both GA and PSO are similar in the sense that they are population-based

search methods and they search for the optimal solution by updating generations.

While GA is inherently discrete, i.e., it encodes the design variables into bits of

0’s and 1’s, and therefore easily handles discrete design variables, PSO is inherently

continuous and must be modified to handle discrete design variables.

4.2.1 Problem Formulation

Objective Function

53
 The objective function consists of two parts: one is to maximize SAIDI savings

that we can obtain by applying a new set of switches to the power system, and the other is

to minimize the cost related to those new switches. This can be expressed by

=
max f ( x) ( SAIDI _ Savings − Switch _ Cost )
x

Where x is a switch placement. The following steps describe how each term is derived:

1. Switch Cost is calculated based on the number of switches used and the number of

phases in each switch. It is assumed that the cost of a single-phase switch is $2,500,

the cost of a three-phase switch is $5,500, and the cost of a two-phase switch is twice

that of a single-phase switch, $5,000.

2. Total Possible SAIDI Minutes is calculated by dividing Total Sustained Fault Time in

Profile times NT by NT , where NT is the total number of customers served in an area.

3. Total SAIDI minutes saved is calculated by subtracting total SAIDI minutes

experienced from total possible SAIDI minutes.

4. Annual Average % SAIDI Reduction is calculated by dividing SAIDI Minutes Saved

by Total Possible SAIDI Minutes.

5. Annual Average % SAIDI Reduction is converted to a Total $ Savings amount by

using the estimated fault cost per SAIDI minute calculated from the Lawrence

Berkeley National Labs corresponding interest rate taken from current US t-bond rate

by varying the payback period, the weighting of switch cost vs. SAIDI benefits

changes.

54
6. The objective function is calculated as: Total $ Savings from SAIDI Reduction –

Total Switch Cost, and this value is to be maximized.

Representation of Switches

Basically, the IEEE Node-123 Test Feeder has 120 buses where a smart switch

can be added. Those 120 buses are represented by a 1×120 binary string, of which each

bit represents whether a switch is attached to the corresponding bus. If a bit is set to 0,

there is no additional switch on the bus; if a bit is set to 1, then a switch is added to the

bus. This binary string can be directly applied to a GA because it is inherently a discrete

optimization technique; however, since PSO is a continuous optimization technique, the

particles of a PSO need to be expressed in a binary string.

4.2.2 Particle Swarm Optimization

Particle swarm optimization (PSO) is one of the evolutionary computation

techniques motivated by the simulation of social behavior. PSO was developed by

Kennedy and Eberhart in the mid 1990s. It is a population-based search algorithm and is

initialized with a population of random solutions, called particles. PSO is similar to

genetic algorithm (GA) in that the system is initialized with a population of random

solutions. It is unlike GA, however, in that each potential solution is also assigned a

randomized velocity, and each particle flies over the search space at velocity dynamically

adjusted according to the historical behaviors of the particle and its companions.

Key Features of Particle Swarm Optimization [8]:

 It uses function values only.

55
 It requires relatively little fine tuning.

 It is easy to program.

 It requires a very large number of function evaluations. However, if approximation

techniques are used, this may be offset by the fact that function evaluations are very

fast.

4.2.2.1 PSO Algorithm

PSO makes use of a velocity vector to update the current position of each particle

in the swarm. The position of each particle is updated based on the social behavior that a

population of individuals, the swarm in the case of PSO, adapts to its environment by

returning to promising regions that were previously discovered [9]. The process is

stochastic in nature and makes use of the memory of each particle, as well as the

knowledge gained by the swarm as a whole. The outline of a basic PSO algorithm is as

follows:

1. Start with an initial set of particles, typically randomly distributed throughout the

design space.

2. Calculate a velocity vector for each particle in the swarm.

3. Update the position of each particle, using its previous position and the updated

velocity vector.

4. Go to step 2 and repeat until convergence.

The scheme for updating the position of each particle is xki +=

1 xki + vki +1

Where xki +1 is the position of the particle i at the iteration k+1 and vki +1 is the

56
corresponding velocity vector.

The scheme for updating the velocity vector of each particle depends on the

particular PSO algorithm under consideration. A commonly used scheme, introduced by

Shi and Eberhart [10], is

vki +1 = wvki + c1r1 ( p i − xki ) + c2 r2 ( pkg − xki )

Where r1 and r2 are independent random numbers between 0 and 1, p i is the best

position found by particle i so far, and pkg is the best position in the swarm at time k.

There are three problem-dependent parameters, the inertia of the particle w and two

‘trust’ parameters c1 and c2, which regulate the relative velocity toward local and global

best, respectively. The inertia controls the exploration properties of the algorithm, with

larger values facilitating a more global behavior and smaller values facilitating a more

local behavior. It is indicated in [10] that decreasing the inertia weight from about 0.9 to

0.4 during an optimization run provides improved performance in a number of

applications. The trust parameters indicate how much confidence the current particle has

in itself, c1, and how much confidence it has in the swarm, c2. In [10] it is proposed that

using c1=c2=2.05 for multimodal functions will provide improved optimization

performance.

Clerc in [11] indicates that use of a constriction factor may be necessary to insure

convergence of the particle swarm algorithm. A simplified method of incorporating a

constriction factor appears in the equation

57
vk=
i
+1 K  wv i
 k 11 + c r ( p i
− xk
i
) + c r
2 2 ( pk
g
− x i

k )

Where K is the a function of c1 and c2 as in the equation:

2
K= , ϕ =+
c1 c2 , ϕ > 4
2 − ϕ − ϕ − 4ϕ
2

A maximum allowable velocity vector Vmax clamps velocities of particles on each

dimension. If the acceleration causes the velocity on a dimension to exceed Vmax specified

by the user, then the velocity on that dimension is limited to Vmax .

Eberhart and Shi compared the performance of PSO with its different versions,

and concluded that the best approach is to use the constriction factor while limiting the

maximum velocity Vmax to the dynamic range of the variable X max on each dimension

[10]. In this work, the PSO algorithm described in [9] and [11] is implemented.

4.2.2.2 Implementation of PSO

Optimization Parameters

The inertia weight (w) is set to 0.9 at the beginning of the run, and made to

decrease linearly to 0.4 at the maximum number of iterations. Vmax is set to maximum

range X max . Each of the two (p-x) terms is multiplied by acceleration constants , c1 and

c2, of 2.05 (times a random number between 0 and 1). Also, Clerc’s constriction method

is used. ϕ is set to 4.1 and the constant multiplier K is thus 0.729, and each of the two (p-

x) terms being multiplied by 0.729×2.05=1.49445 (times a random number between 0

58
and 1). From [10], Vmax is set to be equal to X max . This significantly improves results

when using the constriction approach, creates the most consistent way to obtain good

results, and almost always creates the fastest optimization [10].

Initial Swarm

The initial swarm is generally created such that the particles are randomly

distributed throughout the design space, each with a random initial velocity vector. The

following equations are used to obtain the random initial position and velocity vectors:

x0i =X min + r3 ( X max − X min )

v0i =X min + r4 ( X max − X min )

where x0i is the initial position vector, v0i is the initial velocity vector of particle i, r3 and

r4 are independent random numbers between 0 and 1, X min is the vector of lower bounds,

and X max is the vector of upper bounds for the design variables.

The initial swarm distribution has an influence on the effectiveness of the PSO

algorithm; however, the influence of the initial swarm distribution is not important, as

long as it is fairly well distributed throughout the design space [9]. In this work, all initial

swarms were randomly distributed.

4.2.2.3 Procedures of the PSO

The following is the procedure of the PSO we implemented, from [12]:

1. Initialize a population of N particles. For the ith particle, its location xi in the search

59
space is randomly placed. Its velocity vector is vi in which the velocity in the dth

dimension is=
vid rand × Vmax , where is the random number in the range of (-1, 1).

2. Assign ϕ . Calculate K, and assign c=

1 c=
2 ϕ / 2.

3. Set the number of iteration k=1, and evaluate the fitness function for each particle. Let

pbest equal the fitness value of each particle. Let gbest equal the index of the particle which

gives best fitness value.

4. Compare the evaluated fitness value of each particle with its pbest . If the current value

is better than pbest , then set the current location as the pbest location.

5. Furthermore, if the current value is better than gbest , then reset gbest to the current

index in the particle array.

6. Change the velocity and location of the particles according to the equations,

respectively. Check if the velocity and location exceed the constraints.

7. Set k=k+1, repeat Steps 4-5 until the number of iteration is greater than the allowable

maximum iteration number Tmax .

4.2.2.4 Validation of PSO Algorithm

The implemented PSO algorithm was applied to a sample optimization problem

involving a two design variable, unconstrained mathematical function. The objective

function to be minimized is shown below, with the two design variables allowed to vary

between -50 and 10:

F ( x1 , x2 ) =
x12 − 100(cos x1 )2 − 100(cos x12 / 30) + x22 − 100(cos x2 )2 − 100(cos x22 / 30) + 1400

60
 Figure 4.2 Contour plot of sample objective function

This function has many local minima, and the global optimum has a value of 1000

at x=
1 x=
2 0 as shown in Figure 4.2 above. A swarm sizes of 20 particles, was

considered with c1=1.5 and c2=2.5, and the craziness operator was applied. The

optimization was repeated 50 times and the best, worst, mean, and standard deviation of

the best objective function and the number of function evaluations to convergence were

calculated for each of 50 repetitions. In addition, the reliability of the algorithm,

measured in terms of the success rate by counting the number of optimization runs that

found an objective function value within 1% of the best objective function value from all

50 runs, was also calculated. The results are summarized in Table 4.1 and Table 4.2.

Table4.1 - Number of function evaluation for continuous case

Mean Standard Deviation Best Worst

1966 240 300 2000

61
 Table 4.2 - Objective function for continuous case

Mean Standard Best Worst Success

Deviation Rate, %
1000.35 2.47 1000.00 1017.43 98.00

The results of Table 4.1 and Table 4.2 clearly show that the implemented PSO

algorithm is capable of solving the unconstrained, continuous example problem very

accurately. It is also able to deal with the many local minima present in the objective

function and finds the global optimum reliably.

Furthermore, to see if the implemented PSO can deal with integer design

variables like in our formulation of the switch placement optimization, the design

variables were converted to integers and the PSO algorithm was applied to the integer

design variable optimization problem. The simulation results are summarized in Table

4.3 and Table 4.4.

Table 4.3 - Number of function evaluation for integer case

Mean Standard Deviation Best Worst

1660 686 240 2000

Table 4.4 - Objective function for integer case

Mean Standard Best Worst Success

Deviation Rate, %
1005.24 18.02 1000.00 1123.03 80.00

As shown in Table 4.3 and Table 4.4, the implemented PSO can also deal with the

integer design variable reliably like as in the continuous case. Therefore, the implemented

PSO algorithm is validated to deal with both continuous and integer design variable

62
problems. Figure 4.3 shows how particles fly over the design space to get to the optimum

value at x=
1 x=
2 0 . The blue circles represent the initial particles and the red crosses

represent the final position of the particles. It can be seen that particles, which are

initially scattered in the design space, flock together at the optimum point through the

implemented PSO algorithm.

Figure 4.3 Validation of PSO Algorithm

4.2.2.5 Verification of the PSO Algorithm for Switch Placement Optimization

In order to verify the implemented PSO algorithm, two simple cases are

optimized; the first involves adding 1 switch to the system and the second involves

adding 2 switches in the system. To see if the PSO algorithm is capable of finding a

“near-optimal” switch placement, the theoretical optimal solution is required. Since there

are only 120 possible ways to add 1 switch and 7140 ways to add 2 switches, to find a

theoretical optimum for both cases is quite simple by brute force, and thus the switch

placement optimizations of 1 switch and 2 switches are considered to verify the

63
implemented algorithm. The other cases are not considered because due to exponentially

increasing simulation times; for example, if 3 switches are to be added, then there are

over 280 thousand cases to be run to find a theoretical optimum solution, as calculated

below.

120!
= 280840
(120 − 3)!3!

For the purpose of finding the theoretical optimum, the objective function, SAIDI

reduction, and cost for all the possible switch location were first simulated by adding a

switch or two switches to all possible buses of the IEEE Node-123 Test Feeder and

running the power model. Table 4.5 compares both cases.

Table 4.5 - Comparison of 1-switch and 2-switch simulations

No. of No. of possible No. of Fault Payback Interest Switch Cost

switches placement Profiles Applied Period Rate
1 120 119 1 year 0.30% $2,500/1-phase
$5,000/2-phase
2 7140 1 1 year 0.30% $5,500/3-phase

Table 4.6 summarizes the objective function (i.e. FVAL), SAIDI reduction, and

switch cost against the switch location for the 1-switch case (Payback Period: 1 year;

Interest Rate: 0.30%). Note that Table 4.6 shows the benefits, SAIDI reduction, and cost

which can be obtained by adding a switch to all the possible bus location, not

optimization results. As seen in Table 4.6, for the 1-switch case, adding a switch to Bus

58 has the most benefit, that is, Bus 58 is the theoretical, optimum location of a new

switch.

64
 Table 4.6 - Fval, SAIDI Reduction, & Cost of 1-Switch Verification Case

PO FVAL SAIDI(% SAIDI($ COS PO FVAL SAIDI(% SAIDI($ COS

S ) ) T S ) ) T
1 -5458.77 0.0003 41.23 5500 61 -2436.66 0.0005 63.34 2500
2 13137.8 0.1357 18637.83 5500 62 -2435.53 0.0005 64.47 2500
3
3 -2435.53 0.0005 64.47 2500 63 286.75 0.0421 5786.75 5500
4 3057.85 0.0405 5557.85 2500 64 63.33 0.0405 5563.33 5500
5 -2436.65 0.0005 63.35 2500 65 62072.1 0.4919 67572.13 5500
3
6 3149.4 0.0411 5649.4 2500 66 154.88 0.0412 5654.88 5500
7 3195.21 0.0415 5695.21 2500 67 246.42 0.0418 5746.42 5500
8 14032.3 0.1422 19532.36 5500 68 256.84 0.0419 5756.84 5500
6
9 23778.0 0.2131 29278.06 5500 69 276.34 0.042 5776.34 5500
6
10 3057.85 0.0405 5557.85 2500 70 43031.9 0.3533 48531.92 5500
2
11 -2435.53 0.0005 64.47 2500 71 29426.0 0.2542 34926.06 5500
6
12 3149.4 0.0411 5649.4 2500 72 3012.05 0.0401 5512.05 2500
13 52626.6 0.4231 58126.65 5500 73 3057.85 0.0405 5557.85 2500
5
14 63212.4 0.5002 68712.42 5500 74 3149.4 0.0411 5649.4 2500
2
15 3057.85 0.0405 5557.85 2500 75 3195.21 0.0415 5695.21 2500
16 3195.21 0.0415 5695.21 2500 76 39278.6 0.326 44778.63 5500
3
17 -2435.53 0.0005 64.47 2500 77 -2441.16 0.0004 58.84 2500
18 3195.21 0.0415 5695.21 2500 78 -2438.91 0.0004 61.09 2500
19 -2435.53 0.0005 64.47 2500 79 -2436.65 0.0005 63.35 2500
20 22990.4 0.2074 28490.46 5500 80 -262.59 0.0381 5237.41 5500
6
21 19494.2 0.1819 24994.28 5500 81 22186.7 0.2015 27686.75 5500
8 5
22 -2438.91 0.0004 61.09 2500 82 -171.04 0.0388 5328.96 5500
23 -2436.65 0.0005 63.35 2500 83 -5436.65 0.0005 63.35 5500
24 16733.3 0.1618 22233.37 5500 84 -79.5 0.0395 5420.5 5500
7
25 -2436.65 0.0005 63.35 2500 85 12.05 0.0401 5512.05 5500
26 14946.0 0.1488 20446.01 5500 86 -5437.78 0.0005 62.22 5500
1
27 -2436.65 0.0005 63.35 2500 87 -2437.78 0.0005 62.22 2500
28 1000.21 0.0473 6500.21 5500 88 -5486.66 0.0001 13.34 5500
29 6146.1 0.0811 11146.1 5000 89 -2436.65 0.0005 63.35 2500
30 -4936.65 0.0005 63.35 5000 90 16733.3 0.1618 22233.37 5500
7
31 3195.14 0.0415 5695.14 2500 91 13158.6 0.1358 18658.65 5500
5
32 -2436.65 0.0005 63.35 2500 92 -2436.65 0.0005 63.35 2500

65
 Table 4.6 continued

33 -1144.66 0.0317 4355.34 5500 93 11371.29 0.1228 16871.29 5500

34 -3289.53 0.0161 2210.47 5500 94 -2436.65 0.0005 63.35 2500
35 -5434.4 0.0005 65.6 5500 95 4220.42 0.0708 9720.42 5500
36 -2435.53 0.0005 64.47 2500 96 -2436.65 0.0005 63.35 2500
37 3149.4 0.0411 5649.4 2500 97 2254.3 0.0564 7754.3 5500
38 15583.37 0.1535 21083.37 5500 98 -2436.65 0.0005 63.35 2500
39 603.59 0.0408 5603.59 5000 99 3240.95 0.0418 5740.95 2500
40 -2436.65 0.0005 63.35 2500 100 1000.21 0.0473 6500.21 5500
41 -2436.66 0.0005 63.34 2500 101 9567.63 0.1097 15067.63 5500
42 -2435.53 0.0005 64.47 2500 102 -1144.66 0.0317 4355.34 5500
43 10461.16 0.1162 15961.16 5500 103 -3289.53 0.0161 2210.47 5500
44 3240.95 0.0418 5740.95 2500 104 -5434.4 0.0005 65.6 5500
45 -2436.65 0.0005 63.35 2500 105 5013.1 0.0765 10513.1 5500
46 10866.46 0.1191 16366.46 5500 106 -2440.04 0.0004 59.96 2500
47 5611.41 0.0809 11111.41 5500 107 -2438.91 0.0004 61.09 2500
48 3195.14 0.0415 5695.14 2500 108 -2436.65 0.0005 63.35 2500
49 3240.95 0.0418 5740.95 2500 109 -33.76 0.0398 5466.24 5500
50 -5435.12 0.0005 64.88 5500 110 3103.66 0.0408 5603.66 2500
51 5698.07 0.0815 11198.07 5500 111 -2436.65 0.0005 63.35 2500
52 5736.58 0.0818 11236.58 5500 112 -5434.4 0.0005 65.6 5500
53 5917.42 0.0831 11417.42 5500 113 2966.24 0.0398 5466.24 2500
54 286.75 0.0421 5786.75 5500 114 3057.79 0.0405 5557.79 2500
55 63212.61 0.5002 68712.61 5500 115 3240.95 0.0418 5740.95 2500
56 63212.8 0.5002 68712.8 5500 116 -2438.91 0.0004 61.09 2500
57 195.14 0.0415 5695.14 5500 117 -2437.78 0.0005 62.22 2500
58 63537.04 0.5026 69037.04 5500 118 -2435.53 0.0005 64.47 2500
59 240.95 0.0418 5740.95 5500 119 -2434.4 0.0005 65.6 2500
60 63358.35 0.5013 68858.35 5500 120 -5434.4 0.0005 65.6 5500

For the 2-switch case fifty-six different 2-switch combinations show the maximal

economic benefits, which is about $120K under the assumption of 1-year payback period

and an interest rate of 0.30%. Due to its enormous size, the simulation result for finding a

theoretical optimum of the 2-switch case is not included in this report. In order to reduce

the optimization run time, a simplified set of fault scenarios was applied to the 2-switch

case while a comprehensive set of 119 fault profiles was applied to the 1-switch case. For

both cases, 100 optimization runs were performed.

Optimization Results of the 1-Switch Case

66
 The verification results are summarized in Table 4.7 (Payback Period: 1 year;

Interest Rate: 0.30%) and Table 4.8 (Payback Period: 1 year; Interest Rate: 0.30%). Of

100 runs, 98 runs yield Bus 58 as an optimum location of switch. Two runs found Bus 60

as an optimum location for a new switch; in other words, the success rate is 98%.

However, Bus 60 is a near-optimum location as indicated by the standard deviation of

$25.14 in Table 4.8.

Table 4.7 - Number of Function Evaluations of 1-Switch PSO Optimization

Mean Standard Deviation Best Worst

679,728 152,199.75 523,600 1,142,400

Table 4.8 - Objective Function of 1-Switch PSO Optimization

Mean Standard Best Worst Success Rate, % Optimum

Deviation Location
$63,533.47 $25.14 $63,537.04 $63,358.35 98.00 Bus#58

Since 119 fault profiles were applied to the SAIDI reduction calculation of each

particle, 119 function evaluations are required to evaluate a particle’s objective function.

Hence, 2380 (=119*20) function evaluations are required for 20 particles. This is why the

number of function evaluations is so large. Based upon these results, even though the

standard deviation of function evaluations is quite large, the implemented PSO algorithm

is capable of finding an optimum switch placement with great accuracy.

Optimization Results of the 2-Switch Case

The optimization results of the 2-switch case are summarized in Table 4.9

(Payback Period: 1 year; Interest Rate: 0.30%) and Table 4.10 (Payback Period: 1 year;

67
Interest Rate: 0.30%).

Table 4.9 - Number of Function Evaluations of 2-Switch PSO Optimization

Mean Standard Deviation Best Worst

511.2 81.7150 440 920

Table 4.10 - Objective Function of 2-Switch PSO Optimization

Mean Standard Best Worst Success Rate, % Optimum
Deviation Location
$120,760.00 $184.00 $120,792.20 $ 119,718.99 97.00 56 Unique
Switch
Combinations

Like the 1-switch case, the implemented PSO algorithm finds an optimum

solution very accurately. Of 100 optimization runs, only 3 runs found very different

locations which, however, are near-optimum solutions as indicated by the standard

deviation of $184.00 in Table 4.10.

4.2.2.6 Conclusion of the PSO Algorithm Verification

The IEEE 123 Node Test Feeder is similar to the simple mathematical function in

that it has multiple local maxima; however, it differs from the simple mathematical

function in the sense that it doesn’t have near-optimal solutions around its local optima.

Therefore, it was not possible to observe that each particle representing switch locations

flies towards the optimum solution, and it is also not possible to display simulation

results in a figurative way as shown in Figure 4.4 (red line: optimal location, black circles:

best locations at every iteration, blue star: particles’ positions at every iteration). The

implemented PSO algorithm shows a behavior similar to one of a general random search.

That’s the reason in this work the inertia coefficient w was changed to 1.4, not

68
dynamically decreased from 0.9 to 0.4, in order to ensure a global search.

Figure 4.4 Verification of PSO Algorithm with 1 Switch

Initial Optimization Results of PSO Algorithm

First, for simplicity, only 10 random, but pre-defined, fault profiles were applied

and the switch placement was optimized using the implemented PSO algorithm under the

assumption of 2-year payback period and the interest rate of 0.61%.

To find the optimum switch placement, the PSO algorithm optimized the switch

placement for each number of switches varying from 2 to 15. After the optimum switch

placement for each number of switches had been obtained, the optimum number of

switches which provides the highest benefits was determined off-line among those 14

69
optimization results. Table 4.11 (Payback Period: 2 years; Interest Rate: 0.61%)

summarizes the optimization results with 10 fault profiles.

Table 4.11 - Optimization Results of PSO with 10 Fault Profiles

No. of Switches Benefits ($) SAIDI Reduction

(%)
2 159,290 62.07
3 169,010 67.62
4 171,890 69.58
5 179,450 73.25
6 182,620 77.51
7 184,720 79.18
8 176,960 77.08
9 179,180 80.08
10 189,190 84.46
11 171,590 79.13
12 169,860 81.60
13 163,680 79.17
14 172,340 83.23
15 172,990 88.57

The PSO algorithm yielded 10 switches as an optimum number of switches that

gives the highest benefit, $189,190; however, 15 switches gives better SAIDI reduction.

The reason why 15-switch provides the more SAIDI reduction with the lower benefits is

that the increase in the number of switches added to the power system not only reduces

the SAIDI more but also increase the switch cost. It is obvious as shown in Table 4.12:

the addition of 120 switches dramatically reduces the SAIDI score, which is almost 98%

reduction, but the switch cost also dramatically increases so as to have an adverse effect

on the benefits.

70
 Table 4.12: Objective Function, SAIDI Reduction, and Cost of 120-Switch Case

Payback Interest Rate Benefits ($) SAIDI SAIDI Switch Cost

Period Reduction Reduction ($)
(%) ($)
1 Year 0.30% -356,032.94 97.88 134,467.06 490.500.00
2 Years 0.61% -223,207.28 97.88 267,292.72 490,500.00

Figure 4.5 (Payback Period; 2 years; Interest Rate: 0.61%) and Figure 4.6

(Payback Period; 2 years; Interest Rate: 0.61%) below show the benefits (objective

function) and SAIDI reduction against the number of switches with 10 fault profiles

applied, respectively. It can be also seen in Figure 4.6 that the more switches are added to

the system, the more reduction in SAIDI scores can be achieved.

Figure 4.5 Objective Function of 1-Switch PSO Optimization

71
 Figure 4.6 SAIDI Reduction of 1-Switch PSO Optimization

4.2.3 Genetic Algorithm

Genetic algorithms are an optimization method that employs a search process

imitated from the mechanism of biological selection and biological genetics. They

combine survival of the fittest among those feasible solutions in the form of string

structures (or genes: in binary form), and a randomized formation exchange to form a

search algorithm [13]. In every generation, a new set of solutions is created from the

fittest of the previous set of solutions.

The control variables have to be represented as strings. During the search process,

solutions are initially randomized. Then, the fitness of each solution is computed. The

solution with higher fitness has higher probability to be chosen for generating a new

72
generation. This procedure is called reproduction or selection. A crossover is used for

innovating of solutions, and a mutation can help solutions have a wider range of feasible

solutions. After these three genetic operations, the new generation is obtained and it starts

the genetic operations again and again until the optimum solution is found.

Key Features of Genetic Algorithms [9]:

 It uses function values only.

 It naturally handles discrete variable.

 It is easy to program.

 It requires a very large number of function evaluations. However, if approximation

techniques are used, this may be offset by the fact that function evaluations are very

fast.

4.2.3.1 Implementation of a GA

Representation

The representation scheme determines how the problem is structured in a GA and

also determines the genetic operators that can be used. For this study, a binary

representation is used. IEEE Node-123 Test Feeder model has 120 line segments on

which a smart switch is possibly added. We represent these 120 segments by a 120×1

vector of which each bit represents whether a switch is added on a line segment or not. If

a bit is set to 1, then the corresponding line segment will have a new switch; otherwise, it

has no switch.

Initialization of a Population

73
 A GA must be provided with an initial population. In the initial stage of the GA, a

matrix of random switch locations with the number of rows equal to the population size

and the number of columns equal to the number of the line segments on which a new

switch can be added is created. A total population size of 80, which was recommended in

[14], is used.

Selection

The selection determines which of the individuals will survive and continue on to

the next generation. The GA carries out the selection each generation after all the new

children have been evaluated to create a new population from the old one. In this work,

two selection techniques are considered: one is the Elitest technique and the other is to

use a Roulette Wheel. The Elitest technique is to choose the predefined number of

individuals that have the best fitness values. On the other hand, the Roulette Wheel

technique uses a roulette wheel, each area of which is proportional to each individual’s

fitness, selects individuals on the portions of which randomly generated numbers fall.

Thus, it prevents the algorithm from being trapped in local optima. Figure 4.7 illustrates

how differently these two techniques work.

Figure 4.7 Selection - Elitest vs. Roulette Wheel

74
Reproduction

To produce a new generation, two genetic operators, two-point crossover and

non-uniform mutation, are used. From Figure 4.7, it is best in an off-line GA optimization

for 90% of the new generation to be created by crossover and mutation, while 10% be

carried over exactly from the previous generation.

Crossover

Crossover takes two individuals and produces two new as individuals, as

illustrated in Figure 4.8, where and are randomly generated numbers. From Figure 4.7,

45% of the members of the next generation should be created by crossover in an off-line

GA optimization.

Figure 4.8 Two-Point Crossover

Termination

The termination determines when to stop the simulated evolution and return the

resulting switch placement. A maximum generation criteria is used to stop a simulation.

A maximum generation of 50 is used.

Mutation

Mutation alters one individual to produce a single new solution as illustrated in

Figure 4.9 below. From Figure 4.7, the probability of any specific bit being mutated is set

75
to 0.01.

Figure 4.9 Mutation

4.2.3.2 Procedures of the GA

The basic genetic algorithm implemented in this work consists of the following

steps [8]:

1. Create a random initial population.

2. Calculate all fitness values.

3. Get the sum of all the fitness values Fsum = ∑ Fi .

4. Construct a roulette wheel, with each binary string occupying an area on the wheel in

proportion to the ratio Fi / Fsum .

5. Send 10% of the current population to the next generation unaltered, using the roulette

wheel to choose, and always including the fittest individual. Note: Here a combination of

the ‘roulette wheel’ and elitist strategies is employed.

6. For the remaining members of the next generation, use a random number 0-1 to pick

pairs on the wheel as ‘mating pairs’ that will reproduce.

7. Perform crossover. Use a weighted coin toss to pick the probability of crossover, with

the probability of crossover being 0.45.

8. If crossover is dictated, pick two integer numbers between 1 and the length of the

76
binary string to establish the starting and ending crossover locations. Exchange values in

the string between parents.

9. Perform the mutation operation on the child, with the probability of mutation being

0.01 for every bit.

10. Repeat the process until convergence is achieved.

4.2.3.3 Comparison of PSO and GA

Table 4.13 summarizes how differently the PSO and the GA optimization

algorithm are implemented in this work. Both of them are applied to the IEEE Node-123

Test Feeder and have the same objective function, maximize the difference between

SAIDI savings and related switch cost. However, the PSO can optimize only switch

locations with the fixed number of switches because, if the number of switches varies, it

can’t update its particles’ velocities and positions during optimization runs while the GA

can optimize both switch locations and the number of switches simultaneously. In other

words, the GA can take into account reliability indices and cost at the same time, but the

PSO can consider cost off-line after all optimizations for a various number of switches

are done. Due to the same reason, the representations of design variables are different.

For the PSO, a 120×1 binary string can have only the fixed number of 1’s. On the other

hand, the GA has a 120×1 binary string which can have varying number of 1’s. The same

convergence criteria are applied to both techniques: an optimization process is terminated

if it seems to converge to an optimal switch placement or if the maximum number of

iterations is reached. The outcome of the optimization process is optimized switch

location for the PSO and optimized number of switches and their locations for the GA.

77
 Table 4.13 - Comparison of PSO and GA

4.2.4 Results and Conclusions

From our genetic algorithm, the result given by the best found switch placement is

shown in Table 4.14 below.

Table 4.14 - Optimal GA Switch Configuration Results

Payback Interest Net Present SAIDI Annual Switch

Period Rate Value of Reduction SAIDI Costs ($)
Benefits ($) (%) Savings ($)
1 Year 0.30% 84,238.50 81.3 111,404.50 27,500

Five switches need to be added to our baseline system, with the locations chosen

by the GA as shown in Figure 4.10 below, in order to achieve to above results.

78
 Figure 4.10 Optimal switch locations for 1 year payback period.

4.3 Automatic Network Reconfiguration

4.3.1 Rule-based Reconfiguration

Basically, the rule-based reconfiguration is done by finding the switches near to a

fault, operating them to isolate the faulted zone and to restore power to unfaulted zone.

To increase the efficiency to search for nearest switches, a data tree structure for the

power model is constructed by applying the object-oriented programming. The tree

structure is updated every iteration of the power model and identifies the network

configuration at an iteration. The search algorithm starts searching from the location of

fault, which is assumed to be provided by measuring devices. Once the switches nearest

79
to the fault are identified, the reconfiguration algorithm operates them to isolate the

faulted zone and restores power to the other zones as the ideal case that is briefly

explained in the previous section.

Rule-based Reconfiguration Algorithm

The very first step of the rule-based reconfiguration algorithm is to construct a

tree data structure. Each node in the network is represented as a struct (short for

“structure”). For example, the struct for the node 13 is shown in Figure 4.11.

Figure 4.11 Data Structure for Node 13

A node struct has 6 fields: ‘id’, ‘parent’, ‘child’, ‘xpos’, ‘ypos’, and ‘nCustomer’.

The feidls ‘xpos’ and ‘ypos’ are used to plot the grid. The field ‘nCustomer’ contains the

number of customers at a node. The 5 switches are denoted as ‘S1’,’S2’,’S3’,’S4’, and

‘S5’, respectively. The node 149, which is connected to the main feeder (node 150), is

considered as a root node for the whole tree data structure. Suppose that there are two

nodes, N1 and N2, and that they are connected to each other. If node N1 is closer to Node

150, then the node N1 is a parent of the node N2, and N2 is a child of N1. In this way, all

parent-child relationship can be determined because the baseline network is radial. Also,

for the radiality to be held, all nodes must have only one parent. If a node has more than

one parent, the network would have a loop.

80
 Once a fault occurs, the algorithm first starts to search for the switches nearest to

the fault. In order to do this, the algorithm invokes the path-finding routine, which

determines the paths of the switches to the fault and also determines the switch

commands. An example of the path-finding algorithm is shown in Figure 4.12, where

switches S1, S2, S3, and S9 are tie switches, and switches S4 to S8 are automatic

switches. The switches are located as shown in Figure 4.13. This example shows that a

fault occurs at node 19. The algorithm searches for an upstream path from the fault and

terminates the searching process at a node with a switch or at the root node, 149. If there

doesn’t exist a path from the fault to the root node, in other words, the searching process

is terminated at a node with a switch, the fault can be assumed to have occurred

somewhere below the switch from the viewpoint of the tree structure. In this example, the

searching process is terminated at the switch S4, which indicates that the fault has

occurred at the lower level than S4 and that the main feeder can be utilized after the fault

is isolated.

Figure 4.12 Paths finding and switch commands determination

81
 Let the path from the fault to Switch S4 denoted by P0. In case of the switch S5, it

has a common node with the path P0, which means that the switch S5 must be opened to

isolate the fault. The remaining switches, S6 to S8, are determined to remain closed

because they don’t have any common node with the path P0. After the switching

commands for the automatic switches are determined, the algorithm determines the

commands for the tie switches. To avoid a looped network, any two feeders cannot be

utilized at the same time. According to the path table, the tie switch S1 has a path to the

root node 149, which indicates that S1 cannot be connected to the network. In the same

way, the tie switch S3 should remain opened. For the tie switch S2, it has a common node

with the path P1 and hence it must also remain opened. For the last tie switch S9, it is

connected to the switch S5. Since the switch S5 needs to be opened, the switch S9 should

be closed in order to provide power to the area containing nodes 35 to 51.

Figure 4.13 Optimal switch locations for automatic reconfiguration

82
 In summary, the algorithm reconfigures the network as follows:

1) Identify an upstream path, P0, from the fault: ex) P0 = {‘19’,‘18’,‘S4’}

2) Set the command for the switch S4 to ‘OPEN’: ex) SwCmd(S4) = 0

3) Search for a path from the switches until it reaches the fault or the root node 149 and

denote it as P1, P2,…, P9, respectively: ex) P7 = {‘S7’,‘67’,‘60’,‘S6’}

4) If the path Pi, where i=1,…,9, has a common node with the path P0, then set the

command for the corresponding switch to ‘OPEN’: ex) SwCmd(S5) = 0

5) If not, set the command to ‘CLOSED’: ex) SwCmd(S6) = 1

6) If a tie switch has a path to the root node, set the command to ‘OPEN’: ex)

SwCmd(S1) = 0

Figure 4.14 illustrates how the algorithm reconfigures the network when a fault

occurs at the node 19. As explained previously, the switches S4 and S5 are opened to

isolate the fault, and the switches S6, S7, and S8 remain closed. The tie switches S1 and

S3 remain open to make the network radial. The tie switch S2 is opened since it is

connected to the faulty area. Finally, the tie switch S9 gets closed to provide power to the

area containing nodes 35 to 51; otherwise, the area loses power so that it degrades the

SAIDI index score.

83
 Figure 4.14 Network Plot @ t=12s with a Fault at Node 19

4.3.2 Simulation Results of the Optimal Reconfiguration

To verify each reconfiguration technique, 119 fault scenarios are applied over

which SAIDI scores are averaged so that SAIDI reductions of each technique can be

compared with that of the baseline model that is the one with no automatic switches

added. And SAIDI can be improved as high as 82.19%. Table below shows the effect of

types of switches on Power system reliability.

Table 4.15 - Effects of Types of Switches on Reliability

No. of Automatic No. of Manual Averaged, Best

Switches Switches SAIDI Reduction
0 5 0%
2 3 69.65%
3 2 75.86%
4 1 78.07%
5 0 82.19%
5 1 82.94%
5 2 83.27%
5 3 83.78%

84
 CHAPTER 5

CONCLUSION AND FUTURE WORK

The work performed in this thesis focused on the reliability enhancement in

power distribution systems. An algorithm was built to realize optimal relay setting for

several test systems. Both overcurrent relays and distance relays and their coordination

were considered in the optimization problem. In addition, fuse saving scheme was

applied to avoid long-time power outage.

After the fault occurrence, another algorithm was developed to realize optimal

network reconfiguration as quickly as possible to restore as much out-of-service loads as

possible. Some optimization technique was applied to find switches to isolate the fault

and restore power to non faulted sections via sectionalizing switches and open tie

switches. Meanwhile, optimal switch placement algorithm was created to decide the

number and location of additional switches to improve the reliability of the system at a

reasonable cost.

Preliminary results were obtained from the IEEE 123-Node Test Feeder. The

optimization technique can also be used in a realistic distribution system with the

combination of AMI and state estimation. For example, a typical distribution system

usually has 40 to 50 miles of circuits, or a whole substation usually has 8 to 12

distribution feeders. Advanced metering infrastructure (AMI) will be used to realize a

high fidelity real-time monitoring system and a three-phase state estimation algorithm

which filters the measurements to give the real-time model of the system. The results of

the state estimator are used as an input to the optimization algorithm responsible for

generating real-time control signals through which all the available devices will be

coordinated to achieve system level optimal operating conditions. In this case, nominal

85
current can be received from state estimation. Those values will then be transformed into

constraints in the optimization problem. According to the weather information, the

potential area that might be affected by the severe weather will be determined. Since the

location and status of the relays and switches in the area are known in advance, the

optimization algorithm will then decide the optimal relay setting and optimal network

reconfiguration. With the results from the optimization algorithm, the settings of the

relays in the distribution system are updated accordingly and the load break switches and

tie switches are operated in sequence through remote control system so that the system
can react to the fault the fastest and restore power to as many customers as possible.

86
 APPENDIX A

BASIC RELIABILITY CALCULATIONS

System average interruption duration index (SAIDI) is the average duration of all

interruptions per utility customer served. It is determined by dividing the sum of all

customer interruption durations in a year by the number of customers served.

total duration of sustained outages

SAIDI =
number of customers served

System Average Interruption Frequency Index (SAIFI) is the average number of

times that a system customer is interrupted during a time period. It is determined by

dividing the sum of all customer interruption durations during a year by the number of

customers served.

total number of sustained outages

SAIFI =
number of customers served

Customer Average Interruption Duration Index (CAIDI) is the average duration of

interruptions among those customers experiencing interruptions during a year. It is

determined by dividing the sum of all customer interruption durations by the total number

of interruptions over 1-year period.

total duration of sustained outages

CAIDI =
number of sustained outages

87
 Momentary Average Interruption Frequency Index (MAIFI) is the average

number of momentary interruptions that a customer experiences during a year. Electric

power utilities may define momentary interruptions differently, with some considering a

momentary interruption to be an outage of less than 1 minute in duration while others

may consider a momentary interruption to be an outage within several seconds.

number of momentary interruptions

MAIFI =
number of customers served

88
 APPENDIX B

DESCRIPTION OF THE TEST SYSTEM

Figure B.1 shows the 13.8 KV 123 node test system. The ratings of each of the

units are given in the following tables and figures.

Figure B.1 IEEE 123-Node Test Feeder

89
 Table B.1 Line segment data
Node A Node B Length (ft.) Config. 42 44 200 1
1 2 175 10 44 45 200 9
1 3 250 11 44 47 250 1
1 7 300 1 45 46 300 9
3 4 200 11 47 48 150 4
3 5 325 11 47 49 250 4
5 6 250 11 49 50 250 4
7 8 200 1 50 51 250 4
8 12 225 10 52 53 200 1
8 9 225 9 53 54 125 1
8 13 300 1 54 55 275 1
9 14 425 9 54 57 350 3
13 34 150 11 55 56 275 1
13 18 825 2 57 58 250 10
14 11 250 9 57 60 750 3
14 10 250 9 58 59 250 10
15 16 375 11 60 61 550 5
15 17 350 11 60 62 250 12
18 19 250 9 62 63 175 12
18 21 300 2 63 64 350 12
19 20 325 9 64 65 425 12
21 22 525 10 65 66 325 12
21 23 250 2 67 68 200 9
23 24 550 11 67 72 275 3
23 25 275 2 67 97 250 3
25 26 350 7 68 69 275 9
25 28 200 2 69 70 325 9
26 27 275 7 70 71 275 9
26 31 225 11 72 73 275 11
27 33 500 9 72 76 200 3
28 29 300 2 73 74 350 11
29 30 350 2 74 75 400 11
30 250 200 2 76 77 400 6
31 32 300 11 76 86 700 3
34 15 100 11 77 78 100 6
35 36 650 8 78 79 225 6
35 40 250 1 78 80 475 6
36 37 300 9 80 81 475 6
36 38 250 10 81 82 250 6
38 39 325 10 81 84 675 11
40 41 325 11 82 83 250 6
40 42 250 1 84 85 475 11
42 43 500 10 86 87 450 6
42 44 200 1 87 88 175 9
44 45 200 9 87 89 275 6

90
 Table B.1 continued
89 90 225 10 Table B.3 – Overhead Line Configuratioions
89 91 225 6
91 92 300 11 Overhead Line Configurations (Config.)
91 93 225 6
93 94 275 9 Config. Phasing Phase Cond. Neutral Cond. Spacing
93 95 300 6 ACSR ACSR ID
95 96 200 10 1 A B C N 336,400 26/7 4/0 6/1 500
97 98 275 3 2 C A B N 336,400 26/7 4/0 6/1 500
98 99 550 3 3 B C A N 336,400 26/7 4/0 6/1 500
99 100 300 3 4 C B A N 336,400 26/7 4/0 6/1 500
100 450 800 3 5 B A C N 336,400 26/7 4/0 6/1 500
101 102 225 11 6 A C B N 336,400 26/7 4/0 6/1 500
101 105 275 3 7 A CN 336,400 26/7 4/0 6/1 505
102 103 325 11 8 A BN 336,400 26/7 4/0 6/1 505
103 104 700 11 9 AN 1/0 1/0 510
105 106 225 10 10 BN 1/0 1/0 510
105 108 325 3 11 CN 1/0 1/0 510
106 107 575 10
108 109 450 9
108 300 1000 3
109 110 300 9 Table B.4 – Underground Line Configuration
110 111 575 9
110 112 125 9 Underground Line Configuration (Config.)
112 113 525 9
113 114 325 9 Config. Phasing Cable Spacing ID
135 35 375 4 12 ABC 1/0 AA, CN 515
149 1 400 1
152 52 400 1
160 67 350 6
Table B.5 – Shunt Capacitors data
197 101 250 3
Shunt Capacitors

Node Ph-A Ph-B Ph-C

Table B.2 – Three Phase Switches kVAr kVAr kVAr
83 200 200 200
Three Phase Switches
88 50
90 50
Node A Node B Normal
92 50
13 152 closed
Total 250 250 250
18 135 closed
60 160 closed
61 610 closed
97 197 closed
150 149 closed
250 251 open
450 451 open
54 94 open
151 300 open
300 350 open

91
 The spacing ID numbers and type for overhead lines are summarized in Table B.6,

and Figure B.2 shows the spacing distances between the phase conductors and the neutral

conductor for Spacing ID numbers used for the overhead lines.

Table B.6 Overhead Line Spacing

Figure B.2 Overhead Line Spacing

The spacing ID numbers and type for underground lines are summarized in Table

B.7, and Figure B.3 shows the spacing distances between cables for underground lines:

Table B.7 Underground Line Spacing

Figure B.3 Underground Line Spacing

92
 Table B.8 Load data
Spot Loads 58 Y-I 0 0 20 10 0 0
59 Y-PQ 0 0 20 10 0 0
Node Load Ph-1 Ph-1 Ph-2 Ph-2 Ph-3 Ph-3 60 Y-PQ 20 10 0 0 0 0
Model kW kVAr kW kVAr kW kVAr 62 Y-Z 0 0 0 0 40 20
1 Y-PQ 40 20 0 0 0 0 63 Y-PQ 40 20 0 0 0 0
2 Y-PQ 0 0 20 10 0 0 64 Y-I 0 0 75 35 0 0
4 Y-PQ 0 0 0 0 40 20 65 D-Z 35 25 35 25 70 50
5 Y-I 0 0 0 0 20 10 66 Y-PQ 0 0 0 0 75 35
6 Y-Z 0 0 0 0 40 20 68 Y-PQ 20 10 0 0 0 0
7 Y-PQ 20 10 0 0 0 0 69 Y-PQ 40 20 0 0 0 0
9 Y-PQ 40 20 0 0 0 0 70 Y-PQ 20 10 0 0 0 0
10 Y-I 20 10 0 0 0 0 71 Y-PQ 40 20 0 0 0 0
11 Y-Z 40 20 0 0 0 0 73 Y-PQ 0 0 0 0 40 20
12 Y-PQ 0 0 20 10 0 0 74 Y-Z 0 0 0 0 40 20
16 Y-PQ 0 0 0 0 40 20 75 Y-PQ 0 0 0 0 40 20
17 Y-PQ 0 0 0 0 20 10 76 D-I 105 80 70 50 70 50
19 Y-PQ 40 20 0 0 0 0 77 Y-PQ 0 0 40 20 0 0
20 Y-I 40 20 0 0 0 0 79 Y-Z 40 20 0 0 0 0
22 Y-Z 0 0 40 20 0 0 80 Y-PQ 0 0 40 20 0 0
24 Y-PQ 0 0 0 0 40 20 82 Y-PQ 40 20 0 0 0 0
28 Y-I 40 20 0 0 0 0 83 Y-PQ 0 0 0 0 20 10
29 Y-Z 40 20 0 0 0 0 84 Y-PQ 0 0 0 0 20 10
30 Y-PQ 0 0 0 0 40 20 85 Y-PQ 0 0 0 0 40 20
31 Y-PQ 0 0 0 0 20 10 86 Y-PQ 0 0 20 10 0 0
32 Y-PQ 0 0 0 0 20 10 87 Y-PQ 0 0 40 20 0 0
33 Y-I 40 20 0 0 0 0 88 Y-PQ 40 20 0 0 0 0
34 Y-Z 0 0 0 0 40 20 90 Y-I 0 0 40 20 0 0
35 D-PQ 40 20 0 0 0 0 92 Y-PQ 0 0 0 0 40 20
37 Y-Z 40 20 0 0 0 0 94 Y-PQ 40 20 0 0 0 0
38 Y-I 0 0 20 10 0 0 95 Y-PQ 0 0 20 10 0 0
39 Y-PQ 0 0 20 10 0 0 96 Y-PQ 0 0 20 10 0 0
41 Y-PQ 0 0 0 0 20 10 98 Y-PQ 40 20 0 0 0 0
42 Y-PQ 20 10 0 0 0 0 99 Y-PQ 0 0 40 20 0 0
43 Y-Z 0 0 40 20 0 0 100 Y-Z 0 0 0 0 40 20
45 Y-I 20 10 0 0 0 0 102 Y-PQ 0 0 0 0 20 10
46 Y-PQ 20 10 0 0 0 0 103 Y-PQ 0 0 0 0 40 20
47 Y-I 35 25 35 25 35 25 104 Y-PQ 0 0 0 0 40 20
48 Y-Z 70 50 70 50 70 50 106 Y-PQ 0 0 40 20 0 0
49 Y-PQ 35 25 70 50 35 20 107 Y-PQ 0 0 40 20 0 0
50 Y-PQ 0 0 0 0 40 20 109 Y-PQ 40 20 0 0 0 0
51 Y-PQ 20 10 0 0 0 0 111 Y-PQ 20 10 0 0 0 0
52 Y-PQ 40 20 0 0 0 0 112 Y-I 20 10 0 0 0 0
53 Y-PQ 40 20 0 0 0 0 113 Y-Z 40 20 0 0 0 0
55 Y-Z 20 10 0 0 0 0 114 Y-PQ 20 10 0 0 0 0
56 Y-PQ 0 0 20 10 0 0 Total 1420 775 915 515 1155 635

93
 Figure B.4 shows the 123 node test system built in WinIGS-F. Besides the

parameters listed above, for adaptive relay programming, there are additional equipments

in the system such as relays and fuses. Detailed information about those equipments are

shown in Table B.9.

Figure B.4 IEEE 123-Node Test Feeder in WinIGS - F

94
 Table B.9 Switch type and location in the test system

Node A Node B Type Norminal Fuse

current Rating
1 S1 Fuse 14.1 25
8 S8 Fuse 13.5 25
13 S13 Fuse 13.9 25
18 S18 Fuse 10.7 25
18 135 Recloser 30.3
21 S21 Fuse 5.7 12
23 S23 Fuse 5.5 12
25 S25 Fuse 5.5 12
35 S35 Fuse 5.6 12
40 S40 Fuse 2.8 6
42 S42 Fuse 5.7 12
44 S44 Fuse 5.4 12
57 S57 Fuse 5.6 12
60 160 Recloser 65
67 S67 Fuse 15.7 25
72 S72 Fuse 16.3 25
76 S76 Fuse 11 25
87 S87 Fuse 5.2 12
89 S89 Fuse 5.5 12
91 S91 Fuse 5.5 12
93 S93 Fuse 5.2 12
95 S95 Fuse 2.8 6
97 S97 Fuse 5.5 12
101 S101 Fuse 14 25
105 S105 Fuse 11.1 25
108 S108 Fuse 18.1 25
149 150 Breaker 160

In order to calculate SAIDI and MAIFI scores as accurately as possible, each load in

the system needed to be assigned a set of customers. Looking at the customer data for the

U.S. given in the EIA’s (8), the following observations can be seen:

1) 37%, 36%, and 27% of the U.S. electricity load belongs to residential, commercial,

and industrial customers, respectively

95
2) 87%, 12%, and 1% of the U.S. electricity customers are residential, commercial, and

industrial customers, respectively

3) The average industrial load is about 120 times larger than the average residential load

4) The average commercial load is about 7 times larger than the average residential load

Additionally, the following assumptions were made regarding the IEEE 123 node

model:

1) Industrial customers can only be located at 3-phase lines, and are generally

segregated from regular residential/commercial sectors.

2) The average residential load is 1.67 kW and 0.83 kVar

Using these guidelines, a unique customer profile of 878 customers was created for

our 123 node model, with 764 of the 878 customers being residential, 109 being

commercial, and 5 being industrial. Detailed information about the load profile is shown

in Table B.10.

Table B.10 Load data in the test system

Node Residential Commercial Industrial Number of

load load load customer
1 12 1.71 0 13.71
2 6 0.86 0 6.86
4 12 1.71 0 13.71
5 6 0.86 0 6.86
6 12 1.71 0 13.71
7 6 0.86 0 6.86
9 12 1.71 0 13.71
10 6 0.86 0 6.86
11 12 1.71 0 13.71
12 6 0.86 0 6.86
16 12 1.71 0 13.71
17 6 0.86 0 6.86
19 12 1.71 0 13.71

96
 Table B.10 Continued

20 12 1.71 0 13.71
22 12 1.71 0 13.71
24 12 1.71 0 13.71
28 12 1.71 0 13.71
29 12 1.71 0 13.71
30 12 1.71 0 13.71
31 6 0.86 0 6.86
32 6 0.86 0 6.86
33 12 1.71 0 13.71
34 12 1.71 0 13.71
35 12 1.71 0 13.71
37 12 1.71 0 13.71
38 6 0.86 0 6.86
39 6 0.86 0 6.86
41 6 0.86 0 6.86
42 6 0.86 0 6.86
43 12 1.71 0 13.71
45 6 0.86 0 6.86
46 6 0.86 0 6.86
47 1.48 0.21 0.5 2.19
48 2.95 0.42 1 4.38
49 1.97 0.28 0.67 2.92
50 12 1.71 0 13.71
51 6 0.86 0 6.86
52 12 1.71 0 13.71
53 12 1.71 0 13.71
55 6 0.86 0 6.86
56 6 0.86 0 6.86
58 6 0.86 0 6.86
59 6 0.86 0 6.86
60 6 0.86 0 6.86
62 12 1.71 0 13.71
63 12 1.71 0 13.71
64 1.05 0.15 0.36 1.56
65 1.97 0.28 0.67 2.92
66 1.05 0.15 0.36 1.56
68 6 0.86 0 6.86
69 12 1.71 0 13.71
70 6 0.86 0 6.86
71 12 1.71 0 13.71
73 12 1.71 0 13.71
74 12 1.71 0 13.71
75 12 1.71 0 13.71

97
 Table B.10 Continued

76 3.45 0.49 1.17 5.11

77 12 1.71 0 13.71
79 12 1.71 0 13.71
80 12 1.71 0 13.71
82 12 1.71 0 13.71
83 6 0.86 0 6.86
84 6 0.86 0 6.86
85 12 1.71 0 13.71
86 6 0.86 0 6.86
87 12 1.71 0 13.71
88 12 1.71 0 13.71
90 12 1.71 0 13.71
92 12 1.71 0 13.71
94 12 1.71 0 13.71
95 6 0.86 0 6.86
96 6 0.86 0 6.86
98 12 1.71 0 13.71
99 12 1.71 0 13.71
100 12 1.71 0 13.71
102 6 0.86 0 6.86
103 12 1.71 0 13.71
104 12 1.71 0 13.71
106 12 1.71 0 13.71
107 12 1.71 0 13.71
109 12 1.71 0 13.71
111 6 0.86 0 6.86
112 6 0.86 0 6.86
113 12 1.71 0 13.71
114 6 0.86 0 6.86

98
 REFERENCES

[1] Richard E. Brown, “Impact of Smart Grid on distribution system design,” Power and Energy
Society General Meeting, pp. 1-4, July 2008

[2] A. J. Urdaneta, R. Nadira, and L. G. Perez, “Optimal Coordination of Directional

Overcurrent Relays in Interconnected power Systems,”IEEE Transactions on Power
Delivery, vol. 3, pp. 903-911, July 1988.

[3] A.P. Sakis Meliopoulos and George J. Cokkinides, Power System Relaying, Theory and
Applications,1996

[4] B. Chattopadhyay, M. S. Sachdev, and T. S. Sidhu, “An on-line relay coordination algorithm
for adaptive protection using linear programming technique,” IEEE Trans. Power Delivery,
vol. 11, pp. 165–173, Jan. 1996.

[5] A. J. Urdaneta, H. Restrepo, S. Marquez, and J. Sanchez, “Coordination of Directional

Overcurrent Relay Timing Using Linear Programming technique,” IEEE Transactions on
Power Delivery, vol. 11, pp. 122-129, Jan. 1996.

[6] Juan M. Gers and Edward J. Holmes, Protection of Electricity Distribution Networks, The
Institution of Electrical Engineers, London, United Kingdom, 2004, pp.70-73

[7] Perez, L. G., Urdaneta, A.J., “Optimal Computation of Distance Relays Second Zone Timing
in a Mixed Protection Scheme with Directional Overcurrent Relays,” IEEE Trans. On Power
Delivery, vol. 16, NO. 3, pp. 385–388, July 2001

[8] Vanderplaats, Garret N. Numerical Optimization Techniques for Engineering Design.

Colorado Springs : VR&D, 2005..

[9] Venter, Gerhard and Sobieski, Jaroslaw., Particle Swarm Optimization, 2003, AIAA, pp.
1583-1589.

[10] R. C. and Shi, Y., Comparing Inertia Weights and Constriction Factors in Particle Swarm
Optimization. Eberhart, IEEE. pp. 84-88, 2000.

[11] Clerc, M., The Swarm and the Queen: Towards a Deterministic and Adaptive Particle
Swarm Optimization. Washington, D.C., Proc. 1999 ICEC. pp. 1951-1957.

[12] Zhang, Li-ping, Yu, Huan-jun and Hu, Shang-xu., Optimal Choice of Parameters for Particle
Swarm Optimization, Journal of Zhejiang University Science, pp. 528-534, 2005

[13] Mithulananthan, N., Oo, Than and Phu, Le Van., Distributed Generator Placement in Power
Distribution System Using Genetic Algorithm to Reduce Losses., Thammasat Int. J. Sc.
Tech., Vol. 9, No. 3. pp. 55-62, 2004.

99
[14] Grefenstette, John J., Optimization of Control Parameters for Genetic Algorithms, IEEE
Trans. on Systems, Man, and Cybernetics, Vol. SMC-16, No. 1. pp. 122-128,1986.

100