Include constants perturbation for passivation.
@@ -1047,11 +1047,16 @@ class VectorFitting:
A, B, C, D, E = self._get_ABCDE()
dim_A = np.shape(A)[0]
C_t = C
+ D_t = D
delta = 0.999 # predefined tolerance parameter (users should not need to change this)
# calculate coefficient matrix
A_freq = np.linalg.inv(2j * np.pi * freqs_eval[:, None, None] * np.identity(dim_A)[None, :, :] - A[None, :, :])
- coeffs = np.matmul(A_freq, B[None, :, :])
+
+ # construct coefficient matrix with an extra column for the constants
+ coeffs = np.empty((len(freqs_eval), np.shape(B)[0] + 1, np.shape(B)[1]), dtype=complex)
+ coeffs[:, :-1, :] = np.matmul(A_freq, B[None, :, :])
+ coeffs[:, -1, :] = 1
# iterative compensation of passivity violations
t = 0
@@ -1060,7 +1065,7 @@ class VectorFitting:
logging.info('Passivity enforcement; Iteration {}'.format(t + 1))
# calculate S-matrix at this frequency (shape fxNxN)
- s_eval = self._get_s_from_ABCDE(freqs_eval, A, B, C_t, D, E)
+ s_eval = self._get_s_from_ABCDE(freqs_eval, A, B, C_t, D_t, E)
# singular value decomposition
u, sigma, vh = np.linalg.svd(s_eval)
@@ -1096,9 +1101,12 @@ class VectorFitting:
np.imag(s_viol[:, j, i]))),
rcond=None)
- # perturbe residues by subtracting respective row and column in C_t
+ # perturb residues by subtracting respective row and column in C_t
# one half of the solution will always be 0 due to construction of A and B
- C_t[j, :] = C_t[j, :] - x
+ C_t[j, :] = C_t[j, :] - x[:-1]
+
+ # also perturb constants
+ D_t[j, i] = D_t[j, i] - x[-1]
t += 1
self.history_max_sigma.append(sigma_max)
@@ -1131,6 +1139,7 @@ class VectorFitting:
self.residues[i_response, z] = C_t[i, k] + 1j * C_t[i, k + 1]
k += 2
z += 1
+ self.constant_coeff[i_response] = D_t[i, j]
# run final passivity test to make sure passivation was successful
violation_bands = self.passivity_test()
