From 433c1369876ec23a57e386a765995f59f98ac1c7 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Tue, 2 Nov 2021 23:26:56 +0100 Subject: [PATCH 1/2] return frequency response for 0 and 1-state systems directly --- control/statesp.py | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/control/statesp.py b/control/statesp.py index 6b3a1dff3..85eadb506 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -851,7 +851,7 @@ def slycot_laub(self, x): # transformed state matrices, at, bt, ct. # Start at the second frequency, already have the first. - for kk, x_kk in enumerate(x_arr[1:len(x_arr)]): + for kk, x_kk in enumerate(x_arr[1:]): result = tb05ad(n, m, p, x_kk, at, bt, ct, job='NH') # When job='NH', result = (g_i, hinvb, info) @@ -885,15 +885,27 @@ def horner(self, x, warn_infinite=True): Attempts to use Laub's method from Slycot library, with a fall-back to python code. """ + # Make sure the argument is a 1D array of complex numbers + x_arr = np.atleast_1d(x).astype(complex, copy=False) + + # return fast on systems with 0 or 1 state + if self.nstates == 0: + return self.D[:, :, np.newaxis] \ + * np.ones_like(x_arr, dtype=complex) + if self.nstates == 1: + with np.errstate(divide='ignore', invalid='ignore'): + out = (self.C[:, :, np.newaxis] + * (self.B[:, :, np.newaxis] / (x_arr - self.A[0, 0])) + + self.D[:, :, np.newaxis]) + out[np.isnan(out)] = complex(np.inf, np.nan) + return out + try: - out = self.slycot_laub(x) + out = self.slycot_laub(x_arr) except (ImportError, Exception): # Fall back because either Slycot unavailable or cannot handle # certain cases. - # Make sure the argument is a 1D array of complex numbers - x_arr = np.atleast_1d(x).astype(complex, copy=False) - # Make sure that we are operating on a simple list if len(x_arr.shape) > 1: raise ValueError("input list must be 1D") From 5355b025f90ed0838f3e540f8983aa37ce504a43 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Tue, 2 Nov 2021 23:45:20 +0100 Subject: [PATCH 2/2] only return fast in array, not matrix --- control/statesp.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/control/statesp.py b/control/statesp.py index 85eadb506..f04bc55b3 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -889,16 +889,18 @@ def horner(self, x, warn_infinite=True): x_arr = np.atleast_1d(x).astype(complex, copy=False) # return fast on systems with 0 or 1 state - if self.nstates == 0: - return self.D[:, :, np.newaxis] \ - * np.ones_like(x_arr, dtype=complex) - if self.nstates == 1: - with np.errstate(divide='ignore', invalid='ignore'): - out = (self.C[:, :, np.newaxis] - * (self.B[:, :, np.newaxis] / (x_arr - self.A[0, 0])) - + self.D[:, :, np.newaxis]) - out[np.isnan(out)] = complex(np.inf, np.nan) - return out + if not config.defaults['statesp.use_numpy_matrix']: + if self.nstates == 0: + return self.D[:, :, np.newaxis] \ + * np.ones_like(x_arr, dtype=complex) + if self.nstates == 1: + with np.errstate(divide='ignore', invalid='ignore'): + out = self.C[:, :, np.newaxis] \ + / (x_arr - self.A[0, 0]) \ + * self.B[:, :, np.newaxis] \ + + self.D[:, :, np.newaxis] + out[np.isnan(out)] = complex(np.inf, np.nan) + return out try: out = self.slycot_laub(x_arr)