Uv CBPSD lab
YuvrajSingh292
813 views
51 slides
Mar 09, 2018
Slide 1 of 51
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
About This Presentation
CBPSD Lab (electrical engineering 8th sem, RTU)
Slide Content
GLOBAL INSTITUTE OF TECHNOLOGY
JAIPUR
DEPARTMENT OF ELECTRICAL ENGINEERING
COMPUTER BASED POWER SYSTEM
LAB (8EE5A)
Submitted to: Submitted By:
Mr. SHIVAM SHARMA Name-Yuvraj Singh
Assistant professor (EE Dept.) Roll No.- 14EGCEE037
JAIPUR
DEPARTMENT OF ELECTRICAL ENGINEERING
COMPUTER BASED POWER SYSTEM
LAB (8EE5A)
Submitted to: Submitted By:
Mr. SHIVAM SHARMA Name-Yuvraj Singh
Assistant professor (EE Dept.) Roll No.- 14EGCEE037
14EGCEE037 Page 1
EXPERIMENT NO. 1
AIM : Fault Analysis for a 3bus system with verification of results using MATLAB for the
following cases –
1. LG Fault 2. LLG Fault
3. LL Fault 4. 3-phase Fault
Overview of MATLAB Simulink Environment
SIMULATION : Simulation is the virtual realization of some real thing, state of affairs, or process.
Simulation generally entails representing certain key characteristics or behavior of a selected physical or
abstract system.
Simulation is used in many contexts, including the modeling of natural systems or human
systems in order to gain insight into their functioning. Other contexts include simulation of
technology for performance optimization, safety engineering, testing, training and education.
SIMULINK: It is the library browser which contains all the block sets. It models, simulates, and
analyzes dynamic systems. It enables you to pose a question about a system, model the system,
and see what happens.
Simulink supports linear and nonlinear systems, modeled in continuous time, sampled time, or a
hybrid of the two. Systems can also be multirate — having different parts that are sampled or
updated at different rates. Sim Power Systems and Sim Mechanics of the Physical Modeling
product family work together with Simulink to model electrical, mechanical, and control
systems etc.
SIMULATION MODEL FOR FAULTS: The circuit shown below is designed for the
simulation of various types of faults in the transmission lines in MATLAB Simulink model file.
The respective blocks are copied in the model file from the Simulink Library browser.
EXPERIMENT NO. 1
AIM : Fault Analysis for a 3bus system with verification of results using MATLAB for the
following cases –
1. LG Fault 2. LLG Fault
3. LL Fault 4. 3-phase Fault
Overview of MATLAB Simulink Environment
SIMULATION : Simulation is the virtual realization of some real thing, state of affairs, or process.
Simulation generally entails representing certain key characteristics or behavior of a selected physical or
abstract system.
Simulation is used in many contexts, including the modeling of natural systems or human
systems in order to gain insight into their functioning. Other contexts include simulation of
technology for performance optimization, safety engineering, testing, training and education.
SIMULINK: It is the library browser which contains all the block sets. It models, simulates, and
analyzes dynamic systems. It enables you to pose a question about a system, model the system,
and see what happens.
Simulink supports linear and nonlinear systems, modeled in continuous time, sampled time, or a
hybrid of the two. Systems can also be multirate — having different parts that are sampled or
updated at different rates. Sim Power Systems and Sim Mechanics of the Physical Modeling
product family work together with Simulink to model electrical, mechanical, and control
systems etc.
SIMULATION MODEL FOR FAULTS: The circuit shown below is designed for the
simulation of various types of faults in the transmission lines in MATLAB Simulink model file.
The respective blocks are copied in the model file from the Simulink Library browser.
14EGCEE037 Page 2
SIMULATION GRAPHS :
The graphs obtained after the simulation of the various faults are shown for the following faults :
1) THREE PHASE TO GROUND FAULT
2) SINGLE LINE TO GROUND FAULT
3) DOUBLE LINE TO GROUND FAULT
4) LINE TO LINE FAULT
The graphs are shown respectively as follows :
SIMULATION GRAPHS :
The graphs obtained after the simulation of the various faults are shown for the following faults :
1) THREE PHASE TO GROUND FAULT
2) SINGLE LINE TO GROUND FAULT
3) DOUBLE LINE TO GROUND FAULT
4) LINE TO LINE FAULT
The graphs are shown respectively as follows :
14EGCEE037 Page 3
DESCRIPTION OF VARIOUS BLOCKS
Simplified Synchronous Machine
Models the dynamics of a simplified three-phase synchronous machine.
Description : The Simplified Synchronous Machine block models both the electrical and mechanical
characteristics of a simple synchronous machine. The electrical system for each phase consists of a
voltage source in series with an RL impedance, which implements the internal impedance of the
machine. The value of R can be zero but the value of L must be positive.
Three-Phase Series RLC Load
Implements a three-phase series RLC load with selectable connection.
DESCRIPTION OF VARIOUS BLOCKS
Simplified Synchronous Machine
Models the dynamics of a simplified three-phase synchronous machine.
Description : The Simplified Synchronous Machine block models both the electrical and mechanical
characteristics of a simple synchronous machine. The electrical system for each phase consists of a
voltage source in series with an RL impedance, which implements the internal impedance of the
machine. The value of R can be zero but the value of L must be positive.
Three-Phase Series RLC Load
Implements a three-phase series RLC load with selectable connection.
14EGCEE037 Page 4
Description : The Three-Phase Series RLC Load block implements a three-phase balanced load as a
series combination of RLC elements. At the specified frequency, the load exhibits a constant
impedance. The active and reactive powers absorbed by the load are proportional to the square of the
applied voltage.
Three-Phase Transformer (Two Windings)
Implements a three-phase transformer with configurable winding connections.
Description : The three-Phase Transformer (2 Windings) block implements a 3-phase Transformer using
3 single-phase transformers. You can simulate the saturable core or not simply by setting the appropriate
check box in the parameter menu of the block.
Three-Phase Breaker
Implements a three-phase circuit breaker opening at the current zero crossing
Description : Three-Phase Breaker block implements a three-phase circuit breaker where
the opening and closing times can be controlled either from an external Simulink signal
(external control mode), or from an internal control timer (internal control mode). The
Three-Phase Breaker block uses three Breaker blocks connected between the inputs
and the outputs of the block. If the Three-Phase Breaker block is set in external control mode,
a control input appears in the block icon.
Distributed Parameter Line
Implementing an N-phase distributed parameter transmission line model with lumped losses.
Description:- The Distributed Parameter Line block implements an N-phase distributed parameter
line model with lumped losses. The model is based on the Bergeron's traveling wave method used
by the Electromagnetic Transient Program (EMTP) . In this model, the lossless distributed LC line
is characterized by two values (for a single-phase line): the surge impedance Zc = (L/C) and the
phase velocity v= 1/√(LC).
The model uses the fact that the quantity e+Zi (where e is line voltage and i is line current) entering
one end of the line must arrive unchanged at the other end after a transport delay of τ= d/v.
Description : The Three-Phase Series RLC Load block implements a three-phase balanced load as a
series combination of RLC elements. At the specified frequency, the load exhibits a constant
impedance. The active and reactive powers absorbed by the load are proportional to the square of the
applied voltage.
Three-Phase Transformer (Two Windings)
Implements a three-phase transformer with configurable winding connections.
Description : The three-Phase Transformer (2 Windings) block implements a 3-phase Transformer using
3 single-phase transformers. You can simulate the saturable core or not simply by setting the appropriate
check box in the parameter menu of the block.
Three-Phase Breaker
Implements a three-phase circuit breaker opening at the current zero crossing
Description : Three-Phase Breaker block implements a three-phase circuit breaker where
the opening and closing times can be controlled either from an external Simulink signal
(external control mode), or from an internal control timer (internal control mode). The
Three-Phase Breaker block uses three Breaker blocks connected between the inputs
and the outputs of the block. If the Three-Phase Breaker block is set in external control mode,
a control input appears in the block icon.
Distributed Parameter Line
Implementing an N-phase distributed parameter transmission line model with lumped losses.
Description:- The Distributed Parameter Line block implements an N-phase distributed parameter
line model with lumped losses. The model is based on the Bergeron's traveling wave method used
by the Electromagnetic Transient Program (EMTP) . In this model, the lossless distributed LC line
is characterized by two values (for a single-phase line): the surge impedance Zc = (L/C) and the
phase velocity v= 1/√(LC).
The model uses the fact that the quantity e+Zi (where e is line voltage and i is line current) entering
one end of the line must arrive unchanged at the other end after a transport delay of τ= d/v.
14EGCEE037 Page 5
Three-Phase V-I Measurement
Measures three-phase currents and voltages in a circuit.
Description : The Three-PhaseV-I Measurement block is used to measure three-phase voltages and
currents in a circuit. When connected in series with three-phase elements, it returns the three phase-
to-ground or phase-to-phase voltages and the three line currents. The block can output the voltages
and currents in per unit (p.u.) values or in volts and amperes.
Three-Phase Sequence Analyzer
Measures the positive-, negative-, and zero-sequence components of a three-phase signal
Description : The Three-Phase Sequence Analyzer block outputs the magnitude and phase of the
positive- (denoted by the index 1), negative- (index 2), and zero-sequence (index 0) components
of a set of three balanced or unbalanced signals. The signals can contain harmonics or not.
Scope
Displays signals generated during a simulation.
Description : The Scope block displays its input with respect to simulation time. The Scope block can
have multiple axes (one per port); all axes have a common time range with independent y-axes. The
Scope allows. The Scope provides toolbar buttons that enables to zoom in on displayed data, display all
the data input to the Scope, preserve axis settings from one simulation to the next, limit data displayed,
and save data to the workspace.
Three-Phase Fault
Implementing a programmable phase-to-phase and phase-to-ground fault breaker system.
Description : The Three-Phase Fault block implements a three-phase circuit breaker where the opening
and closing times can be controlled either from an external Simulink signal (external control mode),
or from an internal control timer (internal control mode). The Three-Phase Fault block uses three Breaker
blocks that can be individually switched on and off to program phase-to-phase faults, phase-to-ground
faults, or a combination of phase-to-phase and ground faults.
Three-Phase V-I Measurement
Measures three-phase currents and voltages in a circuit.
Description : The Three-PhaseV-I Measurement block is used to measure three-phase voltages and
currents in a circuit. When connected in series with three-phase elements, it returns the three phase-
to-ground or phase-to-phase voltages and the three line currents. The block can output the voltages
and currents in per unit (p.u.) values or in volts and amperes.
Three-Phase Sequence Analyzer
Measures the positive-, negative-, and zero-sequence components of a three-phase signal
Description : The Three-Phase Sequence Analyzer block outputs the magnitude and phase of the
positive- (denoted by the index 1), negative- (index 2), and zero-sequence (index 0) components
of a set of three balanced or unbalanced signals. The signals can contain harmonics or not.
Scope
Displays signals generated during a simulation.
Description : The Scope block displays its input with respect to simulation time. The Scope block can
have multiple axes (one per port); all axes have a common time range with independent y-axes. The
Scope allows. The Scope provides toolbar buttons that enables to zoom in on displayed data, display all
the data input to the Scope, preserve axis settings from one simulation to the next, limit data displayed,
and save data to the workspace.
Three-Phase Fault
Implementing a programmable phase-to-phase and phase-to-ground fault breaker system.
Description : The Three-Phase Fault block implements a three-phase circuit breaker where the opening
and closing times can be controlled either from an external Simulink signal (external control mode),
or from an internal control timer (internal control mode). The Three-Phase Fault block uses three Breaker
blocks that can be individually switched on and off to program phase-to-phase faults, phase-to-ground
faults, or a combination of phase-to-phase and ground faults.
14EGCEE037 Page 6
CONCLUSION
After the MATLAB simulation for faults, it was observed that the voltage and current
waveforms were transient in nature in the initial period after the occurrence of faults.
During the initial part of short circuit, the short circuit current was limited by subtransient
reactance of synchronous machine and impedance of transmission line between the machine
and point of fault. After that, it was limited by transient reactance of synchronous machine and
impedance of line.
Finally, the short circuit current settled down to steady state short circuit value limited by
synchronous reactance of the machine and line impedance. The negative and zero sequence
components were present initially only and they disappeared after the circuit breaker cleared
the fault.
CONCLUSION
After the MATLAB simulation for faults, it was observed that the voltage and current
waveforms were transient in nature in the initial period after the occurrence of faults.
During the initial part of short circuit, the short circuit current was limited by subtransient
reactance of synchronous machine and impedance of transmission line between the machine
and point of fault. After that, it was limited by transient reactance of synchronous machine and
impedance of line.
Finally, the short circuit current settled down to steady state short circuit value limited by
synchronous reactance of the machine and line impedance. The negative and zero sequence
components were present initially only and they disappeared after the circuit breaker cleared
the fault.
14EGCEE037 Page 7
EXPERIMENT NO. 2
AIM : Load Flow Analysis for a 3 to 6 bus system, in MATLAB, using the following methods :
(i) Gauss Seidal Method (ii.) Newton Raphson
(iii.) Fast Decoupled Method.
THEORY
The Newton Raphson method of load flow analysis is an iterative method which approximates
the set of non-linear simultaneous equations to a set of linear simultaneous equations using
Taylor‘s series expansion and the terms are limited to first order approximation. The load flow
equations for Newton Raphson method are non-linear equations in terms of real and imaginary
part of bus voltages.
Where, ep = Real part of Vp
fp = Imaginary part of Vp
Gpq, Bpq = Conductance and Susceptance of admittance Ypq respectively.
Bus system find the load flow using Newton raphson method
DATA
basemva = 100; accuracy = 0.001; accel = 1.8; maxiter = 100
% IEEE 30-BUS TEST SYSTEM (American Electric Power)
% Bus Bus Voltage Angle ---Load---- -------Generator----- Static Mvar
% No code Mag. Degree MW Mvar MW Mvar Qmin Qmax +Qc/-Ql
BUS DATA
1. 1 1.06 0 0 0 0 0 0 0 0
2. 2 1.04 0 21.7 12.7 40 0 -40 50 0
3. 0 1.0 0 2.4 1.2 0 0 0 0 0
4. 0 1.06 0 7.6 1.6 0 0 0 0 0
5. 2 1.01 0 94.2 19 0 0 -40 40 0
6. 0 1 0 0 0 0 0 0 0 0
7. 0 1 0 22.8 10.9 0 0 0 0 0
8. 2 1.01 0 30 30 0 0 -30 40 0
9. 0 1 0 0 0 0 0 0 0 0
EXPERIMENT NO. 2
AIM : Load Flow Analysis for a 3 to 6 bus system, in MATLAB, using the following methods :
(i) Gauss Seidal Method (ii.) Newton Raphson
(iii.) Fast Decoupled Method.
THEORY
The Newton Raphson method of load flow analysis is an iterative method which approximates
the set of non-linear simultaneous equations to a set of linear simultaneous equations using
Taylor‘s series expansion and the terms are limited to first order approximation. The load flow
equations for Newton Raphson method are non-linear equations in terms of real and imaginary
part of bus voltages.
Where, ep = Real part of Vp
fp = Imaginary part of Vp
Gpq, Bpq = Conductance and Susceptance of admittance Ypq respectively.
Bus system find the load flow using Newton raphson method
DATA
basemva = 100; accuracy = 0.001; accel = 1.8; maxiter = 100
% IEEE 30-BUS TEST SYSTEM (American Electric Power)
% Bus Bus Voltage Angle ---Load---- -------Generator----- Static Mvar
% No code Mag. Degree MW Mvar MW Mvar Qmin Qmax +Qc/-Ql
BUS DATA
1. 1 1.06 0 0 0 0 0 0 0 0
2. 2 1.04 0 21.7 12.7 40 0 -40 50 0
3. 0 1.0 0 2.4 1.2 0 0 0 0 0
4. 0 1.06 0 7.6 1.6 0 0 0 0 0
5. 2 1.01 0 94.2 19 0 0 -40 40 0
6. 0 1 0 0 0 0 0 0 0 0
7. 0 1 0 22.8 10.9 0 0 0 0 0
8. 2 1.01 0 30 30 0 0 -30 40 0
9. 0 1 0 0 0 0 0 0 0 0
14EGCEE037 Page 8
10. 0 1 0 5.8 2 0 0 -6 24 19
11. 2 1.08 0 0 0 0 0 0 0 0
12. 0 1 0 11.2 7.5 0 0 0 0 0
13. 2 1.07 0 0 0 0 0 -6 24 0
14. 0 1 0 6.2 1.6 0 0 0 0 0
15. 0 1 0 8.2 2.5 0 0 0 0 0
16. 0 1 0 3.5 1.8 0 0 0 0 0
17. 0 1 0 95.8 0 0 0 0 0 0
18. 0 1 0 3.2 0.9 0 0 0 0 0
19. 0 1 0 9.5 3.4 0 0 0 0 0
20. 0 1 0 2.2 0.7 0 0 0 0 0
DATA
% Line code
% Bus bus R X 1/2 B = 1 for lines
% nl nr p.u. p.u. p.u. > 1 or < 1 tr. tap at bus nl
LINE DATA
1. 1 2 0.0192 0.0575 0.0264 1
2. 1 3 0.0452 0.1852 0.0204 1
3. 2 4 0.057 0.1737 0.0184 1
4. 3 4 0.0132 0.0379 0.0042 1
5. 2 5 0.0472 0.1983 0.0209 1
6. 2 6 0.058 0.1763 0.0187 1
7. 4 6 0.0119 0.0414 0.0045 1
8. 5 7 0.046 0.116 0.0102 1
9. 6 7 0.0267 0.082 0.00850 1
10. 6 8 0.012 0.042 0.00450 1
11. 6 9 0 0.208 0.978 1
12. 6 10 0 0.556 0 0.969
13. 9 11 0 0.2080 0 1
14. 9 10 0 0.11 0 1
15. 4 12 0 0.256 0 0.932
16. 12 13 0 0.14 0 1
17. 12 14 0.123 0.2559 0 1
18. 12 15 0.066 0.1304 0 1
19. 12 16 0.0945 0.1987 0 1
20. 14 15 0.2210 0.1997 0 1
21. 16 17 0.0824 0.1923 0 1
22. 15 18 0.1073 0.2185 0 1
23. 18 19 0.063 0.1292 0 1
24. 19 20 0.034 0.068 0 1
10. 0 1 0 5.8 2 0 0 -6 24 19
11. 2 1.08 0 0 0 0 0 0 0 0
12. 0 1 0 11.2 7.5 0 0 0 0 0
13. 2 1.07 0 0 0 0 0 -6 24 0
14. 0 1 0 6.2 1.6 0 0 0 0 0
15. 0 1 0 8.2 2.5 0 0 0 0 0
16. 0 1 0 3.5 1.8 0 0 0 0 0
17. 0 1 0 95.8 0 0 0 0 0 0
18. 0 1 0 3.2 0.9 0 0 0 0 0
19. 0 1 0 9.5 3.4 0 0 0 0 0
20. 0 1 0 2.2 0.7 0 0 0 0 0
DATA
% Line code
% Bus bus R X 1/2 B = 1 for lines
% nl nr p.u. p.u. p.u. > 1 or < 1 tr. tap at bus nl
LINE DATA
1. 1 2 0.0192 0.0575 0.0264 1
2. 1 3 0.0452 0.1852 0.0204 1
3. 2 4 0.057 0.1737 0.0184 1
4. 3 4 0.0132 0.0379 0.0042 1
5. 2 5 0.0472 0.1983 0.0209 1
6. 2 6 0.058 0.1763 0.0187 1
7. 4 6 0.0119 0.0414 0.0045 1
8. 5 7 0.046 0.116 0.0102 1
9. 6 7 0.0267 0.082 0.00850 1
10. 6 8 0.012 0.042 0.00450 1
11. 6 9 0 0.208 0.978 1
12. 6 10 0 0.556 0 0.969
13. 9 11 0 0.2080 0 1
14. 9 10 0 0.11 0 1
15. 4 12 0 0.256 0 0.932
16. 12 13 0 0.14 0 1
17. 12 14 0.123 0.2559 0 1
18. 12 15 0.066 0.1304 0 1
19. 12 16 0.0945 0.1987 0 1
20. 14 15 0.2210 0.1997 0 1
21. 16 17 0.0824 0.1923 0 1
22. 15 18 0.1073 0.2185 0 1
23. 18 19 0.063 0.1292 0 1
24. 19 20 0.034 0.068 0 1
14EGCEE037 Page 9
PROGRAM
%NEWTON RAPHSON METHOD
clc;
gbus = [1 2.0 1.0 0.0 0.0
2 0.0 0.0 0.5 1.0
3 1.5 0.6 0.0 0.0];
ybus = [5.882-j*23.528 -2.941+j*11.764 -2.941+j*11.764
-2.941+j*11.764 5.882-j*23.528 -2.941+j*11.764
-2.941+j*11.764 -2.941+j*11.764 5.882-j*23.528];
t= 0.001
v1=1.04+j*0;
v2=1+j*0;
v3=1.04+j*0;
del3=angle(v3);
del1=angle(v1);
del2=angle(v2);
%abs(ybus(2,1))
%abs(v2)
for i=1:10
p2=(abs(v2)*abs(v1)*abs(ybus(2,1))*cos((angle(ybus(2,1)))+del1-del2))
+abs(v2)*
abs(v2)*abs(ybus(2,2))*cos((angle(ybus(2,2))))+(abs(v2)*abs(v3)*
abs(ybus(2,3))*cos((angle(ybus(2,3))+del3-del2));
q2=-(abs(v2)*abs(v1)*abs(ybus(2,1))*sin((angle(ybus(2,1)))+del1-del2))-
abs(v2)*abs(v2)*abs(ybus(2,2))*sin((angle(ybus(2,2))))-(abs(v2)*abs(v3)*
abs(ybus(2,3))*sin((angle(ybus(2,3)))+del3-del2));
p3=(abs(v3)*abs(v1)*abs(ybus(3,1))*cos((angle(ybus(3,1)))+del1-
del3))+abs(v3)*abs(v3)*abs(ybus(3,3))*cos((angle(ybus(3,3))))+
(abs(v2)*abs(v3)*
abs(ybus(3,2))*cos((angle(ybus(3,2)))+del2-del3));
delp20=gbus(2,4)-gbus(2,2)-p2;
delp30=gbus(3,4)-gbus(3,2)-p3;
delq20=gbus(2,5)-gbus(2,3)-q2;
J(1,1)=(abs(v2)*abs(v1)*abs(ybus(2,1))*sin((angle(ybus(2,1)))+del1-
del2))+(abs(v2)*abs(v3)*abs(ybus(2,3))*sin((angle(ybus(2,3)))+del3-
del2));
J(1,2)=-(abs(v2)*abs(v3)*abs(ybus(2,3))*sin((angle(ybus(2,3)))+del3-
del2));
PROGRAM
%NEWTON RAPHSON METHOD
clc;
gbus = [1 2.0 1.0 0.0 0.0
2 0.0 0.0 0.5 1.0
3 1.5 0.6 0.0 0.0];
ybus = [5.882-j*23.528 -2.941+j*11.764 -2.941+j*11.764
-2.941+j*11.764 5.882-j*23.528 -2.941+j*11.764
-2.941+j*11.764 -2.941+j*11.764 5.882-j*23.528];
t= 0.001
v1=1.04+j*0;
v2=1+j*0;
v3=1.04+j*0;
del3=angle(v3);
del1=angle(v1);
del2=angle(v2);
%abs(ybus(2,1))
%abs(v2)
for i=1:10
p2=(abs(v2)*abs(v1)*abs(ybus(2,1))*cos((angle(ybus(2,1)))+del1-del2))
+abs(v2)*
abs(v2)*abs(ybus(2,2))*cos((angle(ybus(2,2))))+(abs(v2)*abs(v3)*
abs(ybus(2,3))*cos((angle(ybus(2,3))+del3-del2));
q2=-(abs(v2)*abs(v1)*abs(ybus(2,1))*sin((angle(ybus(2,1)))+del1-del2))-
abs(v2)*abs(v2)*abs(ybus(2,2))*sin((angle(ybus(2,2))))-(abs(v2)*abs(v3)*
abs(ybus(2,3))*sin((angle(ybus(2,3)))+del3-del2));
p3=(abs(v3)*abs(v1)*abs(ybus(3,1))*cos((angle(ybus(3,1)))+del1-
del3))+abs(v3)*abs(v3)*abs(ybus(3,3))*cos((angle(ybus(3,3))))+
(abs(v2)*abs(v3)*
abs(ybus(3,2))*cos((angle(ybus(3,2)))+del2-del3));
delp20=gbus(2,4)-gbus(2,2)-p2;
delp30=gbus(3,4)-gbus(3,2)-p3;
delq20=gbus(2,5)-gbus(2,3)-q2;
J(1,1)=(abs(v2)*abs(v1)*abs(ybus(2,1))*sin((angle(ybus(2,1)))+del1-
del2))+(abs(v2)*abs(v3)*abs(ybus(2,3))*sin((angle(ybus(2,3)))+del3-
del2));
J(1,2)=-(abs(v2)*abs(v3)*abs(ybus(2,3))*sin((angle(ybus(2,3)))+del3-
del2));
14EGCEE037 Page 10
J(1,3)=(abs(v1)*abs(ybus(2,1))*cos((angle(ybus(2,1)))+del1-
del2))+2*(abs(v2)*abs(ybus(2,2))*cos((angle(ybus(2,2))))+
(abs(v3)*abs(ybus(2,3))*
cos((angle(ybus(2,3)))+del3-del2));
J(2,1)=-(abs(v3)*abs(v2)*abs(ybus(3,2))*sin((angle(ybus(3,2)))+del2-
del3));
J(2,2)=(abs(v3)*abs(v1)*abs(ybus(3,1))*sin((angdel3))+(abs(v3)*abs(v2)*abs(ybus(3,2))*sin((a
ngle(ybus(3,2)))+del2-
del3));
J(2,3)=(abs(v3)*abs(ybus(3,2))*cos((angle(ybus(3,2)))+del2-del3));
J(3,1)=(abs(v2)*abs(v1)*abs(ybus(2,1))*cos((angle(ybus(2,1)))+del1-
del2))-(abs(v2)*abs(v3)*abs(ybus(2,3))*cos((angle(ybus(2,3)))+del2-
del3));
J(3,2)=(abs(v2)*abs(v3)*abs(ybus(2,3))*cos((angle(ybus(2,3)))+del2-
del3));
J(3,3)=-(abs(v2)*abs(ybus(2,1))*sin((angle(ybus(2,1)))+del1-
del2))-2*(abs(v2)*abs(ybus(2,2))*sin((angle(ybus(2,2))))-
(abs(v3)*abs(ybus(2,3))*
sin((angle(ybus(2,3)))+del3-del2));
end
J
inv(J);
A=[del2;del3;abs(v2)];
delA0=[delp20;delp30;delq20];
delA1=inv(J)*delA0;
delA1;
b0=abs(v2);
A1=[del2;del3;b0]+delA1;
A1-delA0;
if((A1-delA0)<=t)
break;
del2=A1(1,1);
del3=A1(2,1);
abs(v2)=A1(3,1);
end
A1
RESULT:- Line Flow and Losses
--Line-- Power at bus & line flow --Line loss-- Transformer
from to MW Mvar MVA MW Mvar tap
J(1,3)=(abs(v1)*abs(ybus(2,1))*cos((angle(ybus(2,1)))+del1-
del2))+2*(abs(v2)*abs(ybus(2,2))*cos((angle(ybus(2,2))))+
(abs(v3)*abs(ybus(2,3))*
cos((angle(ybus(2,3)))+del3-del2));
J(2,1)=-(abs(v3)*abs(v2)*abs(ybus(3,2))*sin((angle(ybus(3,2)))+del2-
del3));
J(2,2)=(abs(v3)*abs(v1)*abs(ybus(3,1))*sin((angdel3))+(abs(v3)*abs(v2)*abs(ybus(3,2))*sin((a
ngle(ybus(3,2)))+del2-
del3));
J(2,3)=(abs(v3)*abs(ybus(3,2))*cos((angle(ybus(3,2)))+del2-del3));
J(3,1)=(abs(v2)*abs(v1)*abs(ybus(2,1))*cos((angle(ybus(2,1)))+del1-
del2))-(abs(v2)*abs(v3)*abs(ybus(2,3))*cos((angle(ybus(2,3)))+del2-
del3));
J(3,2)=(abs(v2)*abs(v3)*abs(ybus(2,3))*cos((angle(ybus(2,3)))+del2-
del3));
J(3,3)=-(abs(v2)*abs(ybus(2,1))*sin((angle(ybus(2,1)))+del1-
del2))-2*(abs(v2)*abs(ybus(2,2))*sin((angle(ybus(2,2))))-
(abs(v3)*abs(ybus(2,3))*
sin((angle(ybus(2,3)))+del3-del2));
end
J
inv(J);
A=[del2;del3;abs(v2)];
delA0=[delp20;delp30;delq20];
delA1=inv(J)*delA0;
delA1;
b0=abs(v2);
A1=[del2;del3;b0]+delA1;
A1-delA0;
if((A1-delA0)<=t)
break;
del2=A1(1,1);
del3=A1(2,1);
abs(v2)=A1(3,1);
end
A1
RESULT:- Line Flow and Losses
--Line-- Power at bus & line flow --Line loss-- Transformer
from to MW Mvar MVA MW Mvar tap
14EGCEE037 Page 11
1 260.950 -17.010 261.504
2 177.743 -22.140 179.117 5.461 10.517
3 83.197 5.125 83.354 2.807 7.079
2 18.300 36.126 40.497
1 -172.282 32.657 175.350 5.461 10.517
4 45.702 2.720 45.783 1.106 -0.519
5 82.990 1.704 83.008 2.995 8.178
6 61.905 -0.966 61.913 2.047 2.263
3 -2.400 -1.200 2.683
1 -80.390 1.954 80.414 2.807 7.079
4 78.034 -3.087 78.095 0.771 1.345
4 -7.600 -1.600 7.767
2 -44.596 -3.239 44.713 1.106 -0.519
3 -77.263 4.432 77.390 0.771 1.345
6 70.132 -17.624 72.313 0.605 1.181
12 44.131 14.627 46.492 0.000 4.686 0.932
5 -94.200 16.995 95.721
2 -79.995 6.474 80.256 2.995 8.178
7 -14.210 10.467 17.649 0.151 -1.687
6 0.000 0.000 0.000
2 -59.858 3.229 59.945 2.047 2.263
4 -69.527 18.805 72.026 0.605 1.181
7 37.537 -1.915 37.586 0.368 -0.598
8 29.534 -3.712 29.766 0.103 -0.558
9 27.687 -7.318 28.638 0.000 1.593 0.978
10 15.828 0.656 15.842 -0.000 1.279 0.969
28 18.840 -9.575 21.134 0.060 -13.085
7 -22.800 -10.900 25.272
5 14.361 -12.154 18.814 0.151 -1.687
6 -37.170 1.317 37.193 0.368 -0.598
8 -30.000 0.759 30.010
6 -29.431 3.154 29.599 0.103 -0.558
28 -0.570 -2.366 2.433 0.000 -4.368
1 260.950 -17.010 261.504
2 177.743 -22.140 179.117 5.461 10.517
3 83.197 5.125 83.354 2.807 7.079
2 18.300 36.126 40.497
1 -172.282 32.657 175.350 5.461 10.517
4 45.702 2.720 45.783 1.106 -0.519
5 82.990 1.704 83.008 2.995 8.178
6 61.905 -0.966 61.913 2.047 2.263
3 -2.400 -1.200 2.683
1 -80.390 1.954 80.414 2.807 7.079
4 78.034 -3.087 78.095 0.771 1.345
4 -7.600 -1.600 7.767
2 -44.596 -3.239 44.713 1.106 -0.519
3 -77.263 4.432 77.390 0.771 1.345
6 70.132 -17.624 72.313 0.605 1.181
12 44.131 14.627 46.492 0.000 4.686 0.932
5 -94.200 16.995 95.721
2 -79.995 6.474 80.256 2.995 8.178
7 -14.210 10.467 17.649 0.151 -1.687
6 0.000 0.000 0.000
2 -59.858 3.229 59.945 2.047 2.263
4 -69.527 18.805 72.026 0.605 1.181
7 37.537 -1.915 37.586 0.368 -0.598
8 29.534 -3.712 29.766 0.103 -0.558
9 27.687 -7.318 28.638 0.000 1.593 0.978
10 15.828 0.656 15.842 -0.000 1.279 0.969
28 18.840 -9.575 21.134 0.060 -13.085
7 -22.800 -10.900 25.272
5 14.361 -12.154 18.814 0.151 -1.687
6 -37.170 1.317 37.193 0.368 -0.598
8 -30.000 0.759 30.010
6 -29.431 3.154 29.599 0.103 -0.558
28 -0.570 -2.366 2.433 0.000 -4.368
14EGCEE037 Page 12
9 0.000 0.000 0.000
6 -27.687 8.911 29.086 0.000 1.593
11 0.003 -15.653 15.653 0.000 0.461
10 27.731 6.747 28.540 -0.000 0.811
10 -5.800 17.000 17.962
6 -15.828 0.623 15.840 -0.000 1.279
9 -27.731 -5.936 28.359 -0.000 0.811
20 9.018 3.569 9.698 0.081 0.180
17 5.347 4.393 6.920 0.014 0.037
21 15.723 9.846 18.551 0.110 0.236
22 7.582 4.487 8.811 0.052 0.107
11 0.000 16.113 16.113
9 -0.003 16.114 16.114 0.000 0.461
12 -11.200 -7.500 13.479
4 -44.131 -9.941 45.237 0.000 4.686
13 -0.021 -10.274 10.274 0.000 0.132
14 7.852 2.428 8.219 0.074 0.155
15 17.852 6.968 19.164 0.217 0.428
16 7.206 3.370 7.955 0.053 0.112
13 0.000 10.406 10.406
12 0.021 10.406 10.406 0.000 0.132
14 -6.200 -1.600 6.403
12 -7.778 -2.273 8.103 0.074 0.155
15 1.592 0.708 1.742 0.006 0.006
15 -8.200 -2.500 8.573
12 -17.634 -6.540 18.808 0.217 0.428
14 -1.586 -0.702 1.734 0.006 0.006
18 6.009 1.741 6.256 0.039 0.079
23 5.004 2.963 5.815 0.031 0.063
16 -3.500 -1.800 3.936
12 -7.152 -3.257 7.859 0.053 0.112
17 3.658 1.440 3.931 0.012 0.027
9 0.000 0.000 0.000
6 -27.687 8.911 29.086 0.000 1.593
11 0.003 -15.653 15.653 0.000 0.461
10 27.731 6.747 28.540 -0.000 0.811
10 -5.800 17.000 17.962
6 -15.828 0.623 15.840 -0.000 1.279
9 -27.731 -5.936 28.359 -0.000 0.811
20 9.018 3.569 9.698 0.081 0.180
17 5.347 4.393 6.920 0.014 0.037
21 15.723 9.846 18.551 0.110 0.236
22 7.582 4.487 8.811 0.052 0.107
11 0.000 16.113 16.113
9 -0.003 16.114 16.114 0.000 0.461
12 -11.200 -7.500 13.479
4 -44.131 -9.941 45.237 0.000 4.686
13 -0.021 -10.274 10.274 0.000 0.132
14 7.852 2.428 8.219 0.074 0.155
15 17.852 6.968 19.164 0.217 0.428
16 7.206 3.370 7.955 0.053 0.112
13 0.000 10.406 10.406
12 0.021 10.406 10.406 0.000 0.132
14 -6.200 -1.600 6.403
12 -7.778 -2.273 8.103 0.074 0.155
15 1.592 0.708 1.742 0.006 0.006
15 -8.200 -2.500 8.573
12 -17.634 -6.540 18.808 0.217 0.428
14 -1.586 -0.702 1.734 0.006 0.006
18 6.009 1.741 6.256 0.039 0.079
23 5.004 2.963 5.815 0.031 0.063
16 -3.500 -1.800 3.936
12 -7.152 -3.257 7.859 0.053 0.112
17 3.658 1.440 3.931 0.012 0.027
14EGCEE037 Page 13
17 -9.000 -5.800 10.707
16 -3.646 -1.413 3.910 0.012 0.027
10 -5.332 -4.355 6.885 0.014 0.037
18 -3.200 -0.900 3.324
15 -5.970 -1.661 6.197 0.039 0.079
19 2.779 0.787 2.888 0.005 0.010
19 -9.500 -3.400 10.090
18 -2.774 -0.777 2.881 0.005 0.010
20 -6.703 -2.675 7.217 0.017 0.034
20 -2.200 -0.700 2.309
19 6.720 2.709 7.245 0.017 0.034
10 -8.937 -3.389 9.558 0.081 0.180
21 -17.500 -11.200 20.777
10 -15.613 -9.609 18.333 0.110 0.236
22 -1.849 -1.627 2.463 0.001 0.001
22 0.000 0.000 0.000
10 -7.531 -4.380 8.712 0.052 0.107
21 1.850 1.628 2.464 0.001 0.001
24 5.643 2.795 6.297 0.043 0.067
23 -3.200 -1.600 3.578
15 -4.972 -2.900 5.756 0.031 0.063
24 1.771 1.282 2.186 0.006 0.012
24 -8.700 -2.400 9.025
22 -5.601 -2.728 6.230 0.043 0.067
23 -1.765 -1.270 2.174 0.006 0.012
25 -1.322 1.604 2.079 0.008 0.014
25 0.000 0.000 0.000
24 1.330 -1.590 2.073 0.008 0.014
26 3.520 2.372 4.244 0.044 0.066
27 -4.866 -0.786 4.929 0.026 0.049
26 -3.500 -2.300 4.188
25 -3.476 -2.306 4.171 0.044 0.066
17 -9.000 -5.800 10.707
16 -3.646 -1.413 3.910 0.012 0.027
10 -5.332 -4.355 6.885 0.014 0.037
18 -3.200 -0.900 3.324
15 -5.970 -1.661 6.197 0.039 0.079
19 2.779 0.787 2.888 0.005 0.010
19 -9.500 -3.400 10.090
18 -2.774 -0.777 2.881 0.005 0.010
20 -6.703 -2.675 7.217 0.017 0.034
20 -2.200 -0.700 2.309
19 6.720 2.709 7.245 0.017 0.034
10 -8.937 -3.389 9.558 0.081 0.180
21 -17.500 -11.200 20.777
10 -15.613 -9.609 18.333 0.110 0.236
22 -1.849 -1.627 2.463 0.001 0.001
22 0.000 0.000 0.000
10 -7.531 -4.380 8.712 0.052 0.107
21 1.850 1.628 2.464 0.001 0.001
24 5.643 2.795 6.297 0.043 0.067
23 -3.200 -1.600 3.578
15 -4.972 -2.900 5.756 0.031 0.063
24 1.771 1.282 2.186 0.006 0.012
24 -8.700 -2.400 9.025
22 -5.601 -2.728 6.230 0.043 0.067
23 -1.765 -1.270 2.174 0.006 0.012
25 -1.322 1.604 2.079 0.008 0.014
25 0.000 0.000 0.000
24 1.330 -1.590 2.073 0.008 0.014
26 3.520 2.372 4.244 0.044 0.066
27 -4.866 -0.786 4.929 0.026 0.049
26 -3.500 -2.300 4.188
25 -3.476 -2.306 4.171 0.044 0.066
14EGCEE037 Page 14
27 0.000 0.000 0.000
25 4.892 0.835 4.963 0.026 0.049
28 -18.192 -4.152 18.660 -0.000 1.310
29 6.178 1.675 6.401 0.086 0.162
30 7.093 1.663 7.286 0.162 0.304
28 0.000 0.000 0.000
27 18.192 5.463 18.994 -0.000 1.310 0.968
8 0.570 -2.003 2.082 0.000 -4.368
6 -18.780 -3.510 19.106 0.060 -13.085
29 -2.400 -0.900 2.563
27 -6.093 -1.513 6.278 0.086 0.162
30 3.716 0.601 3.764 0.034 0.063
30 -10.600 -1.900 10.769
27 -6.932 -1.359 7.064 0.162 0.304
29 -3.683 -0.537 3.722 0.034 0.063
Total loss 17.594 22.233
Result:- we have successfully analysis and study of Load Flow Analysis for a 3 to 6 bus system,
in MATLAB, using the following methods :
(i) Gauss Seidal Method (ii.) Newton Raphson
(iii.) Fast Decoupled Method.
27 0.000 0.000 0.000
25 4.892 0.835 4.963 0.026 0.049
28 -18.192 -4.152 18.660 -0.000 1.310
29 6.178 1.675 6.401 0.086 0.162
30 7.093 1.663 7.286 0.162 0.304
28 0.000 0.000 0.000
27 18.192 5.463 18.994 -0.000 1.310 0.968
8 0.570 -2.003 2.082 0.000 -4.368
6 -18.780 -3.510 19.106 0.060 -13.085
29 -2.400 -0.900 2.563
27 -6.093 -1.513 6.278 0.086 0.162
30 3.716 0.601 3.764 0.034 0.063
30 -10.600 -1.900 10.769
27 -6.932 -1.359 7.064 0.162 0.304
29 -3.683 -0.537 3.722 0.034 0.063
Total loss 17.594 22.233
Result:- we have successfully analysis and study of Load Flow Analysis for a 3 to 6 bus system,
in MATLAB, using the following methods :
(i) Gauss Seidal Method (ii.) Newton Raphson
(iii.) Fast Decoupled Method.
14EGCEE037 Page 15
EXPERIMENT NO. 3
AIM: Study of voltage security in power system Real-Time Voltage Security Assessment
(RTVSA)
THEORY:
Voltage stability is the ability of a power system to maintain acceptable voltages at all buses in
the system under normal operating conditions and after being subjected to a disturbance. A
system enters a state of voltage instability when a disturbance, increase in load demand, or
change in system condition cause a progressive and uncontrollable decline in voltage. The main
factor causing voltage instability is the inability of the power system to meet the demand for
reactive power. Voltage collapse is the process or sequence of events accompanying voltage
instability which leads to a low unacceptable voltage profile in a significant part of the system.
For maintaining voltage security within the system, the following need to be monitored from
time to time :
1. Available voltage security margin.
2. The most dangerous stresses in the system leading to voltage collapse.
3. Worst-case contingencies resulting in voltage collapse and/or contingencie with
insufficient voltage stability margin.
4. Contingency ranking according to a severity index for voltage stability related
system problems.
5. Weakest elements within the grid and the regions most affected by potential voltage
problems.
6. Controls to increase the available stability margin and avoid instability.
7. Information about voltage problems at the look-ahead operating conditions and for
the worst-case contingencies (contingencies with large severity ranks) that may
appear in the future.
8. A real-time dispatcher‘s situational awareness-type wide area graphic and geographic
displays.
It is known that voltage magnitudes alone are poor indicators of voltage stability or security.
Voltages can be near normal with generators, synchronous condensers, and Static VAR
compensators (SVCs) near current limiting levels, thus resulting in a possible voltage collapse.
However, as a security problem distinct from voltage collapse, it is also desirable that the system
voltage magnitudes remain within limits, and some of the control actions to
maintain voltage magnitudes may also be of benefit in avoiding voltage instability. Sufficient
reactive power reserves at generators and SVCs contribute strongly to maintaining voltage
stability, but do not measure the ability of the transmission system to transmit reactive power.
Both voltage magnitudes and reactive load margins are useful indicators; however, the voltage
stability margin is the more accurate and complete metric for the proximity to voltage collapse.
EXPERIMENT NO. 3
AIM: Study of voltage security in power system Real-Time Voltage Security Assessment
(RTVSA)
THEORY:
Voltage stability is the ability of a power system to maintain acceptable voltages at all buses in
the system under normal operating conditions and after being subjected to a disturbance. A
system enters a state of voltage instability when a disturbance, increase in load demand, or
change in system condition cause a progressive and uncontrollable decline in voltage. The main
factor causing voltage instability is the inability of the power system to meet the demand for
reactive power. Voltage collapse is the process or sequence of events accompanying voltage
instability which leads to a low unacceptable voltage profile in a significant part of the system.
For maintaining voltage security within the system, the following need to be monitored from
time to time :
1. Available voltage security margin.
2. The most dangerous stresses in the system leading to voltage collapse.
3. Worst-case contingencies resulting in voltage collapse and/or contingencie with
insufficient voltage stability margin.
4. Contingency ranking according to a severity index for voltage stability related
system problems.
5. Weakest elements within the grid and the regions most affected by potential voltage
problems.
6. Controls to increase the available stability margin and avoid instability.
7. Information about voltage problems at the look-ahead operating conditions and for
the worst-case contingencies (contingencies with large severity ranks) that may
appear in the future.
8. A real-time dispatcher‘s situational awareness-type wide area graphic and geographic
displays.
It is known that voltage magnitudes alone are poor indicators of voltage stability or security.
Voltages can be near normal with generators, synchronous condensers, and Static VAR
compensators (SVCs) near current limiting levels, thus resulting in a possible voltage collapse.
However, as a security problem distinct from voltage collapse, it is also desirable that the system
voltage magnitudes remain within limits, and some of the control actions to
maintain voltage magnitudes may also be of benefit in avoiding voltage instability. Sufficient
reactive power reserves at generators and SVCs contribute strongly to maintaining voltage
stability, but do not measure the ability of the transmission system to transmit reactive power.
Both voltage magnitudes and reactive load margins are useful indicators; however, the voltage
stability margin is the more accurate and complete metric for the proximity to voltage collapse.
14EGCEE037 Page 16
MODES OF OPERATION
1) Real-Time Modes - Real-time operations mode
- Real-time look-ahead mode
Under the ‗Real Time Operations Mode‘, a real time assessment of the most current state
estimation is done.
On the other hand, in the ‗Real Time Look-Ahead Mode‘ we perform a 2-hour ―look-ahead‖
predictive assessment by applying planned outage information available and load forecast over
the next 2 hours.
This is done by means of simulation software.
2) Study Mode - Study mode offers off-line analysis capabilities on either the real-time data or
on modified version of real-time solved cases.
Such study cases are:
(i) Real time RTVSA solved cases archived overtime within the Flat Files Storage (under
Central Server)
(ii) Modified versions of the above mentioned real-time solved cases to study hypothetical
scenarios. For instance, a study mode user may extract a previously archived RTVSA solved
case from the Flat Files Storage, remove one or more transmission lines, manually specify
stressing directions, resolve using the RTVSA simulation capabilities and perform a complete
voltage security assessment, and export this as a new ―study case‖ to the central server if so
desired.
RTVSA CAPABILITIES
The RTVSA application shall offer the following categories of functional capabilities:
1) Contingency screening and ranking with respect to voltage limit violations or loading margins
associated with known stressing direction – The application should perform such contingency
analysis under all N-1 conditions and some user defined N-2 conditions within each 5 minute
real time cycle. A directional stressing, representative of the actual system loading conditions
based on the real time dispatch schedule and load forecast, will be used for this analysis and
the most binding contingency shall be identified.
2) Wide area monitoring capabilities offering real time situational awareness to the operators on
key indicators that are closely associated with voltage security – These include voltage
profiles at select buses, real or reactive reserves at key generators both under base case and the
most binding contingency within geographic visualization. It also includes animated power
flow visuals at the higher voltage levels (e.g. 500 kV, 230 kV, and 138 kV). The application
will also have the capability of sending real time alarms to the end-users on voltage violations
and insufficient real or reactive loading margins.
3) Real time voltage stability analysis with known stressing direction – The application shall
present the loading margins (real or reactive) to the point of collapse under the base case and
the most binding contingency, allowing for an additional 2.5% and 5% (user configurable)
safety margins for N-1 and N-2 contingencies, respectively.
MODES OF OPERATION
1) Real-Time Modes - Real-time operations mode
- Real-time look-ahead mode
Under the ‗Real Time Operations Mode‘, a real time assessment of the most current state
estimation is done.
On the other hand, in the ‗Real Time Look-Ahead Mode‘ we perform a 2-hour ―look-ahead‖
predictive assessment by applying planned outage information available and load forecast over
the next 2 hours.
This is done by means of simulation software.
2) Study Mode - Study mode offers off-line analysis capabilities on either the real-time data or
on modified version of real-time solved cases.
Such study cases are:
(i) Real time RTVSA solved cases archived overtime within the Flat Files Storage (under
Central Server)
(ii) Modified versions of the above mentioned real-time solved cases to study hypothetical
scenarios. For instance, a study mode user may extract a previously archived RTVSA solved
case from the Flat Files Storage, remove one or more transmission lines, manually specify
stressing directions, resolve using the RTVSA simulation capabilities and perform a complete
voltage security assessment, and export this as a new ―study case‖ to the central server if so
desired.
RTVSA CAPABILITIES
The RTVSA application shall offer the following categories of functional capabilities:
1) Contingency screening and ranking with respect to voltage limit violations or loading margins
associated with known stressing direction – The application should perform such contingency
analysis under all N-1 conditions and some user defined N-2 conditions within each 5 minute
real time cycle. A directional stressing, representative of the actual system loading conditions
based on the real time dispatch schedule and load forecast, will be used for this analysis and
the most binding contingency shall be identified.
2) Wide area monitoring capabilities offering real time situational awareness to the operators on
key indicators that are closely associated with voltage security – These include voltage
profiles at select buses, real or reactive reserves at key generators both under base case and the
most binding contingency within geographic visualization. It also includes animated power
flow visuals at the higher voltage levels (e.g. 500 kV, 230 kV, and 138 kV). The application
will also have the capability of sending real time alarms to the end-users on voltage violations
and insufficient real or reactive loading margins.
3) Real time voltage stability analysis with known stressing direction – The application shall
present the loading margins (real or reactive) to the point of collapse under the base case and
the most binding contingency, allowing for an additional 2.5% and 5% (user configurable)
safety margins for N-1 and N-2 contingencies, respectively.
14EGCEE037 Page 17
4) Quantify the efficacy of reactive power support at the most effective buses in terms of their
sensitivities (These sensitivities translate to a linear constraint and is representative of the
voltage stability limit associated with the unidirectional stressing which can be incorporated
into Security Constrained Unit Commitment (SCUC) and Security Constrained Economic
Dispatch (SCED) applications in the future).
5) Rank available corrective controls based on their effectiveness – These actions may include
enhancement controls that optimally increase the loading margin with respect to the stressing
direction, or remedial controls in the situation that a contingency may lead the system state
into an insecure region.
6) Identify the weak elements within the system associated with the one-dimensional stressing –
These are buses/regions with the grid that experience severe degradation in their voltage
profile at the voltage collapse caused by the additional stressing. The proportions by which the
voltage magnitudes will fall at these buses shall be presented. Comprehensive Voltage
Security Assessment under Multi-Directional Stressing This is generalization of the above
mentioned capabilities to a multi-directional stressing situation presenting the interaction and
tradeoffs between different stressing directions, and the associated interpretation of the safe-
operating region as a 2-D or 3-D (or higher dimensional) nomogram. The application shall:
a) Develop and update voltage security regions offline on demand based on a set of
predefined stressing directions – The boundaries of these regions shall be expressed as piece-
wise linear approximations (i.e., hyperplanes) in coordinates of key descriptive parameters
(such as MW transfers, total MW generation, total MW loading, etc) associated with the
stressing directions. As with the unidirectional stressing case, these security region boundaries
too shall be representative of the most binding contingency in the various stressing directions.
(These hyperplanes are representative of the voltage stability limits associated with various
stressing scenarios which can ultimately be embedded into SCUC and SCED applications).
b) Real time voltage security assessment with respect to the multidirectional stressing – The
voltage stability margins between the most current base case operating condition and the
security region boundaries shall be evaluated within each 5 minute real time cycle.
c) Suggest appropriate controls to enhance margin to the boundary – While the current
operating point is within the security region, the application should also suggest appropriate
control actions to optimally steer away from the closest boundary.
DATA DESCRIPTION
The following are details on the required list of data :
1. Detailed Network Model : Contains information in a volume sufficient for detailed power
flow simulations, under the CA ISO standards, i.e., branch information (connectivity data, line
impedance), breaker status, etc.
4) Quantify the efficacy of reactive power support at the most effective buses in terms of their
sensitivities (These sensitivities translate to a linear constraint and is representative of the
voltage stability limit associated with the unidirectional stressing which can be incorporated
into Security Constrained Unit Commitment (SCUC) and Security Constrained Economic
Dispatch (SCED) applications in the future).
5) Rank available corrective controls based on their effectiveness – These actions may include
enhancement controls that optimally increase the loading margin with respect to the stressing
direction, or remedial controls in the situation that a contingency may lead the system state
into an insecure region.
6) Identify the weak elements within the system associated with the one-dimensional stressing –
These are buses/regions with the grid that experience severe degradation in their voltage
profile at the voltage collapse caused by the additional stressing. The proportions by which the
voltage magnitudes will fall at these buses shall be presented. Comprehensive Voltage
Security Assessment under Multi-Directional Stressing This is generalization of the above
mentioned capabilities to a multi-directional stressing situation presenting the interaction and
tradeoffs between different stressing directions, and the associated interpretation of the safe-
operating region as a 2-D or 3-D (or higher dimensional) nomogram. The application shall:
a) Develop and update voltage security regions offline on demand based on a set of
predefined stressing directions – The boundaries of these regions shall be expressed as piece-
wise linear approximations (i.e., hyperplanes) in coordinates of key descriptive parameters
(such as MW transfers, total MW generation, total MW loading, etc) associated with the
stressing directions. As with the unidirectional stressing case, these security region boundaries
too shall be representative of the most binding contingency in the various stressing directions.
(These hyperplanes are representative of the voltage stability limits associated with various
stressing scenarios which can ultimately be embedded into SCUC and SCED applications).
b) Real time voltage security assessment with respect to the multidirectional stressing – The
voltage stability margins between the most current base case operating condition and the
security region boundaries shall be evaluated within each 5 minute real time cycle.
c) Suggest appropriate controls to enhance margin to the boundary – While the current
operating point is within the security region, the application should also suggest appropriate
control actions to optimally steer away from the closest boundary.
DATA DESCRIPTION
The following are details on the required list of data :
1. Detailed Network Model : Contains information in a volume sufficient for detailed power
flow simulations, under the CA ISO standards, i.e., branch information (connectivity data, line
impedance), breaker status, etc.
14EGCEE037 Page 18
2. System Component Status Information : Includes current status of generators, transmission
circuits, transformers, switching devices, and other components.
3. Available Power System Controls and their priorities : These include information of -
- Tap Changers
- Static VAR Compensator (SVC)
- Fixed and Controllable Shunt
- Generator Redispatch, etc.
4. Limits (Voltage, Thermal, MVar, Others) : Consists of operational limits of system
facilities/components that are to be specified in appropriate units, e.g. transformer limits in
MVA, line limits in Amps, etc
5. Generator Model : Required information for generator modeling, such as:
- MVA ratings
- Qmax, Qmin values
- Leading and lagging power factor
6. Distributed Slack Bus Information : Required for governor power flow simulations
7. Low Voltage Load Models : These models (static characteristics) should cover the low
voltage load behavior and voltage collapse situations.
8. HVDC Models & Control Schemes
9. Contingency List : Consists of -
- All (N-1) and some (N-2) contingencies, or
- User specified contingency list
- Any Remedial Action Schemes (RASs) associated with these contingencies
10. Stressing Directions & Descriptor Variables : Contains -
- Generator dispatch sequence & pattern
- Load stress pattern
(Should feature the capability to assign participation factors to loads on an individual, area or
zonal basis)
Descriptor variables are parameters that influence the voltage stability margin in certain parts of
the system (voltage stability problem areas). Examples of descriptor variables are: total area
load, power flows in certain transmission paths, total area generation, and so on. The operating
engineers‘ should be able to define/modify these variables for the known voltage problem areas
in the course of offline studies.
SPECIAL PROTECTION SCHEMES/REMEDIAL ACTION SCHEMES
During the system stressing process and contingency analysis, it is required for the RTVSA tool
to automatically trigger Remedial Action Schemes (RAS) or Special Protection Schemes (SPS)
to provide realistic voltage stability margins.
2. System Component Status Information : Includes current status of generators, transmission
circuits, transformers, switching devices, and other components.
3. Available Power System Controls and their priorities : These include information of -
- Tap Changers
- Static VAR Compensator (SVC)
- Fixed and Controllable Shunt
- Generator Redispatch, etc.
4. Limits (Voltage, Thermal, MVar, Others) : Consists of operational limits of system
facilities/components that are to be specified in appropriate units, e.g. transformer limits in
MVA, line limits in Amps, etc
5. Generator Model : Required information for generator modeling, such as:
- MVA ratings
- Qmax, Qmin values
- Leading and lagging power factor
6. Distributed Slack Bus Information : Required for governor power flow simulations
7. Low Voltage Load Models : These models (static characteristics) should cover the low
voltage load behavior and voltage collapse situations.
8. HVDC Models & Control Schemes
9. Contingency List : Consists of -
- All (N-1) and some (N-2) contingencies, or
- User specified contingency list
- Any Remedial Action Schemes (RASs) associated with these contingencies
10. Stressing Directions & Descriptor Variables : Contains -
- Generator dispatch sequence & pattern
- Load stress pattern
(Should feature the capability to assign participation factors to loads on an individual, area or
zonal basis)
Descriptor variables are parameters that influence the voltage stability margin in certain parts of
the system (voltage stability problem areas). Examples of descriptor variables are: total area
load, power flows in certain transmission paths, total area generation, and so on. The operating
engineers‘ should be able to define/modify these variables for the known voltage problem areas
in the course of offline studies.
SPECIAL PROTECTION SCHEMES/REMEDIAL ACTION SCHEMES
During the system stressing process and contingency analysis, it is required for the RTVSA tool
to automatically trigger Remedial Action Schemes (RAS) or Special Protection Schemes (SPS)
to provide realistic voltage stability margins.
14EGCEE037 Page 19
Modeling Details
Accurate modeling of voltage stability conditions and parameters that influence them is a must
for the RTVSA application. This includes the following requirements:
(1) Voltage stability conditions simulated using full power flow Jacobian singularity conditions.
(2) The algorithms used must converge accurately to the power system equilibrium in all cases in
which that equilibrium exists, including cases at and nearly at voltage collapse.
(3) Low voltage/voltage stability load models including the models reflecting the OLTC action
(e.g., constant active and reactive power for the OLTC regulation range), static
characteristics representing load behavior outside the regulation range of the OLTC, and
static characteristics approximately reflecting load behavior at the low voltage conditions
(4) Special Protection Schemes (SPS), Under-Voltage Protection schemes, and Remedial actions
schemes (including remote RAS)
(5) Consistent treatment of the discrete event sequences, for example, the switching sequence of
capacitors (non-uniqueness of these sequences for a given stressing path is not acceptable)
(6) Distributed slack bus/post-transient power flow (governor) model
(7) Generation dispatch options reflecting California ISO models and practices (e.g. generators
maximum and minimum active power output, reliability must-run units, emission-induced
constraints, etc)
(8) Multi-area power flow
(9) Adequate modeling of the reduced (equivalent) parts of the system, especially, voltage and
governor responses of the reduced part of the system.
The minimum requirement for the data that is required to correctly describe the system
equipments have been briefly mentioned here.
1. Bus data
- Consisting of all bus types: swing/slack, PQ, PV, HVDC
- Representation with breaker information and status
2. Transmission line data
- Consisting of: out-of-service, in-service, bypassed, and HVDC lines
- Representation with lossy model
3. Transformers and tap control data
- Model types: 2 & 3 winding transformers
- Control types: Fixed impedance with no control, voltage, MW, and MVAR control
4. Generator data
- Generator remote regulation
- Reactive power limits as Qmin/Qmax
5. Load data
- Static model as described under Modeling Details (3) above
6. Fixed Shunt data
7. Controllable shunt and static VAR devices (SVD) data
- SVD control types: locked, stepwise control, continuous control, stepwise control with
deadband, on/off control with deadband
- Models any controllable capacitive/inductive devices, such as:
Static VAR compensators (SVC)
Mechanically Switched Capacitors (MSC)
Synchronous Condensers
Modeling Details
Accurate modeling of voltage stability conditions and parameters that influence them is a must
for the RTVSA application. This includes the following requirements:
(1) Voltage stability conditions simulated using full power flow Jacobian singularity conditions.
(2) The algorithms used must converge accurately to the power system equilibrium in all cases in
which that equilibrium exists, including cases at and nearly at voltage collapse.
(3) Low voltage/voltage stability load models including the models reflecting the OLTC action
(e.g., constant active and reactive power for the OLTC regulation range), static
characteristics representing load behavior outside the regulation range of the OLTC, and
static characteristics approximately reflecting load behavior at the low voltage conditions
(4) Special Protection Schemes (SPS), Under-Voltage Protection schemes, and Remedial actions
schemes (including remote RAS)
(5) Consistent treatment of the discrete event sequences, for example, the switching sequence of
capacitors (non-uniqueness of these sequences for a given stressing path is not acceptable)
(6) Distributed slack bus/post-transient power flow (governor) model
(7) Generation dispatch options reflecting California ISO models and practices (e.g. generators
maximum and minimum active power output, reliability must-run units, emission-induced
constraints, etc)
(8) Multi-area power flow
(9) Adequate modeling of the reduced (equivalent) parts of the system, especially, voltage and
governor responses of the reduced part of the system.
The minimum requirement for the data that is required to correctly describe the system
equipments have been briefly mentioned here.
1. Bus data
- Consisting of all bus types: swing/slack, PQ, PV, HVDC
- Representation with breaker information and status
2. Transmission line data
- Consisting of: out-of-service, in-service, bypassed, and HVDC lines
- Representation with lossy model
3. Transformers and tap control data
- Model types: 2 & 3 winding transformers
- Control types: Fixed impedance with no control, voltage, MW, and MVAR control
4. Generator data
- Generator remote regulation
- Reactive power limits as Qmin/Qmax
5. Load data
- Static model as described under Modeling Details (3) above
6. Fixed Shunt data
7. Controllable shunt and static VAR devices (SVD) data
- SVD control types: locked, stepwise control, continuous control, stepwise control with
deadband, on/off control with deadband
- Models any controllable capacitive/inductive devices, such as:
Static VAR compensators (SVC)
Mechanically Switched Capacitors (MSC)
Synchronous Condensers
14EGCEE037 Page 20
VOLTAGE SECURITY ASSESSMENT
The display capabilities under this category demonstrate results of the Voltage Security
Assessment tool under the look-ahead scenario with respect to key stressing direction(s).
Such scenarios may be based on current operating conditions or under the worst case
contingency. These illustrate voltage security conditions and metrics that help users study
voltage stability and take decisions to prevent adverse situations. These capabilities include, but
are not limited to:
- Real and reactive loading margins
Margin at base case to point of collapse (POC)
Margin under worst case contingency base case to POC
- Contingency ranking based on severity index (voltage margin, loading margin, etc.)
- Operating nomograms i.e. the chart representing numerical relationships.
- Distance to instability
- Weak elements information
- Corrective actions (preventive control, enhancement control)
Result:- We have successfully Study of voltage security in power system Real-Time Voltage
Security Assessment (RTVSA).
VOLTAGE SECURITY ASSESSMENT
The display capabilities under this category demonstrate results of the Voltage Security
Assessment tool under the look-ahead scenario with respect to key stressing direction(s).
Such scenarios may be based on current operating conditions or under the worst case
contingency. These illustrate voltage security conditions and metrics that help users study
voltage stability and take decisions to prevent adverse situations. These capabilities include, but
are not limited to:
- Real and reactive loading margins
Margin at base case to point of collapse (POC)
Margin under worst case contingency base case to POC
- Contingency ranking based on severity index (voltage margin, loading margin, etc.)
- Operating nomograms i.e. the chart representing numerical relationships.
- Distance to instability
- Weak elements information
- Corrective actions (preventive control, enhancement control)
Result:- We have successfully Study of voltage security in power system Real-Time Voltage
Security Assessment (RTVSA).
14EGCEE037 Page 21
EXPERIMENT NO. 4
AIM : Study of overload security analysis in power system.
THEORY
The recent evolution of the electric power industry has brought about new needs in terms of
assessing the reliability of the transmission system. Perhaps the most important of these include
its accurate assessment and the need to integrate reliability into economic decision making.
These needs exist at the operational level. In this paper, we address them in terms of the one year
planning problem. There are today a number of commercial software packages that include the
influence of circuit overload in a reliability assessment scheme. All programs develop
probabilistic indices characterizing the power system reliability level, although some use
analytical approaches, sometimes called contingency enumeration, while others use Monte Carlo
simulation. Some of the most well known in North America include TRELSS [5], TPLAN,
PROCOSE [6], and CREAM. The approaches for assessing circuit overload for planning
purposes used in these and other programs have rested on two main assumptions.
1) The circuit overload reliability level is indicated by a measurement of the amount of load shed
necessary to avoid circuit overload; loss of load probability (LOLP) and expected unserved
energy (EUE) are two of the most common measurements used.
2) Measurements taken on one or a limited number of selected base cases are sufficient to
indicate the reliability of the system.
In this paper, we report on an approach that avoids these two assumptions. Our approach is
predicated on a desire to provide a measure of risk, a probabilistic expectation, as the product of
probability and monetary consequence of each outcome, summed over all possible outcomes. An
important attribute of using this measure is that it reflects event likelihood and consequence, the
two factors which, according to industry developed disturbance-performance criteria, determine
reliability level.
In measuring risk, it is essential that we distinguish between an outcome and a decision. One
important distinction is that an outcome is an unavoidable result of a decision. Based on this
distinction, we categorize transmission reinforcement, unit commitment, economic dispatch, and
load interruption as decisions.
In the context of overload assessment, the outcome of these decisions are the effects on the
circuits. These effects, which include equipment damage and equipment unavailability, are
random because they are heavily dependent on weather and on loadings, the randomness of the
latter caused mainly by uncertainties in demand and equipment outages. In contrast to the
commonly made assumption 1, our assessment is in terms of the probability and monetary
impact of these effects, given a decision, where the decision includes transmission reinforcement
and policies on unit commitment, redispatch and load interruption.
Calculation of the risk index has the distinct advantage of providing a uniform basis of
comparing various decisions.
EXPERIMENT NO. 4
AIM : Study of overload security analysis in power system.
THEORY
The recent evolution of the electric power industry has brought about new needs in terms of
assessing the reliability of the transmission system. Perhaps the most important of these include
its accurate assessment and the need to integrate reliability into economic decision making.
These needs exist at the operational level. In this paper, we address them in terms of the one year
planning problem. There are today a number of commercial software packages that include the
influence of circuit overload in a reliability assessment scheme. All programs develop
probabilistic indices characterizing the power system reliability level, although some use
analytical approaches, sometimes called contingency enumeration, while others use Monte Carlo
simulation. Some of the most well known in North America include TRELSS [5], TPLAN,
PROCOSE [6], and CREAM. The approaches for assessing circuit overload for planning
purposes used in these and other programs have rested on two main assumptions.
1) The circuit overload reliability level is indicated by a measurement of the amount of load shed
necessary to avoid circuit overload; loss of load probability (LOLP) and expected unserved
energy (EUE) are two of the most common measurements used.
2) Measurements taken on one or a limited number of selected base cases are sufficient to
indicate the reliability of the system.
In this paper, we report on an approach that avoids these two assumptions. Our approach is
predicated on a desire to provide a measure of risk, a probabilistic expectation, as the product of
probability and monetary consequence of each outcome, summed over all possible outcomes. An
important attribute of using this measure is that it reflects event likelihood and consequence, the
two factors which, according to industry developed disturbance-performance criteria, determine
reliability level.
In measuring risk, it is essential that we distinguish between an outcome and a decision. One
important distinction is that an outcome is an unavoidable result of a decision. Based on this
distinction, we categorize transmission reinforcement, unit commitment, economic dispatch, and
load interruption as decisions.
In the context of overload assessment, the outcome of these decisions are the effects on the
circuits. These effects, which include equipment damage and equipment unavailability, are
random because they are heavily dependent on weather and on loadings, the randomness of the
latter caused mainly by uncertainties in demand and equipment outages. In contrast to the
commonly made assumption 1, our assessment is in terms of the probability and monetary
impact of these effects, given a decision, where the decision includes transmission reinforcement
and policies on unit commitment, redispatch and load interruption.
Calculation of the risk index has the distinct advantage of providing a uniform basis of
comparing various decisions.
14EGCEE037 Page 22
I. SMV MODEL DESCRIPTION
The SMV model first uses the expected annual load curve, sampled hourly,1 to arrange the
maintenance and unit commitment
schedules, then employs time invariant variances to represent normally distributed load
uncertainties. The expected annual load curve can be obtained from load forecasting, or it can be
obtained from the load curve of the previous year, with an appropriate scaling to account for load
growth. There are various methods to identify the maintenance, unit commitment schedule and
load forecasting error.
For load forecasting error identification, we first employ time series analysis to identify the
structure and parameters of an ARIMA (autoregressive integrated moving average) model used
to represent the load series. This provides a load value for each hour. We assume that each
hourly load value used in our trajectory has associated with it some error.
Let‘s consider a single hour,h, a single contingency state,s, and a single branch,b, denoted by
and if we have a function which gives the expected monetary impact of each
flow Ib on branch b, the component risk, then we can compute the thermal overload (TOL) risk
for the particular contingency state, s, in hour,h, as
The total risk for this branch in hour over all contingency states is then
From (2), we may sum over all branches to obtain total risk for a particular hour, or we may sum
over all hours to obtain the cumulative risk for a particular branch. These kinds of calculations
reflect the decomposition capability of this approach and attractive for identifying the reasons for
high risk. In addition, we may evaluate total cumulative risk as
These calculations, together with those required to obtain , are referred to as thermal overload
risk assessment. Its use, together with the trajectory development, are illustrated in Fig. 1. From
this figure
I. SMV MODEL DESCRIPTION
The SMV model first uses the expected annual load curve, sampled hourly,1 to arrange the
maintenance and unit commitment
schedules, then employs time invariant variances to represent normally distributed load
uncertainties. The expected annual load curve can be obtained from load forecasting, or it can be
obtained from the load curve of the previous year, with an appropriate scaling to account for load
growth. There are various methods to identify the maintenance, unit commitment schedule and
load forecasting error.
For load forecasting error identification, we first employ time series analysis to identify the
structure and parameters of an ARIMA (autoregressive integrated moving average) model used
to represent the load series. This provides a load value for each hour. We assume that each
hourly load value used in our trajectory has associated with it some error.
Let‘s consider a single hour,h, a single contingency state,s, and a single branch,b, denoted by
and if we have a function which gives the expected monetary impact of each
flow Ib on branch b, the component risk, then we can compute the thermal overload (TOL) risk
for the particular contingency state, s, in hour,h, as
The total risk for this branch in hour over all contingency states is then
From (2), we may sum over all branches to obtain total risk for a particular hour, or we may sum
over all hours to obtain the cumulative risk for a particular branch. These kinds of calculations
reflect the decomposition capability of this approach and attractive for identifying the reasons for
high risk. In addition, we may evaluate total cumulative risk as
These calculations, together with those required to obtain , are referred to as thermal overload
risk assessment. Its use, together with the trajectory development, are illustrated in Fig. 1. From
this figure
14EGCEE037 Page 23
II. COMPONENT RISK
Equation (1) requires , which is the expected monetary impact on branch due to
overload given the flow on branch, b. If branch b is a transmission line, then, depending on the
weather conditions, conductor type, and flow duration, the flow Ib causes conductor heating
which can result in one or both of the following:
• Loss of clearance due to sag: Here, the thermal expansion of the conductor results in sag. In the
worst case, the line can touch an underlying object, resulting in a permanent fault and subsequent
outage.
• Loss of strength due to annealing: Annealing, the recrystallization of metal, is a gradual and
irreversible process when the grain matrix established by cold working is consumed causing loss
of tensile strength. In [12], we have shown how to use weather statistics to obtain f(θ ׀ Ib), the
pdf for conductor temperature, θ. This can be used to obtain the desired risk expression as
Where ImL1(θ) and ImL2(θ) and express the monetary impact on the transmission line of sag and
annealing, respectively, as a function of conductor temperature, also described in [12]. Equation
(4) can be evaluated for a range of flows, resulting in a component risk curve for branch b, as
shown in Fig. 2, where the pdfs for ambient temperature and wind speed are typically chosen.
The same pdf for ambient temperature is also used in transformer risk assessment.
We can then use an expression just like (4) to evaluate the thermal overload risk, except here, θ
represents the hottest spot temperature, and ImL1 and ImL2 represent the monetary impact on the
transformer of failure and loss of life, respectively, as described in [13] and [14]. With these
modifications, we can evaluate eqt. (4) for a range of flows, resulting in a component curve for
branch b, as shown in Fig. 3. Here, 1.0 pu risk equals the cost to rebuild the transformer. It is
chosen to be $1 000 000 in [13] and [14].
III. PROBABILITY DENSITIES FOR CIRCUIT FLOWS
The pdfs of currents can be identified by probabilistic load flow methods. The probabilistic load
flow proposed in [15] uses DC power flow and convolution to deal with load uncertainties. The
stochastic load flow (also called AEP method) proposed
in [16] linearizes the power system around the expected point (which is obtained by iterations in
order to account for the nonlinear nature of the power flow equations), and then applies linear
transformation of Gaussian distributions. Some refinements for these two methods have also
been proposed [17],[18]. In addition, efforts have been made to perform risk assessment for
power system planning in [19]. In this paper, we linearize the system around the operating point
at every hour, then use a convolution method to obtain pdfs.
A. Assumptions
We assume all loads are normally distributed random variables at every hour. Each hourly bus
load is assigned a mean and variance equal to a fixed percentage of the total load forecasted
mean and variance. We also assume the covariance matrix of loads is available. In practice, this
covariance matrix can be estimated by statistical methods. Generator outages are also considered.
We assume that each bus is a single independent generating company.
II. COMPONENT RISK
Equation (1) requires , which is the expected monetary impact on branch due to
overload given the flow on branch, b. If branch b is a transmission line, then, depending on the
weather conditions, conductor type, and flow duration, the flow Ib causes conductor heating
which can result in one or both of the following:
• Loss of clearance due to sag: Here, the thermal expansion of the conductor results in sag. In the
worst case, the line can touch an underlying object, resulting in a permanent fault and subsequent
outage.
• Loss of strength due to annealing: Annealing, the recrystallization of metal, is a gradual and
irreversible process when the grain matrix established by cold working is consumed causing loss
of tensile strength. In [12], we have shown how to use weather statistics to obtain f(θ ׀ Ib), the
pdf for conductor temperature, θ. This can be used to obtain the desired risk expression as
Where ImL1(θ) and ImL2(θ) and express the monetary impact on the transmission line of sag and
annealing, respectively, as a function of conductor temperature, also described in [12]. Equation
(4) can be evaluated for a range of flows, resulting in a component risk curve for branch b, as
shown in Fig. 2, where the pdfs for ambient temperature and wind speed are typically chosen.
The same pdf for ambient temperature is also used in transformer risk assessment.
We can then use an expression just like (4) to evaluate the thermal overload risk, except here, θ
represents the hottest spot temperature, and ImL1 and ImL2 represent the monetary impact on the
transformer of failure and loss of life, respectively, as described in [13] and [14]. With these
modifications, we can evaluate eqt. (4) for a range of flows, resulting in a component curve for
branch b, as shown in Fig. 3. Here, 1.0 pu risk equals the cost to rebuild the transformer. It is
chosen to be $1 000 000 in [13] and [14].
III. PROBABILITY DENSITIES FOR CIRCUIT FLOWS
The pdfs of currents can be identified by probabilistic load flow methods. The probabilistic load
flow proposed in [15] uses DC power flow and convolution to deal with load uncertainties. The
stochastic load flow (also called AEP method) proposed
in [16] linearizes the power system around the expected point (which is obtained by iterations in
order to account for the nonlinear nature of the power flow equations), and then applies linear
transformation of Gaussian distributions. Some refinements for these two methods have also
been proposed [17],[18]. In addition, efforts have been made to perform risk assessment for
power system planning in [19]. In this paper, we linearize the system around the operating point
at every hour, then use a convolution method to obtain pdfs.
A. Assumptions
We assume all loads are normally distributed random variables at every hour. Each hourly bus
load is assigned a mean and variance equal to a fixed percentage of the total load forecasted
mean and variance. We also assume the covariance matrix of loads is available. In practice, this
covariance matrix can be estimated by statistical methods. Generator outages are also considered.
We assume that each bus is a single independent generating company.
14EGCEE037 Page 24
B. Analytical Development
From the DC power flow formulation, we can obtain the following expressions for branch flows
corresponding to any outage state s, and the normal state 1.
where
• PG is the vector of real power generator levels at each bus.
• PD is the vector of real power load levels at each bus.
• Bps is the B-matrix for outage state .
• A(s) is the connection matrix of the network for outage state s, having rows corresponding to
buses (excluding the swing bus) and columns corresponding to branches.
• Pi(s) is the vector of branch power flows for outage state, s.
Subtracting (5) from (6) we get
where .
Since Xi(s)are independent of generation and load level, we can calculate and store them
beforehand to save computation time.If we set PG and PD to be their expected values for the hour,
and use a full AC power flow solution to obtain Pl(1), then (7) provides the expected flows for all
branches at the hour. We will use this below in obtaining the distribution of flows due to
uncertainty in generation. Now we define ∆PG and ∆PD as the vectors of random variables
corresponding to generation and load levels, respectively.
We describe each component of ∆PG with a two state probability mass function. We describe
each element of ∆PD with a normal distribution having a mean equal to 0 and a standard
deviation derived from our load model assumption. The vector of random variables
corresponding to variations in branch flows at outage state are then given as
(8)
and the vector of random variables corresponding to the branch flows
(9)
Substitution of (8) into (9) yields
(10)
Here PLG is the vector of random variables corresponding to branch flows due to variation in
generation, given by
(11)
The mean of PLG(s) is Pl(s) . Also, is the vector of random variables corresponding to branch
flows due to variations in
loads, given by
(12)
The mean of PLD(s) is zero.
The covariance matrix for PLD(s) is obtained from
(13)
and the diagonal elements of (13) are the line flow variances.
B. Analytical Development
From the DC power flow formulation, we can obtain the following expressions for branch flows
corresponding to any outage state s, and the normal state 1.
where
• PG is the vector of real power generator levels at each bus.
• PD is the vector of real power load levels at each bus.
• Bps is the B-matrix for outage state .
• A(s) is the connection matrix of the network for outage state s, having rows corresponding to
buses (excluding the swing bus) and columns corresponding to branches.
• Pi(s) is the vector of branch power flows for outage state, s.
Subtracting (5) from (6) we get
where .
Since Xi(s)are independent of generation and load level, we can calculate and store them
beforehand to save computation time.If we set PG and PD to be their expected values for the hour,
and use a full AC power flow solution to obtain Pl(1), then (7) provides the expected flows for all
branches at the hour. We will use this below in obtaining the distribution of flows due to
uncertainty in generation. Now we define ∆PG and ∆PD as the vectors of random variables
corresponding to generation and load levels, respectively.
We describe each component of ∆PG with a two state probability mass function. We describe
each element of ∆PD with a normal distribution having a mean equal to 0 and a standard
deviation derived from our load model assumption. The vector of random variables
corresponding to variations in branch flows at outage state are then given as
(8)
and the vector of random variables corresponding to the branch flows
(9)
Substitution of (8) into (9) yields
(10)
Here PLG is the vector of random variables corresponding to branch flows due to variation in
generation, given by
(11)
The mean of PLG(s) is Pl(s) . Also, is the vector of random variables corresponding to branch
flows due to variations in
loads, given by
(12)
The mean of PLD(s) is zero.
The covariance matrix for PLD(s) is obtained from
(13)
and the diagonal elements of (13) are the line flow variances.
14EGCEE037 Page 25
The pdfs for the elements of PLG(s) can be obtained by a convolution algorithm. In order to
facilitate speed of this algorithm while maintaining reasonable accuracy, we have used
segmentwise cluster-based convolution. The flowchart for thermal overload risk assessment is
shown in Fig. 4. For the first hour, we form the generation random
variable ∆PG, i.e., we obtain their pdfs. After the first hour, we use convolution (when some unit
is started) and deconvolution
(when some unit is shut down) to update ∆PG. Then for each hour and each credible contingency,
we linearize the system, calculate the pdfs of line flows contributed by generation and by load by
using (11) and (12), respectively. Then we convolve the two pdfs to get the pdfs for active flows,
which is the term . This term can be transformed into based on the
assumption that the reactive flow remains constant.
IV. CONCLUSION
We have proposed the sequential mean variance (SMV) model together with a risk index to
assess power system reliability over a mid-term planning period. We have shown that the SMV
model enables assessment of loading periods and inter-temporal affects that may not be captured
by so-called snapshot models. Yet, it does so with reduced computational requirements relative
to the sequential Monte Carlo model.
The risk index used in the hourly assessment provides a compact evaluation of the hour‘s
reliability level for overload that does not require the representation of the operator‘s load
shedding policy, considered here to be a decision which could be assessed by the risk index. This
risk assessment is performed based on linearization around the operating point and convolution
between random variables.
Result:- We have successfully Study of overload security analysis in power system.
The pdfs for the elements of PLG(s) can be obtained by a convolution algorithm. In order to
facilitate speed of this algorithm while maintaining reasonable accuracy, we have used
segmentwise cluster-based convolution. The flowchart for thermal overload risk assessment is
shown in Fig. 4. For the first hour, we form the generation random
variable ∆PG, i.e., we obtain their pdfs. After the first hour, we use convolution (when some unit
is started) and deconvolution
(when some unit is shut down) to update ∆PG. Then for each hour and each credible contingency,
we linearize the system, calculate the pdfs of line flows contributed by generation and by load by
using (11) and (12), respectively. Then we convolve the two pdfs to get the pdfs for active flows,
which is the term . This term can be transformed into based on the
assumption that the reactive flow remains constant.
IV. CONCLUSION
We have proposed the sequential mean variance (SMV) model together with a risk index to
assess power system reliability over a mid-term planning period. We have shown that the SMV
model enables assessment of loading periods and inter-temporal affects that may not be captured
by so-called snapshot models. Yet, it does so with reduced computational requirements relative
to the sequential Monte Carlo model.
The risk index used in the hourly assessment provides a compact evaluation of the hour‘s
reliability level for overload that does not require the representation of the operator‘s load
shedding policy, considered here to be a decision which could be assessed by the risk index. This
risk assessment is performed based on linearization around the operating point and convolution
between random variables.
Result:- We have successfully Study of overload security analysis in power system.
14EGCEE037 Page 26
EXPERIMENT NO. 5
AIM : Study the various methods of Economic Load Dispatch.
THEORY
Economic load dispatch (ELD) is an important function in power system planning and operation.
ELD solutions are found by solving the conventional load flow equations while at the same time
minimizing fuel costs. The resulting optimization problem has nonlinear constraints from the
load flow nodal equations and simple bound constraints on the variables from the load bus
voltage magnitudes.
Methods of Economic Load Dispatch :
1. Lambda Search
2. Gradient Algorithm
3. Newton‘s Method
For solving the problem with economic dispatch we first need to define or formulate the
problem.
Economic Load Dispatch (ELD) problem formulation
The ELD problem is considered as a general minimization problem with constraints, and can be
written in the following form:
Minimize f(x) (1)
Subject to: g(x) = 0 (2)
h(x) ≤ 0 (3)
f(x) is the objective function, g(x) and h(x) are respectively the set of equality and inequality
constraints. x is the vector of control and state variables. The control variables are generator
active and reactive power outputs, bus voltages, shunt capacitors/reactors and transformers tap-
setting. The state variables are voltage and angle of load buses.
Objective function
The objective function for the ELD reflects the costs associated with generating power in the
system. The quadratic cost model is used. The objective function for the entire power system can
then be written as the sum of the quadratic cost model for each generator:
(4)
Where, ng is the number of thermal units, Pgi is the active power generation at unit i and
ai, bi and ci are the cost coefficients of the i
th
generator.
Equality constraints
The equality constraints g(x) of the ELD problem are represented by the power balance
constraint, where the total power generation must cover the total power demand and the power
loss. This implies solving the load flow problem, which has equality constraints on active and
reactive power at each bus as follows [4]:
EXPERIMENT NO. 5
AIM : Study the various methods of Economic Load Dispatch.
THEORY
Economic load dispatch (ELD) is an important function in power system planning and operation.
ELD solutions are found by solving the conventional load flow equations while at the same time
minimizing fuel costs. The resulting optimization problem has nonlinear constraints from the
load flow nodal equations and simple bound constraints on the variables from the load bus
voltage magnitudes.
Methods of Economic Load Dispatch :
1. Lambda Search
2. Gradient Algorithm
3. Newton‘s Method
For solving the problem with economic dispatch we first need to define or formulate the
problem.
Economic Load Dispatch (ELD) problem formulation
The ELD problem is considered as a general minimization problem with constraints, and can be
written in the following form:
Minimize f(x) (1)
Subject to: g(x) = 0 (2)
h(x) ≤ 0 (3)
f(x) is the objective function, g(x) and h(x) are respectively the set of equality and inequality
constraints. x is the vector of control and state variables. The control variables are generator
active and reactive power outputs, bus voltages, shunt capacitors/reactors and transformers tap-
setting. The state variables are voltage and angle of load buses.
Objective function
The objective function for the ELD reflects the costs associated with generating power in the
system. The quadratic cost model is used. The objective function for the entire power system can
then be written as the sum of the quadratic cost model for each generator:
(4)
Where, ng is the number of thermal units, Pgi is the active power generation at unit i and
ai, bi and ci are the cost coefficients of the i
th
generator.
Equality constraints
The equality constraints g(x) of the ELD problem are represented by the power balance
constraint, where the total power generation must cover the total power demand and the power
loss. This implies solving the load flow problem, which has equality constraints on active and
reactive power at each bus as follows [4]:
14EGCEE037 Page 27
where: i=1,2,..., n and θij = θi - θj
Pi, Qi: injected active and reactive power at bus I
Pdi, Qdi: active and reactive power demand at bus i
Vi, θi: bus voltage magnitude and angle at bus i
Gij, Bij: conductance and susceptance of the (i,j) element in the admittance matrix.
Inequality constraints
The inequality constraints h(x) reflect the limits on physical devices in the power system as well
as the limits created to ensure system security:
• Upper and lower bounds on the active and reactive generations:
• Upper and lower bounds on the tap ratio (t) and phase shifting (α) of variable transformers:
• Upper limit on the active power flow (Pij) of line i-j:
• Upper and lower bounds on the bus voltage magnitude:
The ELD problem formulated above can be solved by any of the above listed three methods.
This is explained as follows :
1) LAMBDA SEARCH METHOD
This method can be described for 2 cases – neglecting or including losses.
A. ELD Neglecting Losses
Let us consider a system of N thermal-generating units connected to a single bus-bar serving
a received electrical load Pload. The input to each unit, shown as Fi, represents the cost
rate of the unit. The output of each unit, Pi, is the electrical power generated by that
particular unit. The total cost rate of this system is, of course, the sum of the costs of
each of the individual units. The essential constraint on the operation of this system is that
the sum of the output powers must equal the load demand. That is, an objective function, FT,
is equal to the total cost for supplying the indicated load. The problem is to minimize FT
subject to the constraints specified above.
Note that any transmission losses are neglected and any operating limits are not explicitly
stated when formulating this problem.
where: i=1,2,..., n and θij = θi - θj
Pi, Qi: injected active and reactive power at bus I
Pdi, Qdi: active and reactive power demand at bus i
Vi, θi: bus voltage magnitude and angle at bus i
Gij, Bij: conductance and susceptance of the (i,j) element in the admittance matrix.
Inequality constraints
The inequality constraints h(x) reflect the limits on physical devices in the power system as well
as the limits created to ensure system security:
• Upper and lower bounds on the active and reactive generations:
• Upper and lower bounds on the tap ratio (t) and phase shifting (α) of variable transformers:
• Upper limit on the active power flow (Pij) of line i-j:
• Upper and lower bounds on the bus voltage magnitude:
The ELD problem formulated above can be solved by any of the above listed three methods.
This is explained as follows :
1) LAMBDA SEARCH METHOD
This method can be described for 2 cases – neglecting or including losses.
A. ELD Neglecting Losses
Let us consider a system of N thermal-generating units connected to a single bus-bar serving
a received electrical load Pload. The input to each unit, shown as Fi, represents the cost
rate of the unit. The output of each unit, Pi, is the electrical power generated by that
particular unit. The total cost rate of this system is, of course, the sum of the costs of
each of the individual units. The essential constraint on the operation of this system is that
the sum of the output powers must equal the load demand. That is, an objective function, FT,
is equal to the total cost for supplying the indicated load. The problem is to minimize FT
subject to the constraints specified above.
Note that any transmission losses are neglected and any operating limits are not explicitly
stated when formulating this problem.
14EGCEE037 Page 28
This is solved by Lagrangian multiplier method. Let λ be the Lagrangian multiplier, then the
function becomes,
Finding the minimum value of the above function gives the following equation,
After solving the above equation it gives following condition for economic dispatch neglecting
losses
When we recognize the inequality constraints, then the necessary conditions may be
expanded slightly as shown in the following equations :
B. ELD Considering Losses
In this case the power generated will be equal to the sum of power demand and power losses,
Ploss. Thus the following equation follows :
(11)
Now by Lagrange‘s multiplier method we get the following equation
(12)
(13)
Equations (11) to (13) are known as coordination equations.
The procedure for solution by this method is defined below :
Step 1: Pick a set of starting values for PI, P2, and P3 that sum to the load.
Step 2: Calculate the incremental losses dPL/dPi, as well as the total losses.
Step 3: Assume the value of λ greater than the largest intercept of the cost functions. Evaluate
values of Pl, P2…..Pn by (11), (12) and (13).
This is solved by Lagrangian multiplier method. Let λ be the Lagrangian multiplier, then the
function becomes,
Finding the minimum value of the above function gives the following equation,
After solving the above equation it gives following condition for economic dispatch neglecting
losses
When we recognize the inequality constraints, then the necessary conditions may be
expanded slightly as shown in the following equations :
B. ELD Considering Losses
In this case the power generated will be equal to the sum of power demand and power losses,
Ploss. Thus the following equation follows :
(11)
Now by Lagrange‘s multiplier method we get the following equation
(12)
(13)
Equations (11) to (13) are known as coordination equations.
The procedure for solution by this method is defined below :
Step 1: Pick a set of starting values for PI, P2, and P3 that sum to the load.
Step 2: Calculate the incremental losses dPL/dPi, as well as the total losses.
Step 3: Assume the value of λ greater than the largest intercept of the cost functions. Evaluate
values of Pl, P2…..Pn by (11), (12) and (13).
14EGCEE037 Page 29
Step 4: Check whether generations P1, P2….Pn are within the prescribed limits.
If Pi>Pmax, then set Pi=Pmax, and distribute the remaining load between the remaining units by
equation (11), or
If Pi<Pmin, then set Pi=Pmin, and distribute the remaining load between the remaining units by
equation (11
Step 5 : Check if value of equation (11) < € (defined tolerance), then increase value of λ slightly
by λ+Δλ, otherwise decrease it by λ-Δλ, and go to step 2.
If the solution of (11) is within specified tolerance limit, then go to step 6.
Step 6 : Print the results of Pi and Ploss, and stop.
2) GRADIENT ALGORITHM METHOD
Note that the lambda search technique always requires that one be able to find the power
output of a generator, given an incremental cost for that generator. In the case of a quadratic
function for the cost function, or in the case where the incremental cost function is
represented by a piecewise linear function, this is possible. However, it is often the case
that the cost function is much more complex, such as the one below:
In this case, we shall propose that a more basic method of solution for the optimum be
found. This method works on the principle that the minimum of a function, f(x), can be
found by a series of steps that always take us in a downward direction. From any starting point,
x
o
, we may find the direction of ―steepest descent‖ by noting that the gradient of f, i.e.,
always points in the direction of maximum ascent. Therefore, if we want to move in the
direction of maximum descent, we negate the gradient. Then we should go from x
o
to x
1
using:
Where α is a scalar to allow us to guarantee that the process converges. The best value of
α must be determined by experiment.
Economic Dispatch by Gradient Search
In the case of power system economic dispatch this becomes:
and the object is to drive the function to its minimum. However, we have to be concerned
with the constraint function:
Step 4: Check whether generations P1, P2….Pn are within the prescribed limits.
If Pi>Pmax, then set Pi=Pmax, and distribute the remaining load between the remaining units by
equation (11), or
If Pi<Pmin, then set Pi=Pmin, and distribute the remaining load between the remaining units by
equation (11
Step 5 : Check if value of equation (11) < € (defined tolerance), then increase value of λ slightly
by λ+Δλ, otherwise decrease it by λ-Δλ, and go to step 2.
If the solution of (11) is within specified tolerance limit, then go to step 6.
Step 6 : Print the results of Pi and Ploss, and stop.
2) GRADIENT ALGORITHM METHOD
Note that the lambda search technique always requires that one be able to find the power
output of a generator, given an incremental cost for that generator. In the case of a quadratic
function for the cost function, or in the case where the incremental cost function is
represented by a piecewise linear function, this is possible. However, it is often the case
that the cost function is much more complex, such as the one below:
In this case, we shall propose that a more basic method of solution for the optimum be
found. This method works on the principle that the minimum of a function, f(x), can be
found by a series of steps that always take us in a downward direction. From any starting point,
x
o
, we may find the direction of ―steepest descent‖ by noting that the gradient of f, i.e.,
always points in the direction of maximum ascent. Therefore, if we want to move in the
direction of maximum descent, we negate the gradient. Then we should go from x
o
to x
1
using:
Where α is a scalar to allow us to guarantee that the process converges. The best value of
α must be determined by experiment.
Economic Dispatch by Gradient Search
In the case of power system economic dispatch this becomes:
and the object is to drive the function to its minimum. However, we have to be concerned
with the constraint function:
14EGCEE037 Page 30
To solve the economic dispatch problem which involves minimizing the objective function
and keeping the equality constraint, we must apply the gradient technique directly to the
Lagrange function itself. The Lagrange function is:
and the gradient of this function is:
The problem with this formulation is the lack of a guarantee that the new points generated each
step will lie on the surface φ. We shall see that this can be overcome by a simple
variation of the gradient method. The economic dispatch algorithm requires a starting λ
value and starting values for Pl, P2, and P3. The gradient for is calculated as above and
the new values of λ, Pl, P2, and P3, etc., are found from:
where the vector x is:
3) NEWTON'S METHOD
We may wish to go a further step beyond the simple gradient method and try to solve the
economic dispatch by observing that the aim is to always drive :
Since this is a vector function, we can formulate the problem as one of finding the correction that
exactly drives the gradient to zero (i.e., to a vector, all of whose elements are zero). We know
how to find this, however, since we can
use Newton's method. Newton's method for a function of more than one variable is developed as
follows.
Suppose we wish to drive the function g(x) to zero. The function g is a vector and the unknowns,
x, are also vectors. Then, to use Newton's method, we observe:
If we let the function be defined as:
To solve the economic dispatch problem which involves minimizing the objective function
and keeping the equality constraint, we must apply the gradient technique directly to the
Lagrange function itself. The Lagrange function is:
and the gradient of this function is:
The problem with this formulation is the lack of a guarantee that the new points generated each
step will lie on the surface φ. We shall see that this can be overcome by a simple
variation of the gradient method. The economic dispatch algorithm requires a starting λ
value and starting values for Pl, P2, and P3. The gradient for is calculated as above and
the new values of λ, Pl, P2, and P3, etc., are found from:
where the vector x is:
3) NEWTON'S METHOD
We may wish to go a further step beyond the simple gradient method and try to solve the
economic dispatch by observing that the aim is to always drive :
Since this is a vector function, we can formulate the problem as one of finding the correction that
exactly drives the gradient to zero (i.e., to a vector, all of whose elements are zero). We know
how to find this, however, since we can
use Newton's method. Newton's method for a function of more than one variable is developed as
follows.
Suppose we wish to drive the function g(x) to zero. The function g is a vector and the unknowns,
x, are also vectors. Then, to use Newton's method, we observe:
If we let the function be defined as:
14EGCEE037 Page 31
Then
which is the familiar Jacobian matrix. The adjustment at each step is then:
Now, if we let the g function be the gradient vector we get:
For our economic dispatch problem this takes the form:
and is as it was defined before. The Jacobian matrix now becomes one made up of second
derivatives and is called the Hessian matrix:
Generally, Newton's method will solve for the correction that is much closer to the minimum
generation cost in one step than would the gradient method.
Result:- We have successfully Study the various methods of Economic Load Dispatch.
Then
which is the familiar Jacobian matrix. The adjustment at each step is then:
Now, if we let the g function be the gradient vector we get:
For our economic dispatch problem this takes the form:
and is as it was defined before. The Jacobian matrix now becomes one made up of second
derivatives and is called the Hessian matrix:
Generally, Newton's method will solve for the correction that is much closer to the minimum
generation cost in one step than would the gradient method.
Result:- We have successfully Study the various methods of Economic Load Dispatch.
14EGCEE037 Page 32
EXPERIMENT NO. 6
AIM:- Study of transient stability analysis using MATLAB/ETAP Software.
THEORY
Stability: Stability problem is concerned with the behavior of power system when it is subjected
to disturbance and is classified into small signal stability problem if the disturbances are small
and transient stability problem when the disturbances are large.
Transient stability: When a power system is under steady state, the load plus transmission loss
equals to the generation in the system. The generating units run at synchronous speed and system
frequency, voltage, current and power flows are steady. When a large disturbance such as three
phase fault, loss of load, loss of generation etc., occurs the power balance is upset and the
generating units rotors experience either acceleration or deceleration. The system may come
back to a steady state condition maintaining synchronism or it may break into subsystems or one
or more machines may pull out of synchronism.
Small signal stability: When a power system is under steady state, normal operating condition,
the system may be subjected to small disturbances such as variation in load and generation,
change in field voltage, change in mechanical toque etc., the nature of system response to small
disturbance depends on the operating conditions, the transmission system strength, types of
controllers etc.
Instability that may result from small disturbance may be of two forms
(i) Steady increase in rotor angle due to lack of synchronizing torque.
(ii) (ii) Rotor oscillations of increasing magnitude due to lack of sufficient damping torque.
FORMULA
Reactive power Qe = sin(cos-1(p.f))
Reactive power
Stator Current
Voltage behind transient condition
Voltage of infinite bus
Where,
EXPERIMENT NO. 6
AIM:- Study of transient stability analysis using MATLAB/ETAP Software.
THEORY
Stability: Stability problem is concerned with the behavior of power system when it is subjected
to disturbance and is classified into small signal stability problem if the disturbances are small
and transient stability problem when the disturbances are large.
Transient stability: When a power system is under steady state, the load plus transmission loss
equals to the generation in the system. The generating units run at synchronous speed and system
frequency, voltage, current and power flows are steady. When a large disturbance such as three
phase fault, loss of load, loss of generation etc., occurs the power balance is upset and the
generating units rotors experience either acceleration or deceleration. The system may come
back to a steady state condition maintaining synchronism or it may break into subsystems or one
or more machines may pull out of synchronism.
Small signal stability: When a power system is under steady state, normal operating condition,
the system may be subjected to small disturbances such as variation in load and generation,
change in field voltage, change in mechanical toque etc., the nature of system response to small
disturbance depends on the operating conditions, the transmission system strength, types of
controllers etc.
Instability that may result from small disturbance may be of two forms
(i) Steady increase in rotor angle due to lack of synchronizing torque.
(ii) (ii) Rotor oscillations of increasing magnitude due to lack of sufficient damping torque.
FORMULA
Reactive power Qe = sin(cos-1(p.f))
Reactive power
Stator Current
Voltage behind transient condition
Voltage of infinite bus
Where,
14EGCEE037 Page 33
Angular separation between and
Pre fault operation:
Power
During fault condition:
Critical clearing angle:
Critical clearing time:
Secs
PROCEDURE
1. Enter the command window of the MATLAB.
2. Create a new M – file by selecting File - New – M – File
3. Type and save the program.
4. Execute the program by pressing Tools – Run
5. View the results.
Angular separation between and
Pre fault operation:
Power
During fault condition:
Critical clearing angle:
Critical clearing time:
Secs
PROCEDURE
1. Enter the command window of the MATLAB.
2. Create a new M – file by selecting File - New – M – File
3. Type and save the program.
4. Execute the program by pressing Tools – Run
5. View the results.
14EGCEE037 Page 34
EXERCISE
1. Transient stability analysis of a 9-bus, 3-machine, 60 Hz power system with the following
system modeling requirements:
i.) Classical model for all synchronous machines, models for excitation and speed governing
systems not included.
(a) Simulate a three-phase fault at the end of the line from bus 5 to bus 7 near bus 7 at time = 0.0
sec. Assume that the fault is cleared successfully by opening the line 5-7 after 5 cycles ( 0.083
sec) . Observe the system for 2.0 seconds
(b) Obtain the following time domain plots:
- Relative angles of machines 2 and 3 with respect to machine 1
- Angular speed deviations of machines 1, 2 and 3 from synchronous speed
- Active power variation of machines 1, 2 and 3.
(c) Determine the critical clearing time by progressively increasing the fault clearing time.
PROGRAM
alpha = [500; 400; 200];
beta = [5.3; 5.5; 5.8]; gamma = [0.004; 0.006; 0.009];
PD = 800;
DelP = 10;
lamda = input('Enter estimated value of Lamda = ');
fprintf(' ')
disp(['Lamda P1 P2 P3 DP'...
' grad Delamda'])
iter = 0;
while abs(DelP) >= 0.001
EXERCISE
1. Transient stability analysis of a 9-bus, 3-machine, 60 Hz power system with the following
system modeling requirements:
i.) Classical model for all synchronous machines, models for excitation and speed governing
systems not included.
(a) Simulate a three-phase fault at the end of the line from bus 5 to bus 7 near bus 7 at time = 0.0
sec. Assume that the fault is cleared successfully by opening the line 5-7 after 5 cycles ( 0.083
sec) . Observe the system for 2.0 seconds
(b) Obtain the following time domain plots:
- Relative angles of machines 2 and 3 with respect to machine 1
- Angular speed deviations of machines 1, 2 and 3 from synchronous speed
- Active power variation of machines 1, 2 and 3.
(c) Determine the critical clearing time by progressively increasing the fault clearing time.
PROGRAM
alpha = [500; 400; 200];
beta = [5.3; 5.5; 5.8]; gamma = [0.004; 0.006; 0.009];
PD = 800;
DelP = 10;
lamda = input('Enter estimated value of Lamda = ');
fprintf(' ')
disp(['Lamda P1 P2 P3 DP'...
' grad Delamda'])
iter = 0;
while abs(DelP) >= 0.001
14EGCEE037 Page 35
iter = iter + 1;
P = (lamda - beta)./(2*gamma);
DelP = PD - sum(P);
J = sum(ones(length(gamma),1)./(2*gamma));
Delamda = DelP/J;
disp([lamda,P(1),P(2),P(3),DelP,J,Delamda])
lamda = lamda + Delamda;
end
totalcost = sum(alpha + beta.*P + gamma.*P.^2)
----------------------------------------------------------------------------
Pm = 0.8; E = 1.17; V = 1.0;
X1 = 0.65; X2 = inf; X3 = 0.65;
eacfault (Pm, E, V, X1, X2, X3)
For b)
Pm = 0.8; E = 1.17; V = 1.0;
X1 = 0.65; X2 = 1.8; X3 = 0.8;
eacfault (Pm, E, V, X1, X2, X3)
------------------------------------------------------------------------------
function eacpower(P0, E, V, X)
if exist('P0')~=1
P0 = input('Generator initial power in p.u. P0 = '); else, end
if exist('E')~=1
E = input('Generator e.m.f. in p.u. E = '); else, end
if exist('V')~=1
V = input('Infinite bus-bar voltage in p.u. V = '); else, end
if exist('X')~=1
X = input('Reactance between internal emf and infinite bus in p.u. X = '); else, end
Pemax= E*V/X;
if P0 >= Pemax
fprintf('\nP0 must be less than the peak electrical power Pemax = %5.3f p.u. \n', Pemax)
fprintf('Try again. \n\n')
return, end
d0=asin(P0/Pemax);
delta = 0:.01:pi;
Pe = Pemax*sin(delta);
dmax=pi;
Ddmax=1;
while abs(Ddmax) > 0.00001
iter = iter + 1;
P = (lamda - beta)./(2*gamma);
DelP = PD - sum(P);
J = sum(ones(length(gamma),1)./(2*gamma));
Delamda = DelP/J;
disp([lamda,P(1),P(2),P(3),DelP,J,Delamda])
lamda = lamda + Delamda;
end
totalcost = sum(alpha + beta.*P + gamma.*P.^2)
----------------------------------------------------------------------------
Pm = 0.8; E = 1.17; V = 1.0;
X1 = 0.65; X2 = inf; X3 = 0.65;
eacfault (Pm, E, V, X1, X2, X3)
For b)
Pm = 0.8; E = 1.17; V = 1.0;
X1 = 0.65; X2 = 1.8; X3 = 0.8;
eacfault (Pm, E, V, X1, X2, X3)
------------------------------------------------------------------------------
function eacpower(P0, E, V, X)
if exist('P0')~=1
P0 = input('Generator initial power in p.u. P0 = '); else, end
if exist('E')~=1
E = input('Generator e.m.f. in p.u. E = '); else, end
if exist('V')~=1
V = input('Infinite bus-bar voltage in p.u. V = '); else, end
if exist('X')~=1
X = input('Reactance between internal emf and infinite bus in p.u. X = '); else, end
Pemax= E*V/X;
if P0 >= Pemax
fprintf('\nP0 must be less than the peak electrical power Pemax = %5.3f p.u. \n', Pemax)
fprintf('Try again. \n\n')
return, end
d0=asin(P0/Pemax);
delta = 0:.01:pi;
Pe = Pemax*sin(delta);
dmax=pi;
Ddmax=1;
while abs(Ddmax) > 0.00001
14EGCEE037 Page 36
Df = cos(d0) - (sin(dmax)*(dmax-d0)+cos(dmax));
J=cos(dmax)*(dmax-d0);
Ddmax=Df/J;
dmax=dmax+Ddmax;
end
dc=pi-dmax;
Pm2=Pemax*sin(dc);
Pmx =[0 pi-d0]*180/pi; Pmy=[P0 P0];
Pm2x=[0 dmax]*180/pi; Pm2y=[Pm2 Pm2];
x0=[d0 d0]*180/pi; y0=[0 Pm2]; xc=[dc dc]*180/pi; yc=[0 Pemax*sin(dc)];
xm=[dmax dmax]*180/pi; ym=[0 Pemax*sin(dmax)];
d0=d0*180/pi; dmax=dmax*180/pi; dc=dc*180/pi;
x=(d0:.1:dc);
y=Pemax*sin(x*pi/180);
%y1=Pe2max*sin(d0*pi/180);
%y2=Pe2max*sin(dc*pi/180);
x=[d0 x dc];
y=[Pm2 y Pm2];
xx=dc:.1:dmax;
h=Pemax*sin(xx*pi/180);
xx=[dc xx dmax];
hh=[Pm2 h Pm2];
delta=delta*180/pi;
%clc
fprintf('\nInitial power =%7.3f p.u.\n', P0)
fprintf('Initial power angle =%7.3f degrees \n', d0)
fprintf('Sudden additional power =%7.3f p.u.\n', Pm2-P0)
fprintf('Total power for critical stability =%7.3f p.u.\n', Pm2)
fprintf('Maximum angle swing =%7.3f degrees \n', dmax)
fprintf('New operating angle =%7.3f degrees \n\n\n', dc)
fill(x,y,'m')
hold;
fill(xx,hh,'c')
plot(delta, Pe,'-', Pmx, Pmy,'g', Pm2x,Pm2y,'g', x0,y0,'c', xc,yc, xm,ym,'r'), grid
Title('Equal-area criterion applied to the sudden change in power')
xlabel('Power angle, degree'), ylabel(' Power, per unit')
axis([0 180 0 1.1*Pemax])
hold off;
Df = cos(d0) - (sin(dmax)*(dmax-d0)+cos(dmax));
J=cos(dmax)*(dmax-d0);
Ddmax=Df/J;
dmax=dmax+Ddmax;
end
dc=pi-dmax;
Pm2=Pemax*sin(dc);
Pmx =[0 pi-d0]*180/pi; Pmy=[P0 P0];
Pm2x=[0 dmax]*180/pi; Pm2y=[Pm2 Pm2];
x0=[d0 d0]*180/pi; y0=[0 Pm2]; xc=[dc dc]*180/pi; yc=[0 Pemax*sin(dc)];
xm=[dmax dmax]*180/pi; ym=[0 Pemax*sin(dmax)];
d0=d0*180/pi; dmax=dmax*180/pi; dc=dc*180/pi;
x=(d0:.1:dc);
y=Pemax*sin(x*pi/180);
%y1=Pe2max*sin(d0*pi/180);
%y2=Pe2max*sin(dc*pi/180);
x=[d0 x dc];
y=[Pm2 y Pm2];
xx=dc:.1:dmax;
h=Pemax*sin(xx*pi/180);
xx=[dc xx dmax];
hh=[Pm2 h Pm2];
delta=delta*180/pi;
%clc
fprintf('\nInitial power =%7.3f p.u.\n', P0)
fprintf('Initial power angle =%7.3f degrees \n', d0)
fprintf('Sudden additional power =%7.3f p.u.\n', Pm2-P0)
fprintf('Total power for critical stability =%7.3f p.u.\n', Pm2)
fprintf('Maximum angle swing =%7.3f degrees \n', dmax)
fprintf('New operating angle =%7.3f degrees \n\n\n', dc)
fill(x,y,'m')
hold;
fill(xx,hh,'c')
plot(delta, Pe,'-', Pmx, Pmy,'g', Pm2x,Pm2y,'g', x0,y0,'c', xc,yc, xm,ym,'r'), grid
Title('Equal-area criterion applied to the sudden change in power')
xlabel('Power angle, degree'), ylabel(' Power, per unit')
axis([0 180 0 1.1*Pemax])
hold off;
14EGCEE037 Page 37
RESULT
RESULT
14EGCEE037 Page 38
(Beyond the Syllabus)
EXPERIMENT 7
AIM : To develop a program to obtain bus impedance matrix z-bus of the given power system.
THEORY:
Z-bus matrix is an important matrix used in different kinds of power system study such as short
circuit study, load flow study etc.
In short circuit analysis the generator uses transformer impedance must be taken into account. In
quality analysis the two-short element are neglected by forming the z-bus matrix which is used to
compute the voltage distribution factor. This can be largely obtained by reversing the y-bus
formed by resection method or by analytical method.
Taking inverse of the y-bus for large system in time conditioning managing modification in the
system requires whole process to be repeated to voltage changes in the system. In such cases z-
bus computed to z-bus solving algorithm.
ALGORITHM:
Step 1: Read the values such as number of lines, number of buses and line data, generator data
and transformer data.
Step 2: Initialize y-bus matrix y-bus[i] [j] =complex.
Step 3: Compute y-bus matrix by considering only line data.
Step 4: Modifies the y-bus matrix by adding the transformer and the generator admittance to the
respected diagonal elements for y-bus matrix.
Step 5: Compute the z-bus matrix by the modified y-bus matrix.
Step 6: Check the iteration by multiplying modified y-bus matrix to z-bus and check whether the
resulting matrix is a unit matrix or not.
Step 7: Print the z-bus matrix.
(Beyond the Syllabus)
EXPERIMENT 7
AIM : To develop a program to obtain bus impedance matrix z-bus of the given power system.
THEORY:
Z-bus matrix is an important matrix used in different kinds of power system study such as short
circuit study, load flow study etc.
In short circuit analysis the generator uses transformer impedance must be taken into account. In
quality analysis the two-short element are neglected by forming the z-bus matrix which is used to
compute the voltage distribution factor. This can be largely obtained by reversing the y-bus
formed by resection method or by analytical method.
Taking inverse of the y-bus for large system in time conditioning managing modification in the
system requires whole process to be repeated to voltage changes in the system. In such cases z-
bus computed to z-bus solving algorithm.
ALGORITHM:
Step 1: Read the values such as number of lines, number of buses and line data, generator data
and transformer data.
Step 2: Initialize y-bus matrix y-bus[i] [j] =complex.
Step 3: Compute y-bus matrix by considering only line data.
Step 4: Modifies the y-bus matrix by adding the transformer and the generator admittance to the
respected diagonal elements for y-bus matrix.
Step 5: Compute the z-bus matrix by the modified y-bus matrix.
Step 6: Check the iteration by multiplying modified y-bus matrix to z-bus and check whether the
resulting matrix is a unit matrix or not.
Step 7: Print the z-bus matrix.
14EGCEE037 Page 39
HAND CALCULATION
Step1: Read the line data.
Step2: Initialize y-bus data. y-bus[i] [j] =complex (0.0, 0.0).
Step3: Compute y-bus matrix.
Y11 = (1/Z12) + (1/Z14) + Yc12 + Yc14
= (1/ (0.1+0.4j)) + (1/ (0.1+0.25j)) + 0.0075j + 0.006j
= 1.3430 - 4.981j
Y22 = (1/Z23) + (1/Z24) + (1/Z21) + Yc23 + Yc24 + Y21
= (1/ (0.15+0.6j)) + (1/ (0.25+0.7j)) + (1/ (0.1+0.4j)) + 0.01j + 0.015j + 0.0075j
= 1.4329 - 5.1562j
Y33 = (1/Z34) + (1/Z32) + Y34 + Y32
= (1/ (0.018+0.55j)) + (1/ (0.15+0.6j)) + 0.009j + 0.01j
= 0.926 - 3.1919j
Y13 = Y31= 0.0006 + 0.00j
Y23 = -1/ (0.15 + 0.6j)
= -0.8922 + 1.5686j
Y14 = Y41 = -1/ (0.1 + 0.35j)
= 0.7547 + 2.64j
Y24 = Y42 = -1/ (0.2 + 0.7j)
= -0.425 + 1.2670j
Y34 = Y43 = -1/ (0.18 + 0.55j)
= -0.5375 + 1.6423j
Y44 = (1/Z34) + (1/Z41) + (1/Z42) + Y34 + Y41 + Y42
= 1/ (0.18 + 0.55j) + 1/ (0.1 + 0.359j) + 1/ (0.25 + 0.1j) + 0.089j + 0.006j + 0.015j
Step4:
Modified Y-matrix
Adding the transformation and the generator admittance to the respective diagonal elements Y-bus matrix is
formed as follows.
HAND CALCULATION
Step1: Read the line data.
Step2: Initialize y-bus data. y-bus[i] [j] =complex (0.0, 0.0).
Step3: Compute y-bus matrix.
Y11 = (1/Z12) + (1/Z14) + Yc12 + Yc14
= (1/ (0.1+0.4j)) + (1/ (0.1+0.25j)) + 0.0075j + 0.006j
= 1.3430 - 4.981j
Y22 = (1/Z23) + (1/Z24) + (1/Z21) + Yc23 + Yc24 + Y21
= (1/ (0.15+0.6j)) + (1/ (0.25+0.7j)) + (1/ (0.1+0.4j)) + 0.01j + 0.015j + 0.0075j
= 1.4329 - 5.1562j
Y33 = (1/Z34) + (1/Z32) + Y34 + Y32
= (1/ (0.018+0.55j)) + (1/ (0.15+0.6j)) + 0.009j + 0.01j
= 0.926 - 3.1919j
Y13 = Y31= 0.0006 + 0.00j
Y23 = -1/ (0.15 + 0.6j)
= -0.8922 + 1.5686j
Y14 = Y41 = -1/ (0.1 + 0.35j)
= 0.7547 + 2.64j
Y24 = Y42 = -1/ (0.2 + 0.7j)
= -0.425 + 1.2670j
Y34 = Y43 = -1/ (0.18 + 0.55j)
= -0.5375 + 1.6423j
Y44 = (1/Z34) + (1/Z41) + (1/Z42) + Y34 + Y41 + Y42
= 1/ (0.18 + 0.55j) + 1/ (0.1 + 0.359j) + 1/ (0.25 + 0.1j) + 0.089j + 0.006j + 0.015j
Step4:
Modified Y-matrix
Adding the transformation and the generator admittance to the respective diagonal elements Y-bus matrix is
formed as follows.
14EGCEE037 Page 40
Y11
mod
= Y11
old
+ 1/ (ZT4 + ZT1)
= 1.3430 - 4.981j + 1/ (0.25 + 0.08j)
= 1.3430 – 8.5524j
Y22
mod
= Y22
old
+1/(ZT2+ZT3)
=1.4329-8.0132j
Y33
mod
=Y33
old
+0.9296-3.1919j
Y44
mod
=Y44
old
+1.7447-5.5208j
The modified Y-bus matrix is as follows
1.3430-8.5524j -0.588+2.35j 0.00+0.00j -0.75+2.64j
Ybus = -0.58+2.35j 1.432-8.013j -0.39+1.58j -0.45+1.26j
0.00+0.00j -0.89+5.68j 0.929+3.19j -0.53+1.64j
-0.75+2.6145j -0.452+1.267j - 0.537+1.54j 1.744+5.250j
Step5:
0.0069+0.195j -0.008+0.011j -0.043+0.13j 0.012+0.138j
Zbus= -0.0084+0.11j 0.11+0.217j 0.0064+0.88j -0.007+0.159j
-0.0084+0.11j 0.006+0.1881j 0.1008+0.517j 0.282+0.263j
0.012+0.158j -0.007+0.159j 0.282+0.26j 0.060+3.55j
Step6:
Verify check matrix
Check matrix=[YBUS
MOD
]*[ZBUS]
1.00+4.16j 1.38+4.16e
-17
j 0.00-0.378e
-19
j 0.00+0.00j
2.77e
-16
+2.77e-17 1.00+2.77e
-17
5.511e
-17
-1.38e
-17
-1.1
e-16
-0.00j
2.22e
-16
+5.5e
-17
j 2.22e
-16
+4.16e
-17
j 1.00-5.5e
-17
j 0.00-5.5e
-17
j
2.22e
-16
+5.55e
-17
j 2.20e
-16
+0.00j 0.00+0.00j 1.00+0.00j
Y11
mod
= Y11
old
+ 1/ (ZT4 + ZT1)
= 1.3430 - 4.981j + 1/ (0.25 + 0.08j)
= 1.3430 – 8.5524j
Y22
mod
= Y22
old
+1/(ZT2+ZT3)
=1.4329-8.0132j
Y33
mod
=Y33
old
+0.9296-3.1919j
Y44
mod
=Y44
old
+1.7447-5.5208j
The modified Y-bus matrix is as follows
1.3430-8.5524j -0.588+2.35j 0.00+0.00j -0.75+2.64j
Ybus = -0.58+2.35j 1.432-8.013j -0.39+1.58j -0.45+1.26j
0.00+0.00j -0.89+5.68j 0.929+3.19j -0.53+1.64j
-0.75+2.6145j -0.452+1.267j - 0.537+1.54j 1.744+5.250j
Step5:
0.0069+0.195j -0.008+0.011j -0.043+0.13j 0.012+0.138j
Zbus= -0.0084+0.11j 0.11+0.217j 0.0064+0.88j -0.007+0.159j
-0.0084+0.11j 0.006+0.1881j 0.1008+0.517j 0.282+0.263j
0.012+0.158j -0.007+0.159j 0.282+0.26j 0.060+3.55j
Step6:
Verify check matrix
Check matrix=[YBUS
MOD
]*[ZBUS]
1.00+4.16j 1.38+4.16e
-17
j 0.00-0.378e
-19
j 0.00+0.00j
2.77e
-16
+2.77e-17 1.00+2.77e
-17
5.511e
-17
-1.38e
-17
-1.1
e-16
-0.00j
2.22e
-16
+5.5e
-17
j 2.22e
-16
+4.16e
-17
j 1.00-5.5e
-17
j 0.00-5.5e
-17
j
2.22e
-16
+5.55e
-17
j 2.20e
-16
+0.00j 0.00+0.00j 1.00+0.00j
14EGCEE037 Page 41
G1 G2
0.1+0.4j
0.15+0.6j
0.1+0.35j
0.25+0.7j
0.18+0.55j
T1 T2
1
4 3
2
Line
Number
Starting
Bus
Ending
Bus
Series Line
Impedance
Line Changing
Admittance
1 1 2 0.1+0.4j 0.15j
2 2 3 0.15+0.6j 0.02j
3 2 4 0.18+0.55j 0.018j
4 3 4 0.1+0.35j 0.012j
5 4 1 0.25+0.7j 0.03j
G1 G2
0.1+0.4j
0.15+0.6j
0.1+0.35j
0.25+0.7j
0.18+0.55j
T1 T2
1
4 3
2
Line
Number
Starting
Bus
Ending
Bus
Series Line
Impedance
Line Changing
Admittance
1 1 2 0.1+0.4j 0.15j
2 2 3 0.15+0.6j 0.02j
3 2 4 0.18+0.55j 0.018j
4 3 4 0.1+0.35j 0.012j
5 4 1 0.25+0.7j 0.03j
14EGCEE037 Page 42
START
READ No. OF BUSES; No. OF LINES
AND LINE DATAS
INITIALIZE THE Y-BUS MATRIX
Is
I = nl
Stop
Consider line 1 = 1
i = sb(1); I = eb(1)
Y(I,j) = Y(I,i)+Yseries(I)+0.5Yseries(I)
Y(I,j) = Y(I,1)+Yseries(i)+0.5Yseries(I)
Y(I,j) = Yseries(I)
Y(I,1) = Y(I,j)
I = I + 1
Point Y - Bus
Compute Z-bus by inverting
Modified Y-Bus
Print the Result
START
READ No. OF BUSES; No. OF LINES
AND LINE DATAS
INITIALIZE THE Y-BUS MATRIX
Is
I = nl
Stop
Consider line 1 = 1
i = sb(1); I = eb(1)
Y(I,j) = Y(I,i)+Yseries(I)+0.5Yseries(I)
Y(I,j) = Y(I,1)+Yseries(i)+0.5Yseries(I)
Y(I,j) = Yseries(I)
Y(I,1) = Y(I,j)
I = I + 1
Point Y - Bus
Compute Z-bus by inverting
Modified Y-Bus
Print the Result
14EGCEE037 Page 43
PROGRAM
clear all; %CLEARS THE CONTENT OF MEMORY OR WORKSPACE
clc %CLEARS THE SCREEN
% GETTING INPUT DATA FROM THE PROGRAMMER
disp(' THE INPUT DATAS ARE AS FOLLOWS ');
nb=input('the number of buses\n');
nl=input('the number of lines\n');
sb=input('the number of starting buses\n');
eb=input('the number of ending buses\n');
sli=input('enter the details of series line impedance\n');
lca=input('enter the detail of line impedance\n');
% EXECUTION STARTS HERE
for m=1:nb
for n=1:nb
if m==n||m~=n
if sli(m,n)==0
n=n+1;
else
y(m,n)=+sli(m,n)^-1+lca(m,n);
y(n,n)=+sli(m,n)^-1+lca(m,n);
y(m,n)=-sli(m,n)^-1;
y(n,m)=y(m,n);
end
end
end
end
% OUTPUT OF Y-BUS
disp(' THE RESULT OF Y BUS MATRIX IS ');
ybus=y
% OUTPUT OF Z-BUS
PROGRAM
clear all; %CLEARS THE CONTENT OF MEMORY OR WORKSPACE
clc %CLEARS THE SCREEN
% GETTING INPUT DATA FROM THE PROGRAMMER
disp(' THE INPUT DATAS ARE AS FOLLOWS ');
nb=input('the number of buses\n');
nl=input('the number of lines\n');
sb=input('the number of starting buses\n');
eb=input('the number of ending buses\n');
sli=input('enter the details of series line impedance\n');
lca=input('enter the detail of line impedance\n');
% EXECUTION STARTS HERE
for m=1:nb
for n=1:nb
if m==n||m~=n
if sli(m,n)==0
n=n+1;
else
y(m,n)=+sli(m,n)^-1+lca(m,n);
y(n,n)=+sli(m,n)^-1+lca(m,n);
y(m,n)=-sli(m,n)^-1;
y(n,m)=y(m,n);
end
end
end
end
% OUTPUT OF Y-BUS
disp(' THE RESULT OF Y BUS MATRIX IS ');
ybus=y
% OUTPUT OF Z-BUS
14EGCEE037 Page 44
zbus=ybus^-1;
disp(' THE RESULT CHECK MATRIX IS ');
chkmat=zbus*ybus
disp(' THE RESULT OF Z BUS MATRIX IS ');
zbus=zbus
Output:-
THE INPUT DATAS ARE AS FOLLOWS
the number of buses 3
the number of lines 3
the number of starting buses [1 2 3]
the number of ending buses [2 3 1]
enter the details of series line impedance
[0 0.1+0.3j -0.2-0.6j; -0.1-0.3j 0 0.15+0.5j; 0.2+0.6j -0.15-0.5j 0]
enter the detail of line impedance
[0 0.02j -0.028j; -0.02j 0 0.0125j; 0.028j -0.0125j 0]
THE RESULT OF Y BUS MATRIX IS
0.5000 - 1.4720i 1.0000 - 3.0000i -0.5000 + 1.5000i
YBUS = 1.0000 - 3.0000i -0.5505 + 1.8224i 0.5505 - 1.8349i
-0.5000 + 1.5000i 0.5505 - 1.8349i 0.5505 - 1.8224i
THE RESULT CHECK MATRIX IS
1.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i
chkmat = -0.0000 1.0000 + 0.0000i 0.0000 - 0.0000i
0.0000 0 1.0000 + 0.0000i
zbus=ybus^-1;
disp(' THE RESULT CHECK MATRIX IS ');
chkmat=zbus*ybus
disp(' THE RESULT OF Z BUS MATRIX IS ');
zbus=zbus
Output:-
THE INPUT DATAS ARE AS FOLLOWS
the number of buses 3
the number of lines 3
the number of starting buses [1 2 3]
the number of ending buses [2 3 1]
enter the details of series line impedance
[0 0.1+0.3j -0.2-0.6j; -0.1-0.3j 0 0.15+0.5j; 0.2+0.6j -0.15-0.5j 0]
enter the detail of line impedance
[0 0.02j -0.028j; -0.02j 0 0.0125j; 0.028j -0.0125j 0]
THE RESULT OF Y BUS MATRIX IS
0.5000 - 1.4720i 1.0000 - 3.0000i -0.5000 + 1.5000i
YBUS = 1.0000 - 3.0000i -0.5505 + 1.8224i 0.5505 - 1.8349i
-0.5000 + 1.5000i 0.5505 - 1.8349i 0.5505 - 1.8224i
THE RESULT CHECK MATRIX IS
1.0000 - 0.0000i 0.0000 + 0.0000i -0.0000 - 0.0000i
chkmat = -0.0000 1.0000 + 0.0000i 0.0000 - 0.0000i
0.0000 0 1.0000 + 0.0000i
14EGCEE037 Page 45
THE RESULT OF Z BUS MATRIX IS
0.0551 + 0.1533i 0.0627 + 0.1922i -0.0213 - 0.0647i
ZBUS = 0.0627 + 0.1922i -0.0043 - 0.0094i 0.0519 + 0.1705i
-0.0213 - 0.0647i 0.0519 + 0.1705i 0.0832 + 0.2772i
Result:- The program has been developed to compute the Z-bus matrix and thus result are
verified by multiplying the modified Y-bus matrix and Z-bus matrix to get unit matrix.
THE RESULT OF Z BUS MATRIX IS
0.0551 + 0.1533i 0.0627 + 0.1922i -0.0213 - 0.0647i
ZBUS = 0.0627 + 0.1922i -0.0043 - 0.0094i 0.0519 + 0.1705i
-0.0213 - 0.0647i 0.0519 + 0.1705i 0.0832 + 0.2772i
Result:- The program has been developed to compute the Z-bus matrix and thus result are
verified by multiplying the modified Y-bus matrix and Z-bus matrix to get unit matrix.
14EGCEE037 Page 46
EXPERIMENT-8
AIM :- To develop a program for formation of y-bus by the method of inspection
THEORY
Bus admittance matrix is often used in power system studies in most of power system studies it
is necessary to form Y-bus matrix of the system by considering certain power system parameter
depending upon the type of analysis.
For example in load flow analysis it is necessary to form Y-bus matrix without taking into
account the generator impedance and load impedance.
In short circuit analysis the generator transient reactance and transformer impedance taken in
account, in addition to line data.
Y-bus may be computed by inspection method only if there is no natural coupling between the
line shunt admittance are added to the diagonal elements corresponding to the admittance.
FORMULA USED
Yij=∑Yij for i=1 to n
Yij= -Yij= -1/Zij
Yij=Yij
Where Yij= Sum of admittance connected to bus
Yij= Negative admittance between buses
Y11=(1/Z12)+(1/Z13)+Yc12+Yc13
=(1/(0.1+0.3j))+(1/(0.2+0.6j))+0.014j+0.01j
=1.5-4.48j
Y22=(1/Z21)+(1/Z23)+Yc21+Yc23
=(1/0.1+0.3j)+(1/0.15+0.5j)+0.01j+0.025j
=1.55-4.81j
Buses of which there are connected the half diagonal elements are affected. The equivalent
circuit of changing transformer may be considering forming Y-bus matrix.The dimension of Y-
bus matrix is n X n where n is the number of buses in a power system network, each bus is
connected only to two other buses. So the Y-bus of the large network is high is properly not
evident in small systems.
EXPERIMENT-8
AIM :- To develop a program for formation of y-bus by the method of inspection
THEORY
Bus admittance matrix is often used in power system studies in most of power system studies it
is necessary to form Y-bus matrix of the system by considering certain power system parameter
depending upon the type of analysis.
For example in load flow analysis it is necessary to form Y-bus matrix without taking into
account the generator impedance and load impedance.
In short circuit analysis the generator transient reactance and transformer impedance taken in
account, in addition to line data.
Y-bus may be computed by inspection method only if there is no natural coupling between the
line shunt admittance are added to the diagonal elements corresponding to the admittance.
FORMULA USED
Yij=∑Yij for i=1 to n
Yij= -Yij= -1/Zij
Yij=Yij
Where Yij= Sum of admittance connected to bus
Yij= Negative admittance between buses
Y11=(1/Z12)+(1/Z13)+Yc12+Yc13
=(1/(0.1+0.3j))+(1/(0.2+0.6j))+0.014j+0.01j
=1.5-4.48j
Y22=(1/Z21)+(1/Z23)+Yc21+Yc23
=(1/0.1+0.3j)+(1/0.15+0.5j)+0.01j+0.025j
=1.55-4.81j
Buses of which there are connected the half diagonal elements are affected. The equivalent
circuit of changing transformer may be considering forming Y-bus matrix.The dimension of Y-
bus matrix is n X n where n is the number of buses in a power system network, each bus is
connected only to two other buses. So the Y-bus of the large network is high is properly not
evident in small systems.
14EGCEE037 Page 47
Generalized [Y-bus]
Each admittance Yii(i=1 to n) is called the self admittance or driving point admittance
terminating on the particular box.
Each half admittance term Yij (j,i = 1, 2-n, i=j)(I,1) and Yii = Yii.
Yji – Yij, where Yij is admittance connected between buses i and j.
Y33 = (1/Z21)+(1/Z23)+ Yc21 + Yc23
= (1/(0.2+0.6j))-(1/(0.15+0.5j))+0.014j + 0.0125j
= 1.05 – 33j
Y12 = -1/ (0.1+0.3j) = 1-3j
Y23 = -1/ (0.15+0.5j) = 0.55 – 1.83j
ALGORITHM
Step1: Initialize Y-bus matrix that is replaced by all zero Yji = -Yij = half diagonal element.
Step2: Compute Yij = ΣYij = diagonal element.
Line Number Starting Bus Ending Bus Series Line Impedance Line Changing Admittance
1 1 2 0.1+0.3j 0.02j
2 2 3 0.15+0.5j 0.0125j
3 3 1 0.2+0.6j 0.028j G1 G2
0.1+0.3j
0.15+0.5j
0.2+0.6j
T1 T2
0.02j
0.0125j
0.028j
1 2
3
Generalized [Y-bus]
Each admittance Yii(i=1 to n) is called the self admittance or driving point admittance
terminating on the particular box.
Each half admittance term Yij (j,i = 1, 2-n, i=j)(I,1) and Yii = Yii.
Yji – Yij, where Yij is admittance connected between buses i and j.
Y33 = (1/Z21)+(1/Z23)+ Yc21 + Yc23
= (1/(0.2+0.6j))-(1/(0.15+0.5j))+0.014j + 0.0125j
= 1.05 – 33j
Y12 = -1/ (0.1+0.3j) = 1-3j
Y23 = -1/ (0.15+0.5j) = 0.55 – 1.83j
ALGORITHM
Step1: Initialize Y-bus matrix that is replaced by all zero Yji = -Yij = half diagonal element.
Step2: Compute Yij = ΣYij = diagonal element.
Line Number Starting Bus Ending Bus Series Line Impedance Line Changing Admittance
1 1 2 0.1+0.3j 0.02j
2 2 3 0.15+0.5j 0.0125j
3 3 1 0.2+0.6j 0.028j G1 G2
0.1+0.3j
0.15+0.5j
0.2+0.6j
T1 T2
0.02j
0.0125j
0.028j
1 2
3
14EGCEE037 Page 48
START
READ No. OF BUSES; No. OF LINES
AND LINE DATAS
INITIALIZE THE Y-BUS MATRIX
Is
I = nl
STOP
Consider line 1 = 1
i = sb(1); I = eb(1)
Y(I,j) = Y(I,i)+Yseries(I)+0.5Yseries(I)
Y(I,j) = Y(I,1)+Yseries(i)+0.5Yseries(I)
Y(I,j) = Yseries(I)
Y(I,1) = Y(I,j)
I = I + 1
Point Y - Bus
START
READ No. OF BUSES; No. OF LINES
AND LINE DATAS
INITIALIZE THE Y-BUS MATRIX
Is
I = nl
STOP
Consider line 1 = 1
i = sb(1); I = eb(1)
Y(I,j) = Y(I,i)+Yseries(I)+0.5Yseries(I)
Y(I,j) = Y(I,1)+Yseries(i)+0.5Yseries(I)
Y(I,j) = Yseries(I)
Y(I,1) = Y(I,j)
I = I + 1
Point Y - Bus
14EGCEE037 Page 49
PROGRAM:-
clear all; %CLEARS THE CONTENT OF MEMORY OR WORKSPACE
clc %CLEARS THE SCREEN
% GETTING INPUT DATA FROM THE PROGRAMMER
disp(' THE INPUT DATAS ARE AS FOLLOWS ');
nb=input('the number of buses\n');
nl=input('the number of lines\n');
sb=input('the number of starting buses\n');
eb=input('the number of ending buses\n');
sli=input('enter the details of series line impedance\n');
lca=input('enter the details of line changing admittance\n');
% EXECUTION STARTS HERE
for m=1:nb
for n=1:nb
if m==n||m~=n
if sli(m,n)==0
n=n+1;
else
y(m,n)=+sli(m,n)^-1+lca(m,n);
y(n,n)=+sli(m,n)^-1+lca(m,n);
y(m,n)=-sli(m,n)^-1;
y(n,m)=y(m,n);
end
end
end
end
% OUTPUT OF Y-BUS
disp(' THE RESULT OF Y BUS MATRIX IS ');
ybus=y
PROGRAM:-
clear all; %CLEARS THE CONTENT OF MEMORY OR WORKSPACE
clc %CLEARS THE SCREEN
% GETTING INPUT DATA FROM THE PROGRAMMER
disp(' THE INPUT DATAS ARE AS FOLLOWS ');
nb=input('the number of buses\n');
nl=input('the number of lines\n');
sb=input('the number of starting buses\n');
eb=input('the number of ending buses\n');
sli=input('enter the details of series line impedance\n');
lca=input('enter the details of line changing admittance\n');
% EXECUTION STARTS HERE
for m=1:nb
for n=1:nb
if m==n||m~=n
if sli(m,n)==0
n=n+1;
else
y(m,n)=+sli(m,n)^-1+lca(m,n);
y(n,n)=+sli(m,n)^-1+lca(m,n);
y(m,n)=-sli(m,n)^-1;
y(n,m)=y(m,n);
end
end
end
end
% OUTPUT OF Y-BUS
disp(' THE RESULT OF Y BUS MATRIX IS ');
ybus=y
14EGCEE037 Page 50
Output:-
THE INPUT DATAS ARE AS FOLLOWS
the number of buses 3
the number of lines 3
the number of starting buses [1 2 3]
the number of ending buses [2 3 1]
enter the details of series line impedance
[0 0.1+0.3j -0.2-0.6j; -0.1-0.3j 0 0.15+0.5j; 0.2+0.6j -0.15-0.5j 0]
enter the details of line changing admittance
[0 0.02j -0.028j; -0.02j 0 0.0125j; 0.028j -0.0125j 0]
THE RESULT OF Y BUS MATRIX IS
0.5000 - 1.4720i 1.0000 - 3.0000i -0.5000 + 1.5000i
Ybus = 1.0000 - 3.0000i -0.5505 + 1.8224i 0.5505 - 1.8349i
-0.5000 + 1.5000i 0.5505 - 1.8349i 0.5505 - 1.8224i
Result:- The program for the Y-bus formation by the method of inspection is written and the line
Data are given and input and output is verified using hand calculation.
Output:-
THE INPUT DATAS ARE AS FOLLOWS
the number of buses 3
the number of lines 3
the number of starting buses [1 2 3]
the number of ending buses [2 3 1]
enter the details of series line impedance
[0 0.1+0.3j -0.2-0.6j; -0.1-0.3j 0 0.15+0.5j; 0.2+0.6j -0.15-0.5j 0]
enter the details of line changing admittance
[0 0.02j -0.028j; -0.02j 0 0.0125j; 0.028j -0.0125j 0]
THE RESULT OF Y BUS MATRIX IS
0.5000 - 1.4720i 1.0000 - 3.0000i -0.5000 + 1.5000i
Ybus = 1.0000 - 3.0000i -0.5505 + 1.8224i 0.5505 - 1.8349i
-0.5000 + 1.5000i 0.5505 - 1.8349i 0.5505 - 1.8224i
Result:- The program for the Y-bus formation by the method of inspection is written and the line
Data are given and input and output is verified using hand calculation.
Tags
Categories
BusinessDownload
Download Slideshow Get the original presentation fileQuick Actions
Statistics
- Views 813
- Slides 51
- Favorites 2
- Age 2828 days
Related Slideshows
1
DTI BPI Pivot Small Business - BUSINESS START UP PLAN
MeljunCortes
32 views
1
CATHOLIC EDUCATIONAL Corporate Responsibilities
MeljunCortes
35 views
11
Karin Schaupp – Evocation; lançamento: 2000
alfeuRIO
32 views
10
Pillars of Biblical Oneness in the Book of Acts
JanParon
29 views
31
7-10. STP + Branding and Product & Services Strategies.pptx
itsyash298
30 views
44
Business Legislation PPT - UNIT 1 jimllpkggg
slogeshk98
33 views