[PYTHON] Let's write a simple DC current solver

Introduction

This article is the 23rd day of Competitive Programming (2) Advent Calendar 2019. I wrote a DC current solver that calculates how many currents and combined resistances flow through each resistor in a circuit in which several electrical resistors are connected in series or in parallel. (Although it is rarely asked in the current competitive programming, it is okay because it is a graph and calculation)

DC current solver

For example, suppose you want to know the combined resistance of the following resistor circuit made by combining resistors with 1.0Ω resistors. Although it looks simple, it is very difficult to actually calculate it, so the program that calculates this is called the DC current solver in this article. As a premise

--Each resistance value is a known value --The power supply is connected somewhere --I want to calculate the DC current that flows constantly without a time element.

How to solve

By solving simultaneous equations formulated according to Kirchhoff's voltage law and current law, the potential of each node and the current of each resistor can be obtained. Therefore, it is necessary to define the variables appropriately.

Formulate the part in red as a variable. Since we end up with a linear equation, we make the minimum unknown element a variable. In high school physics classes and exams, current calculation for systems with many unknown variables is rarely asked, but if there are many variables, it is necessary to devise an implementation to determine how many linearly independent equations should be formulated. I will. By properly considering the connection from the power supply as a variable, the necessary and sufficient equation can be formulated as follows.

You will get as many expressions as there are defined variables. The point is that the formula obtained without connecting the power supply does not become linearly independent, and the linearly independent formula can be obtained no matter how many nodes connecting the power supply are increased (unless the power supplies are short-circuited). .. (Can be proved inductively (should))

Example

For example, in the case of such a simple series resistor

--Three formulas for the current flowing into the node --Two equations for voltage drop at each resistor --Two formulas for power connection

A total of 7 can be formulated as follows.

Implemented in Python

Since there is a linear solver that can be used easily, I will write it in Python. In writing the solver, the following policy

--The number of nodes is fixed at instantiation --Register the resistance definition in the instance with (from which node to which node, resistance value) --Define by (power supply name, supply potential) so that multiple power supplies can be used. -Register to the instance which node to connect the defined power supply to

• call solve()！

The source is on GitHub. It's almost like this

DCsolver.py

# coding: utf-8
import sys
import time
import numpy as np
from collections import deque

