# -*- coding: iso-8859-1 -*-
# pz.py
# Pole-Zero evaluation methods
# Copyright 2014 Giuseppe Venturini
# This file is part of the ahkab simulator.
#
# Ahkab is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, version 2 of the License.
#
# Ahkab is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License v2
# along with ahkab. If not, see <http://www.gnu.org/licenses/>.
"""
This module offers the functions needed to perform a numeric
pole-zero extraction.
Currently, this module implements the MD algorithm, more may be
added in the future.
A description of the algorithm is found in the following references:
Haley, S.B., "The generalized eigenproblem: pole-zero computation,"
*Proceedings of the IEEE*, vol.76, no.2, pp.103,120, Feb 1988
and:
Raghuram, R.; Divekar, D.; Wang, P., "Implementation of pole-zero
analysis in SPICE based on the MD method," Circuits and Systems, 1991.,
*Proceedings of the 34th Midwest Symposium on*, pp.380,
383 vol.1, 14-17 May 1991
Frequency sweeping -- or shifting -- is performed with a random frequency
kick, currently, hoping not to kick so hard that we end up on the negative side.
A bisection method would be better and hopefully will be implemented soon.
Overview
~~~~~~~~
Two main methods are available in this module:
* :func:`calculate_singularities`, which computes both zeros and poles,
* :func:`calculate_poles`, which only computes the poles.
Currently this module uses dense matrices.
Reference
~~~~~~~~~
"""
from __future__ import (unicode_literals, absolute_import,
division, print_function)
import copy
import numpy as np
from . import circuit
from . import dc_analysis
from . import components
from . import transient
from . import plotting
from . import printing
from . import py3compat
from . import options
from . import results
specs = {'pz': {'tokens': ({
'label': 'output',
'pos': 0,
'type': str,
'needed': False,
'dest': 'output_port',
'default': None
},
{
'label': 'input',
'pos': 1,
'type': str,
'needed': False,
'dest': 'input_source',
'default': None
},
{
'label': 'shift',
'pos': 2,
'type': float,
'needed': False,
'dest': 'shift',
'default': 0.0
})
}
}
def _enlarge_matrix(M):
if M is None:
return np.zeros((1, 1))
else:
return np.vstack((np.hstack((M, np.zeros((M.shape[0], 1)))),
np.zeros((1, M.shape[1] + 1))))
[docs]def calculate_poles(mc, MNA=None, x0=None, outfile=None, verbose=0):
"""Calculate the circuit poles.
**Parameters:**
mc : circuit instance
The circuit to be analyzed.
MNA : ndarray, optional
The Modified Nodal Analysis matrix, if available.
In case the circuit is non-linear, MNA should include the contributes
of the non-linear elements (ie the Jacobian :math:`J`).
x0 : ndarray or op_solution, optional
The linearization point. Only needed for non-linear circuits.
outfile : str or None, optional
The data filename.
verbose : int, optional
Verbosity level, from 0 (silent, default) to 6 (debug).
**Returns:**
pz_sol : pz_solution instance
The PZ solution, with no zeros.
"""
return calculate_singularities(mc, input_source=None, output_port=None,
MNA=None, shift=0)[0]
[docs]def calculate_singularities(mc, input_source=None, output_port=None, MNA=None,
x0=None, shift=0, outfile=None, verbose=0):
"""Calculate poles and zeros.
By default, only poles are calculated, as they need no information
other than the circuit description.
To activate zeros calculation, it is necessary:
* to specify an input source (``input_source``),
* to specify an output port (``output_port``).
**Parameters:**
mc : circuit instance
The circuit to be analyzed.
input_source : string or element, optional
If zeros are to be calculated, set this to the input surce.
output_port : external node (ref. to gnd) or tuple of external nodes, opt
If zeros are to be calculated, set this to the output nodes.
MNA : ndarray, optional
The Modified Nodal Analysis matrix, if available.
In case the circuit is non-linear, MNA should include the contributes
of the non-linear elements (ie the Jacobian :math:`J`).
x0 : ndarray or op_solution, optional
The linearization point. Only needed for non-linear circuits.
shift : float, optional
Shift frequency at which the algorithm should be run.
outfile : str or None, optional
The data filename.
verbose : int, optional
Verbosity level, from 0 (silent, default) to 6 (debug).
**Returns:**
pz_sol : pz_solution instance
The PZ solution
"""
calc_zeros = (input_source is not None) and (output_port is not None)
if calc_zeros:
if type(input_source) not in py3compat.string_types:
input_source = input_source.part_id
if type(output_port) in py3compat.string_types:
output_port = plotting._split_netlist_label(output_port)[0]
output_port = [o[1:] for o in output_port]
output_port = [o.lower() for o in output_port]
if len(output_port) == 1:
# we refer to the ground implicitely
output_port += ['0']
if np.isscalar(output_port):
output_port = (output_port, mc.gnd)
if (type(output_port) == tuple or type(output_port) == list) \
and type(output_port[0]) in py3compat.string_types:
output_port = [mc.ext_node_to_int(o) for o in output_port]
we_got_source = False
for e in mc:
if e.part_id == input_source:
we_got_source = True
break
if not we_got_source:
raise ValueError('Source %s not found in circuit.' % input_source)
RIIN = []
ROUT = []
if MNA is None:
MNA, N = dc_analysis.generate_mna_and_N(mc)
if mc.is_nonlinear():
# setup x0
if x0 is None:
printing.print_warning("PZ: No linearization point provided. Using x0 = 0.")
x0 = np.zeros((MNA.shape[0] - 1, 1))
else:
if isinstance(x0, results.op_solution):
x0 = x0.asarray()
# else
# hopefully x0 is an ndarray!
printing.print_info_line(("Using the supplied op as " +
"linearization point.", 5), verbose)
J, _ = dc_analysis.build_J_and_Tx(x0, MNA.shape[0]-1, mc, time=0., sparse=False)
MNA[1:, 1:] += J
D = transient.generate_D(mc, MNA[1:, 1:].shape)
MNAinv = np.linalg.inv(MNA[1:, 1:] + shift*D[1:, 1:])
nodes_m1 = mc.get_nodes_number() - 1
vde1 = -1
MC = np.zeros((MNA.shape[0] - 1, 1))
TCM = None
dei_source = 0
for e1 in mc:
if circuit.is_elem_voltage_defined(e1):
vde1 += 1
if isinstance(e1, components.Capacitor):
MC[e1.n1 - 1, 0] += 1. if e1.n1 > 0 else 0.
MC[e1.n2 - 1, 0] -= 1. if e1.n2 > 0 else 0.
elif isinstance(e1, components.Inductor):
MC[nodes_m1 + vde1] += -1.
elif calc_zeros and e1.part_id == input_source:
if isinstance(e1, components.sources.VSource):
MC[nodes_m1 + vde1] += -1.
elif isinstance(e1, components.sources.ISource):
MC[e1.n1 - 1, 0] += 1. if e1.n1 > 0 else 0.
MC[e1.n2 - 1, 0] -= 1. if e1.n2 > 0 else 0.
else:
raise Exception("Unknown input source type %s" % input_source)
else:
continue
TV = -1. * np.dot(MNAinv, MC)
dei_victim = 0
vde2 = -1
for e2 in mc:
if circuit.is_elem_voltage_defined(e2):
vde2 += 1
if isinstance(e2, components.Capacitor):
v = 0
if e2.n1:
v += TV[e2.n1 - 1, 0]
if e2.n2:
v -= TV[e2.n2 - 1, 0]
elif isinstance(e2, components.Inductor):
v = TV[nodes_m1 + vde2, 0]
else:
continue
if calc_zeros and e1.part_id == input_source:
RIIN += [v]
else:
if not dei_source:
TCM = _enlarge_matrix(TCM)
TCM[dei_victim, dei_source] = v*e1.value
dei_victim += 1
if calc_zeros and e1.part_id == input_source:
ROUTIN = 0
o1, o2 = output_port
if o1:
ROUTIN += TV[o1 - 1, 0]
if o2:
ROUTIN -= TV[o2 - 1, 0]
else:
dei_source += 1
# reset, get ready to restart
MC[:, :] = 0.
if TCM is not None:
if np.linalg.det(TCM):
poles = 1./(2.*np.pi)*(1./np.linalg.eigvals(TCM) + shift)
else:
return calculate_singularities(mc, input_source, output_port, MNA=MNA,
outfile=outfile,
shift=shift+np.abs(1+np.random.uniform())*1e3)
else:
poles = []
if calc_zeros and TCM is not None:
# re-loop, get the ROUT elements
vde1 = -1
MC = np.zeros((MNA.shape[0] - 1, 1))
ROUT = []
for e1 in mc:
if circuit.is_elem_voltage_defined(e1):
vde1 += 1
if isinstance(e1, components.Capacitor):
MC[e1.n1 - 1, 0] += 1. if e1.n1 > 0 else 0.
MC[e1.n2 - 1, 0] -= 1. if e1.n2 > 0 else 0.
elif isinstance(e1, components.Inductor):
MC[nodes_m1 + vde1, 0] += -1.
else:
continue
TV = -1.*np.dot(MNAinv, MC)
v = 0
o1, o2 = output_port
if o1:
v += TV[o1 - 1, 0]
if o2:
v -= TV[o2 - 1, 0]
ROUT += [v*e1.value]
# reset, get ready to restart
MC[:, :] = 0.
# Reshape the matrices and evaluate the zero correction.
RIIN = np.array(RIIN).reshape((-1, 1))
RIIN = np.tile(RIIN, (1, RIIN.shape[0]))
ROUT = np.diag(np.atleast_1d(np.array(ROUT)))
if ROUT.any():
try:
if not np.all(ROUTIN) or not np.isfinite(ROUTIN):
# immediate catch
raise ValueError("ROUT-IN is either Inf, NaN or null")
ZTCM = TCM - np.dot(RIIN, ROUT)/ROUTIN
if not (np.isfinite(ZTCM).all()):
raise ValueError("Array contains infs, NaNs or both")
# immediate catch
eigvals = np.linalg.eigvals(ZTCM)
if not np.all(eigvals) or not np.isfinite(eigvals).all():
# immediate catch
raise ValueError("ZTCM eigenvalues contain either Inf, NaN or null values")
zeros = 1./(2.*np.pi)*(1./np.linalg.eigvals(ZTCM) + shift)
except ValueError:
return calculate_singularities(mc, input_source, output_port,
MNA=MNA, outfile=outfile,
shift=shift+np.abs(np.random.uniform()+1)*1e3)
elif shift < 1e12:
return calculate_singularities(mc, input_source, output_port,
MNA=MNA, outfile=outfile,
shift=shift*np.abs(np.random.uniform()+1)*10)
else:
zeros = []
else:
zeros = []
poles = np.array([a for a in poles if np.abs(a) < options.pz_max], dtype=np.complex_)
zeros = np.array([a for a in zeros if np.abs(a) < options.pz_max], dtype=np.complex_)
poles = np.sort_complex(poles)
zeros = np.sort_complex(zeros)
res = results.pz_solution(mc, poles, zeros, outfile)
return res
def _check_singularities(res, ref, atol=1e-4, rtol=1e-3):
ref = copy.copy(ref)
for i in res:
success = False
for si, s in enumerate(ref):
if np.allclose(i, s, atol=atol, rtol=rtol):
success = True
break
assert si != len(ref) - 1
ref.pop(si)
assert not len(ref)