Movatterモバイル変換

[0]ホーム
URL:

Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly

 Sign up
Appearance settings

Commitcb25633

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

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(notslycot_check(),"slycot not installed")
499+
500+
deftest_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
deftest_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(notslycot_check(),"slycot not installed")
512560
deftest_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
deftest_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
importscipyassp
60-
fromnumpy.polynomial.polynomialimportpolyfromroots
6160
fromscipy.signalimportlti,tf2zpk,zpk2tf,cont2discrete
6261
fromcopyimportdeepcopy
6362
fromwarningsimportwarn
@@ -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 thetolerance to use in decidingof a pole is complex
844+
#Thetolerance to use in decidingif a pole is complex
844845
if (imag_tolisNone):
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= [[]forjinrange(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
forjinrange(self.inputs):
872873
foriinrange(self.outputs):
873874
currentpoles=poleset[i][j][1]
874875
nothave=ones(currentpoles.shape,dtype=bool)
875876
forip,pinenumerate(poles[j]):
876-
idx,=nonzero(
877-
(abs(currentpoles-p)<epsnm)*nothave)
878-
iflen(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+
ifnp.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
forh,cinzip(nothave,currentpoles):
884889
ifh:
890+
ifabs(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)forpinpoles])
891-
den=zeros((self.inputs,npmax+1),dtype=float)
897+
maxindex=max([len(p)forpinpoles])
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
foriinrange(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
foriinrange(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
foripinchain(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 atnp to line up with what td04ad expects...
926-
num[i,j,np+1-len(numpoly):np+1]=numpoly[::-1]
930+
# ending atmaxindex 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
returnnum,den,denorder
935935

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

0 commit comments

Comments
 (0)
[8]ページ先頭
©2009-2025 Movatter.jp