Determine mesh currents for the circuit and apply the method of Diakoptics by removing the branches $Z_7$ and $Z_8$
$z_1$=$1 \Omega$, $z_2$=$1 \Omega$, $z_3$=$1 \Omega$, $z_4$=$1 \Omega$, $z_5$=$1 \Omega$, $z_7$=$1 \Omega$, $z_8$=$1 \Omega$
$z_6$=$2 \Omega$
$E_1$=$15 V$, $E_2$=$3 V$, $E_3$=$-5 V$, $E_4$=$5 V$
import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(sign=' ')
from IPython.display import display, Math
# Circuit parameters:
z1 =1.0
z2 = 1.0
z3=1.0
z4 = 1.0
z5 = 1.0
z7 = 1.0
z8 = 1.0
z6 = 2.0
E1 = 15.0
E2=3.0
E3=-5.0
E4=5.0
E_orgn = np.array([ [E1],
[E2],
[E3],
[E4]])
Z_orgn = np.array([[(z1+z3+z7), -z3, -z7, 0],
[-z3, (z3+z5+z8), 0, -z8],
[-z7, 0, (z2+z4+z7), -z4],
[0, -z8, -z4, (z4+z6+z8)]])
# Mesh currents:
i_orgn = np.matmul(np.linalg.inv(Z_orgn),E_orgn)
print('Impedance matrix: \n' + str(Z_orgn))
print('\nProblem size: 4x4')
print('\nMesh currents: \n' + str(i_orgn))
# Removed branches
y7 = 1/z7
y8 = 1/z8
Y = np.array([ [y7, 0],
[ 0, y8]])
print('Admittance matrix of the removed network')
display(Math(r'Y = '))
print(str(Y))
# Connection matrix
C = np.array([ [-1, 0, 1, 0],
[ 0, 1, 0, -1]])
print('\nConnection matrix')
display(Math(r'C = '))
print(str(C))
# Impedance matrix for subsystem 1 in the equivalent network
Z11 = np.array([[(z1+z3), -z3],
[-z3, (z3+z5)]])
# Impedance matrix for subsystem 2 in the equivalent network
Z22 = np.array([[(z2+z4), -z4],
[-z4, (z4+z6)]])
# Block-diagonal impedance matrix of the equivalent network
Z = np.block([[Z11, np.zeros((2, 2))],
[np.zeros((2, 2)), Z22]])
print('\nBlock-diagonal impedance matrix of the equivalent network')
display(Math(r'\tilde{Z} = '))
print(str(Z))
E = np.array([[E1],
[E2],
[E3],
[E4]])
print('\nVoltage source vector')
display(Math(r'E = '))
print(str(E))
print('\nBlock-diagonal impedance matrix of the equivalent network: \n' + str(Z))
print('\nInverse of the block-diagonal impedance matrix can be calculated by independently calculating inverse of 2 blocks of size 2x2:')
# Matrix inverse of impedance matrix of the equivalent network
Z_inv = np.block([[np.linalg.inv(Z11), np.zeros((2, 2))],
[np.zeros((2, 2)), np.linalg.inv(Z22)]])
display(Math(r'\tilde{Z}^{-1} = '))
print(str(Z_inv))
# Matrix \tilde{Y} must be inverted, the size is 2x2
print('\nMatrix ')
display(Math(r'\tilde{Y} = Y + C \tilde{Z}^{-1} C^T'))
print('can be calculated:')
display(Math(r'\tilde{Y}='))
Y_psi = Y + np.matmul(C, np.matmul(Z_inv,C.transpose()))
print(str(Y_psi))
print('\nMatrix inversion')
display(Math(r'\tilde{Y}^{-1}='))
Y_psi_inv = np.linalg.inv(Y_psi)
print(str(Y_psi_inv))
# Solution obtained based on the method of diakoptics
print('\nFinally, mesh currents can be determined based on:')
display(Math(r'i= \tilde{Z}^{-1} E - \tilde{Z}^{-1} C^T \tilde{Y}^{-1} C \tilde{Z}^{-1} E'))
i_diak= np.dot(Z_inv,E) - np.matmul( np.matmul(np.matmul(Z_inv,C.transpose()) , np.matmul(Y_psi_inv,C)) , np.matmul(Z_inv,E))
display(Math(r'i='))
print(str(i_diak))
print('\nMesh currents obtained based on standard meshed current analysis: \n' + str(i_diak))
The solution obtained by the method of diakoptics is an exact solution of original problem.
Computational effort:
Standard solution requires 1 matrix inversion of matrix size 4x4.
Solution based on diakoptics requires 3 matrix inversions of matrix size 2x2.
If we assume that matrix inversion scales with the cube of the matrix size, the following is valid:
$3*2^3 < 4^3$