Solving electrical radial lines with Python

Electrical transmission systems are something we all take for granted. They work in a reliable manner ensuring high quality of service and as few minutes lost per year as possible.

Low voltage lines and most of medium voltage lines are radial lines, that is they are like a branch of a tree with a unique power supply location. Radial lines are relatively easy to work with, you can solve a radial line problem (ie you can get currents and voltages) by applying Boucherot’s theorem to each section of the line. This is a good news for a Python enthusiast as myself, since repetitive tasks lends themselves to be automated with programming.

Since Boucherot’s theorem uses the absolute values of electrical quantities, it is useful when you are working with AC lines and you would like to know the magnitude (rms) of currents and voltages while at the same time you do not really care about the phase differences. The rms values are used, for instance, when you need to choose the protection systems to install and what kind of electrical wires to use.

Suppose you are given the following three phase balanced radial line:

You can think of C1, C2 and C3 as industrial motors or any other kind of three phased balanced loads.

 

We can make the following assumptions:
- The system is 3 phase, symmetrical and balanced.
- The loads are R – L loads (most loads are Ohmic inductive loads).
- The lines capacity can be neglected (this assumption is reasonable for overhead medium power lines).
- The load at the end of the line is working at nominal voltage.

These assumptions are reasonable and they mostly hold since three phase systems are designed to be as symmetrical as possible.

Under these assumptions, the following python script applies Bocherot’s theorem to calculate voltage and current rms at each section of the line.

The script contains 3 classes one for each element:
1) Electric load class. This class models an electric RL load assuming that the load is operated at nominal voltage.
2) The transformer class models a real transformer.
3) The transmission line class models a typical short power line.

Once the classes have been written, a few lines are required for solving the problem. The results are shown in the plot below.

As expected the results show a voltage drop across all sections (typical for RL loads) and, lower current on the primary side of the transformer. This is to be expected due to the voltage/current relationship of the transformer. As a next step, we could proceed to choose the appropriate security sistems to install at the beginning of each line and the electrical wires. Maybe this could be the next step in this project.

