ElecSolver formalizes electric systems as linear problems, suitable for both temporal and frequency-domain studies. It focuses on constructing the linear system representation, leaving the actual numerical resolution to the user.
This repository is not a general-purpose electrical system solver. Instead, it acts as a bridge between:
- The graph-based description of an electric network
- The corresponding sparse linear system to solve
Its main goal is to provide a friendly Python interface for simulating analog electric systems. While suitable for small circuit simulations, its strength lies in scalability—handling millions of nodes and components, provided that you possess sufficient computational resources.
Note
Non-linear components are not supported. You must manage event detection and system updates yourself.
For now this package is distributed on pypi and can be installed using pip and conda/mamba
pip install ElecSolver
or
conda install elecsolver
For solving the linear systems we advise using MUMPS through python-mumps on linux that can be installed via conda
conda install python-mumps
This class handles frequency-domain analysis of linear electric systems.
- Supports tension and intensity sources
- Models inductive and resistive mutuals
- Detects and couples multiple subsystems
- Accepts arbitrary complex impedances and mutuals
- Constructs sparse linear systems (COO format)
Tip
Some solvers do not support complex-valued systems. Use the utility function cast_complex_system_in_real_system in utils.py to convert an n-dimensional complex system into a 2n-dimensional real system.
We would like to study the following system:
this can simply be defined in the following manner (We took R=1, L=1 and M=2):
import numpy as np
from scipy.sparse.linalg import spsolve
from ElecSolver import FrequencySystemBuilder
# Complex and sparse impedance matrix
# notice coil impedence between points 0 and 2, and coil impedence between 3 and 4
impedence_coords = np.array([[0,0,1,3],[1,2,2,4]], dtype=int)
impedence_data = np.array([1, 1j, 1, 1j], dtype=complex)
# Mutual inductance or coupling
# The indexes here are the impedence indexes in impedence_data
# The coupling is inductive
mutuals_coords = np.array([[1],[3]], dtype=int)
mutuals_data = np.array([2.j], dtype=complex)
electric_sys = FrequencySystemBuilder(
impedence_coords,
impedence_data,
mutuals_coords,
mutuals_data
)
# Set ground
# 2 values because one for each subsystem
electric_sys.set_ground(0, 3)
# Building system
electric_sys.build_system()
electric_sys.build_second_member_intensity(intensity=10, input_node=2, output_node=0)
# Get and solve the system
sys, b = electric_sys.get_system()
sol = spsolve(sys.tocsr(), b)
frequencial_response = electric_sys.build_intensity_and_voltage_from_vector(sol)
## We see a tension appearing on the lonely coil (between node 3 and 4)
print(frequencial_response.potentials[3]-frequencial_response.potentials[4])We want to add components in parallel with existing components for instance inserting a resistor in parallel with the first inductance (between nodes 0 and 2)
In python, simply add the resistance to the list of impedence in the very first lines of the script:
import numpy as np
from scipy.sparse.linalg import spsolve
from ElecSolver import FrequencySystemBuilder
# We add an additionnal resistance between 0 and 2
impedence_coords = np.array([[0,0,1,3,0],[1,2,2,4,2]], dtype=int)
impedence_data = np.array([1, 1j,1, 1j,1], dtype=complex)
# No need to change the couplings since indexes of the coils did not change
mutuals_coords = np.array([[1],[3]], dtype=int)
mutuals_data = np.array([2.j], dtype=complex)This class models time-dependent systems using resistors, capacitors, coils, and mutuals.
- Supports tension and intensity sources
- Models inductive and resistive mutuals
- Detects and couples multiple subsystems
- Accepts 3 dipole types: resistances, capacities and coils
- Constructs sparse linear systems (COO format)
We would like to study the following system:
with R=1, L=0.1, C=2 this gives:
import numpy as np
from scipy.sparse.linalg import spsolve
from ElecSolver import TemporalSystemBuilder
## Defining resistances
res_coords = np.array([[0,2],[1,3]],dtype=int)
res_data = np.array([1,1],dtype=float)
## Defining coils
coil_coords = np.array([[1,0],[3,2]],dtype=int)
coil_data = np.array([0.1,0.1],dtype=float)
## Defining capacities
capa_coords = np.array([[1,3],[2,0]],dtype=int)
capa_data = np.array([2,2],dtype=float)
## Defining empty mutuals here
mutuals_coords=np.array([[],[]],dtype=int)
mutuals_data = np.array([],dtype=float)
res_mutuals_coords=np.array([[],[]],dtype=int)
res_mutuals_data = np.array([],dtype=float)
## initializing system
elec_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
## Seting ground at point 0
elec_sys.set_ground(0)
## Build second member
elec_sys.build_system()
elec_sys.build_second_member_intensity(10,1,0)
# getting initial condition system
S_i,b = elec_sys.get_init_system()
# initial condition
sol = spsolve(S_i.tocsr(),b)
# get system (S1 is real part, S2 derivative part)
S1,S2,rhs = elec_sys.get_system()
## Solving using euler implicit scheme
dt=0.08
vals_res1 = []
vals_res2 = []
for i in range(50):
temporal_response = elec_sys.build_intensity_and_voltage_from_vector(sol)
vals_res1.append(temporal_response.intensities_res[1])
vals_res2.append(temporal_response.intensities_res[0])
## implicit euler time iterations
sol = spsolve(S2+dt*S1,b*dt+S2@sol)
import matplotlib.pyplot as plt
plt.xlabel("Time")
plt.ylabel("Intensity")
plt.plot(vals_res1,label="intensity res 1")
plt.plot(vals_res2,label="intensity res 2")
plt.legend()
plt.savefig("intensities_res.png")This outputs the following graph that displays the intensity passing through the resistances
- For small or moderately sized systems, the built-in
scipy.sparse.linalg.spsolveis effective. - For large-scale temporal problems, consider using MUMPS (via
pyMUMPS). MUMPS is more efficient when only the second member (b) changes during time-stepping.
Tip
See example tests.test_temporal_system in the tests on how to use pyMUMPS for solving the resulting system efficiently.
This repository can be used as is in order to model the mass flow or thermal flux in respectively Hydraulic networks or Thermal networks where a difference of pressure or a difference of temperature can be assimilated to a tension source. Since electric potentials are always computed relatively to the ground node you might need to rescale the resulting potentials:
We are considering the following hydraulic problem:
Taking R=1 this gives
import numpy as np
from scipy.sparse.linalg import spsolve
from ElecSolver import TemporalSystemBuilder
## Defining resistances
R = 1
res_coords = np.array([[0,2,1,0,1,3],[1,3,3,2,2,0]],dtype=int)
res_data = R*np.array([2,3,1,1,1,1],dtype=float)
## Here we are not using coils, capacities or mutuals we defined them as empty
## Defining 0 coil
coil_coords = np.array([[],[]],dtype=int)
coil_data = np.array([],dtype=float)
## Defining 0 capacity
capa_coords = np.array([[],[]],dtype=int)
capa_data = np.array([],dtype=float)
## Defining no mutual
mutuals_coords=np.array([[],[]],dtype=int)
mutuals_data = np.array([],dtype=float)
res_mutuals_coords=np.array([[],[]],dtype=int)
res_mutuals_data = np.array([],dtype=float)
## initializing system
hydraulic_sys = TemporalSystemBuilder(coil_coords,coil_data,res_coords,res_data,capa_coords,capa_data,mutuals_coords,mutuals_data,res_mutuals_coords,res_mutuals_data)
## Seting ground at point 0
hydraulic_sys.set_ground(0)
## Build second member
hydraulic_sys.build_system()
## enforcing a pressure delta of 10 Pa
hydraulic_sys.build_second_member_tension(10,1,0)
# get system (S1 is real part, S2 derivative part)
# the problem is only resitive thus S2 =0
S1,S2,rhs = hydraulic_sys.get_system()
sol = spsolve(S1.tocsr(),rhs)
solution = hydraulic_sys.build_intensity_and_voltage_from_vector(sol)
# After you computed the solution of the system
pressure_input=10000
pressure_node=0
# Rescaling the potential to the new reference
potentials = solution.potentials - solution.potentials[pressure_node] + pressure_input
print("Pressures in the system:", potentials)
## get the flux passing through the system
print("Debit through the system",solution.intensities_sources[0])A new class, named NetlistParser allows importing passive netlist and building a TemporalSystem instance. Solving the system can then be performed like any other example above.
"""
*test netlist for python solver square.net
Iin 0 1 PWL(0 0 0.000000001 10)
L0 1 3 0.1
L1 2 0 0.1
R1 1 0 1
R2 2 3 1
c2 1 2 2
c3 0 3 2
.tran 0 4 0 0.08
.end
"""
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt
from ElecSolver import NetlistParser
parser = NetlistParser("square.net")
parser.map_netlist()
node_zero = parser.node_map["0"]
node_one = parser.node_map["1"]
elec_sys=parser.generate_temporal_system()
## Seting ground at point 0
elec_sys.set_ground(node_zero)
## Build second member
elec_sys.build_system()
# Set 10 A injection entering in node 1 and exiting in node 0
elec_sys.build_second_member_intensity(10, node_one, node_zero)
# getting initial condition system
S_i,b = elec_sys.get_init_system()
# initial condition
sol = spsolve(S_i.tocsr(),b)
# get system (S1 is real part, S2 derivative part)
S1,S2,rhs = elec_sys.get_system()
## Solving using euler implicit scheme
dt=0.08
steps = 50
vals_res1 = []
vals_res2 = []
vals_L1 = []
voltage_src = []
R1_index = list(parser.resistors.keys()).index("R1")
R2_index = list(parser.resistors.keys()).index("R2")
L1_index = list(parser.inductors.keys()).index("L1")
for i in range(steps):
temporal_response = elec_sys.build_intensity_and_voltage_from_vector(sol)
vals_res1.append(temporal_response.intensities_res[R1_index])
vals_res2.append(temporal_response.intensities_res[R2_index])
vals_L1.append(temporal_response.intensities_coil[L1_index])
voltage_src.append(temporal_response.potentials[node_one]-temporal_response.potentials[node_zero])
## implicit euler time iterations
sol = spsolve(S2+dt*S1,b*dt+S2@sol)
plt.xlabel("Time")
plt.ylabel("Intensity")
plt.plot(arange(steps, dtype=float)*dt, vals_res1, label="intensity res 1")
plt.plot(arange(steps, dtype=float)*dt, vals_res2, label="intensity res 2")
plt.plot(arange(steps, dtype=float)*dt, vals_L1, label="intensity L1")
plt.plot(arange(steps, dtype=float)*dt, voltage_src, label="V(1-0)") # equal to I(R1)
plt.legend()
plt.savefig("intensities_res.png")