Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

TransferFunction _common_den rewrite - issue #194 #206

Merged
merged 13 commits into from
Jul 2, 2018
Merged

TransferFunction _common_den rewrite - issue #194 #206

Changes from 5 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
237597e
partial changes, not complete, to num/den combination
repagh May 29, 2018
5acb42c
have replaced the _common_den function internals. Passes tests
repagh May 29, 2018
05a6934
Merge branch 'master' into tf-combination
repagh May 29, 2018
a0810db
make xferfcn.py python2 compatible again
repagh May 29, 2018
f4c230c
Merge branch 'tf-combination' of https://github.com/repagh/python-con…
repagh May 29, 2018
492635c
somewhere halfway
repagh Jun 4, 2018
324c573
still not figured out
repagh Jun 6, 2018
0722cf3
work in progress
repagh Jun 14, 2018
dfb1507
more work in progress
repagh Jun 15, 2018
50c9bf5
working tf -> ss transformation now, also solved #111
repagh Jun 20, 2018
92a3354
disabled MIMO test when no slycot
repagh Jun 20, 2018
7e594c4
- do not cancel pole/zero pairs before calculating pole() in xferfcn.py
repagh Jun 26, 2018
732c924
- change testModred; that one did state removal on a system
repagh Jul 1, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
198 changes: 70 additions & 128 deletions control/xferfcn.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,12 +55,14 @@
import numpy as np
from numpy import angle, any, array, empty, finfo, insert, ndarray, ones, \
polyadd, polymul, polyval, roots, sort, sqrt, zeros, squeeze, exp, pi, \
where, delete, real, poly, poly1d
where, delete, real, poly, poly1d, nonzero
import scipy as sp
from numpy.polynomial.polynomial import polyfromroots
from scipy.signal import lti, tf2zpk, zpk2tf, cont2discrete
from copy import deepcopy
import warnings
from warnings import warn
from itertools import chain
from .lti import LTI, timebaseEqual, timebase, isdtime

__all__ = ['TransferFunction', 'tf', 'ss2tf', 'tfdata']
Expand Down Expand Up @@ -727,144 +729,84 @@ def _common_den(self, imag_tol=None):
if (imag_tol is None):
imag_tol = 1e-8 # TODO: figure out the right number to use

# A sorted list to keep track of cumulative poles found as we scan
# self.den.
poles = []
# A list to keep track of cumulative poles found as we scan
# self.den[..][..]
poles = [ ]

# A 3-D list to keep track of common denominator poles not present in
# the self.den[i][j].
missingpoles = [[[] for j in range(self.inputs)]
for i in range(self.outputs)]
# RvP, new implementation 180526, issue #194

# pre-calculate the poles for all num, den
# has zeros, poles, gain, list for pole indices not in den,
# number of poles known at the time analyzed
self2 = self.minreal()
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure about this; I think tf([1,1],[1,2,1]).pole() should return [-1, -1].

poleset = []
for i in range(self.outputs):
poleset.append([])
for j in range(self.inputs):
# A sorted array of the poles of this SISO denominator.
currentpoles = sort(roots(self.den[i][j]))

cp_ind = 0 # Index in currentpoles.
p_ind = 0 # Index in poles.

# Crawl along the list of current poles and the list of
# cumulative poles, until one of them reaches the end. Keep in
# mind that both lists are always sorted.
while cp_ind < len(currentpoles) and p_ind < len(poles):
if abs(currentpoles[cp_ind] - poles[p_ind]) < (10 * eps):
# If the current element of both
# lists match, then we're
# good. Move to the next pair of elements.
cp_ind += 1
elif currentpoles[cp_ind] < poles[p_ind]:
# We found a pole in this transfer function that's not
# in the list of cumulative poles. Add it to the list.
poles.insert(p_ind, currentpoles[cp_ind])
# Now mark this pole as "missing" in all previous
# denominators.
for k in range(i):
for m in range(self.inputs):
# All previous rows.
missingpoles[k][m].append(currentpoles[cp_ind])
for m in range(j):
# This row only.
missingpoles[i][m].append(currentpoles[cp_ind])
cp_ind += 1
else:
# There is a pole in the cumulative list of poles that
# is not in our transfer function denominator. Mark
# this pole as "missing", and do not increment cp_ind.
missingpoles[i][j].append(poles[p_ind])
p_ind += 1

if cp_ind == len(currentpoles) and p_ind < len(poles):
# If we finished scanning currentpoles first, then all the
# remaining cumulative poles are missing poles.
missingpoles[i][j].extend(poles[p_ind:])
elif cp_ind < len(currentpoles) and p_ind == len(poles):
# If we finished scanning the cumulative poles first, then
# all the reamining currentpoles need to be added to poles.
poles.extend(currentpoles[cp_ind:])
# Now mark these poles as "missing" in previous
# denominators.
for k in range(i):
for m in range(self.inputs):
# All previous rows.
missingpoles[k][m].extend(currentpoles[cp_ind:])
for m in range(j):
# This row only.
missingpoles[i][m].extend(currentpoles[cp_ind:])

