Skip to content

Commit cb25633

Browse files
murrayrmmurrayrm
authored
Merge pull request #345 from bnavigator/fix-armfail
change root precision tolerance and imaginary detection in xferfcn._common_den()
2 parents cee8bf4 + 712f78a commit cb25633

File tree

2 files changed

+89
-41
lines changed

2 files changed

+89
-41
lines changed

‎control/tests/xferfcn_test.py‎

Lines changed: 53 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -496,8 +496,57 @@ def test_freqresp_mimo(self):
496496
np.testing.assert_array_equal(omega, true_omega)
497497

498498
# Tests for TransferFunction.pole and TransferFunction.zero.
499-
500-
@unittest.skipIf(not slycot_check(), "slycot not installed")
499+
500+
def test_common_den(self):
501+
""" Test the helper function to compute common denomitators."""
502+
503+
# _common_den() computes the common denominator per input/column.
504+
# The testing columns are:
505+
# 0: no common poles
506+
# 1: regular common poles
507+
# 2: poles with multiplicity,
508+
# 3: complex poles
509+
# 4: complex poles below threshold
510+
511+
eps = np.finfo(float).eps
512+
tol_imag = np.sqrt(eps*5*2*2)*0.9
513+
514+
numin = [[[1.], [1.], [1.], [1.], [1.]],
515+
[[1.], [1.], [1.], [1.], [1.]]]
516+
denin = [[[1., 3., 2.], # 0: poles: [-1, -2]
517+
[1., 6., 11., 6.], # 1: poles: [-1, -2, -3]
518+
[1., 6., 11., 6.], # 2: poles: [-1, -2, -3]
519+
[1., 6., 11., 6.], # 3: poles: [-1, -2, -3]
520+
[1., 6., 11., 6.]], # 4: poles: [-1, -2, -3],
521+
[[1., 12., 47., 60.], # 0: poles: [-3, -4, -5]
522+
[1., 9., 26., 24.], # 1: poles: [-2, -3, -4]
523+
[1., 7., 16., 12.], # 2: poles: [-2, -2, -3]
524+
[1., 7., 17., 15.], # 3: poles: [-2+1J, -2-1J, -3],
525+
np.poly([-2 + tol_imag * 1J, -2 - tol_imag * 1J, -3])]]
526+
numref = np.array([
527+
[[0., 0., 1., 12., 47., 60.],
528+
[0., 0., 0., 1., 4., 0.],
529+
[0., 0., 0., 1., 2., 0.],
530+
[0., 0., 0., 1., 4., 5.],
531+
[0., 0., 0., 1., 2., 0.]],
532+
[[0., 0., 0., 1., 3., 2.],
533+
[0., 0., 0., 1., 1., 0.],
534+
[0., 0., 0., 1., 1., 0.],
535+
[0., 0., 0., 1., 3., 2.],
536+
[0., 0., 0., 1., 1., 0.]]])
537+
denref = np.array(
538+
[[1., 15., 85., 225., 274., 120.],
539+
[1., 10., 35., 50., 24., 0.],
540+
[1., 8., 23., 28., 12., 0.],
541+
[1., 10., 40., 80., 79., 30.],
542+
[1., 8., 23., 28., 12., 0.]])
543+
sys = TransferFunction(numin, denin)
544+
num, den, denorder = sys._common_den()
545+
np.testing.assert_array_almost_equal(num[:2, :, :], numref)
546+
np.testing.assert_array_almost_equal(num[2:, :, :],
547+
np.zeros((3, 5, 6)))
548+
np.testing.assert_array_almost_equal(den, denref)
549+
501550
def test_pole_mimo(self):
502551
"""Test for correct MIMO poles."""
503552

@@ -508,13 +557,12 @@ def test_pole_mimo(self):
508557

509558
np.testing.assert_array_almost_equal(p, [-2., -2., -7., -3., -2.])
510559

511-
@unittest.skipIf(not slycot_check(), "slycot not installed")
512560
def test_double_cancelling_poles_siso(self):
513-
561+
514562
H = TransferFunction([1, 1], [1, 2, 1])
515563
p = H.pole()
516564
np.testing.assert_array_almost_equal(p, [-1, -1])
517-
565+
518566
# Tests for TransferFunction.feedback
519567
def test_feedback_siso(self):
520568
"""Test for correct SISO transfer function feedback."""

‎control/xferfcn.py‎

Lines changed: 36 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -57,7 +57,6 @@
5757
polyadd, polymul, polyval, roots, sqrt, zeros, squeeze, exp, pi, \
5858
where, delete, real, poly, nonzero
5959
import scipy as sp
60-
from numpy.polynomial.polynomial import polyfromroots
6160
from scipy.signal import lti, tf2zpk, zpk2tf, cont2discrete
6261
from copy import deepcopy
6362
from warnings import warn
@@ -803,8 +802,6 @@ def _common_den(self, imag_tol=None):
803802
output numerator array num is modified to use the common
804803
denominator for this input/column; the coefficient arrays are also
805804
padded with zeros to be the same size for all num/den.
806-
num is an sys.outputs by sys.inputs
807-
by len(d) array.
808805
809806
Parameters
810807
----------
@@ -815,17 +812,20 @@ def _common_den(self, imag_tol=None):
815812
Returns
816813
-------
817814
num: array
818-
Multi-dimensional array of numerator coefficients. num[i][j]
819-
gives the numerator coefficient array for the ith input and jth
820-
output, also prepared for use in td04ad; matches the denorder
821-
order; highest coefficient starts on the left.
815+
n by n by kd where n = max(sys.outputs,sys.inputs)
816+
kd = max(denorder)+1
817+
Multi-dimensional array of numerator coefficients. num[i,j]
818+
gives the numerator coefficient array for the ith output and jth
819+
input; padded for use in td04ad ('C' option); matches the
820+
denorder order; highest coefficient starts on the left.
822821
823822
den: array
823+
sys.inputs by kd
824824
Multi-dimensional array of coefficients for common denominator
825825
polynomial, one row per input. The array is prepared for use in
826826
slycot td04ad, the first element is the highest-order polynomial
827-
coefficiend of s, matching the order in denorder, if denorder <
828-
number of columns in den, the den is padded with zeros
827+
coefficient of s, matching the order in denorder. If denorder <
828+
number of columns in den, the den is padded with zeros.
829829
830830
denorder: array of int, orders of den, one per input
831831
@@ -839,16 +839,18 @@ def _common_den(self, imag_tol=None):
839839