class DCsolver:

    def __init__(self, nodenum):
        #Solver
        self.linear_solver = np.linalg.solve

        #The number of nodes is fixed. Update the number of resistors and the number of bias connections
        self.n_node = nodenum
        self.n_res = 0
        self.n_bias_probe = 0

        #Resistance: Start point, end point, resistance value
        self.res_from = []
        self.res_to = []
        self.res_value = []

        #Power supply: Access by name for virtual connect
        self.bias_name = []
        self.bias_level = []
        self.bias_index = dict() # dict<string,int>

        # Bias supplied (-1: not biased)
        self.biased = [-1] * self.n_node
        #Determine at calculation which node each bias is probing to
        self.bias_from = None
        self.bias_to = None

        #Linear equation A*X=Solve V
        self._A = None
        self._V = None
        self._X = None 

        # result
        self.node_voltage = None
        self.node_current = None
        self.bias_total_current = None
        self.bias_current_per_node = None
        self.res_current = None

    def add_resistance(self, node_from, node_to, res_value):
        # node_from to node_resistance value res to to_Connect the resistance of value
        #Define the direction of the current in this direction
        assert res_value > 0 , "inhibit res_value <= 0"
        self.res_from.append(node_from)
        self.res_to.append(node_to)
        self.res_value.append(res_value)
        self.n_res += 1

    def define_bias(self, bias_name, bias_level=0.0):
        if bias_name in self.bias_index:
            idx = self.bias_index[bias_name]
            self.bias_level[idx] = bias_level
        else :
            idx = len(self.bias_name)
            self.bias_index[bias_name] = idx
            self.bias_name.append(bias_name)
            self.bias_level.append(bias_level)
    
    def supply_bias_to_node(self, bias_name, node_to):
        #Prevent multiple biases from accessing a node to reflect the latest settings
        assert bias_name in self.bias_index, \
            "{0} is not defined, please define before supply".format(bias_name)
        idx = self.bias_index[bias_name]
        #Warn if bias is already supplied.
        if self.biased[node_to] != -1:
            print("bias on node:{0} is changed: {1} --> {2} ".format(
                node_to, self.bias_name[self.bias_index[self.biased[node_to]]], bias_name
            ))
            self.biased[node_to] = idx
        else :
            self.biased[node_to] = idx
            self.n_bias_probe += 1

    def _create_matrix(self):
        # (A, V)To define
        nv = self.n_node
        nr = self.n_res
        nb = self.n_bias_probe

        #List and index the final bias condition settings
        self.bias_from = []
        self.bias_to = []
        for i in range(nv):
            if self.biased[i] != -1:
                self.bias_from.append(self.biased[i])
                self.bias_to.append(i)
        assert nb == len(self.bias_from)

        #Matrix size=Number of nodes+Number of resistors+Number of bias supply paths
        #Unknown variable
        #  [0...nv-1]:Node i potential
        #  [nv...nv+nr-1] :Resistor current
        #  [nv+nr...n-1] :Bias supply path current
        n = nv + nr + nb
        mat = np.zeros((n,n))
        vec = np.zeros(n)

        # Kirchhoff's Current Law (The sum of outgo and income of each node is 0)
        #i line([0...nv-1])Equation is the sum of the currents for node i
        #The current flowing through the resistor j is from[j]Out go to[j]Count as income
        for j in range(nr):
            mat[self.res_from[j]][nv + j] = 1
            mat[self.res_to[j]][nv + j] = -1

        # Kirchhoff's Voltage Law (The voltage drop of each resistor is a potential difference)
        #  nv+line j([nv...nv+nr-1])Equation is the sum of the currents with respect to the resistor j
        for j in range(nr):
            mat[nv + j][self.res_from[j]] = 1
            mat[nv + j][self.res_to[j]] = -1
            mat[nv + j][nv + j] = -self.res_value[j]
        
        #Bias definition formula
        #  bias_from[i]Potential is bias_level['bias']Is fixed to. (Added to KVL part)
        #  bias_to[i]Current from the power supply flows into(Added to KCL part)
        for j in range(len(self.bias_from)):
            mat[nv + nr + j][self.bias_to[j]] = 1
            vec[nv + nr + j] = self.bias_level[self.bias_from[j]]
            mat[self.bias_to[j]][nv + nr + j] = -1
        
        #Check if there is a node to which Bias is not connected (floating node is not allowed)
        self.check_connention()
        return mat, vec

    def check_connention(self):
        E = [[] for i in range(self.n_node)]
        for i in range(self.n_res):
            E[self.res_from[i]].append(self.res_to[i])
            E[self.res_to[i]].append(self.res_from[i])
        
        q = deque()
        vis = [False] * self.n_node
        for node in self.bias_to:
            q.append(node)
        while(len(q) > 0):
            now = q.popleft()
            vis[now] = True
            for nxt in E[now]:
                if not vis[nxt]:
                    q.append(nxt)
        
        floating_node = []
        for node in range(len(vis)):
            if not vis[node]:
                floating_node.append(node)
        if len(floating_node) != 0:
            print("Some floating node(s) exist:\nnode(s):\n{0}\n--> Aborted".format(floating_node))
            sys.exit(0)

            

    
    def solve(self):
        A, V = self._create_matrix()
        X = self.linear_solver(A, V)
        X = np.array(X)
        self._A = A
        self._V = V
        self._X = X

        self.node_voltage = X[0:self.n_node]
        self.res_current = X[self.n_node: self.n_node + self.n_res]        
        self.bias_current_per_node = dict()
        self.bias_total_current = dict()

        for bname in self.bias_name:
            self.bias_current_per_node[bname] = []
            self.bias_total_current[bname] = 0.0
        for i in range(self.n_bias_probe):
            bname = self.bias_name[self.bias_from[i]]
            self.bias_current_per_node[bname].append( (self.bias_to[i], X[self.n_node + self.n_res + i]) )
            self.bias_total_current[bname] += X[self.n_node + self.n_res + i]
        
        #Node current calculates outgo only ((Σ)|outgo|+Σ|income|)/Calculate with 2)
        self.node_current = np.zeros(self.n_node)
        for i in range(self.n_res):
            self.node_current[self.res_from[i]] += np.abs(self.res_current[i])
            self.node_current[self.res_to[i]] += np.abs(self.res_current[i])
        for bname in self.bias_current_per_node:
            for node, cur in self.bias_current_per_node[bname]:
                self.node_current[node] += np.abs(cur)
        self.node_current /= 2.0

    def print_result_summary(self, showCoef = False):
        if showCoef:
            print('A',self._A)
            print('V',self._V)
            print('X',self._X)
        print('node_voltage\n:{0}\n'.format(self.node_voltage))
        print('res_current\n:{0}\n'.format(self.res_current))
        print('node_cur\n:{0}\n'.format(self.node_current))
        print('bias_cur\n:{0}\n'.format(self.bias_total_current))
        print('bias_cur_per_node\n:{0}\n'.format(self.bias_current_per_node))

The point is that the supply_bias_to_node is checked so that different power supplies are not connected to the same node, and the check_connention is inserted at the timing of creating a matrix so that there are no unsupplied nodes (matrix). The linear solver just throws an error when is degenerate) Since the time complexity of preprocessing is $ O (number of nodes + number of resistors + number of power supply connections) $, the calculation itself should be almost rate-determining by the amount of calculation of the linear solver.

Actually use

I will actually use it. First, a simple example.

This is an example of connecting two resistors in series to three nodes. The combined resistance is the reciprocal of the current output from Vdd, but the Vdd current is 0.0333 ... A, which is certainly 30Ω.

Next, try connecting the second resistor in parallel. Two resistors are defined from node 1 to node 2. The combined resistance is 20Ω, but the Vdd current is 0.05A, which is certain.

A little complicated example

Intuitively, I don't know how many wraparound elements there are, but it seems that the Vdd current is 4Ω at 0.25A. Since V1, V3, and V4 have the same potential, they are virtually degenerate (1/20 + 1/10 + 1 / (5 + 5)) ^ -1.

Mesh resistance

Here are some potential distributions and current distributions when an appropriate potential is supplied to the nodes of the mesh using a solver. It is a mesh consisting of 50 x 50 nodes, and is connected by 1Ω between adjacent nodes.

• sample1

• sample2

• sample3

Although it is a calculation time, it seems that even with the same circuit structure, it differs greatly depending on how the power supply is connected and the number, and sample1 and sample3 could be calculated in less than 1 second, but sample2 took nearly 200 seconds. The number of variables is (50 × 50) + 2 * (50 * 49) + the number of power connections is 7400 or more, so I feel that it is still fast. (I used scipy's sparse matrix library) Since it is calculated by the direct method, it seems that it will be a little faster if it is calculated by the indirect method (iteration method) (not verified this time). I put the code etc. in here.

At the end

Have a fun electric circuit life in the future!