Skip to content

added ability to 'prewarp' the conversion of continuous to discrete-time systems #417

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

Merged
murrayrm merged 2 commits into from
Jul 11, 2020
Merged

added ability to 'prewarp' the conversion of continuous to discrete-time systems #417

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8d30b44
added ability to 'prewarp' the conversion of continuous to discrete-t…
sawyerbfuller Jun 11, 2020
03c04e2
adjusted prewarping to only apply to 'bilinear' cont-to-discrete appr…
sawyerbfuller Jun 11, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 13 additions & 4 deletions control/dtime.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@
__all__ = ['sample_system', 'c2d']

# Sample a continuous time system
def sample_system(sysc, Ts, method='zoh', alpha=None):
def sample_system(sysc, Ts, method='zoh', alpha=None, prewarp_frequency=None):
"""Convert a continuous time system to discrete time

Creates a discrete time system from a continuous time system by
Expand All @@ -67,6 +67,10 @@ def sample_system(sysc, Ts, method='zoh', alpha=None):
method : string
Method to use for conversion: 'matched', 'tustin', 'zoh' (default)

prewarp_frequency : float within [0, infinity)
The frequency [rad/s] at which to match with the input continuous-
time system's magnitude and phase

Returns
-------
sysd : linsys
Expand All @@ -87,10 +91,10 @@ def sample_system(sysc, Ts, method='zoh', alpha=None):
if not isctime(sysc):
raise ValueError("First argument must be continuous time system")

return sysc.sample(Ts, method, alpha)
return sysc.sample(Ts, method, alpha, prewarp_frequency)


def c2d(sysc, Ts, method='zoh'):
def c2d(sysc, Ts, method='zoh', prewarp_frequency=None):
'''
Return a discrete-time system

Expand All @@ -109,9 +113,14 @@ def c2d(sysc, Ts, method='zoh'):
'impulse' Impulse-invariant discretization, currently not implemented
'tustin' Bilinear (Tustin) approximation, only SISO
'matched' Matched pole-zero method, only SISO

prewarp_frequency : float within [0, infinity)
The frequency [rad/s] at which to match with the input continuous-
time system's magnitude and phase

'''
# Call the sample_system() function to do the work
sysd = sample_system(sysc, Ts, method)
sysd = sample_system(sysc, Ts, method, prewarp_frequency)

# TODO: is this check needed? If sysc is StateSpace, sysd is too?
if isinstance(sysc, StateSpace) and not isinstance(sysd, StateSpace):
Expand Down
17 changes: 14 additions & 3 deletions control/statesp.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -841,7 +841,7 @@ def __getitem__(self, indices):
j = indices[1]
return StateSpace(self.A, self.B[:, j], self.C[i, :], self.D[i, j], self.dt)

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

Creates a discrete-time system from a continuous-time system by
Expand All @@ -866,6 +866,12 @@ def sample(self, Ts, method='zoh', alpha=None):
should only be specified with method="gbt", and is ignored
otherwise

prewarp_frequency : float within [0, infinity)
The frequency [rad/s] at which to match with the input continuous-
time system's magnitude and phase (the gain=1 crossover frequency,
for example). Should only be specified with method='bilinear' or
'gbt' with alpha=0.5 and ignored otherwise.

Returns
-------
sysd : StateSpace
Expand All @@ -885,8 +891,13 @@ def sample(self, Ts, method='zoh', alpha=None):
raise ValueError("System must be continuous time system")

sys = (self.A, self.B, self.C, self.D)
Ad, Bd, C, D, dt = cont2discrete(sys, Ts, method, alpha)
return StateSpace(Ad, Bd, C, D, dt)
if (method=='bilinear' or (method=='gbt' and alpha==0.5)) and \
prewarp_frequency is not None:
Twarp = 2*np.tan(prewarp_frequency*Ts/2)/prewarp_frequency
else:
Twarp = Ts
Ad, Bd, C, D, _ = cont2discrete(sys, Twarp, method, alpha)
return StateSpace(Ad, Bd, C, D, Ts)

def dcgain(self):
"""Return the zero-frequency gain
Expand Down
18 changes: 18 additions & 0 deletions control/tests/statesp_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -611,6 +611,24 @@ def test_copy_constructor(self):
linsys.A[0, 0] = -3
np.testing.assert_array_equal(cpysys.A, [[-1]]) # original value

