MSP Exercise - System decoupling based on Diakoptics

Extra Task

Sample Circuit

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$

Circuit and Simulation Setup

In [1]:
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

Standard mesh current analysis

In [2]:
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))
Impedance matrix: 
[[ 3. -1. -1.  0.]
 [-1.  3.  0. -1.]
 [-1.  0.  3. -1.]
 [ 0. -1. -1.  4.]]

Problem size: 4x4

Mesh currents: 
[[ 6.90909091]
 [ 4.1969697 ]
 [ 1.53030303]
 [ 2.68181818]]

Diakoptics

Equivalent network

Removed network

Soluation of the circuit based on Diakoptics

In [3]:
# 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))
Admittance matrix of the removed network
$\displaystyle Y = $
[[ 1.  0.]
 [ 0.  1.]]

Connection matrix
$\displaystyle C = $
[[-1  0  1  0]
 [ 0  1  0 -1]]

Block-diagonal impedance matrix of the equivalent network
$\displaystyle \tilde{Z} = $
[[ 2. -1.  0.  0.]
 [-1.  2.  0.  0.]
 [ 0.  0.  2. -1.]
 [ 0.  0. -1.  3.]]

Voltage source vector
$\displaystyle E = $
[[15.]
 [ 3.]
 [-5.]
 [ 5.]]

Circuit solution

In [4]:
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))
Block-diagonal impedance matrix of the equivalent network: 
[[ 2. -1.  0.  0.]
 [-1.  2.  0.  0.]
 [ 0.  0.  2. -1.]
 [ 0.  0. -1.  3.]]

Inverse of the block-diagonal impedance matrix can be calculated by independently calculating inverse of 2 blocks of size 2x2:
$\displaystyle \tilde{Z}^{-1} = $
[[ 0.66666667  0.33333333  0.          0.        ]
 [ 0.33333333  0.66666667  0.          0.        ]
 [ 0.          0.          0.6         0.2       ]
 [ 0.          0.          0.2         0.4       ]]

Matrix 
$\displaystyle \tilde{Y} = Y + C \tilde{Z}^{-1} C^T$
can be calculated:
$\displaystyle \tilde{Y}=$
[[ 2.26666667 -0.53333333]
 [-0.53333333  2.06666667]]

Matrix inversion
$\displaystyle \tilde{Y}^{-1}=$
[[ 0.46969697  0.12121212]
 [ 0.12121212  0.51515152]]

Finally, mesh currents can be determined based on:
$\displaystyle i= \tilde{Z}^{-1} E - \tilde{Z}^{-1} C^T \tilde{Y}^{-1} C \tilde{Z}^{-1} E$
$\displaystyle i=$
[[ 6.90909091]
 [ 4.1969697 ]
 [ 1.53030303]
 [ 2.68181818]]

Mesh currents obtained based on standard meshed current analysis: 
[[ 6.90909091]
 [ 4.1969697 ]
 [ 1.53030303]
 [ 2.68181818]]

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$

In [ ]: