0% found this document useful (0 votes)
7 views16 pages

PS Exp

The document contains various programs for electrical engineering calculations, including transient stability analysis, load flow calculations using the Newton-Raphson method, and fault analysis in power systems. It includes user inputs for parameters such as the number of insulators, line voltage, and conductor diameters, and outputs results like string efficiency, bus voltages, and fault currents. The programs utilize mathematical computations and matrix operations to analyze power system behavior and performance.

Uploaded by

mbasilasil
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as DOCX, PDF, TXT or read online on Scribd
0% found this document useful (0 votes)
7 views16 pages

PS Exp

The document contains various programs for electrical engineering calculations, including transient stability analysis, load flow calculations using the Newton-Raphson method, and fault analysis in power systems. It includes user inputs for parameters such as the number of insulators, line voltage, and conductor diameters, and outputs results like string efficiency, bus voltages, and fault currents. The programs utilize mathematical computations and matrix operations to analyze power system behavior and performance.

Uploaded by

mbasilasil
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as DOCX, PDF, TXT or read online on Scribd
You are on page 1/ 16

Simulation diagram:

Output:
Transient stability analysis of SMIB system:

Output:
PROGRAM:
clc;
n=input('total number of insulators: ');
Vl=input('line voltage(in kV): ');
k=input('enter the value of k: ');
Vph=Vl/(sqrt(3));
if(n==2)
v1=Vph/((k+1)+1)
v2=(k+1)*v1
stringefficiency=(Vph/(n*v2))*100
end
if(n==3)
v1=Vph/((k+1)+((k*k)+(3*k)+1)+1)
v2=(k+1)*v1
v3=((k*k)+(3*k)+1)*v1
stringefficiency=(Vph/(n*v3))*100
end
if(n==4)
v1=Vph/((k+1)+((k*k)+(3*k)+1)+((k*k*k)+(5*k*k)+(6*k)+1)+1)
v2=(k+1)*v1
v3=((k*k)+(3*k)+1)*v1
v4=((k*k*k)+(5*k*k)+(6*k)+1)*v1
stringefficiency=(Vph/(n*v4))*100
end
OUTPUT:
total number of insulators: 3
line voltage(in kV): 33
enter the value of k: 0.11
v1 =
5.5191
v2 =
6.1262
v3 =
7.4072
stringefficiency =
85.7388

