algorithm/ovf
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 1 Wiki Activity
Files
87845cda822a11f0cde23e6ca8a6dd0b66ba1154
ovf/core/vectorFitting.py

2618 lines
115 KiB
Python
Raw Blame History

from __future__ import annotations
import logging
import os
import warnings
from timeit import default_timer as timer
from typing import TYPE_CHECKING, Any
import numpy as np
from scipy.integrate import trapezoid
from .util import Axes, axes_kwarg
# imports for type hinting
if TYPE_CHECKING:
from skrf import Network
logger = logging.getLogger(__name__)
class VectorFitting:
"""
This class provides a Python implementation of the Vector Fitting algorithm and various functions for the fit
analysis, passivity evaluation and enforcement, and export of SPICE equivalent circuits.
Parameters
----------
network : :class:`skrf.network.Network`
Network instance of the :math:`N`-port holding the frequency responses to be fitted, for example a
scattering, impedance or admittance matrix.
Examples
--------
Load the `Network`, create a `VectorFitting` instance, perform the fit with a given number of real and
complex-conjugate starting poles:
>>> nw_3port = skrf.Network('my3port.s3p')
>>> vf = skrf.VectorFitting(nw_3port)
>>> vf.vector_fit(n_poles_real=1, n_poles_cmplx=4)
Notes
-----
The fitting code is based on the original algorithm [#Gustavsen_vectfit]_ and on two improvements for relaxed pole
relocation [#Gustavsen_relaxed]_ and efficient (fast) solving [#Deschrijver_fast]_. See also the Vector Fitting
website [#vectfit_website]_ for further information and download of the papers listed below. A Matlab implementation
is also available there for reference.
References
----------
.. [#Gustavsen_vectfit] B. Gustavsen, A. Semlyen, "Rational Approximation of Frequency Domain Responses by Vector
Fitting", IEEE Transactions on Power Delivery, vol. 14, no. 3, pp. 1052-1061, July 1999,
DOI: https://doi.org/10.1109/61.772353
.. [#Gustavsen_relaxed] B. Gustavsen, "Improving the Pole Relocating Properties of Vector Fitting", IEEE
Transactions on Power Delivery, vol. 21, no. 3, pp. 1587-1592, July 2006,
DOI: https://doi.org/10.1109/TPWRD.2005.860281
.. [#Deschrijver_fast] D. Deschrijver, M. Mrozowski, T. Dhaene, D. De Zutter, "Marcomodeling of Multiport Systems
Using a Fast Implementation of the Vector Fitting Method", IEEE Microwave and Wireless Components Letters,
vol. 18, no. 6, pp. 383-385, June 2008, DOI: https://doi.org/10.1109/LMWC.2008.922585
.. [#vectfit_website] Vector Fitting website: https://www.sintef.no/projectweb/vectorfitting/
"""
def __init__(self, network: Network):
self.network = network
""" Instance variable holding the Network to be fitted. This is the Network passed during initialization,
which may be changed or set to *None*. """
self.poles = None
""" Instance variable holding the list of fitted poles. Will be initialized by :func:`vector_fit`. """
self.residues = None
""" Instance variable holding the list of fitted residues. Will be initialized by :func:`vector_fit`. """
self.proportional_coeff = None
""" Instance variable holding the list of fitted proportional coefficients. Will be initialized by
:func:`vector_fit`. """
self.constant_coeff = None
""" Instance variable holding the list of fitted constants. Will be initialized by :func:`vector_fit`. """
self.max_iterations = 100
""" Instance variable specifying the maximum number of iterations for the fitting process and for the passivity
enforcement. To be changed by the user before calling :func:`vector_fit` and/or :func:`passivity_enforce`. """
self.max_tol = 1e-6
""" Instance variable specifying the convergence criterion in terms of relative tolerance. To be changed by the
user before calling :func:`vector_fit`. """
self.wall_clock_time = 0
""" Instance variable holding the wall-clock time (in seconds) consumed by the most recent fitting process with
:func:`vector_fit`. Subsequent calls of :func:`vector_fit` will overwrite this value. """
self.d_res_history = []
self.delta_max_history = []
self.history_max_sigma = []
self.history_cond_A = []
self.history_rank_deficiency = []
@staticmethod
def get_spurious(poles: np.ndarray, residues: np.ndarray, n_freqs: int = 101, gamma: float = 0.03) -> np.ndarray:
"""
Classifies fitted pole-residue pairs as spurious or not spurious. The implementation is based on the evaluation
of band-limited energy norms (p=2) of the resonance curves of individual pole-residue pairs, as proposed in
[#Grivet-Talocia]_.
Parameters
----------
poles : ndarray, shape (N)
Array of fitted poles
residues : ndarray, shape (M, N)
Array of fitted residues
n_freqs : int, optional
Number of frequencies for the evaluation. The frequency range is chosen automatically and the default
101 frequencies should be appropriate in most cases.
gamma : float, optional
Sensitivity threshold for the classification. Typical values range from 0.01 to 0.05.
Returns
-------
ndarray, bool, shape (M)
Boolean array having the same shape as :attr:`poles`. `True` marks the respective pole as spurious.
References
----------
.. [#Grivet-Talocia] S. Grivet-Talocia and M. Bandinu, "Improving the convergence of vector fitting for
equivalent circuit extraction from noisy frequency responses," in IEEE Transactions on Electromagnetic
Compatibility, vol. 48, no. 1, pp. 104-120, Feb. 2006, DOI: https://doi.org/10.1109/TEMC.2006.870814
"""
omega_eval = np.linspace(np.min(poles.imag) / 3, np.max(poles.imag) * 3, n_freqs)
h = (residues[:, None, :] / (1j * omega_eval[:, None] - poles)
+ np.conj(residues[:, None, :]) / (1j * omega_eval[:, None] - np.conj(poles)))
norm2 = np.sqrt(trapezoid(h.real ** 2 + h.imag ** 2, omega_eval, axis=1))
spurious = np.all(norm2 / np.mean(norm2) < gamma, axis=0)
return spurious
@staticmethod
def get_model_order(poles: np.ndarray) -> int:
"""
Returns the model order calculated with :math:`N_{real} + 2 N_{complex}` for a given set of poles.
Parameters
----------
poles: ndarray
The poles of the model as a list or NumPy array.
Returns
-------
order: int
"""
# poles.imag != 0 is True(1) for complex poles, False (0) for real poles.
# Adding one to each element gives 2 columns for complex and 1 column for real poles.
return np.sum((poles.imag != 0) + 1)
def vector_fit(self, n_poles_real: int = 2, n_poles_cmplx: int = 2, init_pole_spacing: str = 'lin',
parameter_type: str = 's', fit_constant: bool = True, fit_proportional: bool = False) -> None:
"""
Main work routine performing the vector fit. The results will be stored in the class variables
:attr:`poles`, :attr:`residues`, :attr:`proportional_coeff` and :attr:`constant_coeff`.
Parameters
----------
n_poles_real : int, optional
Number of initial real poles. See notes.
n_poles_cmplx : int, optional
Number of initial complex conjugate poles. See notes.
init_pole_spacing : str, optional
Type of initial pole spacing across the frequency interval of the S-matrix. Either *linear* (`'lin'`),
*logarithmic* (`'log'`), or `custom`. In case of `custom`, the initial poles must be stored in :attr:`poles`
as a NumPy array before calling this method. They will be overwritten by the final poles. The
initialization parameters `n_poles_real` and `n_poles_cmplx` will be ignored in case of `'custom'`.
parameter_type : str, optional
Representation type of the frequency responses to be fitted. Either *scattering* (`'s'` or `'S'`),
*impedance* (`'z'` or `'Z'`) or *admittance* (`'y'` or `'Y'`). It's recommended to perform the fit on the
original S parameters. Otherwise, scikit-rf will convert the responses from S to Z or Y, which might work
for the fit but can cause other issues.
fit_constant : bool, optional
Include a constant term **d** in the fit.
fit_proportional : bool, optional
Include a proportional term **e** in the fit.
Returns
-------
None
No return value.
Notes
-----
The required number of real or complex conjugate starting poles depends on the behaviour of the frequency
responses. To fit a smooth response such as a low-pass characteristic, 1-3 real poles and no complex conjugate
poles is usually sufficient. If resonances or other types of peaks are present in some or all of the responses,
a similar number of complex conjugate poles is required. Be careful not to use too many poles, as excessive
poles will not only increase the computation workload during the fitting and the subsequent use of the model,
but they can also introduce unwanted resonances at frequencies well outside the fit interval.
See Also
--------
auto_fit : Automatic vector fitting routine with pole adding and skimming.
"""
timer_start = timer()
# use normalized frequencies during the iterations (seems to be more stable during least-squares fit)
norm = np.average(self.network.f)
# norm = np.exp(np.mean(np.log(self.network.f)))
freqs_norm = np.array(self.network.f) / norm
# get initial poles
poles = self._init_poles(freqs_norm, n_poles_real, n_poles_cmplx, init_pole_spacing)
# check and normalize custom poles
if poles is None:
if self.poles is not None and len(self.poles) > 0:
poles = self.poles / norm
else:
raise ValueError('Initial poles must be provided in `self.poles` when calling with '
'`init_pole_spacing == \'custom\'`.')
# save initial poles (un-normalize first)
initial_poles = poles * norm
max_singular = 1
logger.info('### Starting pole relocation process.\n')
# select network representation type
if parameter_type.lower() == 's':
nw_responses = self.network.s
elif parameter_type.lower() == 'z':
nw_responses = self.network.z
elif parameter_type.lower() == 'y':
nw_responses = self.network.y
else:
warnings.warn('Invalid choice of matrix parameter type (S, Z, or Y); proceeding with scattering '
'representation.', UserWarning, stacklevel=2)
nw_responses = self.network.s
# stack frequency responses as a single vector
# stacking order (row-major):
# s11, s12, s13, ..., s21, s22, s23, ...
freq_responses = []
for i in range(self.network.nports):
for j in range(self.network.nports):
freq_responses.append(nw_responses[:, i, j])
freq_responses = np.array(freq_responses)
# responses will be weighted according to their norm;
# alternative: equal weights with weight_response = 1.0
# or anti-proportional weights with weight_response = 1 / np.linalg.norm(freq_response)
weights_responses = np.linalg.norm(freq_responses, axis=1)
#weights_responses = np.ones(self.network.nports ** 2)
#weights_responses = 10 / np.exp(np.mean(np.log(np.abs(freq_responses)), axis=1))
# ITERATIVE FITTING OF POLES to the provided frequency responses
# initial set of poles will be replaced with new poles after every iteration
iterations = self.max_iterations
self.d_res_history = []
self.delta_max_history = []
self.history_cond_A = []
self.history_rank_deficiency = []
converged = False
# POLE RELOCATION LOOP
while iterations > 0:
logger.info(f'Iteration {self.max_iterations - iterations + 1}')
poles, d_res, cond, rank_deficiency, residuals, singular_vals = self._pole_relocation(
poles, freqs_norm, freq_responses, weights_responses, fit_constant, fit_proportional)
logger.info(f'Condition number of coefficient matrix is {int(cond)}')
self.history_cond_A.append(cond)
self.history_rank_deficiency.append(rank_deficiency)
logger.info(f'Rank deficiency is {rank_deficiency}.')
self.d_res_history.append(d_res)
logger.info(f'd_res = {d_res}')
# calculate relative changes in the singular values; stop iteration loop once poles have converged
new_max_singular = np.amax(singular_vals)
delta_max = np.abs(1 - new_max_singular / max_singular)
self.delta_max_history.append(delta_max)
logger.info(f'Max. relative change in residues = {delta_max}\n')
max_singular = new_max_singular
stop = False
if delta_max < self.max_tol:
if converged:
# is really converged, finish
logger.info(f'Pole relocation process converged after {self.max_iterations - iterations + 1} '
'iterations.')
stop = True
else:
# might be converged, but do one last run to be sure
converged = True
else:
if converged:
# is not really converged, continue
converged = False
iterations -= 1
if iterations == 0:
# loop ran into iterations limit; trying to assess the issue
max_cond = np.amax(self.history_cond_A)
max_deficiency = np.amax(self.history_rank_deficiency)
if max_cond > 1e10:
hint_illcond = ('\nHint: the linear system was ill-conditioned (max. condition number was '
f'{max_cond}).')
else:
hint_illcond = ''
if max_deficiency < 0:
hint_rank = ('\nHint: the coefficient matrix was rank-deficient (max. rank deficiency was '
f'{max_deficiency}).')
else:
hint_rank = ''
if converged and stop is False:
warnings.warn('Vector Fitting: The pole relocation process barely converged to tolerance. '
f'It took the max. number of iterations (N_max = {self.max_iterations}). '
'The results might not have converged properly.'
+ hint_illcond + hint_rank, RuntimeWarning, stacklevel=2)
else:
warnings.warn('Vector Fitting: The pole relocation process stopped after reaching the '
f'maximum number of iterations (N_max = {self.max_iterations}). '
'The results did not converge properly.'
+ hint_illcond + hint_rank, RuntimeWarning, stacklevel=2)
if stop:
iterations = 0
# ITERATIONS DONE
logger.info('Initial poles before relocation:')
logger.info(initial_poles)
logger.info('Final poles:')
logger.info(poles * norm)
logger.info('\n### Starting residues calculation process.\n')
# finally, solve for the residues with the previously calculated poles
residues, constant_coeff, proportional_coeff, residuals, rank, singular_vals = self._fit_residues(
poles, freqs_norm, freq_responses, fit_constant, fit_proportional)
# save poles, residues, d, e in actual frequencies (un-normalized)
self.poles = poles * norm
self.residues = np.array(residues) * norm
self.constant_coeff = np.array(constant_coeff)
self.proportional_coeff = np.array(proportional_coeff) / norm
timer_stop = timer()
self.wall_clock_time = timer_stop - timer_start
logger.info(f'\n### Vector fitting finished in {self.wall_clock_time} seconds.\n')
# raise a warning if the fitted Network is passive but the fit is not (only without proportional_coeff):
if self.network.is_passive() and not fit_proportional:
if not self.is_passive():
warnings.warn('The fitted network is passive, but the vector fit is not passive. Consider running '
'`passivity_enforce()` to enforce passivity before using this model.',
UserWarning, stacklevel=2)
def auto_fit(self, n_poles_init_real: int = 3, n_poles_init_cmplx: int = 3, n_poles_add: int = 3,
model_order_max: int = 100, iters_start: int = 3, iters_inter: int = 3, iters_final: int = 5,
target_error: float = 1e-2, alpha: float = 0.03, gamma: float = 0.03, nu_samples: float = 1.0,
parameter_type: str = 's') -> (np.ndarray, np.ndarray):
"""
Automatic fitting routine implementing the "vector fitting with adding and skimming" algorithm as proposed in
[#Grivet-Talocia]_. This algorithm is able to provide high quality macromodels with automatic model order
optimization, while improving both the rate of convergence and the fit quality in case of noisy data.
The resulting model parameters will be stored in the class variables :attr:`poles`, :attr:`residues`,
:attr:`proportional_coeff` and :attr:`constant_coeff`.
Parameters
----------
n_poles_init_real: int, optional
Number of real poles in the initial model.
n_poles_init_cmplx: int, optional
Number of complex conjugate poles in the initial model.
n_poles_add: int, optional
Number of new poles allowed to be added in each refinement iteration, if possible. This controls how fast
the model order is allowed to grow. Unnecessary poles will have to be skimmed and removed later. This
parameter has a strong effect on the convergence.
model_order_max: int, optional
Maximum model order as calculated with :math:`N_{real} + 2 N_{complex}`. This parameter provides a stopping
criterion in case the refinement process is not converging.
iters_start: int, optional
Number of initial iterations for pole relocations as in regular vector fitting.
iters_inter: int, optional
Number of intermediate iterations for pole relocations during each iteration of the refinement process.
iters_final: int, optional
Number of final iterations for pole relocations after the refinement process.
target_error: float, optional
Target for the model error to be reached during the refinement process. The actual achievable error is
bound by the noise in the data. If specified with a number greater than the noise floor, this parameter
provides another stopping criterion for the refinement process. It therefore affects both the convergence,
the final error, and the final model order (number of poles used in the model).
alpha: float, optional
Threshold for the error decay to stop the refinement loop in case of error stagnation. This parameter
provides another stopping criterion for cases where the model already has enough poles but the target error
still cannot be reached because of excess noise (target error too small for noise level in the data).
gamma: float, optional
Threshold for the detection of spurious poles.
nu_samples: float, optional
Required and enforced (relative) spacing in terms of frequency samples between existing poles and
relocated or added poles. The number can be a float, it does not have to be an integer.
parameter_type: str, optional
Representation type of the frequency responses to be fitted. Either *scattering* (`'s'` or `'S'`),
*impedance* (`'z'` or `'Z'`) or *admittance* (`'y'` or `'Y'`). It's recommended to perform the fit on the
original S parameters. Otherwise, scikit-rf will convert the responses from S to Z or Y, which might work
for the fit but can cause other issues.
Returns
-------
None
No return value.
See Also
--------
vector_fit : Regular vector fitting routine.
References
----------
.. [#Grivet-Talocia] S. Grivet-Talocia and M. Bandinu, "Improving the convergence of vector fitting for
equivalent circuit extraction from noisy frequency responses," in IEEE Transactions on Electromagnetic
Compatibility, vol. 48, no. 1, pp. 104-120, Feb. 2006, DOI: https://doi.org/10.1109/TEMC.2006.870814
"""
self.d_res_history = []
self.delta_max_history = []
self.history_cond_A = []
self.history_rank_deficiency = []
max_singular = 1
error_peak_history = []
model_order_history = []
timer_start = timer()
# use normalized frequencies during the iterations (seems to be more stable during least-squares fit)
norm = np.average(self.network.f)
# norm = np.exp(np.mean(np.log(self.network.f)))
freqs_norm = np.array(self.network.f) / norm
omega_norm = 2 * np.pi * freqs_norm
nu = (omega_norm[1] - omega_norm[0]) * nu_samples
# get initial poles
poles = self._init_poles(freqs_norm, n_poles_init_real, n_poles_init_cmplx, 'lin')
logger.info('### Starting pole relocation process.\n')
# select network representation type
if parameter_type.lower() == 's':
nw_responses = self.network.s
fit_constant = True
fit_proportional = False
elif parameter_type.lower() == 'z':
nw_responses = self.network.z
fit_constant = True
fit_proportional = True
elif parameter_type.lower() == 'y':
nw_responses = self.network.y
fit_constant = True
fit_proportional = True
else:
warnings.warn('Invalid choice of matrix parameter type (S, Z, or Y); proceeding with scattering '
'representation.', UserWarning, stacklevel=2)
nw_responses = self.network.s
fit_constant = True
fit_proportional = False
# stack frequency responses as a single vector
# stacking order (row-major):
# s11, s12, s13, ..., s21, s22, s23, ...
freq_responses = []
for i in range(self.network.nports):
for j in range(self.network.nports):
freq_responses.append(nw_responses[:, i, j])
freq_responses = np.array(freq_responses)
# responses will be weighted according to their norm;
# alternative: equal weights with weight_response = 1.0
# or anti-proportional weights with weight_response = 1 / np.linalg.norm(freq_response)
weights_responses = np.linalg.norm(freq_responses, axis=1)
# weights_responses = np.ones(self.network.nports ** 2)
# weights_responses = 10 / np.exp(np.mean(np.log(np.abs(freq_responses)), axis=1))
# INITIAL POLE RELOCATION FOR i_start ITERATIONS
for _ in range(iters_start):
poles, d_res, cond, rank_deficiency, residuals, singular_vals = self._pole_relocation(
poles, freqs_norm, freq_responses, weights_responses, fit_constant, fit_proportional)
self.d_res_history.append(d_res)
logger.info(f'Condition number of coefficient matrix is {int(cond)}')
self.history_cond_A.append(cond)
self.history_rank_deficiency.append(rank_deficiency)
logger.info(f'Rank deficiency is {rank_deficiency}.')
new_max_singular = np.amax(singular_vals)
delta_max = np.abs(1 - new_max_singular / max_singular)
self.delta_max_history.append(delta_max)
logger.info(f'Max. relative change in residues = {delta_max}\n')
max_singular = new_max_singular
# RESIDUE FITTING FOR ERROR COMPUTATION
residues, constant_coeff, proportional_coeff, residuals, rank, singular_vals = self._fit_residues(
poles, freqs_norm, freq_responses, fit_constant, fit_proportional, enforce_dc=False)
delta = self._get_delta(poles, residues, constant_coeff, proportional_coeff, freqs_norm, freq_responses,
weights_responses)
error_peak = np.max(delta)
error_peak_history.append(error_peak)
model_order = self.get_model_order(poles)
model_order_history.append(model_order)
delta_eps = 10 * alpha
# POLE SKIMMING AND ADDING LOOP
while error_peak > target_error and model_order < model_order_max and delta_eps > alpha:
# SKIMMING OF SPURIOUS POLES
spurious = self.get_spurious(poles, residues, gamma=gamma)
n_skim = np.sum(spurious)
poles = poles[~spurious]
# REPLACING SPURIOUS POLE AND ADDING NEW POLES
idx_freqs_start, idx_freqs_stop, idx_freqs_max, delta_mean_bands = self._find_error_bands(freqs_norm, delta)
n_bands = len(idx_freqs_max)
if n_bands < n_skim:
n_add = n_bands
elif n_bands < n_skim + n_poles_add:
n_add = n_bands
else:
n_add = n_skim + n_poles_add
for i in range(n_add):
omega_add = omega_norm[idx_freqs_max[i]]
pole_add = (-0.01 + 1j) * omega_add
# compute distance to neighbouring poles
abs_poles_existing = np.abs(poles) - pole_add.imag # (equation 16)
#abs_poles_existing = np.abs(poles - pole_add) # (equation 17)
# avoid forbidden bands (too close to neighbour)
if np.min(abs_poles_existing) < nu or pole_add.imag < nu:
# decide shift direction (towards higher or lower frequencies)
if idx_freqs_max[i] > 0:
delta_below = delta[idx_freqs_max[i] - 1]
else:
delta_below = 0
if idx_freqs_max[i] < len(omega_norm) - 1:
delta_above = delta[idx_freqs_max[i] + 1]
else:
delta_above = 0
if delta_above > delta_below:
# shift to higher frequencies
pole_add += 1j * nu
else:
# shift to lower frequencies
pole_add -= 1j * nu
poles = np.append(poles, [pole_add])
# INTERMEDIATE POLE RELOCATION FOR i_inter ITERATIONS
for _ in range(iters_inter):
poles, d_res, cond, rank_deficiency, residuals, singular_vals = self._pole_relocation(
poles, freqs_norm, freq_responses, weights_responses, fit_constant, fit_proportional)
self.d_res_history.append(d_res)
logger.info(f'Condition number of coefficient matrix is {int(cond)}')
self.history_cond_A.append(cond)
self.history_rank_deficiency.append(rank_deficiency)
logger.info(f'Rank deficiency is {rank_deficiency}.')
new_max_singular = np.amax(singular_vals)
delta_max = np.abs(1 - new_max_singular / max_singular)
self.delta_max_history.append(delta_max)
logger.info(f'Max. relative change in residues = {delta_max}\n')
max_singular = new_max_singular
# RESIDUE FITTING FOR ERROR COMPUTATION
residues, constant_coeff, proportional_coeff, residuals, rank, singular_vals = self._fit_residues(
poles, freqs_norm, freq_responses, fit_constant, fit_proportional, enforce_dc=False)
delta = self._get_delta(poles, residues, constant_coeff, proportional_coeff, freqs_norm, freq_responses,
weights_responses)
error_peak_history.append(np.max(delta))
m = 3
if len(error_peak_history) > m:
delta_eps = np.mean(np.abs(np.diff(error_peak_history[-1-m:-1])))
else:
delta_eps = 1
model_order = self.get_model_order(poles)
model_order_history.append(model_order)
# SKIMMING OF SPURIOUS POLES
spurious = self.get_spurious(poles, residues, gamma=gamma)
poles = poles[~spurious]
# FINAL POLE RELOCATION FOR i_final ITERATIONS
for _ in range(iters_final):
poles, d_res, cond, rank_deficiency, residuals, singular_vals = self._pole_relocation(
poles, freqs_norm, freq_responses, weights_responses, fit_constant, fit_proportional)
self.d_res_history.append(d_res)
logger.info(f'Condition number of coefficient matrix is {int(cond)}')
self.history_cond_A.append(cond)
self.history_rank_deficiency.append(rank_deficiency)
logger.info(f'Rank deficiency is {rank_deficiency}.')
new_max_singular = np.amax(singular_vals)
delta_max = np.abs(1 - new_max_singular / max_singular)
self.delta_max_history.append(delta_max)
logger.info(f'Max. relative change in residues = {delta_max}\n')
max_singular = new_max_singular
# FINAL RESIDUE FITTING
residues, constant_coeff, proportional_coeff, residuals, rank, singular_vals = self._fit_residues(
poles, freqs_norm, freq_responses, fit_constant, fit_proportional, enforce_dc=True)
# save poles, residues, d, e in actual frequencies (un-normalized)
self.poles = poles * norm
self.residues = np.array(residues) * norm
self.constant_coeff = np.array(constant_coeff)
self.proportional_coeff = np.array(proportional_coeff) / norm
timer_stop = timer()
self.wall_clock_time = timer_stop - timer_start
@staticmethod
def _init_poles(freqs: list, n_poles_real: int, n_poles_cmplx: int, init_pole_spacing: str):
# create initial poles and space them across the frequencies in the provided Touchstone file
fmin = np.amin(freqs)
fmax = np.amax(freqs)
# poles cannot be at f=0; hence, f_min for starting pole must be greater than 0
if fmin == 0.0:
# random choice: use 1/1000 of first non-zero frequency
fmin = freqs[1] / 1000
init_pole_spacing = init_pole_spacing.lower()
if init_pole_spacing == 'log':
pole_freqs_real = np.geomspace(fmin, fmax, n_poles_real)
pole_freqs_cmplx = np.geomspace(fmin, fmax, n_poles_cmplx)
elif init_pole_spacing == 'lin':
pole_freqs_real = np.linspace(fmin, fmax, n_poles_real)
pole_freqs_cmplx = np.linspace(fmin, fmax, n_poles_cmplx)
elif init_pole_spacing == 'custom':
pole_freqs_real = None
pole_freqs_cmplx = None
else:
warnings.warn('Invalid choice of initial pole spacing; proceeding with linear spacing.',
UserWarning, stacklevel=2)
pole_freqs_real = np.linspace(fmin, fmax, n_poles_real)
pole_freqs_cmplx = np.linspace(fmin, fmax, n_poles_cmplx)
if pole_freqs_real is not None and pole_freqs_cmplx is not None:
# init poles array of correct length
poles = np.zeros(n_poles_real + n_poles_cmplx, dtype=complex)
# add real poles
for i, f in enumerate(pole_freqs_real):
omega = 2 * np.pi * f
poles[i] = -1 * omega
# add complex-conjugate poles (store only positive imaginary parts)
i_offset = len(pole_freqs_real)
for i, f in enumerate(pole_freqs_cmplx):
omega = 2 * np.pi * f
poles[i_offset + i] = (-0.01 + 1j) * omega
return poles
else:
return None
@staticmethod
def _pole_relocation(poles, freqs, freq_responses, weights_responses, fit_constant, fit_proportional):
n_responses, n_freqs = np.shape(freq_responses)
n_samples = n_responses * n_freqs
omega = 2 * np.pi * freqs
s = 1j * omega
# weight of extra equation to avoid trivial solution
weight_extra = np.linalg.norm(weights_responses[:, None] * freq_responses) / n_samples
# weights w are applied directly to the samples, which get squared during least-squares fitting; hence sqrt(w)
weights_responses = np.sqrt(weights_responses)
weight_extra = np.sqrt(weight_extra)
# count number of rows and columns in final coefficient matrix to solve for (c_res, d_res)
# (ratio #real/#complex poles might change during iterations)
# We need two columns for complex poles and one column for real poles in A matrix.
# This number equals the model order.
n_cols_unused = VectorFitting.get_model_order(poles)
n_cols_used = n_cols_unused
n_cols_used += 1
idx_constant = []
idx_proportional = []
if fit_constant:
idx_constant = [n_cols_unused]
n_cols_unused += 1
if fit_proportional:
idx_proportional = [n_cols_unused]
n_cols_unused += 1
real_mask = poles.imag == 0
# list of indices in 'poles' with real values
idx_poles_real = np.nonzero(real_mask)[0]
# list of indices in 'poles' with complex values
idx_poles_complex = np.nonzero(~real_mask)[0]
# positions (columns) of coefficients for real and complex-conjugate terms in the rows of A determine the
# respective positions of the calculated residues in the results vector.
# to have them ordered properly for the subsequent assembly of the test matrix H for eigenvalue extraction,
# place real poles first, then complex-conjugate poles with their respective real and imaginary parts:
# [r1', r2', ..., (r3', r3''), (r4', r4''), ...]
n_real = len(idx_poles_real)
n_cmplx = len(idx_poles_complex)
idx_res_real = np.arange(n_real)
idx_res_complex_re = n_real + 2 * np.arange(n_cmplx)
idx_res_complex_im = idx_res_complex_re + 1
# complex coefficient matrix of shape [N_responses, N_freqs, n_cols_unused + n_cols_used]
# layout of each row:
# [pole1, pole2, ..., (constant), (proportional), pole1, pole2, ..., constant]
A = np.empty((n_responses, n_freqs, n_cols_unused + n_cols_used), dtype=complex)
# calculate coefficients for real and complex residues in the solution vector
#
# real pole-residue term (r = r', p = p'):
# fractional term is r' / (s - p')
# coefficient for r' is 1 / (s - p')
coeff_real = 1 / (s[:, None] - poles[None, idx_poles_real])
# complex-conjugate pole-residue pair (r = r' + j r'', p = p' + j p''):
# fractional term is r / (s - p) + conj(r) / (s - conj(p))
# = [1 / (s - p) + 1 / (s - conj(p))] * r' + [1j / (s - p) - 1j / (s - conj(p))] * r''
# coefficient for r' is 1 / (s - p) + 1 / (s - conj(p))
# coefficient for r'' is 1j / (s - p) - 1j / (s - conj(p))
coeff_complex_re = (1 / (s[:, None] - poles[None, idx_poles_complex]) +
1 / (s[:, None] - np.conj(poles[None, idx_poles_complex])))
coeff_complex_im = (1j / (s[:, None] - poles[None, idx_poles_complex]) -
1j / (s[:, None] - np.conj(poles[None, idx_poles_complex])))
# part 1: first sum of rational functions (variable c)
A[:, :, idx_res_real] = coeff_real
A[:, :, idx_res_complex_re] = coeff_complex_re
A[:, :, idx_res_complex_im] = coeff_complex_im
# part 2: constant (variable d) and proportional term (variable e)
A[:, :, idx_constant] = 1
A[:, :, idx_proportional] = s[:, None]
# part 3: second sum of rational functions multiplied with frequency response (variable c_res)
A[:, :, n_cols_unused + idx_res_real] = -1 * freq_responses[:, :, None] * coeff_real
A[:, :, n_cols_unused + idx_res_complex_re] = -1 * freq_responses[:, :, None] * coeff_complex_re
A[:, :, n_cols_unused + idx_res_complex_im] = -1 * freq_responses[:, :, None] * coeff_complex_im
# part 4: constant (variable d_res)
A[:, :, -1] = -1 * freq_responses
A_ri = np.hstack((A.real, A.imag))
# calculation of matrix sizes after QR decomposition:
# stacked coefficient matrix (A.real, A.imag) has shape (L, M, N)
# with
# L = n_responses = n_ports ** 2
# M = 2 * n_freqs (because of hstack with 2x n_freqs)
# N = n_cols_unused + n_cols_used
# then
# R has shape (L, K, N) with K = min(M, N)
dim_m = 2 * n_freqs
dim_n = n_cols_unused + n_cols_used
dim_k = min(dim_m, dim_n)
# QR decomposition
# R = np.linalg.qr(A_ri, 'r')
# direct QR of stacked matrices for linalg.qr() only works with numpy>=1.22.0
# workaround for old numpy:
R = np.empty((n_responses, dim_k, dim_n))
for i in range(n_responses):
R[i] = np.linalg.qr(A_ri[i], mode='r')
# only R22 is required to solve for c_res and d_res
# R12 and R22 can have a different number of rows, depending on K
if dim_k == dim_m:
# K = M
n_rows_r12 = n_freqs
n_rows_r22 = n_freqs
else:
# K = N
n_rows_r12 = n_cols_unused
n_rows_r22 = n_cols_used
R22 = R[:, n_rows_r12:, n_cols_unused:]
# weighting
R22 = weights_responses[:, None, None] * R22
# assemble compressed coefficient matrix A_fast by row-stacking individual upper triangular matrices R22
dim0 = n_responses * n_rows_r22 + 1
A_fast = np.empty((dim0, n_cols_used))
A_fast[:-1, :] = R22.reshape((dim0 - 1, n_cols_used))
# extra equation to avoid trivial solution
A_fast[-1, idx_res_real] = np.sum(coeff_real.real, axis=0)
A_fast[-1, idx_res_complex_re] = np.sum(coeff_complex_re.real, axis=0)
A_fast[-1, idx_res_complex_im] = np.sum(coeff_complex_im.real, axis=0)
A_fast[-1, -1] = n_freqs
# weighting
A_fast[-1, :] = weight_extra * A_fast[-1, :]
scaling = 1 / np.linalg.norm(A_fast, axis=0)
A_fast = scaling * A_fast
# right hand side vector (weighted)
b = np.zeros(dim0)
b[-1] = weight_extra * n_samples
# check condition of the linear system
cond = np.linalg.cond(A_fast)
full_rank = np.min(A_fast.shape)
# solve least squares for real parts
x, residuals, rank, singular_vals = np.linalg.lstsq(A_fast, b, rcond=None)
x = scaling * x
# rank deficiency
rank_deficiency = full_rank - rank
# assemble individual result vectors from single LS result x
c_res = x[:-1]
d_res = x[-1]
# check if d_res is suited for zeros calculation
tol_res = 1e-8
if np.abs(d_res) < tol_res:
# d_res is too small, discard solution and proceed the |d_res| = tol_res
logger.info(f'Replacing d_res solution as it was too small ({d_res}).')
d_res = tol_res * (d_res / np.abs(d_res))
# build test matrix H, which will hold the new poles as eigenvalues
H = np.zeros((len(c_res), len(c_res)))
poles_real = poles[np.nonzero(real_mask)]
poles_cplx = poles[np.nonzero(~real_mask)]
H[idx_res_real, idx_res_real] = poles_real.real
H[idx_res_real] -= c_res / d_res
H[idx_res_complex_re, idx_res_complex_re] = poles_cplx.real
H[idx_res_complex_re, idx_res_complex_im] = poles_cplx.imag
H[idx_res_complex_im, idx_res_complex_re] = -1 * poles_cplx.imag
H[idx_res_complex_im, idx_res_complex_im] = poles_cplx.real
H[idx_res_complex_re] -= 2 * c_res / d_res
poles_new = np.linalg.eigvals(H)
# replace poles for next iteration
# complex poles need to come in complex conjugate pairs; append only the positive part
poles = poles_new[np.nonzero(poles_new.imag >= 0)]
# flip real part of unstable poles (real part needs to be negative for stability)
poles.real = -1 * np.abs(poles.real)
return poles, d_res, cond, rank_deficiency, residuals, singular_vals
@staticmethod
def _fit_residues(poles, freqs, freq_responses, fit_constant, fit_proportional, enforce_dc=True):
n_responses, n_freqs = np.shape(freq_responses)
omega = 2 * np.pi * freqs
s = 1j * omega
# We need two columns for complex poles and one column for real poles in A matrix.
# This number equals the model order.
n_cols = VectorFitting.get_model_order(poles)
idx_constant = []
idx_proportional = []
if fit_constant:
idx_constant = [n_cols]
n_cols += 1
if fit_proportional:
idx_proportional = [n_cols]
n_cols += 1
# list of indices in 'poles' with real and with complex values
real_mask = poles.imag == 0
idx_poles_real = np.nonzero(real_mask)[0]
idx_poles_complex = np.nonzero(~real_mask)[0]
# find and save indices of real and complex poles in the poles list
i = 0
idx_res_real = []
idx_res_complex_re = []
idx_res_complex_im = []
for pole in poles:
if pole.imag == 0:
idx_res_real.append(i)
i += 1
else:
idx_res_complex_re.append(i)
idx_res_complex_im.append(i + 1)
i += 2
# complex coefficient matrix of shape [N_freqs, n_cols]
# layout of each row:
# [pole1, pole2, ..., (constant), (proportional)]
A = np.empty((n_freqs, n_cols), dtype=complex)
# calculate coefficients for real and complex residues in the solution vector
#
# real pole-residue term (r = r', p = p'):
# fractional term is r' / (s - p')
# coefficient for r' is 1 / (s - p')
coeff_real = 1 / (s[:, None] - poles[None, idx_poles_real])
# complex-conjugate pole-residue pair (r = r' + j r'', p = p' + j p''):
# fractional term is r / (s - p) + conj(r) / (s - conj(p))
# = [1 / (s - p) + 1 / (s - conj(p))] * r' + [1j / (s - p) - 1j / (s - conj(p))] * r''
# coefficient for r' is 1 / (s - p) + 1 / (s - conj(p))
# coefficient for r'' is 1j / (s - p) - 1j / (s - conj(p))
coeff_complex_re = (1 / (s[:, None] - poles[None, idx_poles_complex]) +
1 / (s[:, None] - np.conj(poles[None, idx_poles_complex])))
coeff_complex_im = (1j / (s[:, None] - poles[None, idx_poles_complex]) -
1j / (s[:, None] - np.conj(poles[None, idx_poles_complex])))
# part 1: first sum of rational functions (variable c)
A[:, idx_res_real] = coeff_real
A[:, idx_res_complex_re] = coeff_complex_re
A[:, idx_res_complex_im] = coeff_complex_im
# part 2: constant (variable d) and proportional term (variable e)
A[:, idx_constant] = 1
A[:, idx_proportional] = s[:, None]
scaling = 1 / np.linalg.norm(A, axis=0)
A = scaling * A
# DC POINT ENFORCEMENT
if enforce_dc and freqs[0] == 0.0:
# data contains the dc point; enforce dc point via linear equality constraint:
# 1: remove one variable from the solution vector (constant term, if possible).
# 2: solve remaining linear system (without data at dc) with regular least-squares, as usual. the size of
# the solution vector, the coefficient matrix, and the right-hand side are reduced by 1
# 3: calculate the removed variable (constant term) with the data from the dc point
#
# linear system: A * x = b
# solution vector x contains the unknown residues
# right-hand side b contains the frequency response to be fitted, sorted by ascending frequency (dc first)
# coefficient matrix A and vector b are split: A = [[A11, A12], [A21, A22]], b = [[b1], [b2]]
# [A11, A12] is the first row used later for dc enforcement
# A21 is a column vector, which is not required anymore
# A22 is the rest of the matrix for usual least-squares fitting
# indexing mask of constrained variable in the columns of matrix A
mask_idx_constrained = np.zeros(n_cols, dtype=bool)
if fit_constant:
# use constant term for constrained
mask_idx_constrained[idx_constant] = True
else:
# constant term not present; arbitrarily use first residue instead
mask_idx_constrained[0] = True
A22 = A[1:, ~mask_idx_constrained]
b2 = freq_responses[:, 1:]
A22_ri = np.vstack((A22.real, A22.imag))
b22_ri = np.hstack((b2.real, b2.imag))
logger.info(f'Condition number of coefficient matrix = {int(np.linalg.cond(A22_ri))}')
# solve least-squares and obtain results as stack of real part vector and imaginary part vector
x2, residuals, rank, singular_vals = np.linalg.lstsq(A22_ri, b22_ri.T, rcond=None)
# solve for x1 using the first row (the dc row):
b1 = freq_responses[:, 0]
A11 = A[0, mask_idx_constrained]
A12 = A[0, ~mask_idx_constrained]
x1 = np.real(1 / A11 * (b1 - np.dot(A12, x2)))
# reassemble x from x1 and x2
x = np.empty((n_cols, n_responses))
x[mask_idx_constrained, :] = x1
x[~mask_idx_constrained, :] = x2
else:
# dc point not included; use and solve the entire linear system with least-squares
A_ri = np.vstack((A.real, A.imag))
b_ri = np.hstack((freq_responses.real, freq_responses.imag))
logger.info(f'Condition number of coefficient matrix = {int(np.linalg.cond(A_ri))}')
# solve least-squares and obtain results as stack of real part vector and imaginary part vector
x, residuals, rank, singular_vals = np.linalg.lstsq(A_ri, b_ri.T, rcond=None)
x = scaling[:, None] * x
# extract residues from solution vector and align them with poles to get matching pole-residue pairs
residues = np.empty((len(freq_responses), len(poles)), dtype=complex)
residues[:, idx_poles_real] = np.transpose(x[idx_res_real])
residues[:, idx_poles_complex] = np.transpose(x[idx_res_complex_re] + 1j * x[idx_res_complex_im])
# extract constant and proportional coefficient, if available
if fit_constant:
constant_coeff = x[idx_constant][0]
else:
constant_coeff = np.zeros(n_responses)
if fit_proportional:
proportional_coeff = x[idx_proportional][0]
else:
proportional_coeff = np.zeros(n_responses)
return residues, constant_coeff, proportional_coeff, residuals, rank, singular_vals
@staticmethod
def _get_delta(poles, residues, constant_coeff, proportional_coeff, freqs, freq_responses, weights_responses):
s = 2j * np.pi * freqs
model = proportional_coeff[:, None] * s + constant_coeff[:, None]
for i, pole in enumerate(poles):
if np.imag(pole) == 0.0:
# real pole
model += residues[:, i, None] / (s - pole)
else:
# complex conjugate pole
model += (residues[:, i, None] / (s - pole) +
np.conjugate(residues[:, i, None]) / (s - np.conjugate(pole)))
# compute weighted error and return global maximum at each frequency across all individual responses
delta = np.abs(model - freq_responses) * weights_responses[:, None]
return np.max(delta, axis=0)
@staticmethod
def _find_error_bands(freqs, delta):
# compute error bands (maximal fit deviation)
delta_mean = np.mean(delta)
error = delta - delta_mean
# find limits of error bands
idx_limits = np.nonzero(np.diff(error > 0))[0]
idx_limits_filtered = idx_limits[np.diff(idx_limits, prepend=0) > 2]
freqs_bands = np.split(freqs, idx_limits_filtered)
error_bands = np.split(error, idx_limits_filtered)
n_bands = len(freqs_bands)
idx_freqs_start = []
idx_freqs_stop = []
idx_freqs_max = []
delta_mean_bands = []
for i_band in range(n_bands):
band_error_mean = np.mean(error_bands[i_band])
if band_error_mean > 0:
# band with excess error;
# find frequency index of error maximum inside this band
i_band_max_error = np.argmax(error_bands[i_band])
i_start = np.nonzero(freqs == freqs_bands[i_band][0])[0][0]
i_stop = np.nonzero(freqs == freqs_bands[i_band][-1])[0][0]
i_max = np.nonzero(freqs == freqs_bands[i_band][i_band_max_error])[0][0]
idx_freqs_start.append(i_start)
idx_freqs_stop.append(i_stop)
idx_freqs_max.append(i_max)
delta_mean_bands.append(np.mean(delta[i_start:i_stop]))
idx_freqs_start = np.array(idx_freqs_start)
idx_freqs_stop = np.array(idx_freqs_stop)
idx_freqs_max = np.array(idx_freqs_max)
delta_mean_bands = np.array(delta_mean_bands)
i_sort = np.flip(np.argsort(delta_mean_bands))
return idx_freqs_start[i_sort], idx_freqs_stop[i_sort], idx_freqs_max[i_sort], delta_mean_bands[i_sort]
def get_rms_error(self, i=-1, j=-1, parameter_type: str = 's'):
r"""
Returns the root-mean-square (rms) error magnitude of the fit, i.e.
:math:`\sqrt{ \mathrm{mean}(|S - S_\mathrm{fit} |^2) }`,
either for an individual response :math:`S_{i+1,j+1}` or for larger slices of the network.
Parameters
----------
i : int, optional
Row indices of the responses to be evaluated. Either a single row selected by an integer
:math:`i \in [0, N_\mathrm{ports}-1]`, or multiple rows selected by a list of integers, or all rows
selected by :math:`i = -1` (*default*).
j : int, optional
Column indices of the responses to be evaluated. Either a single column selected by an integer
:math:`j \in [0, N_\mathrm{ports}-1]`, or multiple columns selected by a list of integers, or all columns
selected by :math:`j = -1` (*default*).
parameter_type: str, optional
Representation type of the fitted frequency responses. Either *scattering* (:attr:`s` or :attr:`S`),
*impedance* (:attr:`z` or :attr:`Z`) or *admittance* (:attr:`y` or :attr:`Y`).
Returns
-------
rms_error : ndarray
The rms error magnitude between the vector fitted model and the original network data.
Raises
------
ValueError
If the specified parameter representation type is not :attr:`s`, :attr:`z`, nor :attr:`y`.
"""
if i == -1:
list_i = range(self.network.nports)
elif isinstance(i, int):
list_i = [i]
else:
list_i = i
if j == -1:
list_j = range(self.network.nports)
elif isinstance(j, int):
list_j = [j]
else:
list_j = j
if parameter_type.lower() == 's':
nw_responses = self.network.s
elif parameter_type.lower() == 'z':
nw_responses = self.network.z
elif parameter_type.lower() == 'y':
nw_responses = self.network.y
else:
raise ValueError(f'Invalid parameter type `{parameter_type}`. Valid options: `s`, `z`, or `y`')
error_mean_squared = 0
for i in list_i:
for j in list_j:
nw_ij = nw_responses[:, i, j]
fit_ij = self.get_model_response(i, j, self.network.f)
error_mean_squared += np.mean(np.square(np.abs(nw_ij - fit_ij)))
return np.sqrt(error_mean_squared)
def _get_ABCDE(self) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""
Private method.
Returns the real-valued system matrices of the state-space representation of the current rational model, as
defined in [#]_.
Returns
-------
A : ndarray
State-space matrix A holding the poles on the diagonal as real values with imaginary parts on the sub-
diagonal
B : ndarray
State-space matrix B holding coefficients (1, 2, or 0), depending on the respective type of pole in A
C : ndarray
State-space matrix C holding the residues
D : ndarray
State-space matrix D holding the constants
E : ndarray
State-space matrix E holding the proportional coefficients (usually 0 in case of fitted S-parameters)
Raises
------
ValueError
If the model parameters have not been initialized (by running :func:`vector_fit()` or :func:`read_npz()`).
References
----------
.. [#] B. Gustavsen and A. Semlyen, "Fast Passivity Assessment for S-Parameter Rational Models Via a Half-Size
Test Matrix," in IEEE Transactions on Microwave Theory and Techniques, vol. 56, no. 12, pp. 2701-2708,
Dec. 2008, DOI: 10.1109/TMTT.2008.2007319.
"""
# initial checks
if self.poles is None:
raise ValueError('self.poles = None; nothing to do. You need to run vector_fit() first.')
if self.residues is None:
raise ValueError('self.residues = None; nothing to do. You need to run vector_fit() first.')
if self.proportional_coeff is None:
raise ValueError('self.proportional_coeff = None; nothing to do. You need to run vector_fit() first.')
if self.constant_coeff is None:
raise ValueError('self.constant_coeff = None; nothing to do. You need to run vector_fit() first.')
# assemble real-valued state-space matrices A, B, C, D, E from fitted complex-valued pole-residue model
# determine size of the matrix system
n_ports = int(np.sqrt(len(self.constant_coeff)))
n_poles_real = 0
n_poles_cplx = 0
for pole in self.poles:
if np.imag(pole) == 0.0:
n_poles_real += 1
else:
n_poles_cplx += 1
n_matrix = (n_poles_real + 2 * n_poles_cplx) * n_ports
# state-space matrix A holds the poles on the diagonal as real values with imaginary parts on the sub-diagonal
# state-space matrix B holds coefficients (1, 2, or 0), depending on the respective type of pole in A
# assemble A = [[poles_real, 0, 0],
# [0, real(poles_cplx), imag(poles_cplx],
# [0, -imag(poles_cplx), real(poles_cplx]]
A = np.identity(n_matrix)
B = np.zeros(shape=(n_matrix, n_ports))
i_A = 0 # index on diagonal of A
for j in range(n_ports):
for pole in self.poles:
if np.imag(pole) == 0.0:
# adding a real pole
A[i_A, i_A] = np.real(pole)
B[i_A, j] = 1
i_A += 1
else:
# adding a complex-conjugate pole
A[i_A, i_A] = np.real(pole)
A[i_A, i_A + 1] = np.imag(pole)
A[i_A + 1, i_A] = -1 * np.imag(pole)
A[i_A + 1, i_A + 1] = np.real(pole)
B[i_A, j] = 2
i_A += 2
# state-space matrix C holds the residues
# assemble C = [[R1.11, R1.12, R1.13, ...], [R2.11, R2.12, R2.13, ...], ...]
C = np.zeros(shape=(n_ports, n_matrix))
for i in range(n_ports):
for j in range(n_ports):
# i: row index
# j: column index
i_response = i * n_ports + j
j_residues = 0
for zero in self.residues[i_response]:
if np.imag(zero) == 0.0:
C[i, j * (n_poles_real + 2 * n_poles_cplx) + j_residues] = np.real(zero)
j_residues += 1
else:
C[i, j * (n_poles_real + 2 * n_poles_cplx) + j_residues] = np.real(zero)
C[i, j * (n_poles_real + 2 * n_poles_cplx) + j_residues + 1] = np.imag(zero)
j_residues += 2
# state-space matrix D holds the constants
# assemble D = [[d11, d12, ...], [d21, d22, ...], ...]
D = np.zeros(shape=(n_ports, n_ports))
for i in range(n_ports):
for j in range(n_ports):
# i: row index
# j: column index
i_response = i * n_ports + j
D[i, j] = self.constant_coeff[i_response]
# state-space matrix E holds the proportional coefficients (usually 0 in case of fitted S-parameters)
# assemble E = [[e11, e12, ...], [e21, e22, ...], ...]
E = np.zeros(shape=(n_ports, n_ports))
for i in range(n_ports):
for j in range(n_ports):
# i: row index
# j: column index
i_response = i * n_ports + j
E[i, j] = self.proportional_coeff[i_response]
return A, B, C, D, E
@staticmethod
def _get_s_from_ABCDE(freqs: np.ndarray,
A: np.ndarray, B: np.ndarray, C: np.ndarray, D: np.ndarray, E: np.ndarray) -> np.ndarray:
"""
Private method.
Returns the S-matrix of the vector fitted model calculated from the real-valued system matrices of the state-
space representation, as provided by `_get_ABCDE()`.
Parameters
----------
freqs : ndarray
Frequencies (in Hz) at which to calculate the S-matrices.
A : ndarray
B : ndarray
C : ndarray
D : ndarray
E : ndarray
Returns
-------
ndarray
Complex-valued S-matrices (fxNxN) calculated at frequencies `freqs`.
"""
dim_A = np.shape(A)[0]
stsp_poles = np.linalg.inv(2j * np.pi * freqs[:, None, None] * np.identity(dim_A)[None, :, :] - A[None, :, :])
stsp_S = np.matmul(np.matmul(C, stsp_poles), B)
stsp_S += D + 2j * np.pi * freqs[:, None, None] * E
return stsp_S
def passivity_test(self, parameter_type: str = 's') -> np.ndarray:
"""
Evaluates the passivity of reciprocal vector fitted models by means of a half-size test matrix [#]_. Any
existing frequency bands of passivity violations will be returned as a sorted list.
Parameters
----------
parameter_type: str, optional
Representation type of the fitted frequency responses. Either *scattering* (:attr:`s` or :attr:`S`),
*impedance* (:attr:`z` or :attr:`Z`) or *admittance* (:attr:`y` or :attr:`Y`). Currently, only scattering
parameters are supported for passivity evaluation.
Raises
------
NotImplementedError
If the function is called for `parameter_type` different than `S` (scattering).
ValueError
If the function is used with a model containing nonzero proportional coefficients.
Returns
-------
violation_bands : ndarray
NumPy array with frequency bands of passivity violation:
`[[f_start_1, f_stop_1], [f_start_2, f_stop_2], ...]`.
See Also
--------
is_passive : Query the model passivity as a boolean value.
passivity_enforce : Enforces the passivity of the vector fitted model, if required.
Examples
--------
Load and fit the `Network`, then evaluate the model passivity:
>>> nw_3port = skrf.Network('my3port.s3p')
>>> vf = skrf.VectorFitting(nw_3port)
>>> vf.vector_fit(n_poles_real=1, n_poles_cmplx=4)
>>> violations = vf.passivity_test()
References
----------
.. [#] B. Gustavsen and A. Semlyen, "Fast Passivity Assessment for S-Parameter Rational Models Via a Half-Size
Test Matrix," in IEEE Transactions on Microwave Theory and Techniques, vol. 56, no. 12, pp. 2701-2708,
Dec. 2008, DOI: 10.1109/TMTT.2008.2007319.
"""
if parameter_type.lower() != 's':
raise NotImplementedError('Passivity testing is currently only supported for scattering (S) parameters.')
if parameter_type.lower() == 's' and len(np.flatnonzero(self.proportional_coeff)) > 0:
raise ValueError('Passivity testing of scattering parameters with nonzero proportional coefficients does '
'not make any sense; you need to run vector_fit() with option `fit_proportional=False` '
'first.')
# # the network needs to be reciprocal for this passivity test method to work: S = transpose(S)
# if not np.allclose(self.residues, np.transpose(self.residues)) or \
# not np.allclose(self.constant_coeff, np.transpose(self.constant_coeff)) or \
# not np.allclose(self.proportional_coeff, np.transpose(self.proportional_coeff)):
# logger.error('Passivity testing with unsymmetrical model parameters is not supported. '
# 'The model needs to be reciprocal.')
# return
# get state-space matrices
A, B, C, D, E = self._get_ABCDE()
n_ports = np.shape(D)[0]
# build half-size test matrix P from state-space matrices A, B, C, D
inv_neg = np.linalg.inv(D - np.identity(n_ports))
inv_pos = np.linalg.inv(D + np.identity(n_ports))
prod_neg = np.matmul(np.matmul(B, inv_neg), C)
prod_pos = np.matmul(np.matmul(B, inv_pos), C)
P = np.matmul(A - prod_neg, A - prod_pos)
# extract eigenvalues of P
P_eigs = np.linalg.eigvals(P)
# purely imaginary square roots of eigenvalues identify frequencies (2*pi*f) of borders of passivity violations
freqs_violation = []
for sqrt_eigenval in np.sqrt(P_eigs):
if np.real(sqrt_eigenval) == 0.0:
freqs_violation.append(np.imag(sqrt_eigenval) / 2 / np.pi)
# include dc (0) unless it's already included
if len(np.nonzero(np.array(freqs_violation) == 0.0)[0]) == 0:
freqs_violation.append(0.0)
# sort the output from lower to higher frequencies
freqs_violation = np.sort(freqs_violation)
# identify frequency bands of passivity violations
# sweep the bands between crossover frequencies and identify bands of passivity violations
violation_bands = []
for i, freq in enumerate(freqs_violation):
if i == len(freqs_violation) - 1:
# last band stops always at infinity
f_start = freq
f_stop = np.inf
f_center = 1.1 * f_start # 1.1 is chosen arbitrarily to have any frequency for evaluation
else:
# intermediate band between this frequency and the previous one
f_start = freq
f_stop = freqs_violation[i + 1]
f_center = 0.5 * (f_start + f_stop)
# calculate singular values at the center frequency between crossover frequencies to identify violations
s_center = self._get_s_from_ABCDE(np.array([f_center]), A, B, C, D, E)
sigma = np.linalg.svd(s_center[0], compute_uv=False)
passive = True
for singval in sigma:
if singval > 1:
# passivity violation in this band
passive = False
if not passive:
# add this band to the list of passivity violations
if violation_bands is None:
violation_bands = [[f_start, f_stop]]
else:
violation_bands.append([f_start, f_stop])
return np.array(violation_bands)
def is_passive(self, parameter_type: str = 's') -> bool:
"""
Returns the passivity status of the model as a boolean value.
Parameters
----------
parameter_type : str, optional
Representation type of the fitted frequency responses. Either *scattering* (:attr:`s` or :attr:`S`),
*impedance* (:attr:`z` or :attr:`Z`) or *admittance* (:attr:`y` or :attr:`Y`). Currently, only scattering
parameters are supported for passivity evaluation.
Returns
-------
passivity : bool
:attr:`True` if model is passive, else :attr:`False`.
See Also
--------
passivity_test : Verbose passivity evaluation routine.
passivity_enforce : Enforces the passivity of the vector fitted model, if required.
Examples
--------
Load and fit the `Network`, then check whether or not the model is passive:
>>> nw_3port = skrf.Network('my3port.s3p')
>>> vf = skrf.VectorFitting(nw_3port)
>>> vf.vector_fit(n_poles_real=1, n_poles_cmplx=4)
>>> vf.is_passive() # returns True or False
"""
viol_bands = self.passivity_test(parameter_type)
if len(viol_bands) == 0:
return True
else:
return False
def passivity_enforce(self, n_samples: int = 200, f_max: float = None, parameter_type: str = 's',
preserve_dc: bool = True) -> None:
"""
Enforces the passivity of the vector fitted model, if required. This is an implementation of the methods
presented in [#]_ and [#]_ using singular value perturbation. To preserve the dc point in the model during
passivity enforcement, only the residues are perturbed, not the constant term.
Parameters
----------
n_samples : int, optional
Number of linearly spaced frequency samples at which passivity will be evaluated and enforced.
(Default: 200). If there are very narrow frequency bands of passivity violations, a sufficiently large
number of frequency samples is required.
f_max : float or None, optional
Highest frequency of interest for the passivity enforcement (in Hz, not rad/s). This limit usually
equals the highest sample frequency of the fitted Network. If None, the highest frequency in
:attr:`self.network` is used, which must not be None is this case. If `f_max` is not None, it overrides the
highest frequency in :attr:`self.network`.
parameter_type : str, optional
Representation type of the fitted frequency responses. Either *scattering* (:attr:`s` or :attr:`S`),
*impedance* (:attr:`z` or :attr:`Z`) or *admittance* (:attr:`y` or :attr:`Y`). Currently, only scattering
parameters are supported for passivity evaluation.
preserve_dc : bool, optional
Enables dc point preservation during passivity enforcement. This only works if the fitted model is already
passive at the dc point, which is not always the case. If it is not passive, dc point preservation is
disabled and passivity is also enforced on the dc point.
Returns
-------
None
Raises
------
NotImplementedError
If the function is called for `parameter_type` different than `S` (scattering).
ValueError
If the function is used with a model containing nonzero proportional coefficients. Or if both `f_max` and
:attr:`self.network` are None.
See Also
--------
is_passive : Returns the passivity status of the model as a boolean value.
passivity_test : Verbose passivity evaluation routine.
plot_passivation : Convergence plot for passivity enforcement iterations.
Examples
--------
Load and fit the `Network`, then enforce the passivity of the model:
>>> nw_3port = skrf.Network('my3port.s3p')
>>> vf = skrf.VectorFitting(nw_3port)
>>> vf.vector_fit(n_poles_real=1, n_poles_cmplx=4)
>>> vf.passivity_enforce() # won't do anything if model is already passive
References
----------
.. [#] T. Dhaene, D. Deschrijver and N. Stevens, "Efficient Algorithm for Passivity Enforcement of S-Parameter-
Based Macromodels," in IEEE Transactions on Microwave Theory and Techniques, vol. 57, no. 2, pp. 415-420,
Feb. 2009, DOI: 10.1109/TMTT.2008.2011201
.. [#] D. Deschrijver and T. Dhaene, "DC-Preserving Passivity Enforcement for S-Parameter Based Macromodels,"
in IEEE Transactions on Microwave Theory and Techniques, vol. 58, no. 4, pp. 923-928, April 2010,
DOI: 10.1109/TMTT.2010.2042556
"""
if parameter_type.lower() != 's':
raise NotImplementedError('Passivity testing is currently only supported for scattering (S) parameters.')
if parameter_type.lower() == 's' and len(np.flatnonzero(self.proportional_coeff)) > 0:
raise ValueError('Passivity testing of scattering parameters with nonzero proportional coefficients does '
'not make any sense; you need to run vector_fit() with option `fit_proportional=False` '
'first.')
# always run passivity test first; this will write 'self.violation_bands'
if self.is_passive():
# model is already passive; do nothing and return
logger.info('Passivity enforcement: The model is already passive. Nothing to do.')
return
# check dc passivity and find the highest relevant frequency; either
# 1) the highest frequency of passivity violation (f_viol_max)
# or
# 2) the highest fitting frequency (f_samples_max)
violation_bands = self.passivity_test()
f_viol_min = violation_bands[0, 0]
f_viol_max = violation_bands[-1, 1]
# check passivity at the dc point; 1) in the model, 2) in the original data, if available
if preserve_dc and f_viol_min == 0.0:
# cannot preserve a non-passive dc point during passivity enforcement
preserve_dc = False
hint = ''
if self.network is not None:
if self.network.f[0] == 0.0 and not self.network.is_passive():
hint = '\nHint: The dc point in the original network data is already non-passive.'
warnings.warn('Passivity enforcement: The dc point in the model is not passive. Cannot '
f'preserve the dc point during passivity enforcement. {hint}', UserWarning, stacklevel=2)
if f_max is None:
if self.network is None:
raise RuntimeError('Both `self.network` and parameter `f_max` are None. One of them is required to '
'specify the frequency band of interest for the passivity enforcement.')
else:
f_samples_max = self.network.f[-1]
else:
f_samples_max = f_max
# deal with unbounded violation interval (f_viol_max == np.inf)
if np.isinf(f_viol_max):
f_viol_max = 1.5 * violation_bands[-1, 0]
warnings.warn(
'Passivity enforcement: The passivity violations of this model are unbounded. '
'Passivity enforcement might still work, but consider re-fitting with a lower number of poles '
'and/or without the constants (`fit_constant=False`) if the results are not satisfactory.',
UserWarning, stacklevel=2)
# the frequency band for the passivity evaluation is from dc to 20% above the highest relevant frequency
if f_viol_max < f_samples_max:
f_eval_max = 1.2 * f_samples_max
else:
f_eval_max = 1.2 * f_viol_max
# let's not automatically adjust n_samples. The calculated number can
# be huge (>100k). Combined with a high number of poles in the model, this can bust the memory.
freqs_eval = np.linspace(0, f_eval_max, n_samples)
# get model state-space matrices
A, B, C_t, D, E = self._get_ABCDE()
dim_A = np.shape(A)[0]
# ASYMPTOTIC PASSIVITY ENFORCEMENT
# check if constant term has been fitted (not zero)
# a model without the constant term is always asymptotically passive
if len(np.nonzero(D)[0]) != 0:
# D was fitted;
# asymptotic passivity needs to be checked and enforced, if violated.
# for dc preservation, the asymptotic passivity violations in D are compensated using C
# D is not touched, because it contains the dc point ( lim s --> {inf S(s)} = D)
u, sigma, vh = np.linalg.svd(D, compute_uv=True)
# find and perturb singular values that cause passivity violations
# sigma_viol = sigma * upsilon - psi with
# upsilon[sigma > delta] = 1
# upsilon[sigma <= delta] = 0
# psi[sigma > delta] = delta
# psi[sigma <= delta] = 0
# (implemented below in a more compact form)
delta = 1
idx_viol = np.nonzero(sigma > delta)
sigma_viol = np.zeros_like(sigma)
sigma_viol[idx_viol] = sigma[idx_viol] - delta
# calculate S_viol from perturbed sigma and previous U and Vh
S_viol = np.dot(u * sigma_viol, vh)
# find new set of residues C_viol by solving underdetermined least-squares problem
# S_viol = C_viol * B
#
# mind the transpose of the system to compensate for the exchanged order of matrix multiplication:
# S_viol = C_viol * B <==> transpose(S_viol) = transpose(B) * transpose(C_viol)
C_viol, residuals, rank, singular_vals = np.linalg.lstsq(np.vstack((B.T.real, B.T.imag)),
np.vstack((S_viol.T.real, S_viol.T.imag)),
rcond=None)
C_t -= C_viol.T
# UNIFORM PASSIVITY ENFORCEMENT
# preparing coefficient matrix; can be reused in every iteration
# S(s_eval) = D_t + s_eval * C_t * inv(s_eval * I - A) * B
# = D_t + s_eval * C_t * A_freq * B
# with
# A_freq = inv(s_eval * I - A)
# s_eval = j * omega_eval = 2j * pi * freqs_eval
A_freq = np.linalg.inv(2j * np.pi * freqs_eval[:, None, None] * np.identity(dim_A)[None, :, :] - A[None, :, :])
# construct coefficient matrix for least-squares residue fitting (C_viol)
coeffs = np.matmul(A_freq, B)
C_viol = np.empty_like(C_t)
n_ports = np.shape(C_viol)[0]
model_order = self.get_model_order(self.poles)
# predefined tolerance parameter (users should not need to change this)
delta_threshold = 0.999
sigma_max = 1.1 # just to enter iteration loop for the first time
# iterative compensation of passivity violations
t = 0
self.history_max_sigma = []
while t < self.max_iterations and sigma_max > 1.0:
logger.info(f'Passivity enforcement; Iteration {t + 1}')
# calculate S-matrix of the model at freqs_eval (shape fxNxN)
#S_eval = self._get_s_from_ABCDE(freqs_eval, A, B, C_t, D, E)
S_eval = D + np.matmul(C_t, coeffs) # much faster!
# singular value decomposition,
# shape(u) = (n_samples, n_ports, n_ports)
# shape(sigma) = (n_samples, n_ports)
# shape(vh) = (n_samples, n_ports, n_ports)
u, sigma, vh = np.linalg.svd(S_eval)
# keep track of the greatest singular value in every iteration step
sigma_max = np.amax(sigma)
self.history_max_sigma.append(sigma_max)
if sigma_max > delta_threshold:
delta = delta_threshold
else:
delta = sigma_max
# find and perturb singular values that cause passivity violations
# sigma_viol = sigma * upsilon - psi with
# upsilon[sigma > delta] = 1
# upsilon[sigma <= delta] = 0
# psi[sigma > delta] = delta
# psi[sigma <= delta] = 0
# (implemented below in a more compact form)
idx_viol = np.nonzero(sigma > delta)
sigma_viol = np.zeros_like(sigma)
sigma_viol[idx_viol] = sigma[idx_viol] - delta
S_viol = np.matmul(u * sigma_viol[:, None, :], vh)
# stack frequency responses as a single vector
# stacking order (row-major):
# s11, s12, s13, ..., s21, s22, s23, ...
S_viol_stacked = []
for i in range(n_ports):
for j in range(n_ports):
S_viol_stacked.append(S_viol[:, i, j])
S_viol_stacked = np.array(S_viol_stacked)
# The existing method _fit_residues() can be use here to fit the violation residues. Enabling `fit_constant`
# in combination with `enforce_dc` removes the dc rows from the linear system and enforces the dc solution
# on the constant term. In case of dc preservation during passivity enforcement, we can ignore that constant
# term entirely and only use the violation residues.
# If dc preservation is disabled, we could also perturb the constant term. This is not currently done. In
# this new method, we always only perturb the residues. Disabling `fit_constant` and `preserve_dc` in this
# case will solve for the residues without the constant term in the linear system.
C_viol_stacked, D_viol_stacked, E_viol_stacked, residuals, rank, singular_vals = self._fit_residues(
self.poles, freqs_eval, S_viol_stacked, fit_constant=preserve_dc, fit_proportional=False,
enforce_dc=preserve_dc)
# reshape C_viol into state-space format: [[R1.11, R2.11, R3.11, ..., R1.1N, R2.1N, R3.1N, ...],
# [R1.21, R2.21, R3.21, ..., R1.2N, R3.2N, R3.2N, ...],
# ...
# [R1.N1, R2.N1, R3.N1, ..., R1.NN, R3.NN, R3.NN, ...]]
for i_port in range(n_ports):
for j_port in range(n_ports):
j_residues = 0
for residue in C_viol_stacked[i_port * n_ports + j_port]:
if np.imag(residue) == 0.0:
C_viol[i_port, j_port * model_order + j_residues] = np.real(residue)
j_residues += 1
else:
C_viol[i_port, j_port * model_order + j_residues] = np.real(residue)
C_viol[i_port, j_port * model_order + j_residues + 1] = np.imag(residue)
j_residues += 2
# perturb residues by subtracting respective row and column in C_t
C_t = C_t - C_viol
t += 1
# PASSIVATION PROCESS DONE; model is either passive or max. number of iterations have been exceeded
if t == self.max_iterations:
warnings.warn('Passivity enforcement: Aborting after the max. number of iterations has been '
'exceeded.', RuntimeWarning, stacklevel=2)
# save/update model parameters (perturbed residues)
self.history_max_sigma = np.array(self.history_max_sigma)
n_ports = np.shape(D)[0]
for i in range(n_ports):
k = 0 # column index in C_t
for j in range(n_ports):
i_response = i * n_ports + j
z = 0 # column index self.residues
for pole in self.poles:
if np.imag(pole) == 0.0:
# real pole --> real residue
self.residues[i_response, z] = C_t[i, k]
k += 1
else:
# complex-conjugate pole --> complex-conjugate residue
self.residues[i_response, z] = C_t[i, k] + 1j * C_t[i, k + 1]
k += 2
z += 1
# run final passivity test to make sure passivation was successful
violation_bands = self.passivity_test()
if len(violation_bands) > 0:
# trying to determine the required number of evaluation samples based on the bandwidth and separation
# distance of the violation bands
violation_band_separation = np.diff(violation_bands.flat)
min_spacing_nonzero = np.amin(violation_band_separation[violation_band_separation != 0.0])
# we should need an absolute minimum of 1 sample in each violating frequency band.
# in practice, the frequency spacing should preferrably be much more dense.
# let's recommend 2 samples per violation band.
n_samples_required = int(f_eval_max / min_spacing_nonzero * 2)
if n_samples_required > n_samples:
hint = f'Consider trying again with n_samples > {n_samples_required}.'
else:
hint = ''
warnings.warn('Passivity enforcement was not successful.\nModel is still non-passive in these '
f'frequency bands: {violation_bands}.\nTry running this routine again with a larger number of'
f' samples (parameter `n_samples`). This run was using n_samples = {n_samples}. {hint}',
RuntimeWarning, stacklevel=2)
def write_npz(self, path: str) -> None:
"""
Writes the model parameters in :attr:`poles`, :attr:`residues`,
:attr:`proportional_coeff` and :attr:`constant_coeff` to a labeled NumPy .npz file.
Parameters
----------
path : str
Target path without filename for the export. The filename will be added automatically based on the network
name in :attr:`network`
Returns
-------
None
See Also
--------
read_npz : Reads all model parameters from a .npz file
Examples
--------
Load and fit the `Network`, then export the model parameters to a .npz file:
>>> nw_3port = skrf.Network('my3port.s3p')
>>> vf = skrf.VectorFitting(nw_3port)
>>> vf.vector_fit(n_poles_real=1, n_poles_cmplx=4)
>>> vf.write_npz('./data/')
The filename depends on the network name stored in `nw_3port.name` and will have the prefix `coefficients_`, for
example `coefficients_my3port.npz`. The coefficients can then be read using NumPy's load() function:
>>> coeffs = numpy.load('./data/coefficients_my3port.npz')
>>> poles = coeffs['poles']
>>> residues = coeffs['residues']
>>> prop_coeffs = coeffs['proportionals']
>>> constants = coeffs['constants']
Alternatively, the coefficients can be read directly into a new instance of `VectorFitting`, see
:func:`read_npz`.
"""
if self.poles is None:
warnings.warn('Nothing to export; Poles have not been fitted.', RuntimeWarning, stacklevel=2)
return
if self.residues is None:
warnings.warn('Nothing to export; Residues have not been fitted.', RuntimeWarning, stacklevel=2)
return
if self.proportional_coeff is None:
warnings.warn('Nothing to export; Proportional coefficients have not been fitted.', RuntimeWarning,
stacklevel=2)
return
if self.constant_coeff is None:
warnings.warn('Nothing to export; Constants have not been fitted.', RuntimeWarning, stacklevel=2)
return
filename = self.network.name
logger.info(f'Exporting results as compressed NumPy array to {path}')
np.savez_compressed(os.path.join(path, f'coefficients_{filename}'),
poles=self.poles, residues=self.residues, proportionals=self.proportional_coeff,
constants=self.constant_coeff)
def read_npz(self, file: str) -> None:
"""
Reads all model parameters :attr:`poles`, :attr:`residues`, :attr:`proportional_coeff` and
:attr:`constant_coeff` from a labeled NumPy .npz file.
Parameters
----------
file : str
NumPy .npz file containing the parameters. See notes.
Returns
-------
None
Raises
------
ValueError
If the shapes of the coefficient arrays in the provided file are not compatible.
Notes
-----
The .npz file needs to include the model parameters as individual NumPy arrays (ndarray) labeled '*poles*',
'*residues*', '*proportionals*' and '*constants*'. The shapes of those arrays need to match the network
properties in :class:`network` (correct number of ports). Preferably, the .npz file was created by
:func:`write_npz`.
See Also
--------
write_npz : Writes all model parameters to a .npz file
Examples
--------
Create an empty `VectorFitting` instance (with or without the fitted `Network`) and load the model parameters:
>>> vf = skrf.VectorFitting(None)
>>> vf.read_npz('./data/coefficients_my3port.npz')
This can be useful to analyze or process a previous vector fit instead of fitting it again, which sometimes
takes a long time. For example, the model passivity can be evaluated and enforced:
>>> vf.passivity_enforce()
"""
with np.load(file) as data:
poles = data['poles']
# legacy support for exported residues
if 'zeros' in data:
# old .npz file from deprecated write_npz() with residues called 'zeros'
residues = data['zeros']
else:
# new .npz file from current write_npz()
residues = data['residues']
proportional_coeff = data['proportionals']
constant_coeff = data['constants']
n_ports = int(np.sqrt(len(constant_coeff)))
n_resp = n_ports ** 2
if np.shape(residues)[0] == np.shape(proportional_coeff)[0] == np.shape(constant_coeff)[0] == n_resp:
self.poles = poles
self.residues = residues
self.proportional_coeff = proportional_coeff
self.constant_coeff = constant_coeff
else:
raise ValueError('The shapes of the provided parameters are not compatible. The coefficient file needs '
'to contain NumPy arrays labled `poles`, `residues`, `proportionals`, and '
'`constants`. Their shapes must match the number of network ports and the number of '
'frequencies.')
def get_model_response(self, i: int, j: int, freqs: Any = None) -> np.ndarray:
"""
Returns one of the frequency responses :math:`H_{i+1,j+1}` of the fitted model :math:`H`.
Parameters
----------
i : int
Row index of the response in the response matrix.
j : int
Column index of the response in the response matrix.
freqs : list of float or ndarray or None, optional
List of frequencies for the response plot. If None, the sample frequencies of the fitted network in
:attr:`network` are used.
Returns
-------
response : ndarray
Model response :math:`H_{i+1,j+1}` at the frequencies specified in `freqs` (complex-valued Numpy array).
Examples
--------
Get fitted S11 at 101 frequencies from 0 Hz to 10 GHz:
>>> import skrf
>>> vf = skrf.VectorFitting(skrf.data.ring_slot)
>>> vf.vector_fit(3, 0)
>>> s11_fit = vf.get_model_response(0, 0, numpy.linspace(0, 10e9, 101))
"""
if self.poles is None:
warnings.warn('Returning a zero-vector; Poles have not been fitted.',
RuntimeWarning, stacklevel=2)
return np.zeros_like(freqs)
if self.residues is None:
warnings.warn('Returning a zero-vector; Residues have not been fitted.',
RuntimeWarning, stacklevel=2)
return np.zeros_like(freqs)
if self.proportional_coeff is None:
warnings.warn('Returning a zero-vector; Proportional coefficients have not been fitted.',
RuntimeWarning, stacklevel=2)
return np.zeros_like(freqs)
if self.constant_coeff is None:
warnings.warn('Returning a zero-vector; Constants have not been fitted.',
RuntimeWarning, stacklevel=2)
return np.zeros_like(freqs)
if freqs is None:
freqs = np.linspace(np.amin(self.network.f), np.amax(self.network.f), 1000)
s = 2j * np.pi * np.array(freqs)
n_ports = int(np.sqrt(len(self.constant_coeff)))
i_response = i * n_ports + j
residues = self.residues[i_response]
resp = self.proportional_coeff[i_response] * s + self.constant_coeff[i_response]
for i, pole in enumerate(self.poles):
if np.imag(pole) == 0.0:
# real pole
resp += residues[i] / (s - pole)
else:
# complex conjugate pole
resp += residues[i] / (s - pole) + np.conjugate(residues[i]) / (s - np.conjugate(pole))
return resp
@axes_kwarg
def plot(self, component: str, i: int = -1, j: int = -1, freqs: Any = None,
parameter: str = 's', *, ax: Axes = None) -> Axes:
"""
Plots the specified component of the parameter :math:`H_{i+1,j+1}` in the fit, where :math:`H` is
either the scattering (:math:`S`), the impedance (:math:`Z`), or the admittance (:math:`H`) response specified
in `parameter`.
Parameters
----------
component : str
The component to be plotted. Must be one of the following items:
['db', 'mag', 'deg', 'deg_unwrap', 're', 'im'].
`db` for magnitude in decibels,
`mag` for magnitude in linear scale,
`deg` for phase in degrees (wrapped),
`deg_unwrap` for phase in degrees (unwrapped/continuous),
`re` for real part in linear scale,
`im` for imaginary part in linear scale.
i : int, optional
Row index of the response. `-1` to plot all rows.
j : int, optional
Column index of the response. `-1` to plot all columns.
freqs : list of float or ndarray or None, optional
List of frequencies for the response plot. If None, the sample frequencies of the fitted network in
:attr:`network` are used. This only works if :attr:`network` is not `None`.
parameter : str, optional
The network representation to be used. This is only relevant for the plot of the original sampled response
in :attr:`network` that is used for comparison with the fit. Must be one of the following items unless
:attr:`network` is `None`: ['s', 'z', 'y'] for *scattering* (default), *impedance*, or *admittance*.
ax : :class:`matplotlib.Axes` object or None
matplotlib axes to draw on. If None, the current axes is fetched with :func:`gca()`.
Returns
-------
:class:`matplotlib.Axes`
matplotlib axes used for drawing. Either the passed :attr:`ax` argument or the one fetch from the current
figure.
Raises
------
ValueError
If the `freqs` parameter is not specified while the Network in :attr:`network` is `None`.
Also if `component` and/or `parameter` are not valid.
"""
components = ['db', 'mag', 'deg', 'deg_unwrap', 're', 'im']
if component.lower() in components:
if self.residues is None or self.poles is None:
raise RuntimeError('Poles and/or residues have not been fitted. Cannot plot the model response.')
n_ports = int(np.sqrt(np.shape(self.residues)[0]))
if i == -1:
list_i = range(n_ports)
elif isinstance(i, int):
list_i = [i]
else:
list_i = i
if j == -1:
list_j = range(n_ports)
elif isinstance(j, int):
list_j = [j]
else:
list_j = j
if self.network is not None:
# plot the original network response at each sample frequency (scatter plot)
if parameter.lower() == 's':
responses = self.network.s
elif parameter.lower() == 'z':
responses = self.network.z
elif parameter.lower() == 'y':
responses = self.network.y
else:
raise ValueError('The network parameter type is not valid, must be `s`, `z`, or `y`, '
f'got `{parameter}`.')
i_samples = 0
for i in list_i:
for j in list_j:
if i_samples == 0:
label = 'Samples'
else:
label = '_nolegend_'
i_samples += 1
y_vals = None
if component.lower() == 'db':
y_vals = 20 * np.log10(np.abs(responses[:, i, j]))
elif component.lower() == 'mag':
y_vals = np.abs(responses[:, i, j])
elif component.lower() == 'deg':
y_vals = np.rad2deg(np.angle(responses[:, i, j]))
elif component.lower() == 'deg_unwrap':
y_vals = np.rad2deg(np.unwrap(np.angle(responses[:, i, j])))
elif component.lower() == 're':
y_vals = np.real(responses[:, i, j])
elif component.lower() == 'im':
y_vals = np.imag(responses[:, i, j])
ax.scatter(self.network.f, y_vals, color='r', label=label)
if freqs is None:
# get frequency array from the network
freqs = self.network.f
if freqs is None:
raise ValueError(
'Neither `freqs` nor `self.network` is specified. Cannot plot model response without any '
'frequency information.')
# plot the fitted responses
y_label = ''
i_fit = 0
for i in list_i:
for j in list_j:
if i_fit == 0:
label = 'Fit'
else:
label = '_nolegend_'
i_fit += 1
y_model = self.get_model_response(i, j, freqs)
y_vals = None
if component.lower() == 'db':
y_vals = 20 * np.log10(np.abs(y_model))
y_label = 'Magnitude (dB)'
elif component.lower() == 'mag':
y_vals = np.abs(y_model)
y_label = 'Magnitude'
elif component.lower() == 'deg':
y_vals = np.rad2deg(np.angle(y_model))
y_label = 'Phase (Degrees)'
elif component.lower() == 'deg_unwrap':
y_vals = np.rad2deg(np.unwrap(np.angle(y_model)))
y_label = 'Phase (Degrees)'
elif component.lower() == 're':
y_vals = np.real(y_model)
y_label = 'Real Part'
elif component.lower() == 'im':
y_vals = np.imag(y_model)
y_label = 'Imaginary Part'
ax.plot(freqs, y_vals, color='k', label=label)
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel(y_label)
ax.legend(loc='best')
# only print title if a single response is shown
if i_fit == 1:
ax.set_title(f'Response i={i}, j={j}')
return ax
else:
raise ValueError(f'The specified component ("{component}") is not valid. Must be in {components}.')
def plot_s_db(self, *args, **kwargs) -> Axes:
"""
Plots the magnitude in dB of the scattering parameter response(s) in the fit.
Parameters
----------
*args : any, optional
Additonal arguments to be passed to :func:`plot`.
**kwargs : dict, optional
Additonal keyword arguments to be passed to :func:`plot`.
Returns
-------
:class:`matplotlib.Axes`
matplotlib axes used for drawing. Either the passed :attr:`ax` argument or the one fetch from the current
figure.
Notes
-----
This simply calls ``plot('db', *args, **kwargs)``.
"""
return self.plot('db', *args, **kwargs)
def plot_s_mag(self, *args, **kwargs) -> Axes:
"""
Plots the magnitude in linear scale of the scattering parameter response(s) in the fit.
Parameters
----------
*args : any, optional
Additonal arguments to be passed to :func:`plot`.
**kwargs : dict, optional
Additonal keyword arguments to be passed to :func:`plot`.
Returns
-------
:class:`matplotlib.Axes`
matplotlib axes used for drawing. Either the passed :attr:`ax` argument or the one fetch from the current
figure.
Notes
-----
This simply calls ``plot('mag', *args, **kwargs)``.
"""
return self.plot('mag', *args, **kwargs)
def plot_s_deg(self, *args, **kwargs) -> Axes:
"""
Plots the phase in degrees of the scattering parameter response(s) in the fit.
Parameters
----------
*args : any, optional
Additonal arguments to be passed to :func:`plot`.
**kwargs : dict, optional
Additonal keyword arguments to be passed to :func:`plot`.
Returns
-------
:class:`matplotlib.Axes`
matplotlib axes used for drawing. Either the passed :attr:`ax` argument or the one fetch from the current
figure.
Notes
-----
This simply calls ``plot('deg', *args, **kwargs)``.
"""
return self.plot('deg', *args, **kwargs)
def plot_s_deg_unwrap(self, *args, **kwargs) -> Axes:
"""
Plots the unwrapped phase in degrees of the scattering parameter response(s) in the fit.
Parameters
----------
*args : any, optional
Additonal arguments to be passed to :func:`plot`.
**kwargs : dict, optional
Additonal keyword arguments to be passed to :func:`plot`.
Returns
-------
:class:`matplotlib.Axes`
matplotlib axes used for drawing. Either the passed :attr:`ax` argument or the one fetch from the current
figure.
Notes
-----
This simply calls ``plot('deg_unwrap', *args, **kwargs)``.
"""
return self.plot('deg_unwrap', *args, **kwargs)
def plot_s_re(self, *args, **kwargs) -> Axes:
"""
Plots the real part of the scattering parameter response(s) in the fit.
Parameters
----------
*args : any, optional
Additonal arguments to be passed to :func:`plot`.
**kwargs : dict, optional
Additonal keyword arguments to be passed to :func:`plot`.
Returns
-------
:class:`matplotlib.Axes`
matplotlib axes used for drawing. Either the passed :attr:`ax` argument or the one fetch from the current
figure.
Notes
-----
This simply calls ``plot('re', *args, **kwargs)``.
"""
return self.plot('re', *args, **kwargs)
def plot_s_im(self, *args, **kwargs) -> Axes:
"""
Plots the imaginary part of the scattering parameter response(s) in the fit.
Parameters
----------
*args : any, optional
Additonal arguments to be passed to :func:`plot`.
**kwargs : dict, optional
Additonal keyword arguments to be passed to :func:`plot`.
Returns
-------
:class:`matplotlib.Axes`
matplotlib axes used for drawing. Either the passed :attr:`ax` argument or the one fetch from the current
figure.
Notes
-----
This simply calls ``plot('im', *args, **kwargs)``.
"""
return self.plot('im', *args, **kwargs)
@axes_kwarg
def plot_s_singular(self, freqs: Any = None, *, ax: Axes = None) -> Axes:
"""
Plots the singular values of the vector fitted S-matrix in linear scale.
Parameters
----------
freqs : list of float or ndarray or None, optional
List of frequencies for the response plot. If None, the sample frequencies of the fitted network in
:attr:`network` are used. This only works if :attr:`network` is not `None`.
ax : :class:`matplotlib.Axes` object or None
matplotlib axes to draw on. If None, the current axes is fetched with :func:`gca()`.
Returns
-------
:class:`matplotlib.Axes`
matplotlib axes used for drawing. Either the passed :attr:`ax` argument or the one fetch from the current
figure.
Raises
------
ValueError
If the `freqs` parameter is not specified while the Network in :attr:`network` is `None`.
"""
if freqs is None:
if self.network is None:
raise ValueError(
'Neither `freqs` nor `self.network` is specified. Cannot plot model response without any '
'frequency information.')
else:
freqs = self.network.f
# get system matrices of state-space representation
A, B, C, D, E = self._get_ABCDE()
n_ports = np.shape(D)[0]
# calculate and save singular values for each frequency
u, sigma, vh = np.linalg.svd(self._get_s_from_ABCDE(freqs, A, B, C, D, E))
# plot the frequency response of each singular value
for n in range(n_ports):
ax.plot(freqs, sigma[:, n], label=fr'$\sigma_{n + 1}$')
ax.set_xlabel('Frequency (Hz)')
ax.set_ylabel('Magnitude')
ax.legend(loc='best')
return ax
@axes_kwarg
def plot_convergence(self, ax: Axes = None) -> Axes:
"""
Plots the history of the model residue parameter **d_res** during the iterative pole relocation process of the
vector fitting, which should eventually converge to a fixed value. Additionally, the relative change of the
maximum singular value of the coefficient matrix **A** are plotted, which serve as a convergence indicator.
Parameters
----------
ax : :class:`matplotlib.Axes` object or None
matplotlib axes to draw on. If None, the current axes is fetched with :func:`gca()`.
Returns
-------
:class:`matplotlib.Axes`
matplotlib axes used for drawing. Either the passed :attr:`ax` argument or the one fetch from the current
figure.
"""
ax.semilogy(np.arange(len(self.delta_max_history)) + 1, self.delta_max_history, color='darkblue')
ax.set_xlabel('Iteration step')
ax.set_ylabel('Max. relative change', color='darkblue')
ax2 = ax.twinx()
ax2.plot(np.arange(len(self.d_res_history)) + 1, self.d_res_history, color='orangered')
ax2.set_ylabel('Residue', color='orangered')
return ax
@axes_kwarg
def plot_passivation(self, ax: Axes = None) -> Axes:
"""
Plots the history of the greatest singular value during the iterative passivity enforcement process, which
should eventually converge to a value slightly lower than 1.0 or stop after reaching the maximum number of
iterations specified in the class variable :attr:`max_iterations`.
Parameters
----------
ax : :class:`matplotlib.Axes` object or None
matplotlib axes to draw on. If None, the current axes is fetched with :func:`gca()`.
Returns
-------
:class:`matplotlib.Axes`
matplotlib axes used for drawing. Either the passed :attr:`ax` argument or the one fetch from the current
figure.
"""
ax.plot(np.arange(len(self.history_max_sigma)) + 1, self.history_max_sigma)
ax.set_xlabel('Iteration step')
ax.set_ylabel('Max. singular value')
return ax
def write_spice_subcircuit_s(self, file: str, fitted_model_name: str = "s_equivalent",
create_reference_pins: bool = False) -> None:
"""
Creates an equivalent N-port subcircuit based on its vector fitted scattering (S) parameter responses
in spice simulator netlist syntax (compatible with LTspice, ngspice, Xyce, ...). The circuit synthesis is based
on a direct implementation of the state-space representation of the vector fitted model [#vf-book]_.
Parameters
----------
file : str
Path and filename including file extension (usually .sp) for the subcircuit file.
fitted_model_name: str
Name of the resulting subcircuit, default "s_equivalent"
create_reference_pins: bool
If set to True, the synthesized subcircuit will have N pin-pairs:
p1 p1_ref p2 p2_ref ... pN pN_ref
If set to False, the synthesized subcircuit will have N pins
p1 p2 ... pN
In this case, the reference nodes will be internally connected
to the global ground net 0.
The default is False
Returns
-------
None
Examples
--------
Load and fit the `Network`, then export the equivalent subcircuit:
>>> nw_3port = skrf.Network('my3port.s3p')
>>> vf = skrf.VectorFitting(nw_3port)
>>> vf.auto_fit()
>>> vf.write_spice_subcircuit_s('/my3port_model.sp')
References
----------
.. [#vf-book] S. Grivet-Talocia and B. Gustavsen, "Passive Macromodeling", Wiley, 2016,
doi: https://doi.org/10.1002/9781119140931
"""
if np.any(self.proportional_coeff):
build_e = True
else:
build_e = False
with open(file, 'w') as f:
# write title line
f.write('* EQUIVALENT CIRCUIT FOR VECTOR FITTED S-MATRIX\n')
f.write('* Created using scikit-rf vectorFitting.py\n')
f.write('*\n')
# Create subcircuit pin string and reference nodes
if create_reference_pins:
str_input_nodes = " ".join(map(lambda x: f'p{x + 1} p{x + 1}_ref', range(self.network.nports)))
else:
str_input_nodes = " ".join(map(lambda x: f'p{x + 1}', range(self.network.nports)))
f.write(f'.SUBCKT {fitted_model_name} {str_input_nodes}\n')
for i in range(self.network.nports):
f.write('*\n')
f.write(f'* Port network for port {i + 1}\n')
if create_reference_pins:
node_ref_i = f'p{i + 1}_ref'
else:
node_ref_i = '0'
# reference impedance (real, i.e. resistance) of port i
z0_i = np.real(self.network.z0[0, i])
# transfer gain of the controlled current sources representing the incident power wave a_i at port i
#
# the gain values result from the definition of the incident power wave:
# a_i = 1 / 2 / sqrt(Z0_i) * (V_i + Z0_i * I_i) = 1 / 2 / sqrt(Z0_i) * V_i + sqrt(Z0_i) / 2 * I_i
gain_vccs_a_i = 1 / 2 / np.sqrt(z0_i)
gain_cccs_a_i = np.sqrt(z0_i) / 2
# transfer gain of the controlled current source representing the reflected power wave b_i at port i
#
# the gain values result from the definition of the reflected power wave:
# b_i = 1 / 2 / sqrt(Z0_i) * (V_i - Z0_i * I_i)
#
# depending on the circuit topology used for the equivalent port network, this can be implemented
# with either controlled current and/or controlled voltage sources. in case of the Norton current
# source used in this implementation, the reflected power wave relates to the source current as:
# b_i = sqrt(Z0_i) / 2 * I_b_i <==> I_b_i = 2 / sqrt(Z0_i) * b_i
gain_b_i = 2 / np.sqrt(z0_i)
# dummy voltage source (v = 0) for port current sensing (I_i)
f.write(f'V{i + 1} p{i + 1} s{i + 1} 0\n')
# adding port reference resistor Ri = Z0_i
f.write(f'R{i + 1} s{i + 1} {node_ref_i} {z0_i}\n')
# transfer of states and inputs from port j to input/output network of port i
for j in range(self.network.nports):
if create_reference_pins:
node_ref_j = f'p{j + 1}_ref'
else:
node_ref_j = '0'
# reference impedance (real, i.e. resistance) of port i
z0_j = np.real(self.network.z0[0, j])
# Stacking order in VectorFitting class variables:
# s11, s12, s13, ..., s21, s22, s23, ...
idx_S_i_j = i * self.network.nports + j
# VCCS and CCCS adding their currents to represent the incident wave a_j
gain_vccs_a_j = 1 / 2 / np.sqrt(z0_j)
gain_cccs_a_j = np.sqrt(z0_j) / 2
d = self.constant_coeff[idx_S_i_j]
e = self.proportional_coeff[idx_S_i_j]
if d != 0.0:
# avoid zero-valued coefficients (in case of fit_constant=False)
# input a_j is scaled by constant term d_i_j and by current gain for b_i
g_ij = gain_b_i * d * gain_vccs_a_j
f_ij = gain_b_i * d * gain_cccs_a_j
f.write(f'Gd{i + 1}_{j + 1} {node_ref_i} s{i + 1} p{j + 1} {node_ref_j} {g_ij}\n')
f.write(f'Fd{i + 1}_{j + 1} {node_ref_i} s{i + 1} V{j + 1} {f_ij}\n')
if build_e and e != 0.0:
# avoid zero-valued coefficients (in case of fit_proportional=False)
# proportional coefficients require an extra node for the differentiation using an inductor
# [Y(s) ~ s * E * U(s)]
# differentiated input a_j is scaled by proportional term e_i_j and by current gain for b_i
g_ij = gain_b_i * e
f.write(f'Ge{i + 1}_{j + 1} {node_ref_i} s{i + 1} e{j + 1} 0 {g_ij}\n')
# each residue rk_i_j at port i is multiplied by its respective state signal xk_j
for k in range(len(self.poles)):
pole = self.poles[k]
residue = self.residues[idx_S_i_j, k]
g_re = gain_b_i * np.real(residue)
g_im = gain_b_i * np.imag(residue)
if np.imag(pole) == 0.0:
# Real pole/residue pair; represented by one state
xkj = f'x{k + 1}_a{j + 1}'
f.write(f'Gr{k + 1}_{i + 1}_{j + 1} {node_ref_i} s{i + 1} {xkj} 0 {g_re}\n')
else:
# Complex-conjugate pole/residue pair; represented by two states
# real part at x_{k + 1}_re_{j + 1}
# imaginary part at x_{k + 1}_im_{j + 1}
xk_re_j = f'x{k + 1}_re_a{j + 1}'
xk_im_j = f'x{k + 1}_im_a{j + 1}'
f.write(f'Gr{k + 1}_re_{i + 1}_{j + 1} {node_ref_i} s{i + 1} {xk_re_j} 0 {g_re}\n')
f.write(f'Gr{k + 1}_im_{i + 1}_{j + 1} {node_ref_i} s{i + 1} {xk_im_j} 0 {g_im}\n')
# create state networks driven by this port i (input variable u = a_i)
f.write('*\n')
f.write(f'* State networks driven by port {i + 1}\n')
for k in range(len(self.poles)):
pole = self.poles[k]
pole_re = np.real(pole)
pole_im = np.imag(pole)
# Transfer of input (a_i) to state networks (node xk_i) using VCCS and CCCS
if pole_im == 0.0:
# Real pole; represented by one state, input a_i is scaled by b = 1
xki = f'x{k + 1}_a{i + 1}'
f.write(f'Cx{k + 1}_a{i + 1} {xki} 0 1.0\n') # 1F capacitor makes math easy
f.write(f'Gx{k + 1}_a{i + 1} 0 {xki} p{i + 1} {node_ref_i} {1 * gain_vccs_a_i}\n')
f.write(f'Fx{k + 1}_a{i + 1} 0 {xki} V{i + 1} {1 * gain_cccs_a_i}\n')
f.write(f'Rp{k + 1}_a{i + 1} 0 {xki} {-1 / pole_re}\n')
else:
# Complex pole of a conjugate pair; represented by two states
# real part at x_{k + 1}_re_{i + 1}, input a_i is scaled by b = 2
xk_re_i = f'x{k + 1}_re_a{i + 1}'
xk_im_i = f'x{k + 1}_im_a{i + 1}'
f.write(f'Cx{k + 1}_re_a{i + 1} {xk_re_i} 0 1.0\n') # 1F capacitor makes math easy
f.write(
f'Gx{k + 1}_re_a{i + 1} 0 {xk_re_i} p{i + 1} {node_ref_i} {2 * gain_vccs_a_i}\n')
f.write(f'Fx{k + 1}_re_a{i + 1} 0 {xk_re_i} V{i + 1} {2 * gain_cccs_a_i}\n')
f.write(f'Rp{k + 1}_re_re_a{i + 1} 0 {xk_re_i} {-1 / pole_re}\n')
f.write(f'Gp{k + 1}_re_im_a{i + 1} 0 {xk_re_i} {xk_im_i} 0 {pole_im}\n')
# imaginary part at x_{k + 1}_im_{i + 1}, input a_i is inactive (b = 0)
f.write(f'Cx{k + 1}_im_a{i + 1} {xk_im_i} 0 1.0\n') # 1F capacitor makes math easy
f.write(f'Gp{k + 1}_im_re_a{i + 1} 0 {xk_im_i} {xk_re_i} 0 {-1 * pole_im}\n')
f.write(f'Rp{k + 1}_im_im_a{i + 1} 0 {xk_im_i} {-1 / pole_re}\n')
# create differentiation network for this port i (input variable u = a_i)
if build_e:
f.write('*\n')
f.write(f'* Network with derivative of input a_{i + 1} for proportional term\n')
# voltage on node 'e{i + 1}' to gnd (0) represents time-derivative of input a_i for terms e_j_i
f.write(f'Le{i + 1} e{i + 1} 0 1.0\n') # 1H inductor makes math easy
f.write(f'Ge{i + 1} 0 e{i + 1} p{i + 1} {node_ref_i} {gain_vccs_a_i}\n')
f.write(f'Fe{i + 1} 0 e{i + 1} V{i + 1} {gain_cccs_a_i}\n')
f.write(f'.ENDS {fitted_model_name}\n')
Reference in New Issue View Git Blame Copy Permalink