直流電気回路の電流の計算

これなに

昔作った電気回路の電流を計算するC#のコードをPythonにしてみました。
キルヒホッフの法則を使っています。

Pythonのコード

import networkx as nx
import pandas as pd
from more_itertools import pairwise

class CycleInfo:
    """サイクルの情報"""
    def __init__(self, g, i, cycle):
        """g.edgesのcidsにインデックスを追加"""
        self.cycle = cycle
        self.volt = self.resist = self.cur = 0
        for e in pairwise(cycle + [cycle[0]]):
            if dc := g.edges.get(e):
                coe = 1
            else:
                coe, dc = -1, g.edges[e[1], e[0]]
            self.volt += coe * dc['volt']
            self.resist += dc['resist']
            dc['cids'].append((coe, i))

def calc_current(df):
    """電流計算。結果はcur列に追加"""
    g = nx.DiGraph()
    for row in df.itertuples():
        g.add_edge(row.node1, row.node2, volt=row.volt,
                   resist=row.resist, cids=[])

    cycles = nx.cycle_basis(g.to_undirected())
    cinfos = [CycleInfo(g, i, cycle) for i, cycle in enumerate(cycles)]

    vv = [[0] * len(cinfos) for _ in cinfos]  # 影響量
    from itertools import combinations
    for dc in g.edges.values():
        for (coe1, cid1), (coe2, cid2) in combinations(dc['cids'], 2):
            vv[cid1][cid2] += -coe1 * coe2 * dc['resist']
            vv[cid2][cid1] += -coe1 * coe2 * dc['resist']

    modified = True
    while modified:
        modified = False
        for i, cinfo in enumerate(cinfos):
            if not cinfo.resist:
                if not cinfo.volt:
                    continue
                raise ValueError('回路ショート')
            cur = cinfo.volt / cinfo.resist
            dif = cur - cinfo.cur
            if abs(dif) < 1e-8:
                continue
            modified  = True
            cinfo.cur = cur
            for j, cinfo2 in enumerate(cinfos):
                cinfo2.volt += dif * vv[i][j]
    df['cur'] = [sum(cinfos[i].cur for _, i in dc['cids'])
                 for _, dc in g.edges.items()]
    return df

サンプル1

circ1.csv

node1,node2,volt,resist
0,1,0,5
1,2,0,0
2,3,1.5,0
3,0,0,0

実行結果（cur列）

df = pd.read_csv('circ1.csv')
calc_current(df)

	node1	node2	volt	resist	cur
0	0	1	0	5	0.3
1	1	2	0	0	0.3
2	2	3	1.5	0	0.3
3	3	0	0	0	0.3

サンプル2

circ2.csv

node1,node2,volt,resist
0,1,0,20
0,2,0,10
1,2,0,20
1,3,0,50
2,3,0,40
3,4,0,0
4,5,30,0.98
5,0,0,0

実行結果（cur列）

df = pd.read_csv('circ2.csv')
calc_current(df)

	node1	node2	volt	resist	cur
0	0	1	0	20	0.372554
1	0	2	0	10	0.509811
2	1	2	0	20	-0.0588243
3	1	3	0	50	0.431378
4	2	3	0	40	0.568635
5	3	4	0	0	1.00001
6	4	5	30	0.98	1.00001
7	5	0	0	0	1.00001

追記

サイクルを使わず、また、連立方程式を解く方法。

from itertools import islice

import networkx as nx
import numpy as np


def calc_flow(g: nx.DiGraph) -> None:
    """キルヒホッフの法則で電流を計算。結果はflowに追加"""
    num_nodes = g.number_of_nodes()
    num_edges = g.number_of_edges()
    node_idx = {node: i for i, node in enumerate(g.nodes)}
    size = num_nodes + num_edges
    a = np.zeros((size, size))
    b = np.zeros(size)
    a[0, 0] = 1
    for i, node in enumerate(islice(g.nodes, 1, num_nodes), 1):
        for j, (u, v) in enumerate(g.edges):
            if u == node:
                a[i, num_nodes + j] = 1
            if v == node:
                a[i, num_nodes + j] = -1

    # V_u - V_v + V_src = I * R  =>  V_u - V_v - I * R = -V_src
    for j, (u, v, data) in enumerate(g.edges(data=True)):
        r = data.get("resist", 0)
        v_src = data.get("volt", 0)
        row = num_nodes + j
        a[row, node_idx[u]] = 1
        a[row, node_idx[v]] = -1
        a[row, num_nodes + j] = -r
        b[row] = -v_src

    x = np.linalg.solve(a, b)
    for j, (u, v) in enumerate(g.edges):
        g[u][v]["flow"] = x[num_nodes + j]


if __name__ == "__main__":
    g = nx.DiGraph()
    g.add_edge("0", "1", volt=0, resist=20)
    g.add_edge("0", "2", volt=0, resist=10)
    g.add_edge("1", "2", volt=0, resist=20)
    g.add_edge("1", "3", volt=0, resist=50)
    g.add_edge("2", "3", volt=0, resist=40)
    g.add_edge("3", "4", volt=0, resist=0)
    g.add_edge("4", "5", volt=30, resist=0.98)
    g.add_edge("5", "0", volt=0, resist=0)
    calc_flow(g)
    for u, v, data in g.edges(data=True):
        print(f"Edge ({u}, {v}): flow={data.get('flow', 'N/A'):.4f} A")