PROGRAM:
#include <stdio.h>
#include<math.h>
void main()
{
flot d,r, Dm,Ds,Dab,Dbc,Dca;
double L,C;

printf("enter the diameter of the conductor (in m): ");


scanf("%f",&d);
r=d/2;

printf("enter the distance between the conductors (in m)");


printf("\nDab = ");
scanf("%f",&Dab);
printf("\nDbc = ");
scanf("%f",&Dbc);
printf("\nDca = ");
scanf("%f",&Dca);

Dm = cbrt(Dab*Dbc*Dca);
printf("\nDm = %.3f",Dm);
Ds = 0.7788*r;
printf("\nDs = %.6f",Ds);
L = (2*0.0000001)*log(Dm/Ds);
C = (5.5631*10*0.000000000001)/(log(Dm/r));
printf("\nL = %.18f\n",L);
printf("\nC = %.18f",C);

}
OUTPUT:
enter the diameter of the conductor (in m): 0.012
enter the distance between the conductors (in m)
Dab = 2.4
Dbc = 2.4
Dca = 2.4
Dm = 2.3999
Ds = 0.004672
L = 0.000001248293110518
C = 0.000000000009285042
PROGRAM:
Z_series = 0.02 + 1i * 0.08;
Y_shunt = 1i * 0.02;
bus_data = [
1, 2.0, 1.0, NaN, NaN, 1.04 + 1i * 0;
2, 0.0, 0.0, 0.5, 1.0, NaN + 1i * NaN;
3, 1.5, 0.6, 0.0, NaN, 1.04 + 1i * NaN;
];
Q3_min = 0;
Q3_max = 1.5;
num_buses = size(bus_data, 1);
Y_bus = zeros(num_buses, num_buses);
lines = [
1, 2;
1, 3;
2, 3;
];
for k = 1:size(lines, 1)
i = lines(k, 1);
j = lines(k, 2);
y_series = 1 / Z_series;
y_shunt_half = Y_shunt / 2;
Y_bus(i, j) = Y_bus(i, j) - y_series;
Y_bus(j, i) = Y_bus(j, i) - y_series;
Y_bus(i, i) = Y_bus(i, i) + y_series + y_shunt_half;
Y_bus(j, j) = Y_bus(j, j) + y_series + y_shunt_half;
end
disp('Y_bus Matrix:');
disp(Y_bus);
V = ones(num_buses, 1);
delta = zeros(num_buses, 1);
slack_bus_idx = 1;
for i = 1:num_buses
if ~isnan(real(bus_data(i, 6)))
V(i) = real(bus_data(i, 6));
end
if ~isnan(imag(bus_data(i, 6)))
delta(i) = imag(bus_data(i, 6));
end
end
pv_bus_idx = find(~isnan(bus_data(:, 4)) & isnan(bus_data(:, 5)) & isnan(real(bus_data(:,
6))));
pq_bus_idx = find(isnan(bus_data(:, 4)) & isnan(bus_data(:, 5)) & isnan(real(bus_data(:,
6))));
slack_bus_idx = 1;
pv_bus_idx = 3;
pq_bus_idx = 2;
non_slack_buses = setdiff(1:num_buses, slack_bus_idx);
delta_unknown_idx = setdiff(1:num_buses, slack_bus_idx);
V_unknown_idx = pq_bus_idx;
Q_gen_PV = bus_data(pv_bus_idx, 5);
P_net = bus_data(:, 4) - bus_data(:, 2);
Q_net = bus_data(:, 5) - bus_data(:, 3);
P_net(pq_bus_idx) = bus_data(pq_bus_idx, 4) - bus_data(pq_bus_idx, 2);
Q_net(pq_bus_idx) = bus_data(pq_bus_idx, 5) - bus_data(pq_bus_idx, 3);
P_net(pv_bus_idx) = bus_data(pv_bus_idx, 4) - bus_data(pv_bus_idx, 2);
V(slack_bus_idx) = real(bus_data(slack_bus_idx, 6));
delta(slack_bus_idx) = imag(bus_data(slack_bus_idx, 6));
V(pv_bus_idx) = real(bus_data(pv_bus_idx, 6));
tolerance = 1e-6;
max_iter = 100;
iteration = 0;
fprintf('\nStarting Newton-Raphson Load Flow Calculation...\n');
while iteration < max_iter
iteration = iteration + 1;
P_calc = zeros(num_buses, 1);
Q_calc = zeros(num_buses, 1);
for i = 1:num_buses
for j = 1:num_buses
theta_ij = delta(i) - delta(j);
P_calc(i) = P_calc(i) + V(i) * V(j) * (real(Y_bus(i, j)) * cos(theta_ij) + imag(Y_bus(i,
j)) * sin(theta_ij));
Q_calc(i) = Q_calc(i) + V(i) * V(j) * (real(Y_bus(i, j)) * sin(theta_ij) - imag(Y_bus(i,
j)) * cos(theta_ij));
end
end
deltaP = P_net(non_slack_buses) - P_calc(non_slack_buses);
deltaQ = Q_net(pq_bus_idx) - Q_calc(pq_bus_idx);
Q_calc_pv = Q_calc(pv_bus_idx);
if Q_calc_pv < Q3_min
warning('Q_G3 at Bus 3 is below minimum limit. Bus 3 will be treated as PQ bus for
next iteration.');
Q_net(pv_bus_idx) = Q3_min - bus_data(pv_bus_idx, 3);
V_unknown_idx = unique([V_unknown_idx; pv_bus_idx]);
pv_bus_idx = setdiff(pv_bus_idx, pv_bus_idx);
deltaP = P_net(non_slack_buses) - P_calc(non_slack_buses);
deltaQ = Q_net(pq_bus_idx) - Q_calc(pq_bus_idx);
deltaQ_new = [deltaQ; Q_net(setdiff(non_slack_buses, slack_bus_idx)) -
Q_calc(setdiff(non_slack_buses, slack_bus_idx))];
if Q_calc_pv < Q3_min
Q_net_clamped_pv = Q3_min;
elseif Q_calc_pv > Q3_max
Q_net_clamped_pv = Q3_max;
else
Q_net_clamped_pv = Q_calc_pv;
end
end
if Q_calc_pv > Q3_max
warning('Q_G3 at Bus 3 is above maximum limit. Bus 3 will be treated as PQ bus for
next iteration.');
Q_net(pv_bus_idx) = Q3_max - bus_data(pv_bus_idx, 3);
V_unknown_idx = unique([V_unknown_idx; pv_bus_idx]);
pv_bus_idx = setdiff(pv_bus_idx, pv_bus_idx);
deltaP = P_net(non_slack_buses) - P_calc(non_slack_buses);
deltaQ = Q_net(pq_bus_idx) - Q_calc(pq_bus_idx);
end
mismatch = [deltaP; deltaQ];
if max(abs(mismatch)) < tolerance
fprintf('\nConverged in %d iterations.\n', iteration);
break;
end
J = zeros(length(non_slack_buses) + length(V_unknown_idx), length(non_slack_buses) +
length(V_unknown_idx));
angle_vars_idx = non_slack_buses;
volt_vars_idx = V_unknown_idx;
for i_idx = 1:length(angle_vars_idx)
i = angle_vars_idx(i_idx);
for j_idx = 1:length(angle_vars_idx)
j = angle_vars_idx(j_idx);
theta_ij = delta(i) - delta(j);
if i == j
H_ij = -Q_calc(i) - V(i)^2 * imag(Y_bus(i, i));
else
H_ij = V(i) * V(j) * (real(Y_bus(i, j)) * sin(theta_ij) - imag(Y_bus(i, j)) *
cos(theta_ij));
end
J(i_idx, j_idx) = H_ij;
end
end
for i_idx = 1:length(angle_vars_idx)
i = angle_vars_idx(i_idx);
for j_idx = 1:length(volt_vars_idx)
j = volt_vars_idx(j_idx);
theta_ij = delta(i) - delta(j);
if i == j
N_ij = (P_calc(i) + V(i)^2 * real(Y_bus(i, i))) / V(i);
else
N_ij = V(i) * (real(Y_bus(i, j)) * cos(theta_ij) + imag(Y_bus(i, j)) * sin(theta_ij));
end
J(i_idx, length(angle_vars_idx) + j_idx) = N_ij; % Corrected assignment for N_ij
end
end
for i_idx = 1:length(volt_vars_idx)
i = volt_vars_idx(i_idx);
for j_idx = 1:length(angle_vars_idx)
j = angle_vars_idx(j_idx);
theta_ij = delta(i) - delta(j);
if i == j
M_ij = P_calc(i) - V(i)^2 * real(Y_bus(i, i));
else
M_ij = -V(i) * V(j) * (real(Y_bus(i, j)) * cos(theta_ij) + imag(Y_bus(i, j)) *
sin(theta_ij));
end
J(length(angle_vars_idx) + i_idx, j_idx) = M_ij;
end
end
for i_idx = 1:length(volt_vars_idx)
i = volt_vars_idx(i_idx);
for j_idx = 1:length(volt_vars_idx)
j = volt_vars_idx(j_idx);
theta_ij = delta(i) - delta(j);
if i == j
L_ij = (Q_calc(i) + V(i)^2 * imag(Y_bus(i, i))) / V(i);
else
L_ij = V(i) * (real(Y_bus(i, j)) * sin(theta_ij) - imag(Y_bus(i, j)) * cos(theta_ij));
end
J(length(angle_vars_idx) + i_idx, length(angle_vars_idx) + j_idx) = L_ij;
end
end
if iteration >= max_iter
warning('Newton-Raphson did not converge within the maximum number of
iterations.');
end
P_final = zeros(num_buses, 1);
Q_final = zeros(num_buses, 1);

for i = 1:num_buses
for j = 1:num_buses
theta_ij = delta(i) - delta(j);
P_final(i) = P_final(i) + V(i) * V(j) * (real(Y_bus(i, j)) * cos(theta_ij) +
imag(Y_bus(i, j)) * sin(theta_ij));
Q_final(i) = Q_final(i) + V(i) * V(j) * (real(Y_bus(i, j)) * sin(theta_ij) -
imag(Y_bus(i, j)) * cos(theta_ij));
end
end
P_slack_gen = P_final(slack_bus_idx) + bus_data(slack_bus_idx, 2);
Q_slack_gen = Q_final(slack_bus_idx) + bus_data(slack_bus_idx, 3);
Q_pv_gen = Q_final(pv_bus_idx) + bus_data(pv_bus_idx, 3);
Q_pv_gen = max(Q3_min, min(Q_pv_gen, Q3_max));
del_x = J \ mismatch;
delta(angle_vars_idx) = delta(angle_vars_idx) + del_x(1:length(angle_vars_idx));
V(volt_vars_idx) = V(volt_vars_idx) + del_x(length(angle_vars_idx) + 1 : end);
fprintf('Iteration %d: Max Mismatch = %f\n', iteration, max(abs(mismatch)));
end
fprintf('\n--- Load Flow Results ---\n');
fprintf('%-5s %-10s %-10s %-10s %-10s %-10s %-10s\n', 'Bus', 'V (p.u.)', 'Angle (deg)',
'P_gen', 'Q_gen', 'P_load', 'Q_load');
fprintf('--------------------------------------------------------------------------------\n');
for i = 1:num_buses
P_gen_display = NaN;
Q_gen_display = NaN;
if i == slack_bus_idx
P_gen_display = P_slack_gen;
Q_gen_display = Q_slack_gen;
elseif i == pv_bus_idx % Bus 3
P_gen_display = bus_data(i, 4);
Q_gen_display = Q_pv_gen;
elseif i == pq_bus_idx % Bus 2
P_gen_display = bus_data(i, 4);
Q_gen_display = bus_data(i, 5);
end
fprintf('%-5d %-10.4f %-10.4f %-10.4f %-10.4f %-10.4f %-10.4f\n', ...
i, V(i), rad2deg(delta(i)), P_gen_display, Q_gen_display, bus_data(i, 2), bus_data(i, 3));
end
fprintf('\nFinal Bus Voltages (Magnitude and Angle):\n');
for i = 1:num_buses
fprintf('Bus %d: V = %.4f p.u., Angle = %.4f degrees\n', i, V(i), rad2deg(delta(i)));
end
OUTPUT:

PROGRAM:
clc;
v=1;zf=0.1j;
z10=0.25j; z12=0.125j; z13=0.15j; z20=0.25j; z23=0.25j;
y10=-1/z10; y12=-1/z12; y13=-1/z13; y20=-1/z20; y23=-1/z23;
y11=-(y10+y12+y13);
y21=y12; y22=-(y20+y21+y23);
y31=y13; y32=y23; y33=-(y31+y32);
ybus=[y11 y12 y13;
y21 y22 y23;
y31 y32 y33];
zbus=inv(ybus);
n=input('faulted bus:');
i=0;
while i<=3
i=i+1;
if i==n
z1=zbus(i,i);
end
end
disp('3 phase fault current')
I=v/(z1+zf)
z2=z1;
z10=0.4j; z12=0.30j; z13=0.35j; z20=0.1j; z23=0.7125j;
y10=-1/z10; y12=-1/z12; y13=-1/z13; y20=-1/z20; y23=-1/z23;
y11=-(y10+y12+y13);
y21=y12; y22=-(y20+y21+y23);
y31=y13; y32=y23; y33=-(y31+y32);
ybus=[y11 y12 y13;
y21 y22 y23;
y31 y32 y33];
zbus=inv(ybus);
i=0;
while i<=3
i=i+1;
if i==n
z0=zbus(i,i);
end
end
disp('line to ground fault current')
Ia1=v/(z1+z2+z0+(3*zf));
Ia=Ia1*3

OUTPUT:
faulted bus:3
3 phase fault current

I=
0.0000 - 3.1250i

line to ground fault current

Ia =
0.0000 - 2.7576i

You might also like