Skip to content

Commit 8366806

Browse files
sawyerbfullersawyerbfuller
authored
added ability to 'prewarp' the conversion of continuous to discrete-time systems (#417)
* added ability to 'prewarp' the conversion of continuous to discrete-time systems (in functions sample, sample_system, and c2d) , so that their gain and phase match at a specific desired frequency.
1 parent 155d8eb commit 8366806

5 files changed

Lines changed: 79 additions & 11 deletions

File tree

‎control/dtime.py‎

Lines changed: 13 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,7 @@
5252
__all__ = ['sample_system', 'c2d']
5353

5454
# Sample a continuous time system
55-
def sample_system(sysc, Ts, method='zoh', alpha=None):
55+
def sample_system(sysc, Ts, method='zoh', alpha=None, prewarp_frequency=None):
5656
"""Convert a continuous time system to discrete time
5757
5858
Creates a discrete time system from a continuous time system by
@@ -67,6 +67,10 @@ def sample_system(sysc, Ts, method='zoh', alpha=None):
6767
method : string
6868
Method to use for conversion: 'matched', 'tustin', 'zoh' (default)
6969
70+
prewarp_frequency : float within [0, infinity)
71+
The frequency [rad/s] at which to match with the input continuous-
72+
time system's magnitude and phase
73+
7074
Returns
7175
-------
7276
sysd : linsys
@@ -87,10 +91,10 @@ def sample_system(sysc, Ts, method='zoh', alpha=None):
8791
if not isctime(sysc):
8892
raise ValueError("First argument must be continuous time system")
8993

90-
return sysc.sample(Ts, method, alpha)
94+
return sysc.sample(Ts, method, alpha, prewarp_frequency)
9195

9296

93-
def c2d(sysc, Ts, method='zoh'):
97+
def c2d(sysc, Ts, method='zoh', prewarp_frequency=None):
9498
'''
9599
Return a discrete-time system
96100
@@ -109,9 +113,14 @@ def c2d(sysc, Ts, method='zoh'):
109113
'impulse' Impulse-invariant discretization, currently not implemented
110114
'tustin' Bilinear (Tustin) approximation, only SISO
111115
'matched' Matched pole-zero method, only SISO
116+
117+
prewarp_frequency : float within [0, infinity)
118+
The frequency [rad/s] at which to match with the input continuous-
119+
time system's magnitude and phase
120+
112121
'''
113122
# Call the sample_system() function to do the work
114-
sysd = sample_system(sysc, Ts, method)
123+
sysd = sample_system(sysc, Ts, method, prewarp_frequency)
115124

116125
# TODO: is this check needed? If sysc is StateSpace, sysd is too?
117126
if isinstance(sysc, StateSpace) and not isinstance(sysd, StateSpace):

‎control/statesp.py‎

Lines changed: 14 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -837,7 +837,7 @@ def __getitem__(self, indices):
837837
j = indices[1]
838838
return StateSpace(self.A, self.B[:, j], self.C[i, :], self.D[i, j], self.dt)
839839

840-
def sample(self, Ts, method='zoh', alpha=None):
840+
def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None):
841841
"""Convert a continuous time system to discrete time
842842
843843
Creates a discrete-time system from a continuous-time system by
@@ -862,6 +862,12 @@ def sample(self, Ts, method='zoh', alpha=None):
862862
should only be specified with method="gbt", and is ignored
863863
otherwise
864864
865+
prewarp_frequency : float within [0, infinity)
866+
The frequency [rad/s] at which to match with the input continuous-
867+
time system's magnitude and phase (the gain=1 crossover frequency,
868+
for example). Should only be specified with method='bilinear' or
869+
'gbt' with alpha=0.5 and ignored otherwise.
870+
865871
Returns
866872
-------
867873
sysd : StateSpace
@@ -881,8 +887,13 @@ def sample(self, Ts, method='zoh', alpha=None):
881887
raise ValueError("System must be continuous time system")
882888

883889
sys = (self.A, self.B, self.C, self.D)
884-
Ad, Bd, C, D, dt = cont2discrete(sys, Ts, method, alpha)
885-
return StateSpace(Ad, Bd, C, D, dt)
890+
if (method=='bilinear' or (method=='gbt' and alpha==0.5)) and \
891+
prewarp_frequency is not None:
892+
Twarp = 2*np.tan(prewarp_frequency*Ts/2)/prewarp_frequency
893+
else:
894+
Twarp = Ts
895+
Ad, Bd, C, D, _ = cont2discrete(sys, Twarp, method, alpha)
896+
return StateSpace(Ad, Bd, C, D, Ts)
886897