840840
# Machine precision for floats.
841841
eps = finfo(float).eps
842+
real_tol = sqrt(eps * self.inputs * self.outputs)
842843

843-
# Decide on the tolerance to use in deciding of a pole is complex
844+
# The tolerance to use in deciding if a pole is complex
844845
if (imag_tol is None):
845-
imag_tol = 1e-8 # TODO: figure out the right number to use
846+
imag_tol = 2 * real_tol
846847

847848
# A list to keep track of cumulative poles found as we scan
848849
# self.den[..][..]
849850
poles = [[] for j in range(self.inputs)]
850851

851852
# RvP, new implementation 180526, issue #194
853+
# BG, modification, issue #343, PR #354
852854

853855
# pre-calculate the poles for all num, den
854856
# has zeros, poles, gain, list for pole indices not in den,
@@ -867,30 +869,36 @@ def _common_den(self, imag_tol=None):
867869
poleset[-1].append([z, p, k, [], 0])
868870

869871
# collect all individual poles
870-
epsnm = eps * self.inputs * self.outputs
871872
for j in range(self.inputs):
872873
for i in range(self.outputs):
873874
currentpoles = poleset[i][j][1]
874875
nothave = ones(currentpoles.shape, dtype=bool)
875876
for ip, p in enumerate(poles[j]):
876-
idx, = nonzero(
877-
(abs(currentpoles - p) < epsnm) * nothave)
878-
if len(idx):
879-
nothave[idx[0]] = False
877+
collect = (np.isclose(currentpoles.real, p.real,
878+
atol=real_tol) &
879+
np.isclose(currentpoles.imag, p.imag,
880+
atol=imag_tol) &
881+
nothave)
882+
if np.any(collect):
883+
# mark first found pole as already collected
884+
nothave[nonzero(collect)[0][0]] = False
880885
else:
881886
# remember id of pole not in tf
882887
poleset[i][j][3].append(ip)
883888
for h, c in zip(nothave, currentpoles):
884889
if h:
890+
if abs(c.imag) < imag_tol:
891+
c = c.real
885892
poles[j].append(c)
886893
# remember how many poles now known
887894
poleset[i][j][4] = len(poles[j])
888895

889896
# figure out maximum number of poles, for sizing the den
890-
npmax = max([len(p) for p in poles])
891-
den = zeros((self.inputs, npmax + 1), dtype=float)
897+
maxindex = max([len(p) for p in poles])
898+
den = zeros((self.inputs, maxindex + 1), dtype=float)
892899
num = zeros((max(1, self.outputs, self.inputs),
893-
max(1, self.outputs, self.inputs), npmax + 1),
900+
max(1, self.outputs, self.inputs),
901+
maxindex + 1),
894902
dtype=float)
895903
denorder = zeros((self.inputs,), dtype=int)
896904

@@ -901,12 +909,11 @@ def _common_den(self, imag_tol=None):
901909
for i in range(self.outputs):
902910
num[i, j, 0] = poleset[i][j][2]
903911
else:
904-
# create the denominator matching this input polyfromroots
905-
# gives coeffs in opposite order from what we use coefficients
906-
# should be padded on right, ending at np
907-
np = len(poles[j])
908-
den[j, np::-1] = polyfromroots(poles[j]).real
909-
denorder[j] = np
912+
# create the denominator matching this input
913+
# coefficients should be padded on right, ending at maxindex
914+
maxindex = len(poles[j])
915+
den[j, :maxindex+1] = poly(poles[j])
916+
denorder[j] = maxindex
910917

911918
# now create the numerator, also padded on the right
912919
for i in range(self.outputs):
@@ -915,22 +922,15 @@ def _common_den(self, imag_tol=None):
915922
# add all poles not found in the original denominator,
916923
# and the ones later added from other denominators
917924
for ip in chain(poleset[i][j][3],
918-
range(poleset[i][j][4], np)):
925+
range(poleset[i][j][4], maxindex)):
919926
nwzeros.append(poles[j][ip])
920927

921-
numpoly = poleset[i][j][2] * polyfromroots(nwzeros).real
922-
# print(numpoly, den[j])
923-
# polyfromroots gives coeffs in opposite order => invert
928+
numpoly = poleset[i][j][2] * np.atleast_1d(poly(nwzeros))
924929
# numerator polynomial should be padded on left and right
925-
# ending at np to line up with what td04ad expects...
926-
num[i, j, np + 1 - len(numpoly):np + 1] = numpoly[::-1]
930+
# ending at maxindex to line up with what td04ad expects.
931+
num[i, j, maxindex+1-len(numpoly):maxindex+1] = numpoly
927932
# print(num[i, j])
928933

929-
if (abs(den.imag) > epsnm).any():
930-
print("Warning: The denominator has a nontrivial imaginary part: "
931-
" %f" % abs(den.imag).max())
932-
den = den.real
933-
934934
return num, den, denorder
935935

936936
def sample(self, Ts, method='zoh', alpha=None):

0 commit comments

Comments
 (0)