def test_sample_system_prewarping(self):
"""test that prewarping works when converting from cont to discrete time system"""
A = np.array([
[ 0.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[-3.81097561e+01, -1.12500000e+00, 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, -1.66356135e+04, -1.34748470e+01]])
B = np.array([
[ 0. ], [ 38.1097561 ],[ 0. ],[16635.61352143]])
C = np.array([[0.90909091, 0. , 0.09090909, 0. ],])
wwarp = 50
Ts = 0.025
plant = StateSpace(A,B,C,0)
plant_d_warped = plant.sample(Ts, 'bilinear', prewarp_frequency=wwarp)
np.testing.assert_array_almost_equal(
evalfr(plant, wwarp*1j),
evalfr(plant_d_warped, np.exp(wwarp*1j*Ts)),
decimal=4)

if __name__ == "__main__":
unittest.main()
21 changes: 20 additions & 1 deletion control/tests/xferfcn_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -853,7 +853,26 @@ def test_latex_repr(self):
r'+ 2.3 ' + expmul + ' 10^{-45}'
r'}' + suffix + '$$')
self.assertEqual(H._repr_latex_(), ref)


def test_sample_system_prewarping(self):
"""test that prewarping works when converting from cont to discrete time system"""
A = np.array([
[ 0.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[-3.81097561e+01, -1.12500000e+00, 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00],
[ 0.00000000e+00, 0.00000000e+00, -1.66356135e+04, -1.34748470e+01]])
B = np.array([
[ 0. ], [ 38.1097561 ],[ 0. ],[16635.61352143]])
C = np.array([[0.90909091, 0. , 0.09090909, 0. ],])
wwarp = 50
Ts = 0.025
plant = StateSpace(A,B,C,0)
plant = ss2tf(plant)
plant_d_warped = plant.sample(Ts, 'bilinear', prewarp_frequency=wwarp)
np.testing.assert_array_almost_equal(
evalfr(plant, wwarp*1j),
evalfr(plant_d_warped, np.exp(wwarp*1j*Ts)),
decimal=4)

if __name__ == "__main__":
unittest.main()
17 changes: 14 additions & 3 deletions control/xferfcn.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -961,7 +961,7 @@ def _common_den(self, imag_tol=None, allow_nonproper=False):

return num, den, denorder

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

Creates a discrete-time system from a continuous-time system by
Expand All @@ -985,6 +985,12 @@ def sample(self, Ts, method='zoh', alpha=None):
The generalized bilinear transformation weighting parameter, which
should only be specified with method="gbt", and is ignored
otherwise.

prewarp_frequency : float within [0, infinity)
The frequency [rad/s] at which to match with the input continuous-
time system's magnitude and phase (the gain=1 crossover frequency,
for example). Should only be specified with method='bilinear' or
'gbt' with alpha=0.5 and ignored otherwise.

Returns
-------
Expand All @@ -1010,8 +1016,13 @@ def sample(self, Ts, method='zoh', alpha=None):
if method == "matched":
return _c2d_matched(self, Ts)
sys = (self.num[0][0], self.den[0][0])
numd, dend, dt = cont2discrete(sys, Ts, method, alpha)
return TransferFunction(numd[0, :], dend, dt)
if (method=='bilinear' or (method=='gbt' and alpha==0.5)) and \
prewarp_frequency is not None:
Twarp = 2*np.tan(prewarp_frequency*Ts/2)/prewarp_frequency
else:
Twarp = Ts
numd, dend, _ = cont2discrete(sys, Twarp, method, alpha)
return TransferFunction(numd[0, :], dend, Ts)

def dcgain(self):
"""Return the zero-frequency (or DC) gain
Expand Down