The code I used is shown below. I tried to make it as modular as possible so that it could be applied to almost every radial line. You can add as many loads, transformers and transmission lines as you like and then compute the values for each section. Originally this is how I checked my assignments but I thought it could be an interesting mini-project to share. I feel like the code is pretty self explanatory but feel free to leave a comment if you are not sure on what the code is doing!

	import matplotlib.pyplot as plt
	import numpy as np
	plt.style.use("ggplot")
	#-------------------------------------------------------------------------------
	# Auxiliary functions
	# Complex power
	def complex_power(p, q):
	return np.sqrt(p**2 + q**2)
	# Current
	def current(s, V):
	return s/(np.sqrt(3)*V)
	# Voltage
	def voltage(s, I):
	return s/(np.sqrt(3)*I)
	# Power (3 phase)
	def power_3f(a, I):
	return 3 * a * I**2
	# Vectorized version of round
	round_ = np.vectorize(round)
	#-------------------------------------------------------------------------------
	# 3 phase electric load
	class electric_load(object):
	def __init__(self, pf, Vn, S = None, P = None):
	# Electric load
	#
	# Power factor pf
	# Nominal voltage Vn [V]
	# Apparent (or complex) power [VA]
	# Power [W]
	#
	if S is None and P is None:
	raise ValueError("Either S or P should be assigned a value...")
	elif S is None:
	self.S = P/pf
	self.P = P
	else:
	self.S = S
	self.P = S*pf
	self.pf = pf
	self.Q = self.P * np.tan(np.arccos(pf))
	self.Vn = Vn
	def current_draw(self):
	# Get the current draw of the electric load
	# Output: I
	#
	return self.S/(np.sqrt(3)*self.Vn)
	def get_electric_load(self):
	# Return the load characteristics
	#
	# Output: V,I,p,q
	#
	p_ = self.P
	q_ = self.Q
	V = self.Vn
	I = self.current_draw()
	return V, I, p_, q_
	def add_load_in_parallel(self, p, q, V):
	# Add a load in parallel in a given line.
	# Get V,I,p,q at the top of the transmission line
	#
	# Input: p,q,V at the bottom of the line before the load
	# Output: V,I,p,q at the top of the line
	#
	p_ = p + self.P
	q_ = q + self.Q
	s_ = complex_power(p_, q_)
	I = current(s_, V)
	return V, I, p_, q_
	#-------------------------------------------------------------------------------
	# Transformer class
	class transformer(object):
	def __init__(self, V1n, V2n, Sn, vcc_pc):
	# Transfomer
	#
	# V1n = nominal voltage at primary [V]
	# V2n = nominal voltage at secondary [V]
	# Sn = nominal apparent power [VA]
	# vcc_pc = short circuit voltage [%]
	#
	self.V1n = V1n
	self.V2n = V2n
	self.Sn = Sn
	self.k = V1n/V2n
	self.vcc_pc = vcc_pc
	self.Rcc = 0
	self.Xcc = (vcc_pc * V2n**2)/(100 * Sn)
	def primary_values(self, p, q, I2):
	# Get the values at the primary side of the transfomer
	#
	# Input: p,q,I at the secondary side of the transfomer
	# Output: V,I,p,q at the primary side of the transfomer
	#
	p_ = p + power_3f(self.Rcc, I2)
	q_ = q + power_3f(self.Xcc, I2)
	s_ = complex_power(p_, q_)
	V2 = voltage(s_, I2)
	V1 = self.k * V2
	I1 = I2 / self.k
	return V1, I1, p_, q_
	#-------------------------------------------------------------------------------
	# Transmission line class
	class transmission_line(object):
	def __init__(self, R, X, L):
	# Transmission line
	#
	# R = resistance [Ohm/km]
	# X = reactance [Ohm/km]
	# L = length [km]
	#
	self.L = L
	self.R = R * L
	self.X = X * L
	self.Z = np.sqrt(R**2 + X**2)
	#self.Z_cmp = R + 1j * X
	def same_current(self, p, q, I):
	# Get V,I,p,q at the top of the transmission line
	#
	# Input: p,q,I at the bottom of the line
	# Output: V,I,p,q at the top of the line
	#
	p_ = p + power_3f(self.R, I)
	q_ = q + power_3f(self.X, I)
	s_ = complex_power(p_, q_)
	V = voltage(s_, I)
	return V, I, p_, q_
	#-------------------------------------------------------------------------------
	# Example of application
	# Loads
	c1 = electric_load(pf=0.90, Vn=15000, P=1000000)
	c2 = electric_load(pf=0.90, Vn=15000, S=1000000)
	c3 = electric_load(pf=0.95, Vn=15000, S=2000000)
	# Transmission lines
	tl_2 = transmission_line(R=0.1, X=0.1, L=10)
	tl_1 = transmission_line(R=0.1, X=0.1, L=5)
	# Transformers
	T_1 = transformer(V1n=132000, V2n=15000, Sn=25000000, vcc_pc=12.5)
	# Solve transmission line
	Va, Ia, Pa, Qa = c1.get_electric_load()
	Vb, Ib, Pb, Qb = tl_2.same_current(Pa, Qa, Ia)
	Vb, Ib, Pb, Qb = c2.add_load_in_parallel(Pb, Qb, Vb)
	Vb, Ib, Pb, Qb = c3.add_load_in_parallel(Pb, Qb, Vb)
	Vc, Ic, Pc, Qc = tl_1.same_current(Pb, Qb, Ib)
	Vd, Id, Pd, Qd = T_1.primary_values(Pc, Qc, Ic)
	# Get ready for printing out results
	sections = ["a", "b", "c", "d"]
	V = round_([Va, Vb, Vc, Vd], 2)
	I = round_([Ia, Ib, Ic, Id], 2)
	dV_ = ["Vdc", "Vcb", "Vba"]
	dV = [Vd/132000 - Vc/15000, Vc/15000 - Vb/15000, Vb/15000 - Va/15000]
	dV = [i * 100 for i in dV]
	# Print results
	print("Section", "\tVoltage [kV]", "\tCurrent [A]")
	for i,k in enumerate(sections):
	print(str(k), "\t\t" + str(round(V[i]/1000, 2)), "\t\t" + str(I[i]))
	print("\nVoltage drop:")
	for i,k in enumerate(dV_):
	print(str(k), "\t" + str(round(dV[i], 2)) + "%")
	#-------------------------------------------------------------------------------
	# Plot
	plt.figure()
	plt.subplot(3,1,1)
	plt.plot(I)
	plt.title("Current [A]")
	plt.xticks(range(len(sections)), sections, size='small')
	plt.subplot(3,1,2)
	plt.plot(V/1000)
	plt.xticks(range(len(sections)), sections, size='small')
	plt.title("Voltage [kV]")
	plt.subplot(3,1,3)
	plt.plot(dV)
	plt.xticks(range(len(dV_)), dV_, size='small')
	plt.title("Voltage drop %")
	plt.show()
	################################################################################
	# OUTPUT
	################################################################################
	#
	#Section Voltage [kV] Current [A]
	#a 15.0 42.77
	#b 15.1 157.12
	#c 15.28 157.12
	#d 135.48 17.85
	#
	#Voltage drop:
	#Vdc 0.8%
	#Vcb 1.18%
	#Vba 0.66%
	#>>>

view raw radial_lines.py hosted with ❤ by GitHub