Power System Analysis of The IEEE 14-Bus Test
System Using PSAT and MATLAB
Submitted By
1. Pratik Roy (ID-104925203)
2. Pronay Kumar Chakrobarty (ID-104907898)
Procedure and Analysis:
1. Power flows using PSAT:
a) The provided MATLAB file entitled “d14bus_LF.m” was loaded into the PSAT program. Then
in general setting, PV reactive limits was checked for considering the generator reactive power
limits.
Fig-1: PSAT setting to enable reactive limits
Then power flow was run on PSAT and the obtained results are summarized on the attached excel
file named “Power Flow with Q Limit”.
2
Fig-2: Running power flow in PSAT considering reactive limits
b) Then in general setting, PV reactive limits was unchecked in order to obtain power flow without
considering the generator reactive power limits. Again, power flow was run on PSAT and the
obtained results are summarized on the attached excel file named “Power Flow without Q Limit”.
3
Fig-2: Running power flow in PSAT without considering reactive limits
For the power flow results without generator Q-limits, all PV generator buses were able to properly
control the bus voltage at the specified voltage magnitude. With the limits enforced, bus 8 reached
the upper limit of 0.24 p.u. of reactive power as seen in Figure 3. This caused the voltage magnitude
to decrease from the desired 1.09 p.u. to 1.0867 p.u. In addition, the remaining reactive power now
has to be compensated by the other generator buses. Lastly, it causes a change in phase and voltage
magnitude at buses that do not have a controlled voltage.
4
Figure-3: Comparison of power flow a) with considering & b) without considering Q limit
c) Shunt capacitors improve a network’s power factor, voltage stability, and reduce network loss.
A shunt capacitor is used to manipulate a bus voltage magnitude in the project. This is done by the
shunt capacitor’s ability to inject reactive power into the bus. To determine the amount of reactive
power required by the shunt capacitor, the subject bus was treated as a PV generator bus that
generated zero reactive power and had a fixed voltage of one per unit which was the desired
outcome. The Q-limits were set to very large values to ensure that a limit was not reached resulting
in an undesired voltage magnitude at the target bus. Another specification of the project is to
increase the P and Q load at all buses by 40% which was done by multiplying each load by 1.4.
The whole modification was done in the provided MATLAB file “d14bus_LF.m”. The overall
modification can be seen in figure -4.
5
Fig-4: Modification of MATLAB file for shunt capacitor at Bus 14
Then the modified MATLAB file was loaded to PSAT and power flow was run again. The results
are summarized in a excel file named “Power Flow with 40% increase and added shunt
capacitor for Bus 14”.
6
Fig-5: Running power flow in PSAT with 40% increase in P & Q load and added shunt capacitor
From power flow results, it can be shown that the shunt capacitor must inject 0.25675 p.u. or
25.675 MVAr to the bus in order to obtain a bus voltage magnitude of 1 p.u.
Fig-6: Results for Reactive Power Generation at Bus 14
7
Now, the actual value of capacitance for the shunt capacitor can be determined using the
calculations in below equations:
13.8kV
2
Vc 2
Xc 7.41
Q 25.675Mvar
1 1
C 357.97uF
2fX c 2*3.14159*60*7.41
The same process is repeated for bus 13 by changing the PV bus number as seen in Figure 6. The
modified MATLAB file “d14bus_LF.m” and the excel file for the power flow results named
“Power Flow with 40% increase and added shunt capacitor for Bus 13” were attached in the
report.
Fig-7: Modification of MATLAB file for shunt capacitor at Bus 13
Fig-8: Results for Reactive Power Generation at Bus 13
8
Fig-9: Running power flow in PSAT for current case
From power flow results (as per figure-8), it can be shown that the shunt capacitor must inject
0.20197 p.u. or 20.197 MVAr to the bus in order to obtain a bus voltage magnitude of 1 p.u. Again,
the actual value of capacitance for the shunt capacitor can be determined using the calculations in
below equations:
13.8kV
2
V2
Xc c 9.43
Q 20.197 Mvar
1 1
C 281.29uF
2fX c 2*3.14159*60*9.43
Comparing the two results in Figure 6 and Figure 8, it can be seen that the shunt capacitor
placement at bus 14 is the better option. It requires a greater value of capacitance and therefore
more production of reactive power by the capacitor bank, but the amount of active power generated
by the entire is smaller and thus, it has a better power factor. In addition, all values of voltage fall
between the allowable 5% fluctuation with the exception of bus one that has a controlled voltage
of 1.06 in both scenarios.
9
2. Short circuits using MATLAB:
a) The YBUS Matrices for all three sequence were made by hand calculation. The two YBUS
matrices were calculated for positive and zero sequence circuit. For negative sequence circuit
YBUS matrix is same as positive sequence circuit. The p.u. V, Pload and Qload were taken from the
power flow results of Q.1(a) to calculate the impedance, Z and later the admittance, Y. The other
admittances were calculated using the values of R, X and B from given table.
Fig-10: Taken values of Pload and Qload from power flow results of Q.1(a)
The impedance, Z and later the admittance, Y for load at each buses was calculated by below
equations:
|𝑉|2
𝑍𝑙𝑜𝑎𝑑 =
𝑆∗
Where,
𝑆 ∗ = 𝑃𝑙𝑜𝑎𝑑 − 𝑄𝑙𝑜𝑎𝑑
Now,
1
𝑌𝑙𝑜𝑎𝑑 =
𝑍𝑙𝑜𝑎𝑑
Also, for the easy calculation, the Y connection between bus 9, 4 and 8 was converted into Δ
connection and the total 13 buses were considered as the point 7 is not a real bus.
10
Fig-11: Hints to ease the calculation of YBUS Matrices
The overall hand calculation was attached as a scanned pdf copy named “YBUS calculation” with
the report.
b) Then in MATLAB those YBUS matrices were put and constructed a below code to determine
fault current at Bus 6 for 1-phase to ground solid faults.
MATLAB Code for determining 1-phase to ground solid faults
a=1*exp(j*120/180*pi);
A=[1 1 1;
1 a*a a;
1 a a*a];
%ZL = 0.911^2/(1*exp(-j*acos(0.9)));
Ybus1 = [6.0259-j*22.474 -5+j*15.26 0 0 -1.0259+j*4.235 0 0 0 0 0 0 0 0;
-5+j*15.26 9.73914-j*34.8074 -1.135+j*4.7819 -1.686+j*5.116 -1.70114+j*5.1939 0 0 0 0 0 0 0 0;
0 0 4.0629-j*14.3515 -1.9859+j*5.0688 0 0 0 0 0 0 0 0 0;
0 0 0 10.9909-j*38.1903 -6.84098+j*21.578 0 j*1.7806 j*4.6493 0 0 0 0 0;
-1.0259+j*4.235 -1.70114+j*5.1939 0 0 9.644-j*34.9428 j*3.9676 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 -1.955+j*4.094 -1.5259+j*3.1759 -3.0989+j*6.1027 0;
0 0 0 0 0 0 0 j*13.4956 j*3.38 0 0 0 0;
0 0 0 j*4.6493 0 0 j*3.38 5.455-j*21.4269 -3.902+j*10.365 0 0 0 -1.424+j*3.029;
0 0 0 0 0 0 0 -3.902+j*10.365 5.872-j*14.823 -1.88+j*4.40 0 0 0;
0 0 0 0 0 -1.955+j*4.094 0 0 -1.88+j*4.40 3.825-j*8.477 0 0 0;
0 0 0 0 0 0 0 0 0 0 4.0759-j*5.44 -2.489+j*2.252 0;
11
0 0 0 0 0 0 0 0 0 0 -2.489+j*2.252 6.8598-j*10.7269 -1.1369+j*2.3149;
0 0 0 0 0 0 0 -1.424+j*3.029 0 0 0 -1.1369+j*2.3149 2.7149-j*5.3959];
Ybus0 = [2.00837-j*6.3464 -1.6664+j*5.0877 0 0 -0.34197+j*1.4117 0 0 0 0 0 0 0 0;
-1.6664+j*5.0877 3.17376-j*9.86618 -0.3783+j*1.5939 -0.56201+j*1.70528 -0.05705+j*1.7313 0
0 0 0 0 0 0 0;
0 0 1.0403-j*3.1659 -0.662+j*1.6896 0 0 0 0 0 0 0 0 0;
0 0 0 3.504-j*16.884 -2.2803+j*7.1929 0 0 0 0 0 0 0 0;
-0.34197+j*1.4117 -0.05705+j*1.7313 0 0 3.18932-j*14.1578 0 0 0 0 0 0 0 0;
0 0 0 0 0 2.193-j*4.458 0 0 0 -0.6517+j*1.3647 -0.5087+j*1.0507 -1.03298+j*2.0343 0;
0 0 0 j*1.7806 0 0 -j*18.4506 0 0 0 0 0 0;
0 0 0 0 0 0 0 1.7754-j*4.4547 -1.3007+j*3.455 0 0 0 -0.9747+j*1.0077;
0 0 0 0 0 0 0 -1.3007+j*3.455 1.9276-j*4.9226 -0.6269+j*1.4676 0 0 0;
0 0 0 0 0 -0.6517+j*1.3647 0 0 -0.6269+j*1.4676 1.2786-j*2.8323 0 0 0;
0 0 0 0 0 -0.5087+j*1.0507 0 0 0 0 1.33837-j*1.80936 -0.82967+j*0.75066 0;
0 0 0 0 0 -1.03298+j*2.0343 0 0 0 0 -0.82967+j*0.75066 2.242-j*3.5567 -0.3789+j*0.7717;
0 0 0 0 0 0 0 -0.9747+j*1.0097 0 0 0 -0.3789+j*0.7717 0.8536-j*1.7814];
%%% Define pre-fault voltages for each sequence
V0_pf=zeros(13,1);
V2_pf=zeros(13,1);
V1_pf=[1;1;1;1;1;1;1;1;1;1;1;1;1];
%%% Define fault admittance matrix:
Yfabc=[999999999 0 0;0 0 0;0 0 0]; % SLG fault
%Yfabc=[999999999 0 0;0 999999999 0;0 0 999999999]; % 3-phase fault
%%% Calculate the fault admittance matrix in 0-1-2 frame
Yf012=(A^-1)*Yfabc*A;
%%% Use LU decomposition to solve for Z3 in each sequence
[L1 U1]=lu(Ybus1);
[L0 U0]=lu(Ybus0);
aux=eye(13);
X1= L1\aux(:,6);
X0= L0\aux(:,6);
Z6_1=U1\X1;
Z6_0=U0\X0;
Z6_2=Z6_1;
%%% Define the Z5_012 matrix (Z5_5 of each sequence in the diagonal)
Z6_012=[Z6_0(6) 0 0;0 Z6_1(6) 0;0 0 Z6_2(6)];
%%% Calculate fault current in 0-1-2 frame
If012=((Yf012*Z6_012+eye(3))^-1)*Yf012*[V0_pf(3);V1_pf(3);V2_pf(3)];
%%% Transform to a-b-c frame
Ifabc=A*If012;
Fault_current = [ abs(Ifabc), angle(Ifabc)*180/pi]
12
The output fault currents are:
Fault_current =
2.9449 178.5741
0.0000 17.9011
0.0000 -115.8976
Again, the MATLAB code was constructed to determine fault current at Bus 6 for 3-phase faults:
MATLAB Code for determining 3-phase faults
a=1*exp(j*120/180*pi);
A=[1 1 1;
1 a*a a;
1 a a*a];
%ZL = 0.911^2/(1*exp(-j*acos(0.9)));
Ybus1 = [6.0259-j*22.474 -5+j*15.26 0 0 -1.0259+j*4.235 0 0 0 0 0 0 0 0;
-5+j*15.26 9.73914-j*34.8074 -1.135+j*4.7819 -1.686+j*5.116 -1.70114+j*5.1939 0 0 0 0 0 0 0 0;
0 0 4.0629-j*14.3515 -1.9859+j*5.0688 0 0 0 0 0 0 0 0 0;
0 0 0 10.9909-j*38.1903 -6.84098+j*21.578 0 j*1.7806 j*4.6493 0 0 0 0 0;
-1.0259+j*4.235 -1.70114+j*5.1939 0 0 9.644-j*34.9428 j*3.9676 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 -1.955+j*4.094 -1.5259+j*3.1759 -3.0989+j*6.1027 0;
0 0 0 0 0 0 0 j*13.4956 j*3.38 0 0 0 0;
0 0 0 j*4.6493 0 0 j*3.38 5.455-j*21.4269 -3.902+j*10.365 0 0 0 -1.424+j*3.029;
0 0 0 0 0 0 0 -3.902+j*10.365 5.872-j*14.823 -1.88+j*4.40 0 0 0;
0 0 0 0 0 -1.955+j*4.094 0 0 -1.88+j*4.40 3.825-j*8.477 0 0 0;
0 0 0 0 0 0 0 0 0 0 4.0759-j*5.44 -2.489+j*2.252 0;
0 0 0 0 0 0 0 0 0 0 -2.489+j*2.252 6.8598-j*10.7269 -1.1369+j*2.3149;
0 0 0 0 0 0 0 -1.424+j*3.029 0 0 0 -1.1369+j*2.3149 2.7149-j*5.3959];
Ybus0 = [2.00837-j*6.3464 -1.6664+j*5.0877 0 0 -0.34197+j*1.4117 0 0 0 0 0 0 0 0;
-1.6664+j*5.0877 3.17376-j*9.86618 -0.3783+j*1.5939 -0.56201+j*1.70528 -0.05705+j*1.7313 0
0 0 0 0 0 0 0;
0 0 1.0403-j*3.1659 -0.662+j*1.6896 0 0 0 0 0 0 0 0 0;
0 0 0 3.504-j*16.884 -2.2803+j*7.1929 0 0 0 0 0 0 0 0;
-0.34197+j*1.4117 -0.05705+j*1.7313 0 0 3.18932-j*14.1578 0 0 0 0 0 0 0 0;
0 0 0 0 0 2.193-j*4.458 0 0 0 -0.6517+j*1.3647 -0.5087+j*1.0507 -1.03298+j*2.0343 0;
0 0 0 j*1.7806 0 0 -j*18.4506 0 0 0 0 0 0;
0 0 0 0 0 0 0 1.7754-j*4.4547 -1.3007+j*3.455 0 0 0 -0.9747+j*1.0077;
0 0 0 0 0 0 0 -1.3007+j*3.455 1.9276-j*4.9226 -0.6269+j*1.4676 0 0 0;
0 0 0 0 0 -0.6517+j*1.3647 0 0 -0.6269+j*1.4676 1.2786-j*2.8323 0 0 0;
0 0 0 0 0 -0.5087+j*1.0507 0 0 0 0 1.33837-j*1.80936 -0.82967+j*0.75066 0;
0 0 0 0 0 -1.03298+j*2.0343 0 0 0 0 -0.82967+j*0.75066 2.242-j*3.5567 -0.3789+j*0.7717;
13
0 0 0 0 0 0 0 -0.9747+j*1.0097 0 0 0 -0.3789+j*0.7717 0.8536-j*1.7814];
%%% Define pre-fault voltages for each sequence
V0_pf=zeros(13,1);
V2_pf=zeros(13,1);
V1_pf=[1;1;1;1;1;1;1;1;1;1;1;1;1];
%%% Define fault admittance matrix:
%Yfabc=[999999999 0 0;0 0 0;0 0 0]; % SLG fault
Yfabc=[999999999 0 0;0 999999999 0;0 0 999999999]; % 3-phase fault
%%% Calculate the fault admittance matrix in 0-1-2 frame
Yf012=(A^-1)*Yfabc*A;
%%% Use LU decomposition to solve for Z3 in each sequence
[L1 U1]=lu(Ybus1);
[L0 U0]=lu(Ybus0);
aux=eye(13);
X1= L1\aux(:,6);
X0= L0\aux(:,6);
Z6_1=U1\X1;
Z6_0=U0\X0;
Z6_2=Z6_1;
%%% Define the Z5_012 matrix (Z5_5 of each sequence in the diagonal)
Z6_012=[Z6_0(6) 0 0;0 Z6_1(6) 0;0 0 Z6_2(6)];
%%% Calculate fault current in 0-1-2 frame
If012=((Yf012*Z6_012+eye(3))^-1)*Yf012*[V0_pf(3);V1_pf(3);V2_pf(3)];
%%% Transform to a-b-c frame
Ifabc=A*If012;
Fault_current = [ abs(Ifabc), angle(Ifabc)*180/pi]
The output fault currents are:
Fault_current =
2.5051 116.8367
2.5051 -3.1633
2.5051 -123.1633
The MATLAB code “shortcircuit.m” was attached with the project report.
14
3. Voltage Stability Analysis using PSAT:
At first the MATLAB file “Base_14bus_VS.m” was loaded in PSAT and in general settings
below options were checked.
Fig-12: General settings for power flow
Then power flow was run and the results are summarized in a excel file named “VS Power Flow
Results”.
15
Fig-13: Running power flow in PSAT
Then CPF settings were checked and changed according to below figure -14 and then CPF was
run.
16
Fig-14: CPF Settings in PSAT
Fig-15: Running CPF in PSAT
17
Then plot the V_Bus 11 for PV curve.
Fig-16: PV Curve at bus 11 by running CPF
Then, the code has been modified as figure 11 for running CPF with line 2-4 tripped, and with
the consider of voltage limits at the generator buses of 0.9-1.1 p.u and 0.95-1.05 p.u at the load
buses. The figure 12 is the result of the modified CPF.
Fig-17: Modification of Code
18
For these case, same thing was again done to plot V_Bus 11. In this case, power flow results are
summarized in an excel file named “Different VS Power Flow Results” and modified MATLAB
file “Base_14bus_VS.m” was also attached with the report.
Fig-18: PV Curve at bus 11 by running CPF
From the power flow result, the Vmin = 0.993395 p.u. was taken and λmax = 1.1411 was calculated
at Vmin. Also the overall active load power and losses were taken from power flow result. Now,
the calculation was done by using below equations:
𝑇𝑇𝐶 = λ max∗ ∑(𝑃𝐿0 + 𝑙𝑜𝑠𝑠𝑒𝑠) = 1.1411 ∗ (2.59 + 0.157929) = 3.877 𝑝. 𝑢.
𝐴𝑇𝐶 = 𝑇𝑇𝐶 − 𝐸𝑇𝐶 − 𝑇𝑅𝑀 = (3.877 ∗ 100 𝑀𝑊) − 3.5 𝑀𝑊 − 5% ∗ (3.877 ∗ 100 𝑀𝑊)
= 364.85 𝑀𝑊
Since line 2-4 trip is the worst contingency, if it has been cut off from the system, the rest system
still could flow successfully. As the result of that, the given system operation condition could meet
the N-1 contingency criteria.
19
4. Angle Stability Analysis using PSAT (Bonus Question):
a) At first the MATLAB file “Base_14bus_lambda.m” was loaded in PSAT. Before that
[Link] and [Link] were activated in that program for dynamic loading.
Fig-19: Activation of [Link] and [Link]
Then in general settings of power flow, PV reactive limit were checked and power flow was run.
After that time domain analysis and eigen value analysis were done and eigen values are plotted
and report was made in a excel file named “Report Base Case” which is attached with the report.
Fig-20: Eigen value analysis for base case
20
Now, in contingency case everything was repeated after modifying the MATLAB code as per
below:
Fig-21: Modification of connection status in code
The found report of eigen value analysis was attached with the report in a excel file named “Report
Contingency”.
Fig-22: Eigen value analysis for contingency case
For both cases, by increasing the load in which both sides of imaginary part coincide in line have
to be noticed. From that eigen value at which a Hopf bifurcation appears can be determined.
b) For this case, connection status was again changed in to 0 to 1 and 8th element was changed to
100 from 200.
21
Fig-23: Modification in code
Then eigen value and time domain value analysis was done to present the eigen, delta and omega
value curves.
Fig-24: Eigen Value Graph
22
Fig-25: Delta vs. Time Graph
Fig-26: Omega vs. Time Graph
23
Those graphs are the indicators that the system becomes unstable due to the existence of a Hopf
bifurcation.
24