From daf176e918df60c4fd590fbd1ce72ac0aacf16e9 Mon Sep 17 00:00:00 2001 From: Folkert Kevelam Date: Sun, 9 Jun 2024 17:52:33 +0200 Subject: [PATCH] Python Modified Nodal Analysis test --- test.py | 329 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 329 insertions(+) create mode 100644 test.py diff --git a/test.py b/test.py new file mode 100644 index 0000000..0c45767 --- /dev/null +++ b/test.py @@ -0,0 +1,329 @@ +import matplotlib.pyplot as plt +import numpy as np +from numpy import linalg +from math import pi, sin + +def Create_MNA(mat_cons): + Y = np.zeros(mat_cons[0]) + B = np.zeros(mat_cons[1]) + C = np.zeros(mat_cons[2]) + D = np.zeros(mat_cons[3]) + + return [Y,B,C,D] + +class Voltage_Supply: + def __init__(self, func, a, b): + self.node_a = a - 1 + self.node_b = b - 1 + self.v = 0 + self.func = func + self.branch_current_node = 0 + self.vp = 0 + self.ip = 0 + + def set_voltages_currents(self, v, c): + if self.node_b > 0 and self.node_a > 0: + self.vp = v[self.node_a] - v[self.node_b] + elif self.node_b <= 0 and self.node_a > 0: + self.vp = v[self.node_a] + elif self.node_b > 0 and self.node_b <= 0: + self.vp = v[self.node_b] + + if self.needs_branches: + self.ip = c[self.branch_current_node] + + def set_time(self, t): + self.v = self.func(t) + + def needs_branches(self): + return True + + def set_branch_current_node(self, a): + self.branch_current_node = a + + def get_matrix(self, mat_cons): + Y,B,C,D = Create_MNA(mat_cons) + + if self.node_b >= 0: + B[self.node_b][self.branch_current_node] = 1 + C[self.branch_current_node][self.node_b] = 1 + if self.node_a >= 0: + B[self.node_a][self.branch_current_node] = -1 + C[self.branch_current_node][self.node_a] = -1 + + M = np.block([[Y,B],[C,D]]) + RHS_Top = np.zeros((mat_cons[0][0], 1)) + RHS_Bot = np.zeros((mat_cons[2][0], 1)) + RHS_Bot[self.branch_current_node] = self.v + + RHS = np.block([[RHS_Top], [RHS_Bot]]) + + return [M, RHS] + +class Resistance: + def __init__(self, func, a, b): + self.node_a = a - 1 + self.node_b = b - 1 + self.r = 0 + self.func = func + self.vp = 0 + self.ip = 0 + + def set_voltages_currents(self, v, c): + if self.node_b > 0 and self.node_a > 0: + self.vp = v[self.node_a] - v[self.node_b] + elif self.node_b <= 0 and self.node_a > 0: + self.vp = v[self.node_a] + elif self.node_b > 0 and self.node_b <= 0: + self.vp = v[self.node_b] + + def set_time(self, t): + self.r = self.func(t) + + def needs_branches(self): + return False + + def get_matrix(self, mat_cons): + Y,B,C,D = Create_MNA(mat_cons) + + a = self.node_a + b = self.node_b + + if b >= 0: + Y[b,b] = self.r + if a >= 0: + Y[a,a] = self.r + + if a >= 0 and b >= 0: + Y[a,b] = -self.r + Y[b,a] = -self.r + + M = np.block([[Y,B], [C,D]]) + RHS = np.zeros((mat_cons[0][0] + mat_cons[2][0], 1)) + + return [M, RHS] + +class Inductor: + def __init__(self, func, a, b): + self.node_a = a - 1 + self.node_b = b - 1 + self.l = 0 + self.func = func + self.previous_t = 0.0 + self.current_t = 0.0 + self.vp = 0 + self.ip = 0 + + def set_voltages_currents(self, v, c): + if self.node_b > 0 and self.node_a > 0: + self.vp = v[self.node_a] - v[self.node_b] + elif self.node_b <= 0 and self.node_a > 0: + self.vp = v[self.node_a] + elif self.node_b > 0 and self.node_b <= 0: + self.vp = -v[self.node_b] + + if self.needs_branches: + self.ip = c[self.branch_current_node] + + def set_time(self, t): + self.previous_t = self.current_t + self.current_t = t + self.l = self.func(t) + + def needs_branches(self): + return True + + def set_branch_current_node(self, a): + self.branch_current_node = a + + def get_matrix(self, mat_cons): + Y,B,C,D = Create_MNA(mat_cons) + + a = self.node_a + b = self.node_b + h = self.current_t - self.previous_t + lh = self.l / h + + if self.node_b >= 0: + B[self.node_b][self.branch_current_node] = -1 + C[self.branch_current_node][self.node_b] = -1 + if self.node_a >= 0: + B[self.node_a][self.branch_current_node] = 1 + C[self.branch_current_node][self.node_a] = 1 + + D[self.branch_current_node][self.branch_current_node] = -lh + + M = np.block([[Y,B],[C,D]]) + RHS_Top = np.zeros((mat_cons[0][0], 1)) + RHS_Bot = np.zeros((mat_cons[2][0], 1)) + RHS_Bot[self.branch_current_node] = -lh * self.ip + + RHS = np.block([[RHS_Top], [RHS_Bot]]) + + return [M, RHS] + +class Capacitor: + def __init__(self, func, a, b): + self.node_a = a - 1 + self.node_b = b - 1 + self.c = 0 + self.func = func + self.previous_t = 0.0 + self.current_t = 0.0 + self.vp = 0 + self.ip = 0 + + def set_voltages_currents(self, v, c): + if self.node_b > 0 and self.node_a > 0: + self.vp = v[self.node_a] - v[self.node_b] + elif self.node_b <= 0 and self.node_a > 0: + self.vp = v[self.node_a] + elif self.node_b > 0 and self.node_b <= 0: + self.vp = -v[self.node_b] + + def set_time(self, t): + self.previous_t = self.current_t + self.current_t = t + self.c = self.func(t) + + def needs_branches(self): + return False + + def set_branch_current_node(self, a): + self.branch_current_node = a + + def get_matrix(self, mat_cons): + Y,B,C,D = Create_MNA(mat_cons) + + a = self.node_a + b = self.node_b + + h = self.current_t - self.previous_t + ch = self.c / h + + if b >= 0: + Y[b,b] = ch + if a >= 0: + Y[a,a] = ch + + if a >= 0 and b >= 0: + Y[a,b] = -ch + Y[b,a] = -ch + + M = np.block([[Y,B],[C,D]]) + RHS_Top = np.zeros((mat_cons[0][0], 1)) + RHS_Bot = np.zeros((mat_cons[2][0], 1)) + + if a > 0: + RHS_Top[a] = ch * self.vp + if b > 0: + RHS_Top[b] = - ch * self.vp + + RHS = np.block([[RHS_Top], [RHS_Bot]]) + + return [M, RHS] + +def run( netlist, t0, t1, points): + counter = 0 + nodes = list() + for j in netlist: + if j.needs_branches(): + j.set_branch_current_node(counter) + counter += 1 + nodes.append(j.node_a) + nodes.append(j.node_b) + + max_node = max(nodes) + 1 + + Y = (max_node,max_node) + B = (max_node,counter) + C = (counter,max_node) + D = (counter,counter) + + M = (max_node+counter, max_node+counter) + RHS = (max_node+counter, 1) + + mat_cons = [Y,B,C,D] + + t_space = np.linspace(t0,t1,points)[1:] + + answer_array = list() + answer_array.append(list()) + + for i in range(0,M[0]): + answer_array.append(list()) + + for t in range(0, len(t_space)): + mat = np.zeros(M) + rhs = np.zeros(RHS) + for i in netlist: + i.set_time(t_space[t]) + m, r = i.get_matrix(mat_cons) + mat += m + rhs += r + + ans = linalg.lstsq(mat,rhs)[0] + + voltages = ans[:max_node] + currents = ans[max_node:] + + voltages = [v[0] for v in voltages] + currents = [c[0] for c in currents] + + for j in netlist: + j.set_voltages_currents(voltages, currents) + + answer_array[0].append(t_space[t]) + for i in range(0, len(ans)): + answer_array[i+1].append(ans[i][0]) + + return answer_array + +def R_func(t): + return 1/1000.0 + +def V1_func(t): + #if t > 0.00001 * 0.5: + # return 5 + #else: + # return 0 + v = 0 + for i in range(1, 10, 2): + v += (10.0/i)*sin(2*pi*i*1000000.0*t) + return v + +def L_func(t): + return 1.0e-9 + +def C_func(t): + return 1000.0e-12 + +def C1_func(t): + return 1.0e-12 + +def V2_func(t): + if t == 0.0: + return -2.5 + else: + return 0 + +netlist = [ + Voltage_Supply(V1_func, 0, 1), + #Inductor(L_func, 1, 2), + Capacitor(C1_func, 1, 2), + Resistance(R_func, 1, 2), + Capacitor(C_func, 2, 0), + Resistance(R_func, 2, 0) + #Voltage_Supply(V2_func, 0, 3) +] + + +ans = run(netlist, 0.0, 0.00001, 100000) + +plt.plot(ans[0],ans[1], 'r') +plt.plot(ans[0],ans[2], 'g') +#plt.plot(t, outs[3]) +plt.show() + +plt.plot(ans[0], ans[3], 'r') +plt.show()