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