tmp2.py

import numpy as np


def run(input_file, step=0.1, iterations=1):
    n = 0  # Number of nodes.
    m = 0  # Number of voltage sources.
    components = []  # List of all components found in the netlist.
    node_res_dic = {}  # Dictionary of resistances where key is the node number and the resistance is the value.
    res_node_dic = {}  # Dictionary of resistances where key is the resistance and the node number is the value.

    # Read file.
    with open(input_file, "r") as lines:
        for line in lines:
            line = line.replace("\n", "")
            line = line.lower()

            # Get component data.
            component = line.split(" ")
            component_type = component[0]
            component_node_1 = int(component[1])
            component_node_2 = int(component[2])

            # Append the components for later use.
            components.append(component)

            # Calculate the n and m to know the size of the matrices.
            n = np.max((n, int(component_node_1), int(component_node_2)))

            if component_type == 'vsrc' or component_type == 'i':
                m += 1

    # Check the matrices sizes.
    print('N =', n, ', M =', m)

    # Create the G, B, C, D and A matrices.
    # Create the X and Z vectors.
    g = np.zeros((n, n))
    b = np.zeros((n, m))
    c = np.zeros((m, n))
    d = np.zeros((m, m))  # In case of not using inductance or capacitance it will remain a zero matrix.
    x = np.zeros((n + m, 1))
    z = np.zeros((n + m, 1))

    cur_m = 0
    cur_m_z = 0
    z_relations = []
    for i in range(n + m):
        z_relations.append([])

    for component in components:
        # Get component data.
        component_type = component[0]
        component_node_1 = int(component[1]) - 1
        component_node_2 = int(component[2]) - 1
        component_value = float(component[3])

        # Update G.
        if component_type == 'r':
            if component_node_1 == -1 and component_node_2 == -1:
                continue

            if component_node_1 == -1:
                g[component_node_2][component_node_2] += (1 / component_value)

            elif component_node_2 == -1:
                g[component_node_1][component_node_1] += (1 / component_value)
            else:
                g[component_node_2][component_node_2] += (1 / component_value)
                g[component_node_1][component_node_1] += (1 / component_value)
                g[component_node_2][component_node_1] -= (1 / component_value)
                g[component_node_1][component_node_2] -= (1 / component_value)

        if component_type == 'c':
            if component_node_1 == -1 and component_node_2 == -1:
                continue

            if component_node_1 == -1:
                g[component_node_2][component_node_2] += (component_value / step)
                z_relations[component_node_2].append((-component_value / step, component_node_2))

            elif component_node_2 == -1:
                g[component_node_1][component_node_1] += (component_value / step)
                z_relations[component_node_1].append((component_value / step, component_node_1))

            else:
                g[component_node_2][component_node_2] += (component_value / step)
                g[component_node_1][component_node_1] += (component_value / step)
                g[component_node_2][component_node_1] -= (component_value / step)
                g[component_node_1][component_node_2] -= (component_value / step)
                # append (c/h, src, dst)
                z_relations[component_node_1].append((component_value / step, component_node_1, component_node_2))
                z_relations[component_node_2].append((-component_value / step, component_node_1, component_node_2))

        # Update b.
        if component_type == 'vsrc':
            # Component node 1 is the +ve and component node 2 is the -ve.
            b[component_node_1] = 1
            b[component_node_2] = -1
            cur_m += 1
            z[n + cur_m_z] = component_value
            cur_m_z += 1
        elif component_type == 'i':
            b[component_node_1] = 1
            b[component_node_2] = -1
            d[cur_m][cur_m] = - component_value / step

            # append (-l/h, current of inductor index)
            z_relations[n + cur_m_z].append((component_value / step, cur_m_z))
            cur_m_z += 1
            cur_m += 1

        # Update z.
        if component_type == 'isrc':
            # Component node 1 is the +ve and component node 2 is the -ve.
            if component_node_1 > -1:
                z_relations[component_node_1].append((component_value))
                z[component_node_1] = component_value
            if component_node_2 > -1:
                z_relations[component_node_2].append((-component_value))
                z[component_node_2] = component_value

    # Get C.
    c = np.transpose(b)

    print('G =\n', g, '\n')
    print('B =\n', b, '\n')
    print('C =\n', c, '\n')
    print('D =\n', d, '\n')

    # Create A.
    a = np.zeros((n + m, n + m))
    a[0:n, 0:n] = g
    a[0:n, n:n + m] = b
    a[n:n + m, 0:n] = c
    a[n:n + m, n:n + m] = d
    print('A =\n', a, '\n')
    print('Z =\n', z, '\n')
    print('X =\n', x, '\n')
    print('ZR =\n', z_relations, '\n')

    # Loop over iterations
    for i in range(iterations):
        x = np.matmul(np.linalg.inv(a), z)
        for j in range(n + m):
            z[j] = 0
            for relation in z_relations[j]:
                if type(relation) is float:
                    z[j] = relation
                elif len(relation) == 1:
                    z[j] += relation[0]
                elif len(relation) == 3:
                    z[j] += relation[0] * (x[int(relation[1])] - x[int(relation[2])])
                elif len(relation) == 2:
                    z[j] += relation[0] * (x[n + int(relation[1])])
        print(x)


run("TestCases/2/input.txt", 0.1, 1)