From 3fb1a516501ce154cd818135d4fd0cf5687ff77e Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Tue, 26 Dec 2023 08:28:36 -0800 Subject: [PATCH 1/3] fix bug in matched transformation + address other issues in #950 --- control/dtime.py | 9 ++++---- control/freqplot.py | 6 ++++-- control/tests/discrete_test.py | 39 ++++++++++++++++++++++++++++++---- control/xferfcn.py | 15 ++++++++----- 4 files changed, 53 insertions(+), 16 deletions(-) diff --git a/control/dtime.py b/control/dtime.py index 2ae482811..9b91eabd3 100644 --- a/control/dtime.py +++ b/control/dtime.py @@ -55,8 +55,7 @@ # Sample a continuous time system def sample_system(sysc, Ts, method='zoh', alpha=None, prewarp_frequency=None, name=None, copy_names=True, **kwargs): - """ - Convert a continuous time system to discrete time by sampling. + """Convert a continuous time system to discrete time by sampling. Parameters ---------- @@ -67,9 +66,9 @@ def sample_system(sysc, Ts, method='zoh', alpha=None, prewarp_frequency=None, method : string Method to use for conversion, e.g. 'bilinear', 'zoh' (default) alpha : float within [0, 1] - The generalized bilinear transformation weighting parameter, which - should only be specified with method="gbt", and is ignored - otherwise. See :func:`scipy.signal.cont2discrete`. + The generalized bilinear transformation weighting parameter, which + should only be specified with method="gbt", and is ignored + otherwise. See :func:`scipy.signal.cont2discrete`. prewarp_frequency : float within [0, infinity) The frequency [rad/s] at which to match with the input continuous- time system's magnitude and phase (only valid for method='bilinear', diff --git a/control/freqplot.py b/control/freqplot.py index 164ec28a2..533515415 100644 --- a/control/freqplot.py +++ b/control/freqplot.py @@ -894,8 +894,10 @@ def gen_zero_centered_series(val_min, val_max, period): # list of systems (e.g., "Step response for sys[1], sys[2]"). # - # Set the initial title for the data (unique system names) - sysnames = list(set([response.sysname for response in data])) + # Set the initial title for the data (unique system names, preserving order) + seen = set() + sysnames = [response.sysname for response in data \ + if not (response.sysname in seen or seen.add(response.sysname))] if title is None: if data[0].title is None: title = "Bode plot for " + ", ".join(sysnames) diff --git a/control/tests/discrete_test.py b/control/tests/discrete_test.py index 96777011e..cccb53708 100644 --- a/control/tests/discrete_test.py +++ b/control/tests/discrete_test.py @@ -5,11 +5,12 @@ import numpy as np import pytest +import cmath -from control import (StateSpace, TransferFunction, bode, common_timebase, - feedback, forced_response, impulse_response, - isctime, isdtime, rss, c2d, sample_system, step_response, - timebase) +import control as ct +from control import StateSpace, TransferFunction, bode, common_timebase, \ + feedback, forced_response, impulse_response, isctime, isdtime, rss, \ + c2d, sample_system, step_response, timebase class TestDiscrete: @@ -526,3 +527,33 @@ def test_signal_names(self, tsys): assert sysd_newnames.find_input('u') is None assert sysd_newnames.find_output('y') == 0 assert sysd_newnames.find_output('x') is None + + +@pytest.mark.parametrize("num, den", [ + ([1], [1, 1]), + ([1, 2], [1, 3]), + ([1, 2], [3, 4, 5]) +]) +@pytest.mark.parametrize("dt", [True, 0.1, 2]) +@pytest.mark.parametrize("method", ['zoh', 'bilinear', 'matched']) +def test_c2d_matched(num, den, dt, method): + sys_ct = ct.tf(num, den) + sys_dt = ct.sample_system(sys_ct, dt, method=method) + assert sys_dt.dt == dt # make sure sampling time is OK + assert cmath.isclose(sys_ct(0), sys_dt(1)) # check zero frequency gain + assert cmath.isclose( + sys_ct.dcgain(), sys_dt.dcgain()) # another way to check + + if method in ['zoh', 'matched']: + # Make sure that poles were properly matched + zpoles = sys_dt.poles() + for cpole in sys_ct.poles(): + zpole = zpoles[(np.abs(zpoles - cmath.exp(cpole * dt))).argmin()] + assert cmath.isclose(cmath.exp(cpole * dt), zpole) + + if method in ['matched']: + # Make sure that zeros were properly matched + zzeros = sys_dt.zeros() + for czero in sys_ct.zeros(): + zzero = zzeros[(np.abs(zzeros - cmath.exp(czero * dt))).argmin()] + assert cmath.isclose(cmath.exp(czero * dt), zzero) diff --git a/control/xferfcn.py b/control/xferfcn.py index b5334b7b8..8f56c96a0 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -1181,7 +1181,9 @@ def sample(self, Ts, method='zoh', alpha=None, prewarp_frequency=None, if not self.issiso(): raise ControlMIMONotImplemented("Not implemented for MIMO systems") if method == "matched": - return _c2d_matched(self, Ts) + if prewarp_frequency is not None: + warn('prewarp_frequency ignored: incompatible conversion') + return _c2d_matched(self, Ts, name=name, **kwargs) sys = (self.num[0][0], self.den[0][0]) if prewarp_frequency is not None: if method in ('bilinear', 'tustin') or \ @@ -1284,9 +1286,12 @@ def _isstatic(self): # c2d function contributed by Benjamin White, Oct 2012 -def _c2d_matched(sysC, Ts): +def _c2d_matched(sysC, Ts, **kwargs): + if sysC.ninputs > 1 or sysC.noutputs > 1: + raise ValueError("MIMO transfer functions not supported") + # Pole-zero match method of continuous to discrete time conversion - szeros, spoles, sgain = tf2zpk(sysC.num[0][0], sysC.den[0][0]) + szeros, spoles, _ = tf2zpk(sysC.num[0][0], sysC.den[0][0]) zzeros = [0] * len(szeros) zpoles = [0] * len(spoles) pregainnum = [0] * len(szeros) @@ -1302,9 +1307,9 @@ def _c2d_matched(sysC, Ts): zpoles[idx] = z pregainden[idx] = 1 - z zgain = np.multiply.reduce(pregainnum) / np.multiply.reduce(pregainden) - gain = sgain / zgain + gain = sysC.dcgain() / zgain sysDnum, sysDden = zpk2tf(zzeros, zpoles, gain) - return TransferFunction(sysDnum, sysDden, Ts) + return TransferFunction(sysDnum, sysDden, Ts, **kwargs) # Utility function to convert a transfer function polynomial to a string From d4525a1dac0de5e7f9c84e1ada6a7e48cfd276af Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Tue, 26 Dec 2023 11:11:31 -0800 Subject: [PATCH 2/3] Update control/xferfcn.py Co-authored-by: Sawyer Fuller <58706249+sawyerbfuller@users.noreply.github.com> --- control/xferfcn.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/control/xferfcn.py b/control/xferfcn.py index 8f56c96a0..59923fff9 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -1287,7 +1287,8 @@ def _isstatic(self): # c2d function contributed by Benjamin White, Oct 2012 def _c2d_matched(sysC, Ts, **kwargs): - if sysC.ninputs > 1 or sysC.noutputs > 1: + if not sysC.issiso(): + raise ControlMIMONotImplemented("Not implemented for MIMO systems") raise ValueError("MIMO transfer functions not supported") # Pole-zero match method of continuous to discrete time conversion From e4bcec01cda450ab12a44b456c8351d8af1d83ce Mon Sep 17 00:00:00 2001 From: Richard Murray Date: Tue, 26 Dec 2023 11:18:44 -0800 Subject: [PATCH 3/3] fixed botched pathc from github --- control/xferfcn.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/control/xferfcn.py b/control/xferfcn.py index 59923fff9..206f0bc92 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -1288,8 +1288,7 @@ def _isstatic(self): # c2d function contributed by Benjamin White, Oct 2012 def _c2d_matched(sysC, Ts, **kwargs): if not sysC.issiso(): - raise ControlMIMONotImplemented("Not implemented for MIMO systems") - raise ValueError("MIMO transfer functions not supported") + raise ControlMIMONotImplemented("Not implemented for MIMO systems") # Pole-zero match method of continuous to discrete time conversion szeros, spoles, _ = tf2zpk(sysC.num[0][0], sysC.den[0][0])