mph- / lcapy

Lcapy is a Python package for symbolic linear circuit analysis and signal processing. It uses SymPy for symbolic mathematics.
GNU Lesser General Public License v2.1
239 stars 46 forks source link

The MNA A matrix is not invertible for time analysis #133

Closed dlee267 closed 2 months ago

dlee267 commented 3 months ago

Python 3.12.2, with lcapy 1.23

Problem

I may not see a voltage source that is short-circuited or all other concerns regarding this error, but I cannot figure out why this configuration is not solvable. Is there any reason why the matrix is not invertible, or is my formatting of my netlist giving the parser trouble? Any help will be appreciated.

Code to Run:

from lcapy import *

cct = Circuit('''
R_top0_3 t0_2 t0_3 0.02; right
R_top0_4 t0_3 t0_4 0.02; right
R_top1_3 t1_2 t1_3 0.02; right
R_top1_4 t1_3 t1_4 0.02; right
R_top2 t2 t2_0 10; right
R_top2_1 t2_0 t2_1 0.02; right
R_top2_2 t2_1 t2_2 0.02; right
R_top2_3 t2_2 t2_3 0.02; right
R_bottom0_2 b0_2 b1_2 0.02; down
R_bottom2_3 b2_3 b3_3 0.02; down
R_bottom3_3 b3_3 b4_3 0.02; down
R_bottom3 b4_3 b3 10; down
R_bottom0_4 b0_4 b1_4 0.02; down
R_middle0_2 t0_2 b0_2 16.291257858276367; rotate=-45
R_middle0_4 t0_4 b0_4 19.355159759521484; rotate=-45
R_middle1_2 t1_2 b1_2 17.680463790893555; rotate=-45
R_middle1_4 t1_4 b1_4 15.43282699584961; rotate=-45
R_middle2_3 t2_3 b2_3 14.204872131347656; rotate=-45
V n_r t2 5
R_ref n_r 0 51.0
W 0 b3
''')
print(cct.unconnected_nodes())
cct["R_ref"].v

Output:

