import numpy as np python…..
# System parameters
n_buses = 14
base_MVA = 1.0 # Base MVA
tolerance = 1e-6 # Convergence tolerance
max_iter = 50 # Maximum iterations
# Build Ybus
Ybus = [Link]((n_buses, n_buses), dtype=complex)
# Line data (a, b, z = R + jX total ohm)
line_data = [
(0, 1, 0.0297 + 0.016335j), # Line 1: LV1-LV2
(0, 4, 0.0396 + 0.02178j), # Line 2: LV1-LV5
(1, 4, 0.0297 + 0.016335j), # Line 3: LV2-LV5
(1, 3, 0.07992 + 0.04356j), # Line 4: LV2-LV4
(3, 4, 0.07992 + 0.04356j), # Line 5: LV4-LV5
(1, 2, 0.07992 + 0.04356j), # Line 6: LV2-LV3
(2, 3, 0.0198 + 0.01089j), # Line 7: LV3-LV4
(6, 8, 0.788 + 0.2336j), # Line 8: MV7-MV9
(5, 10, 2.364 + 0.7008j), # Line 9: MV6-MV11
(5, 11, 2.364 + 0.7008j), # Line 10: MV6-MV12
(5, 12, 1.182 + 0.3504j), # Line 11: MV6-MV13
(9, 10, 2.364 + 0.7008j), # Line 12: MV10-MV11
(12, 13, 1.182 + 0.3504j), # Line 13: MV13-MV14
(8, 13, 0.788 + 0.2336j), # Line 14: MV9-MV14
]
for a, b, z in line_data:
y=1/z
Ybus[a, a] += y
Ybus[b, b] += y
Ybus[a, b] -= y
Ybus[b, a] -= y
# Transformers (assume nominal ratio, add series admittance)
# T1: MV6 (5) - LV5 (4), 1.5 MVA, z_tr_pu = 0.03 + j0.03 on own base
z_t1 = (0.03 + 0.03j) * (base_MVA / 1.5)
y_t1 = 1 / z_t1
Ybus[5, 5] += y_t1
Ybus[4, 4] += y_t1
Ybus[5, 4] -= y_t1
Ybus[4, 5] -= y_t1
# T2: MV14 (13) - LV4 (3), same as T1
y_t2 = y_t1
Ybus[13, 13] += y_t2
Ybus[3, 3] += y_t2
Ybus[13, 3] -= y_t2
Ybus[3, 13] -= y_t2
# TG: Assume between MV7 (6) - MV8 (7), 3.5 MVA, z_tr_pu = 0.015 + j0.015 on own base
z_tg = (0.015 + 0.015j) * (base_MVA / 3.5)
y_tg = 1 / z_tg
Ybus[6, 6] += y_tg
Ybus[7, 7] += y_tg
Ybus[6, 7] -= y_tg
Ybus[7, 6] -= y_tg
# Bus types: 0 = Slack (MV12 index 11), 1 = PV (MV8 index 7), 2 = PQ
bus_types = [Link]([2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 0, 2, 2])
slack_bus = 11
# Specified net P and Q (pu = MW/MVar on base 1 MVA)
P_spec = [Link]([0.04266, 0, -0.06472, 0, 0, 0, 0, 0.69, -0.3273, -0.5724, -0.29028, 1.8102, -0.5862, -
0.11961])
Q_spec = [Link]([0.03045, 0, -0.04011, 0, 0, 0, 0, 0.45, -0.3823, -0.4272, -0.21771, 1.665, -0.4398, -
0.0897])
# Initial flat start
V = [Link](n_buses)
delta = [Link](n_buses)
# Or use table values for faster convergence (optional)
# V = [Link]([0.955, 0.931, 0.930, 0.953, 0.951, 0.966, 0.971, 0.975, 0.966, 0.94, 0.953, 0.957, 0.974,
0.967])
# delta = [Link]([-29.76, -30.76, -31.26, -31.5, -31.25, -30.31, -30.72, -30.84, -30.67, -29.81, -30.01, -
30.11, -30.66, -30.59]) * [Link] / 180
# Set fixed V for PV and slack
V[bus_types == 0] = 1.0 # Slack V
V[bus_types == 1] = 1.0 # PV V
delta[slack_bus] = 0.0 # Slack angle
# Power calculation
def calculate_power(V, delta, Ybus):
V_complex = V * [Link](1j * delta)
I = Ybus @ V_complex
S = V_complex * [Link](I)
return [Link](S), [Link](S)
# Mismatches
def calculate_mismatches(P_calc, Q_calc):
return P_spec - P_calc, Q_spec - Q_calc
# Jacobian
def build_jacobian(V, delta, P_calc, Q_calc, Ybus, bus_types, slack_bus):
n = len(V)
delta_idx = [i for i in range(n) if i != slack_bus]
v_idx = [i for i in range(n) if bus_types[i] == 2]
size_p = len(delta_idx)
size_q = len(v_idx)
J = [Link]((size_p + size_q, size_p + size_q))
for a in range(size_p):
i = delta_idx[a]
for b in range(size_p):
k = delta_idx[b]
G = Ybus[i,k].real
B = Ybus[i,k].imag
theta_ik = delta[i] - delta[k]
if i == k:
J[a, b] = -Q_calc[i] - V[i]**2 * B
else:
J[a, b] = -V[i] * V[k] * (G * [Link](theta_ik) - B * [Link](theta_ik))
for a in range(size_p):
i = delta_idx[a]
for b in range(size_q):
k = v_idx[b]
G = Ybus[i,k].real
B = Ybus[i,k].imag
theta_ik = delta[i] - delta[k]
if i == k:
J[a, size_p + b] = P_calc[i] / V[i] + V[i] * G
else:
J[a, size_p + b] = V[i] * (G * [Link](theta_ik) + B * [Link](theta_ik))
for a in range(size_q):
i = v_idx[a]
for b in range(size_p):
k = delta_idx[b]
G = Ybus[i,k].real
B = Ybus[i,k].imag
theta_ik = delta[i] - delta[k]
if i == k:
J[size_p + a, b] = P_calc[i] - V[i]**2 * G
else:
J[size_p + a, b] = V[i] * V[k] * (G * [Link](theta_ik) + B * [Link](theta_ik))
for a in range(size_q):
i = v_idx[a]
for b in range(size_q):
k = v_idx[b]
G = Ybus[i,k].real
B = Ybus[i,k].imag
theta_ik = delta[i] - delta[k]
if i == k:
J[size_p + a, size_p + b] = Q_calc[i] / V[i] - V[i] * B
else:
J[size_p + a, size_p + b] = V[i] * (G * [Link](theta_ik) - B * [Link](theta_ik))
return J, delta_idx, v_idx
# Newton-Raphson method
iter = 0
converged = False
while not converged and iter < max_iter:
P_calc, Q_calc = calculate_power(V, delta, Ybus)
dP, dQ = calculate_mismatches(P_calc, Q_calc)
mismatch = [Link]((dP[delta_idx], dQ[v_idx]))
if [Link]([Link](mismatch)) < tolerance:
converged = True
break
J, delta_idx, v_idx = build_jacobian(V, delta, P_calc, Q_calc, Ybus, bus_types, slack_bus)
dx = [Link](J, mismatch)
delta[delta_idx] += dx[:len(delta_idx)]
V[v_idx] += dx[len(delta_idx):]
iter += 1
# Results
if converged:
print("Load flow converged in {} iterations.".format(iter))
else:
print("Load flow did not converge!")
print("Bus Voltages (pu): \n", V)
print("Bus Angles (deg): \n", delta * 180 / [Link])
print("Active Power (pu): \n", P_calc)
print("Reactive Power (pu): \n", Q_calc)