# Construct the common denominator.
den = 1.
n = 0
while n < len(poles):
if abs(poles[n].imag) > 10 * eps:
# To prevent buildup of imaginary part error, handle complex
# pole pairs together.
#
# Because we might have repeated real parts of poles
# and the fact that we are using lexigraphical
# ordering, we can't just combine adjacent poles.
# Instead, we have to figure out the multiplicity
# first, then multiple the pairs from the outside in.

# Figure out the multiplicity
m = 1 # multiplicity count
while (n+m < len(poles) and
poles[n].real == poles[n+m].real and
poles[n].imag * poles[n+m].imag > 0):
m += 1

# Multiple pairs from the outside in
for i in range(m):
quad = polymul([1., -poles[n]], [1., -poles[n+2*(m-i)-1]])
assert all(quad.imag < 10 * eps), \
"Quadratic has a nontrivial imaginary part: %g" \
% quad.imag.max()

den = polymul(den, quad.real)
n += 1 # move to next pair
n += m # skip past conjugate pairs
else:
den = polymul(den, [1., -poles[n].real])
n += 1

# Modify the numerators so that they each take the common denominator.
num = deepcopy(self.num)
if isinstance(den, float):
den = array([den])

if abs(self2.num[i][j]).max() <= eps:
poleset[-1].append( [array([], dtype=float),
roots(self2.den[i][j]), 0.0, [], 0 ])
else:
z, p, k = tf2zpk(self2.num[i][j], self2.den[i][j])
poleset[-1].append([ z, p, k, [], 0])

# collect all individual poles
epsnm = eps * self.inputs * self.outputs
for i in range(self.outputs):
for j in range(self.inputs):
# The common denominator has leading coefficient 1. Scale out
# the existing denominator's leading coefficient.
assert self.den[i][j][0], "The i = %i, j = %i denominator has \
a zero leading coefficient." % (i, j)
num[i][j] = num[i][j] / self.den[i][j][0]

# Multiply in the missing poles.
for p in missingpoles[i][j]:
num[i][j] = polymul(num[i][j], [1., -p])

# Pad all numerator polynomials with zeros so that the numerator arrays
# are the same size as the denominator.
currentpoles = poleset[i][j][1]
nothave = ones(currentpoles.shape, dtype=bool)
for ip, p in enumerate(poles):
Copy link
Contributor

@roryyorke roryyorke Jun 2, 2018 •
edited

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is this supposed to be enumerate(poleset)? I can't see any modification to poles after the initial assignment to [].

edit: ah, never mind, I see the poles.append() a bit later.

idx, = nonzero(
(abs(currentpoles - p) < epsnm) * nothave)
if len(idx):
nothave[idx[0]] = False
else:
# remember id of pole not in tf
poleset[i][j][3].append(ip)
for h, c in zip(nothave, currentpoles):
if h:
poles.append(c)
# remember how many poles now known
poleset[i][j][4] = len(poles)

# for only gain systems
if len(poles) == 0:
den = ones((1,), dtype=float)
num = zeros((self.outputs, self.inputs, 1), dtype=float)
for i in range(self.outputs):
for j in range(self.inputs):
num[i,j,0] = poleset[i][j][2]
return num, den

# recreate the denominator
den = polyfromroots(poles)[::-1]
if (abs(den.imag) > epsnm).any():
print("Warning: The denominator has a nontrivial imaginary part: %f"
% abs(den.imag).max())
den = den.real
np = len(poles)

# now supplement numerators with all new poles
num = zeros((self.outputs, self.inputs, len(poles)+1), dtype=float)
for i in range(self.outputs):
for j in range(self.inputs):
pad = len(den) - len(num[i][j])
if (pad > 0):
num[i][j] = insert(
num[i][j], zeros(pad, dtype=int),
zeros(pad))

# Finally, convert the numerator to a 3-D array.
num = array(num)
# Remove trivial imaginary parts.
# Check for nontrivial imaginary parts.
if any(abs(num.imag) > sqrt(eps)):
print ("Warning: The numerator has a nontrivial imaginary part: %g"
% abs(num.imag).max())
num = num.real
# collect as set of zeros
nwzeros = list(poleset[i][j][0])
# add all poles not found in this denominator, and the
# ones later added from other denominators
for ip in chain(poleset[i][j][3],
range(poleset[i][j][4],np)):
nwzeros.append(poles[ip])
m = len(nwzeros) + 1
num[i,j,-m:] = polyfromroots(nwzeros).real[::-1]

# determine tf gain correction
num[i,j] *= poleset[i][j][2]

return num, den


def sample(self, Ts, method='zoh', alpha=None):
"""Convert a continuous-time system to discrete time

Expand Down