887898
def dcgain(self):
888899
"""Return the zero-frequency gain

‎control/tests/statesp_test.py‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -634,6 +634,24 @@ def test_copy_constructor(self):
634634
linsys.A[0, 0] = -3
635635
np.testing.assert_array_equal(cpysys.A, [[-1]]) # original value
636636

637+
def test_sample_system_prewarping(self):
638+
"""test that prewarping works when converting from cont to discrete time system"""
639+
A = np.array([
640+
[ 0.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00],
641+
[-3.81097561e+01, -1.12500000e+00, 0.00000000e+00, 0.00000000e+00],
642+
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00],
643+
[ 0.00000000e+00, 0.00000000e+00, -1.66356135e+04, -1.34748470e+01]])
644+
B = np.array([
645+
[ 0. ], [ 38.1097561 ],[ 0. ],[16635.61352143]])
646+
C = np.array([[0.90909091, 0. , 0.09090909, 0. ],])
647+
wwarp = 50
648+
Ts = 0.025
649+
plant = StateSpace(A,B,C,0)
650+
plant_d_warped = plant.sample(Ts, 'bilinear', prewarp_frequency=wwarp)
651+
np.testing.assert_array_almost_equal(
652+
evalfr(plant, wwarp*1j),
653+
evalfr(plant_d_warped, np.exp(wwarp*1j*Ts)),
654+
decimal=4)
637655

638656
if __name__ == "__main__":
639657
unittest.main()

‎control/tests/xferfcn_test.py‎

Lines changed: 20 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -872,7 +872,26 @@ def test_latex_repr(self):
872872
r'+ 2.3 ' + expmul + ' 10^{-45}'
873873
r'}' + suffix + '$$')
874874
self.assertEqual(H._repr_latex_(), ref)
875-
875+
876+
def test_sample_system_prewarping(self):
877+
"""test that prewarping works when converting from cont to discrete time system"""
878+
A = np.array([
879+
[ 0.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00],
880+
[-3.81097561e+01, -1.12500000e+00, 0.00000000e+00, 0.00000000e+00],
881+
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00],
882+
[ 0.00000000e+00, 0.00000000e+00, -1.66356135e+04, -1.34748470e+01]])
883+
B = np.array([
884+
[ 0. ], [ 38.1097561 ],[ 0. ],[16635.61352143]])
885+
C = np.array([[0.90909091, 0. , 0.09090909, 0. ],])
886+
wwarp = 50
887+
Ts = 0.025
888+
plant = StateSpace(A,B,C,0)
889+
plant = ss2tf(plant)
890+
plant_d_warped = plant.sample(Ts, 'bilinear', prewarp_frequency=wwarp)
891+
np.testing.assert_array_almost_equal(
892+
evalfr(plant, wwarp*1j),
893+
evalfr(plant_d_warped, np.exp(wwarp*1j*Ts)),
894+
decimal=4)
876895

877896
if __name__ == "__main__":
878897
unittest.main()

‎control/xferfcn.py‎

Lines changed: 14 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -983,7 +983,7 @@ def _common_den(self, imag_tol=None, allow_nonproper=False):
983983

984984
return num, den, denorder
985985

986-
def sample(self, Ts, method='zoh', alpha=None):
986+
def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None):
987987
"""Convert a continuous-time system to discrete time
988988
989989
Creates a discrete-time system from a continuous-time system by
@@ -1007,6 +1007,12 @@ def sample(self, Ts, method='zoh', alpha=None):
10071007
The generalized bilinear transformation weighting parameter, which
10081008
should only be specified with method="gbt", and is ignored
10091009
otherwise.
1010+
1011+
prewarp_frequency : float within [0, infinity)
1012+
The frequency [rad/s] at which to match with the input continuous-
1013+
time system's magnitude and phase (the gain=1 crossover frequency,
1014+
for example). Should only be specified with method='bilinear' or
1015+
'gbt' with alpha=0.5 and ignored otherwise.
10101016
10111017
Returns
10121018
-------
@@ -1032,8 +1038,13 @@ def sample(self, Ts, method='zoh', alpha=None):
10321038
if method == "matched":
10331039
return _c2d_matched(self, Ts)
10341040
sys = (self.num[0][0], self.den[0][0])
1035-
numd, dend, dt = cont2discrete(sys, Ts, method, alpha)
1036-
return TransferFunction(numd[0, :], dend, dt)
1041+
if (method=='bilinear' or (method=='gbt' and alpha==0.5)) and \
1042+
prewarp_frequency is not None:
1043+
Twarp = 2*np.tan(prewarp_frequency*Ts/2)/prewarp_frequency
1044+
else:
1045+
Twarp = Ts
1046+
numd, dend, _ = cont2discrete(sys, Twarp, method, alpha)
1047+
return TransferFunction(numd[0, :], dend, Ts)
10371048

10381049
def dcgain(self):
10391050
"""Return the zero-frequency (or DC) gain

0 commit comments

Comments
 (0)