[]

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File [~/.local/lib/python3.12/site-packages/lcapy/matrix.py:344](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/lcapy/matrix.py#line=343), in matrix_solve(M, b, method)
    343 try:
--> 344     sol_num, sol_den = M.to_DM().solve_den(b.to_DM())
    345     x = (sol_num.to_field() [/](http://localhost:8888/) sol_den).to_Matrix()

AttributeError: 'MutableDenseMatrix' object has no attribute 'to_DM'

During handling of the above exception, another exception occurred:

NonInvertibleMatrixError                  Traceback (most recent call last)
File [~/.local/lib/python3.12/site-packages/lcapy/mna.py:205](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/lcapy/mna.py#line=204), in MNA._solve(self)
    204 try:
--> 205     results = matrix_solve(self._A, self._Z, method=solver_method)
    206 except ValueError:

File [~/.local/lib/python3.12/site-packages/lcapy/matrix.py:348](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/lcapy/matrix.py#line=347), in matrix_solve(M, b, method)
    346 except:
    347     # Fallback
--> 348     Minv = matrix_inverse(M, method)
    349     x = Minv * b

File [~/.local/lib/python3.12/site-packages/lcapy/matrix.py:332](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/lcapy/matrix.py#line=331), in matrix_inverse(M, method)
    330             method = matrix_inverse_fallback_method
--> 332 return M.inv(method=method)

File [~/.local/lib/python3.12/site-packages/sympy/matrices/matrices.py:2179](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/sympy/matrices/matrices.py#line=2178), in MatrixBase.inv(self, method, iszerofunc, try_block_diag)
   2178 def inv(self, method=None, iszerofunc=_iszero, try_block_diag=False):
-> 2179     return _inv(self, method=method, iszerofunc=iszerofunc,
   2180             try_block_diag=try_block_diag)

File [~/.local/lib/python3.12/site-packages/sympy/matrices/inverse.py:463](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/sympy/matrices/inverse.py#line=462), in _inv(M, method, iszerofunc, try_block_diag)
    462 elif method == "ADJ":
--> 463     rv = M.inverse_ADJ(iszerofunc=iszerofunc)
    464 elif method == "CH":

File [~/.local/lib/python3.12/site-packages/sympy/matrices/matrices.py:2158](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/sympy/matrices/matrices.py#line=2157), in MatrixBase.inverse_ADJ(self, iszerofunc)
   2157 def inverse_ADJ(self, iszerofunc=_iszero):
-> 2158     return _inv_ADJ(self, iszerofunc=iszerofunc)

File [~/.local/lib/python3.12/site-packages/sympy/matrices/inverse.py:219](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/sympy/matrices/inverse.py#line=218), in _inv_ADJ(M, iszerofunc)
    207 """Calculates the inverse using the adjugate matrix and a determinant.
    208 
    209 See Also
   (...)
    216 inverse_LDL
    217 """
--> 219 d = _verify_invertible(M, iszerofunc=iszerofunc)
    221 return M.adjugate() [/](http://localhost:8888/) d

File [~/.local/lib/python3.12/site-packages/sympy/matrices/inverse.py:202](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/sympy/matrices/inverse.py#line=201), in _verify_invertible(M, iszerofunc)
    201 if zero:
--> 202     raise NonInvertibleMatrixError("Matrix det == 0; not invertible.")
    204 return d

NonInvertibleMatrixError: Matrix det == 0; not invertible.

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
Cell In[36], line 2
      1 print(cct.unconnected_nodes())
----> 2 cct["R_ref"].v

File [~/.local/lib/python3.12/site-packages/lcapy/mnacpts.py:693](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/lcapy/mnacpts.py#line=692), in Cpt.v(self)
    689 @property
    690 def v(self):
    691     """Time-domain voltage drop across component."""
--> 693     return self.cct.get_vd(self.nodes[0].name, self.nodes[1].name)

File [~/.local/lib/python3.12/site-packages/lcapy/netlist.py:324](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/lcapy/netlist.py#line=323), in Netlist.get_vd(self, Np, Nm)
    321 def get_vd(self, Np, Nm=None):
    322     """Time-domain voltage drop between nodes"""
--> 324     return self.get_Vd(Np, Nm).time()

File [~/.local/lib/python3.12/site-packages/lcapy/netlist.py:319](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/lcapy/netlist.py#line=318), in Netlist.get_Vd(self, Np, Nm, **kwargs)
    316 """Voltage drop between nodes (time-domain)"""
    318 Np, Nm = self._check_nodes(Np, Nm)
--> 319 return self._get_Vd(Np, Nm, **kwargs)

File [~/.local/lib/python3.12/site-packages/lcapy/netlist.py:310](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/lcapy/netlist.py#line=309), in Netlist._get_Vd(self, Np, Nm, nowarn)
    308 result = SuperpositionVoltage()
    309 for sub in subs.values():
--> 310     Vd = sub.get_Vd(Np, Nm)
    311     result.add(Vd)
    312 result = result.canonical()

File [~/.local/lib/python3.12/site-packages/lcapy/subnetlist.py:60](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/lcapy/subnetlist.py#line=59), in SubNetlist.get_Vd(self, Np, Nm)
     57 def get_Vd(self, Np, Nm=None):
     58     """Voltage drop between nodes"""
---> 60     return (self.mna.Vdict[Np] - self.mna.Vdict[Nm]).canonical()

File [~/.local/lib/python3.12/site-packages/lcapy/mna.py:340](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/lcapy/mna.py#line=339), in MNA.Vdict(self)
    335 @property
    336 def Vdict(self):
    337     """Return dictionary of transform domain node voltages indexed by node
    338     name"""
--> 340     self._solve()
    341     return self._Vdict

File [~/.local/lib/python3.12/site-packages/lcapy/mna.py:208](http://localhost:8888/home/segatrooper/.local/lib/python3.12/site-packages/lcapy/mna.py#line=207), in MNA._solve(self)
    206 except ValueError:
    207     message = self._failure_reasons()
--> 208     raise ValueError(message)
    210 results = results.subs(cct.context.symbols)
    212 # Handle capacitors at DC by assuming an infinite resistance
    213 # in parallel.

ValueError: The MNA A matrix is not invertible for time analysis: Not all nodes are connected.  Use cct.unconnected_nodes() to find them.

Edit

Since my post, I tried debugging with nodal analysis. Here are the results:

cct = Circuit('''
R_top0_3 t0_2 t0_3 0.02; right
R_top0_4 t0_3 t0_4 0.02; right
R_top1_3 t1_2 t1_3 0.02; right
R_top1_4 t1_3 t1_4 0.02; right
R_top2 t2 t2_0 10; right
R_top2_1 t2_0 t2_1 0.02; right
R_top2_2 t2_1 t2_2 0.02; right
R_top2_3 t2_2 t2_3 0.02; right
R_bottom0_2 b0_2 b1_2 0.02; down
R_bottom2_3 b2_3 b3_3 0.02; down
R_bottom3_3 b3_3 b4_3 0.02; down
R_bottom3 b4_3 b3 10; down
R_bottom0_4 b0_4 b1_4 0.02; down
R_middle0_2 t0_2 b0_2 16.291257858276367; rotate=-45
R_middle0_4 t0_4 b0_4 19.355159759521484; rotate=-45
R_middle1_2 t1_2 b1_2 17.680463790893555; rotate=-45
R_middle1_4 t1_4 b1_4 15.43282699584961; rotate=-45
R_middle2_3 t2_3 b2_3 14.204872131347656; rotate=-45
V n_r t2 5
R_ref n_r 0 51.0
W 0 b3
''')
l = cct.nodal_analysis()
l.nodal_equations()["t2"]; l.nodal_equations()["n_r"]

Output: image image

They seem to only be dependent on each other, at least according to the nodal analysis. However, the netlist shows that there are other dependencies that affect voltage of n_r and t2; am I wrong in this assessment?

mph- commented 3 months ago

There seems to be two independent circuits. You can see this with cct.cg.draw()

dlee267 commented 3 months ago

Thank you for the suggestion @mph ! I knew that I could draw the circuit on latex, but i didn't know you could draw the nodes programmatically! I made an algorithm that creates random circuits for a project, and it accidentally made two disconnected loops.

So, I added this code snippet I just made to only keep the loop that contains the voltage source:

main_graph = set(cct.super_nodes[0])
main_graph_has_changed = True
cct_other = Circuit(str(cct))
cct_main = Circuit()
while main_graph_has_changed:
    print(main_graph)
    cct_other_new = Circuit()
    graph_new = main_graph.copy()
    # print(main_graph)
    main_graph_has_changed = False
    for netlist_line in str(cct_other).split("\n"):
        # print(netlist_line)
        netlist_split = netlist_line.strip().split(" ")
        for i in main_graph:
            # print(i, netlist_split, i == netlist_split[1], netlist_split[2] not in graph_new)
            if (i == netlist_split[1]) and (netlist_split[2] not in graph_new):
                graph_new.add(netlist_split[2])
                cct_main.add(netlist_line)
                main_graph_has_changed = True
                # print(main_graph_has_changed)
                break

            elif (i == netlist_split[2]) and (netlist_split[1] not in graph_new):
                graph_new.add(netlist_split[1])
                cct_main.add(netlist_line)
                main_graph_has_changed = True
                break
            elif netlist_split[1] in graph_new and netlist_split[2] in graph_new:
                cct_main.add(netlist_line)
                break

        else:
            cct_other_new.add(netlist_line)
        # print(main_graph_has_changed)
    main_graph = graph_new
    cct_other = cct_other_new
print(cct_main)
# cct_main.cg.draw()
cct_main["R_ref"].v.evaluate()

I am going to put this in my codebase real quick to see if it works out; if so, I will close the comment. Thank you for the help!

dlee267 commented 3 months ago

On a similar note, I couldn't find a method that removes unnecessary loops given only one voltage source in the Circuit class. I don't want to make a pull request right away in case it is not in the scope of the library, but do you think it is possible to add something like my code snippet as a method, or is it too niche of a method to add to the lcapy library?

mph- commented 3 months ago

There is no Lcapy method to do this. It would be useful to add if it was generic. Another useful method would prune any component from the netlist that doesn't have a path to the ground node.