From 9658e0f55759b5c511bb152941329df47c978975 Mon Sep 17 00:00:00 2001 From: Sawyer Fuller Date: Mon, 17 Aug 2020 15:49:16 -0700 Subject: [PATCH 1/7] beginnings of code starts with a docstring --- control/timeresp.py | 52 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 51 insertions(+), 1 deletion(-) diff --git a/control/timeresp.py b/control/timeresp.py index fa4ced2bd..f46d23569 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -826,7 +826,57 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=None, return T, yout # utility function to find time period and time increment using pole locations -def _ideal_tfinal_and_dt(sys): +def _ideal_tfinal_and_dt(sys, is_step=True): + """helper function to compute ideal simulation duration tfinal and dt, the + time increment. Usually called by _default_time_vector, whose job it is to + bring reality into the picture. + + For discrete-time models, dt is inherent and only tfinal is computed. + + Parameters + ---------- + sys : StateSpace or TransferFunction + The system whose time response is to be computed + is_step : bool + Scales the dc value by the magnitude of the nonzero mode since + integrating the impulse response gives ∫exp(-λt) = -exp(-λt)/λ. + Default is True. + Returns + ------- + tfinal : float + The final time instance for which the simulation will be performed. + dt : float + The estimated sampling period for the simulation. + + Notes + ----- + Just by evaluating the fastest mode for dt and slowest for tfinal often + leads to unnecessary, bloated sampling (e.g., Transfer(1,[1,1001,1000])) + since dt will be very small and tfinal will be too large though the fast + mode hardly ever contributes. Similarly, change the numerator to [1, 2, 0] + and the simulation would be unnecessarily long and the plot is virtually + an L shape since the decay is so fast. + + Instead, a modal decomposition in time domain hence a truncated ZIR and ZSR + can be used such that only the modes that have significant effect on the + time response are taken. But the sensitivity of the eigenvalues complicate + the matter since dλ = with = 1. Hence we can only work + with simple poles with this formulation. See Golub, Van Loan Section 7.2.2 + for simple eigenvalue sensitivity about the nonunity of . The size of + the response is dependent on the size of the eigenshapes rather than the + eigenvalues themselves. + + By Ilhan Polat, with modifications by Sawyer Fuller to integrate into + python-control 2020.08.17 + """ + + sqrt_eps = np.sqrt(np.spacing(1.)) + default_tfinal = 5 # Default simulation horizon + total_cycles = 5 # number of cycles for oscillating modes + pts_per_cycle = 25 # Number of points divide a period of oscillation + log_decay_percent = np.log(100) # Factor of reduction for real pole decays + + if sys.is_static_gain(): pass constant = 7.0 tolerance = 1e-10 A = ssdata(sys)[0] From afdc5e34857a79c8bef438e136075cb53034188d Mon Sep 17 00:00:00 2001 From: bnavigator Date: Tue, 18 Aug 2020 13:37:06 +0200 Subject: [PATCH 2/7] first commit that adds ilayn's default time for step code --- control/statesp.py | 4 + control/tests/sisotool_test.py | 14 +- control/tests/timeresp_test.py | 65 ++++++--- control/timeresp.py | 260 ++++++++++++++++++++++----------- control/xferfcn.py | 12 +- 5 files changed, 237 insertions(+), 118 deletions(-) diff --git a/control/statesp.py b/control/statesp.py index ca68fc22b..dd0ea6f5e 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -944,6 +944,10 @@ def dcgain(self): gain = np.tile(np.nan, (self.outputs, self.inputs)) return np.squeeze(gain) + def is_static_gain(self): + """True if and only if the system has no dynamics, that is, + if A and B are zero. """ + return not np.any(self.A) and not np.any(self.B) # TODO: add discrete time check def _convertToStateSpace(sys, **kw): diff --git a/control/tests/sisotool_test.py b/control/tests/sisotool_test.py index f93de54f8..5b627c22d 100644 --- a/control/tests/sisotool_test.py +++ b/control/tests/sisotool_test.py @@ -32,9 +32,12 @@ def test_sisotool(self): initial_point_2, 4) # Check the step response before moving the point + # new array needed because change in compute step response default time step_response_original = np.array( - [0., 0.0217, 0.1281, 0.3237, 0.5797, 0.8566, 1.116, - 1.3261, 1.4659, 1.526]) + [0. , 0.0069, 0.0448, 0.124 , 0.2427, 0.3933, 0.5653, 0.7473, + 0.928 , 1.0969]) + #old: np.array([0., 0.0217, 0.1281, 0.3237, 0.5797, 0.8566, 1.116, + # 1.3261, 1.4659, 1.526]) assert_array_almost_equal( ax_step.lines[0].get_data()[1][:10], step_response_original, 4) @@ -77,9 +80,12 @@ def test_sisotool(self): bode_mag_moved, 4) # Check if the step response has changed + # new array needed because change in compute step response default time step_response_moved = np.array( - [0., 0.0239, 0.161 , 0.4547, 0.8903, 1.407, - 1.9121, 2.2989, 2.4686, 2.353]) + [0. , 0.0072, 0.0516, 0.1554, 0.3281, 0.5681, 0.8646, 1.1987, + 1.5452, 1.875 ]) + #old: array([0., 0.0239, 0.161 , 0.4547, 0.8903, 1.407, + # 1.9121, 2.2989, 2.4686, 2.353]) assert_array_almost_equal( ax_step.lines[0].get_data()[1][:10], step_response_moved, 4) diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index 5549b2a88..b33dd5969 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -29,6 +29,12 @@ def setUp(self): # Create some transfer functions self.siso_tf1 = TransferFunction([1], [1, 2, 1]) self.siso_tf2 = _convert_to_transfer_function(self.siso_ss1) + + # tests for pole cancellation + self.pole_cancellation = TransferFunction([1.067e+05, 5.791e+04], + [10.67, 1.067e+05, 5.791e+04]) + self.no_pole_cancellation = TransferFunction([1.881e+06], + [188.1, 1.881e+06]) # Create MIMO system, contains ``siso_ss1`` twice A = np.matrix("1. -2. 0. 0.;" @@ -167,6 +173,14 @@ def test_step_info(self): 2.50, rtol=rtol) + # confirm that pole-zero cancellation doesn't perturb results + # https://github.com/python-control/python-control/issues/440 + step_info_no_cancellation = step_info(self.no_pole_cancellation) + step_info_cancellation = step_info(self.pole_cancellation) + for key in step_info_no_cancellation: + np.testing.assert_allclose(step_info_no_cancellation[key], + step_info_cancellation[key], rtol=1e-4) + def test_impulse_response(self): # Test SISO system sys = self.siso_ss1 @@ -348,33 +362,41 @@ def test_step_robustness(self): sys2 = TransferFunction(num, den2) # Compute step response from input 1 to output 1, 2 - t1, y1 = step_response(sys1, input=0, T_num=100) - t2, y2 = step_response(sys2, input=0, T_num=100) + t1, y1 = step_response(sys1, input=0, T=2, T_num=100) + t2, y2 = step_response(sys2, input=0, T=2, T_num=100) np.testing.assert_array_almost_equal(y1, y2) def test_auto_generated_time_vector(self): - # confirm a TF with a pole at p simulates for 7.0/p seconds + # confirm a TF with a pole at p simulates for ratio/p seconds p = 0.5 + ratio = 9.21034*p # taken from code + ratio2 = 25*p np.testing.assert_array_almost_equal( _ideal_tfinal_and_dt(TransferFunction(1, [1, .5]))[0], - (7/p)) + (ratio/p)) np.testing.assert_array_almost_equal( _ideal_tfinal_and_dt(TransferFunction(1, [1, .5]).sample(.1))[0], - (7/p)) - # confirm a TF with poles at 0 and p simulates for 7.0/p seconds + (ratio2/p)) + # confirm a TF with poles at 0 and p simulates for ratio/p seconds np.testing.assert_array_almost_equal( _ideal_tfinal_and_dt(TransferFunction(1, [1, .5, 0]))[0], - (7/p)) + (ratio2/p)) + # confirm a TF with a natural frequency of wn rad/s gets a - # dt of 1/(7.0*wn) + # dt of 1/(ratio*wn) wn = 10 + ratio_dt = 1/(0.025133 * ratio * wn) np.testing.assert_array_almost_equal( _ideal_tfinal_and_dt(TransferFunction(1, [1, 0, wn**2]))[1], - 1/(7.0*wn)) + 1/(ratio_dt*ratio*wn)) + wn = 100 + np.testing.assert_array_almost_equal( + _ideal_tfinal_and_dt(TransferFunction(1, [1, 0, wn**2]))[1], + 1/(ratio_dt*ratio*wn)) zeta = .1 np.testing.assert_array_almost_equal( _ideal_tfinal_and_dt(TransferFunction(1, [1, 2*zeta*wn, wn**2]))[1], - 1/(7.0*wn)) + 1/(ratio_dt*ratio*wn)) # but a smapled one keeps its dt np.testing.assert_array_almost_equal( _ideal_tfinal_and_dt(TransferFunction(1, [1, 2*zeta*wn, wn**2]).sample(.1))[1], @@ -384,37 +406,32 @@ def test_auto_generated_time_vector(self): .1) np.testing.assert_array_almost_equal( _ideal_tfinal_and_dt(TransferFunction(1, [1, 2*zeta*wn, wn**2]))[1], - 1/(7.0*wn)) + 1/(ratio_dt*ratio*wn)) + + # TF with fast oscillations simulates only 5000 time steps even with long tfinal self.assertEqual(5000, len(_default_time_vector(TransferFunction(1, [1, 0, wn**2]),tfinal=100))) - # and simulates for 7.0/dt time steps - self.assertEqual( - len(_default_time_vector(TransferFunction(1, [1, 0, wn**2]))), - int(7.0/(1/(7.0*wn)))) sys = TransferFunction(1, [1, .5, 0]) sysdt = TransferFunction(1, [1, .5, 0], .1) # test impose number of time steps self.assertEqual(10, len(step_response(sys, T_num=10)[0])) - self.assertEqual(10, len(step_response(sysdt, T_num=10)[0])) + # test that discrete ignores T_num + self.assertNotEqual(15, len(step_response(sysdt, T_num=15)[0])) # test impose final time np.testing.assert_array_almost_equal( 100, - step_response(sys, 100)[0][-1], - decimal=.5) + np.ceil(step_response(sys, 100)[0][-1])) np.testing.assert_array_almost_equal( 100, - step_response(sysdt, 100)[0][-1], - decimal=.5) + np.ceil(step_response(sysdt, 100)[0][-1])) np.testing.assert_array_almost_equal( 100, - impulse_response(sys, 100)[0][-1], - decimal=.5) + np.ceil(impulse_response(sys, 100)[0][-1])) np.testing.assert_array_almost_equal( 100, - initial_response(sys, 100)[0][-1], - decimal=.5) + np.ceil(initial_response(sys, 100)[0][-1])) def test_time_vector(self): diff --git a/control/timeresp.py b/control/timeresp.py index f46d23569..42e4cf125 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -71,6 +71,9 @@ # Libraries that we make use of import scipy as sp # SciPy library (used all over) import numpy as np # NumPy library +from scipy.linalg import eig, eigvals, matrix_balance, norm +from numpy import (einsum, maximum, minimum, + atleast_1d) import warnings from .lti import LTI # base class of StateSpace, TransferFunction from .statesp import _convertToStateSpace, _mimo2simo, _mimo2siso, ssdata @@ -84,7 +87,7 @@ def _check_convert_array(in_obj, legal_shapes, err_msg_start, squeeze=False, transpose=False): """ - Helper function for checking array-like parameters. + Helper function for checking array_like parameters. * Check type and shape of ``in_obj``. * Convert ``in_obj`` to an array if necessary. @@ -201,20 +204,20 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, Parameters ---------- - sys: LTI (StateSpace, or TransferFunction) + sys: LTI (StateSpace or TransferFunction) LTI system to simulate - T: array-like, optional for discrete LTI `sys` + T: array_like, optional for discrete LTI `sys` Time steps at which the input is defined; values must be evenly spaced. - U: array-like or number, optional + U: array_like or float, optional Input array giving input at each time `T` (default = 0). If `U` is ``None`` or ``0``, a special algorithm is used. This special algorithm is faster than the general algorithm, which is used otherwise. - X0: array-like or number, optional + X0: array_like or float, optional Initial condition (default = 0). transpose: bool, optional (default=False) @@ -459,20 +462,21 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None, Parameters ---------- - sys: StateSpace, or TransferFunction + sys: StateSpace or TransferFunction LTI system to simulate - T: array-like or number, optional + T: array_like or float, optional Time vector, or simulation time duration if a number. If T is not provided, an attempt is made to create it automatically from the dynamics of sys. If sys is continuous-time, the time increment dt is chosen small enough to show the fastest mode, and the simulation time period tfinal long enough to show the slowest mode, excluding - poles at the origin. If this results in too many time steps (>5000), - dt is reduced. If sys is discrete-time, only tfinal is computed, and - tfinal is reduced if it requires too many simulation steps. + poles at the origin and pole-zero cancellations. If this results in + too many time steps (>5000), dt is reduced. If sys is discrete-time, + only tfinal is computed, and final is reduced if it requires too + many simulation steps. - X0: array-like or number, optional + X0: array_like or float, optional Initial condition (default = 0) Numbers are converted to constant arrays with the correct shape. @@ -484,7 +488,7 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None, Index of the output that will be used in this simulation. Set to None to not trim outputs - T_num: number, optional + T_num: int, optional Number of time steps to use in simulation if T is not provided as an array (autocomputed if not given); ignored if sys is discrete-time. @@ -527,7 +531,7 @@ def step_response(sys, T=None, X0=0., input=None, output=None, T_num=None, """ sys = _get_ss_simo(sys, input, output) if T is None or np.asarray(T).size == 1: - T = _default_time_vector(sys, N=T_num, tfinal=T) + T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=True) U = np.ones_like(T) T, yout, xout = forced_response(sys, T, U, X0, transpose=transpose, @@ -546,21 +550,21 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, Parameters ---------- - sys: StateSpace, or TransferFunction + sys : StateSpace or TransferFunction LTI system to simulate - T: array-like or number, optional + T : array_like or float, optional Time vector, or simulation time duration if a number (time vector is autocomputed if not given, see :func:`step_response` for more detail) - T_num: number, optional + T_num : int, optional Number of time steps to use in simulation if T is not provided as an array (autocomputed if not given); ignored if sys is discrete-time. - SettlingTimeThreshold: float value, optional + SettlingTimeThreshold : float value, optional Defines the error to compute settling time (default = 0.02) - RiseTimeLimits: tuple (lower_threshold, upper_theshold) + RiseTimeLimits : tuple (lower_threshold, upper_theshold) Defines the lower and upper threshold for RiseTime computation Returns @@ -587,7 +591,7 @@ def step_info(sys, T=None, T_num=None, SettlingTimeThreshold=0.02, ''' sys = _get_ss_simo(sys) if T is None or np.asarray(T).size == 1: - T = _default_time_vector(sys, N=T_num, tfinal=T) + T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=True) T, yout = step_response(sys, T) @@ -636,49 +640,49 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None, Parameters ---------- - sys: StateSpace, or TransferFunction + sys : StateSpace or TransferFunction LTI system to simulate - T: array-like or number, optional + T : array_like or float, optional Time vector, or simulation time duration if a number (time vector is autocomputed if not given; see :func:`step_response` for more detail) - X0: array-like or number, optional + X0 : array_like or float, optional Initial condition (default = 0) Numbers are converted to constant arrays with the correct shape. - input: int + input : int Ignored, has no meaning in initial condition calculation. Parameter ensures compatibility with step_response and impulse_response - output: int + output : int Index of the output that will be used in this simulation. Set to None to not trim outputs - T_num: number, optional + T_num : int, optional Number of time steps to use in simulation if T is not provided as an array (autocomputed if not given); ignored if sys is discrete-time. - transpose: bool + transpose : bool If True, transpose all input and output arrays (for backward compatibility with MATLAB and :func:`scipy.signal.lsim`) - return_x: bool + return_x : bool If True, return the state vector (default = False). - squeeze: bool, optional (default=True) + squeeze : bool, optional (default=True) If True, remove single-dimensional entries from the shape of the output. For single output systems, this converts the output response to a 1D array. Returns ------- - T: array + T : array Time values of the output - yout: array + yout : array Response of the system - xout: array + xout : array Individual response of each x variable See Also @@ -699,7 +703,7 @@ def initial_response(sys, T=None, X0=0., input=0, output=None, T_num=None, # Create time and input vectors; checking is done in forced_response(...) # The initial vector X0 is created in forced_response(...) if necessary if T is None or np.asarray(T).size == 1: - T = _default_time_vector(sys, N=T_num, tfinal=T) + T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=False) U = np.zeros_like(T) T, yout, _xout = forced_response(sys, T, U, X0, transpose=transpose, @@ -726,48 +730,48 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=None, Parameters ---------- - sys: StateSpace, TransferFunction + sys : StateSpace, TransferFunction LTI system to simulate - T: array-like or number, optional - Time vector, or simulation time duration if a number (time vector is + T : array_like or float, optional + Time vector, or simulation time duration if a scalar (time vector is autocomputed if not given; see :func:`step_response` for more detail) - X0: array-like or number, optional + X0 : array_like or float, optional Initial condition (default = 0) Numbers are converted to constant arrays with the correct shape. - input: int + input : int Index of the input that will be used in this simulation. - output: int + output : int Index of the output that will be used in this simulation. Set to None to not trim outputs - T_num: number, optional + T_num : int, optional Number of time steps to use in simulation if T is not provided as an array (autocomputed if not given); ignored if sys is discrete-time. - transpose: bool + transpose : bool If True, transpose all input and output arrays (for backward compatibility with MATLAB and :func:`scipy.signal.lsim`) - return_x: bool + return_x : bool If True, return the state vector (default = False). - squeeze: bool, optional (default=True) + squeeze : bool, optional (default=True) If True, remove single-dimensional entries from the shape of the output. For single output systems, this converts the output response to a 1D array. Returns ------- - T: array + T : array Time values of the output - yout: array + yout : array Response of the system - xout: array + xout : array Individual response of each x variable See Also @@ -803,7 +807,7 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=None, # Compute T and U, no checks necessary, will be checked in forced_response if T is None or np.asarray(T).size == 1: - T = _default_time_vector(sys, N=T_num, tfinal=T) + T = _default_time_vector(sys, N=T_num, tfinal=T, is_step=False) U = np.zeros_like(T) # Compute new X0 that contains the impulse @@ -827,10 +831,10 @@ def impulse_response(sys, T=None, X0=0., input=0, output=None, T_num=None, # utility function to find time period and time increment using pole locations def _ideal_tfinal_and_dt(sys, is_step=True): - """helper function to compute ideal simulation duration tfinal and dt, the - time increment. Usually called by _default_time_vector, whose job it is to - bring reality into the picture. - + """helper function to compute ideal simulation duration tfinal and dt, the + time increment. Usually called by _default_time_vector, whose job it is to + choose a realistic time vector. Considers both poles and zeros. + For discrete-time models, dt is inherent and only tfinal is computed. Parameters @@ -841,22 +845,23 @@ def _ideal_tfinal_and_dt(sys, is_step=True): Scales the dc value by the magnitude of the nonzero mode since integrating the impulse response gives ∫exp(-λt) = -exp(-λt)/λ. Default is True. + Returns ------- tfinal : float The final time instance for which the simulation will be performed. dt : float The estimated sampling period for the simulation. - + Notes - ----- + ----- Just by evaluating the fastest mode for dt and slowest for tfinal often leads to unnecessary, bloated sampling (e.g., Transfer(1,[1,1001,1000])) since dt will be very small and tfinal will be too large though the fast mode hardly ever contributes. Similarly, change the numerator to [1, 2, 0] and the simulation would be unnecessarily long and the plot is virtually an L shape since the decay is so fast. - + Instead, a modal decomposition in time domain hence a truncated ZIR and ZSR can be used such that only the modes that have significant effect on the time response are taken. But the sensitivity of the eigenvalues complicate @@ -865,65 +870,142 @@ def _ideal_tfinal_and_dt(sys, is_step=True): for simple eigenvalue sensitivity about the nonunity of . The size of the response is dependent on the size of the eigenshapes rather than the eigenvalues themselves. - - By Ilhan Polat, with modifications by Sawyer Fuller to integrate into + + By Ilhan Polat, with modifications by Sawyer Fuller to integrate into python-control 2020.08.17 """ sqrt_eps = np.sqrt(np.spacing(1.)) default_tfinal = 5 # Default simulation horizon + default_dt = 0.1 total_cycles = 5 # number of cycles for oscillating modes pts_per_cycle = 25 # Number of points divide a period of oscillation log_decay_percent = np.log(100) # Factor of reduction for real pole decays - if sys.is_static_gain(): pass - constant = 7.0 - tolerance = 1e-10 - A = ssdata(sys)[0] - if A.shape == (0,0): - # no dynamics - tfinal = constant * 1.0 - dt = sys.dt if isdtime(sys, strict=True) else 1.0 - else: - poles = sp.linalg.eigvals(A) - # calculate ideal dt - if isdtime(sys, strict=True): - # z-poles to s-plane using s=(lnz)/dt, no ln(0) - poles = np.log(poles[abs(poles) > 0])/sys.dt - dt = sys.dt - else: - fastest_natural_frequency = max(abs(poles)) - dt = 1/constant / fastest_natural_frequency - - # calculate ideal tfinal - poles = poles[abs(poles.real) > tolerance] # ignore poles near im axis - if poles.size == 0: - slowest_decay_rate = 1.0 - else: - slowest_decay_rate = min(abs(poles.real)) - tfinal = constant / slowest_decay_rate + if sys.is_static_gain(): + tfinal = default_tfinal + dt = sys.dt if isdtime(sys, strict=True) else default_dt + elif isdtime(sys, strict=True): + dt = sys.dt + A = _convertToStateSpace(sys).A + tfinal = default_tfinal + p = eigvals(A) + # Array Masks + # unstable + m_u = (np.abs(p) >= 1 + sqrt_eps) + p_u, p = p[m_u], p[~m_u] + if p_u.size > 0: + m_u = (p_u.real < 0) & (np.abs(p_u.imag) < sqrt_eps) + t_emp = np.max(log_decay_percent / np.abs(np.log(p_u[~m_u])/dt)) + tfinal = max(tfinal, t_emp) + + # zero - negligible effect on tfinal + m_z = np.abs(p) < sqrt_eps + p = p[~m_z] + # Negative reals- treated as oscillary mode + m_nr = (p.real < 0) & (np.abs(p.imag) < sqrt_eps) + p_nr, p = p[m_nr], p[~m_nr] + if p_nr.size > 0: + t_emp = np.max(log_decay_percent / np.abs((np.log(p_nr)/dt).real)) + tfinal = max(tfinal, t_emp) + # discrete integrators + m_int = (p.real - 1 < sqrt_eps) & (np.abs(p.imag) < sqrt_eps) + p_int, p = p[m_int], p[~m_int] + # pure oscillatory modes + m_w = (np.abs(np.abs(p) - 1) < sqrt_eps) + p_w, p = p[m_w], p[~m_w] + if p_w.size > 0: + t_emp = total_cycles * 2 * np.pi / np.abs(np.log(p_w)/dt).min() + tfinal = max(tfinal, t_emp) + + if p.size > 0: + t_emp = log_decay_percent / np.abs((np.log(p)/dt).real).min() + tfinal = max(tfinal, t_emp) + + if p_int.size > 0: + tfinal = tfinal * 5 + else: # cont time + sys_ss = _convertToStateSpace(sys) + # Improve conditioning via balancing and zeroing tiny entries + # See for [[1,2,0], [9,1,0.01], [1,2,10*np.pi]] before/after balance + b, (sca, perm) = matrix_balance(sys_ss.A, separate=True) + p, l, r = eig(b, left=True, right=True) + # Reciprocal of inner product for each λ, (bound the ~infs by 1e12) + # G = Transfer([1], [1,0,1]) gives zero sensitivity (bound by 1e-12) + eig_sens = np.reciprocal(maximum(1e-12, einsum('ij,ij->j', l, r).real)) + eig_sens = minimum(1e12, eig_sens) + # Tolerances + p[np.abs(p) < np.spacing(eig_sens * norm(b, 1))] = 0. + # Incorporate balancing to outer factors + l[perm, :] *= np.reciprocal(sca)[:, None] + r[perm, :] *= sca[:, None] + w, v = sys_ss.C @ r, l.T.conj() @ sys_ss.B + + origin = False + # Computing the "size" of the response of each simple mode + wn = np.abs(p) + if np.any(wn == 0.): + origin = True + + dc = np.zeros_like(p, dtype=float) + # well-conditioned nonzero poles, np.abs just in case + ok = np.abs(eig_sens) <= 1/sqrt_eps + # the averaged t→∞ response of each simple λ on each i/o channel + # See, A = [[-1, k], [0, -2]], response sizes are k-dependent (that is + # R/L eigenvector dependent) + dc[ok] = norm(v[ok, :], axis=1)*norm(w[:, ok], axis=0)*eig_sens[ok] + dc[wn != 0.] /= wn[wn != 0] if is_step else 1. + dc[wn == 0.] = 0. + # double the oscillating mode magnitude for the conjugate + dc[p.imag != 0.] *= 2 + + # Now get rid of noncontributing integrators and simple modes if any + relevance = (dc > 0.1*dc.max()) | ~ok + psub = p[relevance] + wnsub = wn[relevance] + + tfinal, dt = [], [] + ints = wnsub == 0. + iw = (psub.imag != 0.) & (np.abs(psub.real) <= sqrt_eps) + + # Pure imaginary? + if np.any(iw): + tfinal += (total_cycles * 2 * np.pi / wnsub[iw]).tolist() + dt += (2 * np.pi / pts_per_cycle / wnsub[iw]).tolist() + # The rest ~ts = log(%ss value) / exp(Re(λ)t) + texp_mode = log_decay_percent / np.abs(psub[~iw & ~ints].real) + tfinal += texp_mode.tolist() + dt += minimum(texp_mode / 50, + (2 * np.pi / pts_per_cycle / wnsub[~iw & ~ints])).tolist() + + # All integrators? + if len(tfinal) == 0: + return default_tfinal*5, default_dt*5 + + tfinal = np.max(tfinal)*(5 if origin else 1) + dt = np.min(dt) return tfinal, dt -# test below: ct with pole at the origin is 7 seconds, ct with pole at .5 is 14 s long, -def _default_time_vector(sys, N=None, tfinal=None): - """Returns a time vector suitable for observing the response of the - both the slowest poles and fastest resonant modes. if system is - discrete-time, N is ignored """ +def _default_time_vector(sys, N=None, tfinal=None, is_step=True): + """Returns a time vector that has a reasonable number of points. + if system is discrete-time, N is ignored """ N_max = 5000 - N_min_ct = 100 - N_min_dt = 7 # more common to see just a few samples in discrete-time + N_min_ct = 100 # min points for cont time systems + N_min_dt = 20 # more common to see just a few samples in discrete-time - ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(sys) + ideal_tfinal, ideal_dt = _ideal_tfinal_and_dt(sys, is_step=is_step) if isdtime(sys, strict=True): + # only need to use default_tfinal if not given; N is ignored. if tfinal is None: # for discrete time, change from ideal_tfinal if N too large/small N = int(np.clip(ideal_tfinal/sys.dt, N_min_dt, N_max))# [N_min, N_max] tfinal = sys.dt * N else: N = int(tfinal/sys.dt) + tfinal = N * sys.dt # make tfinal an integer multiple of sys.dt else: if tfinal is None: # for continuous time, simulate to ideal_tfinal but limit N diff --git a/control/xferfcn.py b/control/xferfcn.py index b00edc7d8..51e1e732f 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -1091,7 +1091,17 @@ def _dcgain_cont(self): gain[i][j] = np.nan return np.squeeze(gain) - + def is_static_gain(self): + """returns True if and only if all of the numerator and denominator + polynomials of the (possibly MIMO) transfer funnction are zeroth order, + that is, if the system has no dynamics. """ + for list_of_polys in self.num, self.den: + for row in list_of_polys: + for poly in row: + if len(poly) > 1: + return False + return True + # c2d function contributed by Benjamin White, Oct 2012 def _c2d_matched(sysC, Ts): # Pole-zero match method of continuous to discrete time conversion From a6cac86a4252be4e22800c7c9c00c798918821b0 Mon Sep 17 00:00:00 2001 From: sawyerbfuller <58706249+sawyerbfuller@users.noreply.github.com> Date: Wed, 19 Aug 2020 11:11:59 -0700 Subject: [PATCH 3/7] Update control/timeresp.py removed unicode in code to make python 2 happy Co-authored-by: Ben Greiner --- control/timeresp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/control/timeresp.py b/control/timeresp.py index 42e4cf125..b677efabb 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -843,7 +843,7 @@ def _ideal_tfinal_and_dt(sys, is_step=True): The system whose time response is to be computed is_step : bool Scales the dc value by the magnitude of the nonzero mode since - integrating the impulse response gives ∫exp(-λt) = -exp(-λt)/λ. + integrating the impulse response gives :math:`\int e^{-\lambda t} = -e^{-\lambda t}/ \lambda` Default is True. Returns From 65ee5ed6130f9a279ef5d76e6eabbf3544166ba1 Mon Sep 17 00:00:00 2001 From: Sawyer Fuller Date: Wed, 19 Aug 2020 11:34:24 -0700 Subject: [PATCH 4/7] removed one more unicode symbol --- control/timeresp.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/control/timeresp.py b/control/timeresp.py index b677efabb..de282cd2d 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -843,7 +843,8 @@ def _ideal_tfinal_and_dt(sys, is_step=True): The system whose time response is to be computed is_step : bool Scales the dc value by the magnitude of the nonzero mode since - integrating the impulse response gives :math:`\int e^{-\lambda t} = -e^{-\lambda t}/ \lambda` + integrating the impulse response gives + :math:`\int e^{-\lambda t} = -e^{-\lambda t}/ \lambda` Default is True. Returns @@ -865,7 +866,7 @@ def _ideal_tfinal_and_dt(sys, is_step=True): Instead, a modal decomposition in time domain hence a truncated ZIR and ZSR can be used such that only the modes that have significant effect on the time response are taken. But the sensitivity of the eigenvalues complicate - the matter since dλ = with = 1. Hence we can only work + the matter since dlambda = with = 1. Hence we can only work with simple poles with this formulation. See Golub, Van Loan Section 7.2.2 for simple eigenvalue sensitivity about the nonunity of . The size of the response is dependent on the size of the eigenshapes rather than the From 3d25b8b06d2e42da3b78da67f139ead92c8d1ab9 Mon Sep 17 00:00:00 2001 From: Sawyer Fuller Date: Wed, 19 Aug 2020 12:59:29 -0700 Subject: [PATCH 5/7] more unicode expulsion --- control/timeresp.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/control/timeresp.py b/control/timeresp.py index de282cd2d..5df51d0dc 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -931,7 +931,7 @@ def _ideal_tfinal_and_dt(sys, is_step=True): # See for [[1,2,0], [9,1,0.01], [1,2,10*np.pi]] before/after balance b, (sca, perm) = matrix_balance(sys_ss.A, separate=True) p, l, r = eig(b, left=True, right=True) - # Reciprocal of inner product for each λ, (bound the ~infs by 1e12) + # Reciprocal of inner product for each eigval, (bound the ~infs by 1e12) # G = Transfer([1], [1,0,1]) gives zero sensitivity (bound by 1e-12) eig_sens = np.reciprocal(maximum(1e-12, einsum('ij,ij->j', l, r).real)) eig_sens = minimum(1e12, eig_sens) @@ -951,7 +951,7 @@ def _ideal_tfinal_and_dt(sys, is_step=True): dc = np.zeros_like(p, dtype=float) # well-conditioned nonzero poles, np.abs just in case ok = np.abs(eig_sens) <= 1/sqrt_eps - # the averaged t→∞ response of each simple λ on each i/o channel + # the averaged t->inf response of each simple eigval on each i/o channel # See, A = [[-1, k], [0, -2]], response sizes are k-dependent (that is # R/L eigenvector dependent) dc[ok] = norm(v[ok, :], axis=1)*norm(w[:, ok], axis=0)*eig_sens[ok] @@ -973,7 +973,7 @@ def _ideal_tfinal_and_dt(sys, is_step=True): if np.any(iw): tfinal += (total_cycles * 2 * np.pi / wnsub[iw]).tolist() dt += (2 * np.pi / pts_per_cycle / wnsub[iw]).tolist() - # The rest ~ts = log(%ss value) / exp(Re(λ)t) + # The rest ~ts = log(%ss value) / exp(Re(eigval)t) texp_mode = log_decay_percent / np.abs(psub[~iw & ~ints].real) tfinal += texp_mode.tolist() dt += minimum(texp_mode / 50, From 23d90d652113c15fd8e28b6ade1b37f0b3cd1e8f Mon Sep 17 00:00:00 2001 From: Sawyer Fuller Date: Wed, 19 Aug 2020 13:19:49 -0700 Subject: [PATCH 6/7] curse you python 2. removed '@' sign --- control/timeresp.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/control/timeresp.py b/control/timeresp.py index 5df51d0dc..07a37c65d 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -65,6 +65,9 @@ capability and better automatic time vector creation Date: June 2020 +Modified by Ilhan Polat to improve automatic time vector creation +Date: August 17, 2020 + $Id$ """ @@ -940,7 +943,7 @@ def _ideal_tfinal_and_dt(sys, is_step=True): # Incorporate balancing to outer factors l[perm, :] *= np.reciprocal(sca)[:, None] r[perm, :] *= sca[:, None] - w, v = sys_ss.C @ r, l.T.conj() @ sys_ss.B + w, v = sys_ss.C.dot(r), l.T.conj().dot(sys_ss.B) origin = False # Computing the "size" of the response of each simple mode From c588cdcc674a8b7bb42db22a3b4939120d22e272 Mon Sep 17 00:00:00 2001 From: sawyerbfuller <58706249+sawyerbfuller@users.noreply.github.com> Date: Mon, 24 Aug 2020 14:23:43 -0700 Subject: [PATCH 7/7] Update control/xferfcn.py Spelling error in docstring Co-authored-by: Naman Gera --- control/xferfcn.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/control/xferfcn.py b/control/xferfcn.py index 51e1e732f..1cba50bd7 100644 --- a/control/xferfcn.py +++ b/control/xferfcn.py @@ -1093,7 +1093,7 @@ def _dcgain_cont(self): def is_static_gain(self): """returns True if and only if all of the numerator and denominator - polynomials of the (possibly MIMO) transfer funnction are zeroth order, + polynomials of the (possibly MIMO) transfer function are zeroth order, that is, if the system has no dynamics. """ for list_of_polys in self.num, self.den: for row in list_of_polys: