diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py index 70607419e..fd474f9d0 100644 --- a/control/tests/modelsimp_test.py +++ b/control/tests/modelsimp_test.py @@ -65,8 +65,6 @@ def testMarkovSignature(self, matarrayout, matarrayin): markov(Y, U, m) # Make sure markov() returns the right answer - # forced response can return wrong shape until gh-488 is merged - @pytest.mark.xfail @pytest.mark.parametrize("k, m, n", [(2, 2, 2), (2, 5, 5), @@ -81,7 +79,7 @@ def testMarkovResults(self, k, m, n): # m = number of Markov parameters # n = size of the data vector # - # Values should match exactly for n = m, otherewise you get a + # Values *should* match exactly for n = m, otherewise you get a # close match but errors due to the assumption that C A^k B = # 0 for k > m-2 (see modelsimp.py). # @@ -108,7 +106,10 @@ def testMarkovResults(self, k, m, n): Mcomp = markov(Y, U, m) # Compare to results from markov() - np.testing.assert_array_almost_equal(Mtrue, Mcomp) + # experimentally determined probability to get non matching results + # with rtot=1e-6 and atol=1e-8 due to numerical errors + # for k=5, m=n=10: 0.015 % + np.testing.assert_allclose(Mtrue, Mcomp, rtol=1e-6, atol=1e-8) def testModredMatchDC(self, matarrayin): #balanced realization computed in matlab for the transfer function: @@ -219,4 +220,3 @@ def testBalredMatchDC(self, matarrayin): np.testing.assert_array_almost_equal(rsys.B, Brtrue, decimal=4) np.testing.assert_array_almost_equal(rsys.C, Crtrue, decimal=4) np.testing.assert_array_almost_equal(rsys.D, Drtrue, decimal=4) - diff --git a/control/tests/timeresp_test.py b/control/tests/timeresp_test.py index 1c97b9385..a91507a83 100644 --- a/control/tests/timeresp_test.py +++ b/control/tests/timeresp_test.py @@ -1,12 +1,4 @@ -"""timeresp_test.py - test time response functions - -RMM, 17 Jun 2011 (based on TestMatlab from v0.4c) - -This test suite just goes through and calls all of the MATLAB -functions using different systems and arguments to make sure that -nothing crashes. It doesn't test actual functionality; the module -specific unit tests will do that. -""" +"""timeresp_test.py - test time response functions""" from copy import copy from distutils.version import StrictVersion @@ -38,50 +30,46 @@ def __repr__(self): class TestTimeresp: @pytest.fixture - def siso_ss1(self): + def tsystem(self, request): + """Define some test systems""" + """continuous""" A = np.array([[1., -2.], [3., -4.]]) B = np.array([[5.], [7.]]) C = np.array([[6., 8.]]) D = np.array([[9.]]) - T = TSys(StateSpace(A, B, C, D, 0)) - - T.t = np.linspace(0, 1, 10) - T.ystep = np.array([9., 17.6457, 24.7072, 30.4855, 35.2234, - 39.1165, 42.3227, 44.9694, 47.1599, 48.9776]) - - T.yinitial = np.array([11., 8.1494, 5.9361, 4.2258, 2.9118, - 1.9092, 1.1508, 0.5833, 0.1645, -0.1391]) - - return T - - @pytest.fixture - def siso_ss2(self, siso_ss1): - """System siso_ss2 with D=0""" + siso_ss1 = TSys(StateSpace(A, B, C, D, 0)) + siso_ss1.t = np.linspace(0, 1, 10) + siso_ss1.ystep = np.array([9., 17.6457, 24.7072, 30.4855, 35.2234, + 39.1165, 42.3227, 44.9694, 47.1599, + 48.9776]) + siso_ss1.X0 = np.array([[.5], [1.]]) + siso_ss1.yinitial = np.array([11., 8.1494, 5.9361, 4.2258, 2.9118, + 1.9092, 1.1508, 0.5833, 0.1645, -0.1391]) ss1 = siso_ss1.sys - T = TSys(StateSpace(ss1.A, ss1.B, ss1.C, 0, 0)) - T.t = siso_ss1.t - T.ystep = siso_ss1.ystep - 9 - T.initial = siso_ss1.yinitial - 9 - T.yimpulse = np.array([86., 70.1808, 57.3753, 46.9975, 38.5766, - 31.7344, 26.1668, 21.6292, 17.9245, 14.8945]) - return T + """D=0, continuous""" + siso_ss2 = TSys(StateSpace(ss1.A, ss1.B, ss1.C, 0, 0)) + siso_ss2.t = siso_ss1.t + siso_ss2.ystep = siso_ss1.ystep - 9 + siso_ss2.initial = siso_ss1.yinitial - 9 + siso_ss2.yimpulse = np.array([86., 70.1808, 57.3753, 46.9975, 38.5766, + 31.7344, 26.1668, 21.6292, 17.9245, + 14.8945]) - @pytest.fixture - def siso_tf1(self): - # Create some transfer functions - return TSys(TransferFunction([1], [1, 2, 1], 0)) + """System with unspecified timebase""" + siso_ss2_dtnone = TSys(StateSpace(ss1.A, ss1.B, ss1.C, 0, None)) + siso_ss2_dtnone.t = np.arange(0, 10, 1.) + siso_ss2_dtnone.ystep = np.array([0., 86., -72., 230., -360., 806., + -1512., 3110., -6120., 12326.]) - @pytest.fixture - def siso_tf2(self, siso_ss1): - T = copy(siso_ss1) - T.sys = ss2tf(siso_ss1.sys) - return T + siso_tf1 = TSys(TransferFunction([1], [1, 2, 1], 0)) - @pytest.fixture - def mimo_ss1(self, siso_ss1): - # Create MIMO system, contains ``siso_ss1`` twice + siso_tf2 = copy(siso_ss1) + siso_tf2.sys = ss2tf(siso_ss1.sys) + + """MIMO system, contains ``siso_ss1`` twice""" + mimo_ss1 = copy(siso_ss1) A = np.zeros((4, 4)) A[:2, :2] = siso_ss1.sys.A A[2:, 2:] = siso_ss1.sys.A @@ -94,13 +82,10 @@ def mimo_ss1(self, siso_ss1): D = np.zeros((2, 2)) D[:1, :1] = siso_ss1.sys.D D[1:, 1:] = siso_ss1.sys.D - T = copy(siso_ss1) - T.sys = StateSpace(A, B, C, D) - return T + mimo_ss1.sys = StateSpace(A, B, C, D) - @pytest.fixture - def mimo_ss2(self, siso_ss2): - # Create MIMO system, contains ``siso_ss2`` twice + """MIMO system, contains ``siso_ss2`` twice""" + mimo_ss2 = copy(siso_ss2) A = np.zeros((4, 4)) A[:2, :2] = siso_ss2.sys.A A[2:, 2:] = siso_ss2.sys.A @@ -111,93 +96,84 @@ def mimo_ss2(self, siso_ss2): C[:1, :2] = siso_ss2.sys.C C[1:, 2:] = siso_ss2.sys.C D = np.zeros((2, 2)) - T = copy(siso_ss2) - T.sys = StateSpace(A, B, C, D, 0) - return T - - # Create discrete time systems - - @pytest.fixture - def siso_dtf0(self): - T = TSys(TransferFunction([1.], [1., 0.], 1.)) - T.t = np.arange(4) - T.yimpulse = [0., 1., 0., 0.] - return T - - @pytest.fixture - def siso_dtf1(self): - T = TSys(TransferFunction([1], [1, 1, 0.25], True)) - T.t = np.arange(0, 5, 1) - return T - - @pytest.fixture - def siso_dtf2(self): - T = TSys(TransferFunction([1], [1, 1, 0.25], 0.2)) - T.t = np.arange(0, 5, 0.2) - return T - - @pytest.fixture - def siso_dss1(self, siso_dtf1): - T = copy(siso_dtf1) - T.sys = tf2ss(siso_dtf1.sys) - return T - - @pytest.fixture - def siso_dss2(self, siso_dtf2): - T = copy(siso_dtf2) - T.sys = tf2ss(siso_dtf2.sys) - return T - - @pytest.fixture - def mimo_dss1(self, mimo_ss1): - ss1 = mimo_ss1.sys - T = TSys( - StateSpace(ss1.A, ss1.B, ss1.C, ss1.D, True)) - T.t = np.arange(0, 5, 0.2) - return T - - @pytest.fixture - def mimo_dss2(self, mimo_ss1): - T = copy(mimo_ss1) - T.sys = c2d(mimo_ss1.sys, T.t[1]-T.t[0]) - return T - - @pytest.fixture - def mimo_tf2(self, siso_ss2, mimo_ss2): - T = copy(mimo_ss2) - # construct from siso to avoid slycot during fixture setup + mimo_ss2.sys = StateSpace(A, B, C, D, 0) + + """discrete""" + siso_dtf0 = TSys(TransferFunction([1.], [1., 0.], 1.)) + siso_dtf0.t = np.arange(4) + siso_dtf0.yimpulse = [0., 1., 0., 0.] + + siso_dtf1 = TSys(TransferFunction([1], [1, 1, 0.25], True)) + siso_dtf1.t = np.arange(0, 5, 1) + siso_dtf1.ystep = np.array([0. , 0. , 1. , 0. , 0.75]) + + siso_dtf2 = TSys(TransferFunction([1], [1, 1, 0.25], 0.2)) + siso_dtf2.t = np.arange(0, 5, 0.2) + siso_dtf2.ystep = np.array( + [0. , 0. , 1. , 0. , 0.75 , 0.25 , + 0.5625, 0.375 , 0.4844, 0.4219, 0.457 , 0.4375, + 0.4482, 0.4424, 0.4456, 0.4438, 0.4448, 0.4443, + 0.4445, 0.4444, 0.4445, 0.4444, 0.4445, 0.4444, + 0.4444]) + + """Time step which leads to rounding errors for time vector length""" + num = [-0.10966442, 0.12431949] + den = [1., -1.86789511, 0.88255018] + dt = 0.12493963338370018 + siso_dtf3 = TSys(TransferFunction(num, den, dt)) + siso_dtf3.t = np.linspace(0, 9*dt, 10) + siso_dtf3.ystep = np.array( + [ 0. , -0.1097, -0.1902, -0.2438, -0.2729, + -0.2799, -0.2674, -0.2377, -0.1934, -0.1368]) + + """dtf1 converted statically, because Slycot and Scipy produce + different realizations, wich means different initial condtions,""" + siso_dss1 = copy(siso_dtf1) + siso_dss1.sys = StateSpace([[-1., -0.25], + [ 1., 0.]], + [[1.], + [0.]], + [[0., 1.]], + [[0.]], + True) + siso_dss1.X0 = [0.5, 1.] + siso_dss1.yinitial = np.array([1., 0.5, -0.75, 0.625, -0.4375]) + + siso_dss2 = copy(siso_dtf2) + siso_dss2.sys = tf2ss(siso_dtf2.sys) + + mimo_dss1 = TSys(StateSpace(ss1.A, ss1.B, ss1.C, ss1.D, True)) + mimo_dss1.t = np.arange(0, 5, 0.2) + + mimo_dss2 = copy(mimo_ss1) + mimo_dss2.sys = c2d(mimo_ss1.sys, mimo_ss1.t[1]-mimo_ss1.t[0]) + + mimo_tf2 = copy(mimo_ss2) tf_ = ss2tf(siso_ss2.sys) - T.sys = TransferFunction([[tf_.num[0][0], [0]], [[0], tf_.num[0][0]]], - [[tf_.den[0][0], [1]], [[1], tf_.den[0][0]]], - 0) - return T + mimo_tf2.sys = TransferFunction( + [[tf_.num[0][0], [0]], [[0], tf_.num[0][0]]], + [[tf_.den[0][0], [1]], [[1], tf_.den[0][0]]], + 0) - @pytest.fixture - def mimo_dtf1(self, siso_dtf1): - T = copy(siso_dtf1) - # construct from siso to avoid slycot during fixture setup + mimo_dtf1 = copy(siso_dtf1) tf_ = siso_dtf1.sys - T.sys = TransferFunction([[tf_.num[0][0], [0]], [[0], tf_.num[0][0]]], - [[tf_.den[0][0], [1]], [[1], tf_.den[0][0]]], - True) - return T + mimo_dtf1.sys = TransferFunction( + [[tf_.num[0][0], [0]], [[0], tf_.num[0][0]]], + [[tf_.den[0][0], [1]], [[1], tf_.den[0][0]]], + True) - @pytest.fixture - def pole_cancellation(self): # for pole cancellation tests - return TransferFunction([1.067e+05, 5.791e+04], - [10.67, 1.067e+05, 5.791e+04]) + pole_cancellation = TSys(TransferFunction( + [1.067e+05, 5.791e+04], + [10.67, 1.067e+05, 5.791e+04])) - @pytest.fixture - def no_pole_cancellation(self): - return TransferFunction([1.881e+06], - [188.1, 1.881e+06]) + no_pole_cancellation = TSys(TransferFunction( + [1.881e+06], + [188.1, 1.881e+06])) - @pytest.fixture - def siso_tf_type1(self): # System Type 1 - Step response not stationary: G(s)=1/s(s+1) - T = TSys(TransferFunction(1, [1, 1, 0])) - T.step_info = { + siso_tf_type1 = TSys(TransferFunction(1, [1, 1, 0])) + siso_tf_type1.step_info = { 'RiseTime': np.NaN, 'SettlingTime': np.NaN, 'SettlingMin': np.NaN, @@ -207,14 +183,11 @@ def siso_tf_type1(self): 'Peak': np.Inf, 'PeakTime': np.Inf, 'SteadyStateValue': np.NaN} - return T - @pytest.fixture - def siso_tf_kpos(self): # SISO under shoot response and positive final value # G(s)=(-s+1)/(s²+s+1) - T = TSys(TransferFunction([-1, 1], [1, 1, 1])) - T.step_info = { + siso_tf_kpos = TSys(TransferFunction([-1, 1], [1, 1, 1])) + siso_tf_kpos.step_info = { 'RiseTime': 1.242, 'SettlingTime': 9.110, 'SettlingMin': 0.90, @@ -224,14 +197,11 @@ def siso_tf_kpos(self): 'Peak': 1.208, 'PeakTime': 4.282, 'SteadyStateValue': 1.0} - return T - @pytest.fixture - def siso_tf_kneg(self): # SISO under shoot response and negative final value # k=-1 G(s)=-(-s+1)/(s²+s+1) - T = TSys(TransferFunction([1, -1], [1, 1, 1])) - T.step_info = { + siso_tf_kneg = TSys(TransferFunction([1, -1], [1, 1, 1])) + siso_tf_kneg.step_info = { 'RiseTime': 1.242, 'SettlingTime': 9.110, 'SettlingMin': -1.208, @@ -241,13 +211,9 @@ def siso_tf_kneg(self): 'Peak': 1.208, 'PeakTime': 4.282, 'SteadyStateValue': -1.0} - return T - @pytest.fixture - def siso_tf_asymptotic_from_neg1(self): - # Peak_value = Undershoot = y_final(y(t=inf)) - T = TSys(TransferFunction([-1, 1], [1, 1])) - T.step_info = { + siso_tf_asymptotic_from_neg1 = TSys(TransferFunction([-1, 1], [1, 1])) + siso_tf_asymptotic_from_neg1.step_info = { 'RiseTime': 2.197, 'SettlingTime': 4.605, 'SettlingMin': 0.9, @@ -257,15 +223,14 @@ def siso_tf_asymptotic_from_neg1(self): 'Peak': 1.0, 'PeakTime': 0.0, 'SteadyStateValue': 1.0} - T.kwargs = {'step_info': {'T': np.arange(0, 5, 1e-3)}} - return T + siso_tf_asymptotic_from_neg1.kwargs = { + 'step_info': {'T': np.arange(0, 5, 1e-3)}} - @pytest.fixture - def siso_tf_step_matlab(self): # example from matlab online help # https://www.mathworks.com/help/control/ref/stepinfo.html - T = TSys(TransferFunction([1, 5, 5], [1, 1.65, 5, 6.5, 2])) - T.step_info = { + siso_tf_step_matlab = TSys(TransferFunction([1, 5, 5], + [1, 1.65, 5, 6.5, 2])) + siso_tf_step_matlab.step_info = { 'RiseTime': 3.8456, 'SettlingTime': 27.9762, 'SettlingMin': 2.0689, @@ -275,10 +240,7 @@ def siso_tf_step_matlab(self): 'Peak': 2.6873, 'PeakTime': 8.0530, 'SteadyStateValue': 2.5} - return T - @pytest.fixture - def mimo_ss_step_matlab(self): A = [[0.68, -0.34], [0.34, 0.68]] B = [[0.18, -0.05], @@ -287,113 +249,70 @@ def mimo_ss_step_matlab(self): [-1.12, -1.10]] D = [[0, 0], [0.06, -0.37]] - T = TSys(StateSpace(A, B, C, D, 0.2)) - T.kwargs['step_info'] = {'T': 4.6} - T.step_info = [[{'RiseTime': 0.6000, - 'SettlingTime': 3.0000, - 'SettlingMin': -0.5999, - 'SettlingMax': -0.4689, - 'Overshoot': 15.5072, - 'Undershoot': 0., - 'Peak': 0.5999, - 'PeakTime': 1.4000, - 'SteadyStateValue': -0.5193}, - {'RiseTime': 0., - 'SettlingTime': 3.6000, - 'SettlingMin': -0.2797, - 'SettlingMax': -0.1043, - 'Overshoot': 118.9918, - 'Undershoot': 0, - 'Peak': 0.2797, - 'PeakTime': .6000, - 'SteadyStateValue': -0.1277}], - [{'RiseTime': 0.4000, - 'SettlingTime': 2.8000, - 'SettlingMin': -0.6724, - 'SettlingMax': -0.5188, - 'Overshoot': 24.6476, - 'Undershoot': 11.1224, - 'Peak': 0.6724, - 'PeakTime': 1, - 'SteadyStateValue': -0.5394}, - {'RiseTime': 0.0000, # (*) - 'SettlingTime': 3.4000, - 'SettlingMin': -0.4350, # (*) - 'SettlingMax': -0.1485, - 'Overshoot': 132.0170, - 'Undershoot': 0., - 'Peak': 0.4350, - 'PeakTime': .2, - 'SteadyStateValue': -0.1875}]] - # (*): MATLAB gives 0.4 for RiseTime and -0.1034 for - # SettlingMin, but it is unclear what 10% and 90% of - # the steady state response mean, when the step for - # this channel does not start a 0. - return T - - @pytest.fixture - def siso_ss_step_matlab(self, mimo_ss_step_matlab): - T = copy(mimo_ss_step_matlab) - T.sys = T.sys[1, 0] - T.step_info = T.step_info[1][0] - return T + mimo_ss_step_matlab = TSys(StateSpace(A, B, C, D, 0.2)) + mimo_ss_step_matlab.kwargs['step_info'] = {'T': 4.6} + mimo_ss_step_matlab.step_info = [[ + {'RiseTime': 0.6000, + 'SettlingTime': 3.0000, + 'SettlingMin': -0.5999, + 'SettlingMax': -0.4689, + 'Overshoot': 15.5072, + 'Undershoot': 0., + 'Peak': 0.5999, + 'PeakTime': 1.4000, + 'SteadyStateValue': -0.5193}, + {'RiseTime': 0., + 'SettlingTime': 3.6000, + 'SettlingMin': -0.2797, + 'SettlingMax': -0.1043, + 'Overshoot': 118.9918, + 'Undershoot': 0, + 'Peak': 0.2797, + 'PeakTime': .6000, + 'SteadyStateValue': -0.1277}], + [{'RiseTime': 0.4000, + 'SettlingTime': 2.8000, + 'SettlingMin': -0.6724, + 'SettlingMax': -0.5188, + 'Overshoot': 24.6476, + 'Undershoot': 11.1224, + 'Peak': 0.6724, + 'PeakTime': 1, + 'SteadyStateValue': -0.5394}, + {'RiseTime': 0.0000, # (*) + 'SettlingTime': 3.4000, + 'SettlingMin': -0.4350, # (*) + 'SettlingMax': -0.1485, + 'Overshoot': 132.0170, + 'Undershoot': 0., + 'Peak': 0.4350, + 'PeakTime': .2, + 'SteadyStateValue': -0.1875}]] + # (*): MATLAB gives 0.4 for RiseTime and -0.1034 for + # SettlingMin, but it is unclear what 10% and 90% of + # the steady state response mean, when the step for + # this channel does not start a 0. + + siso_ss_step_matlab = copy(mimo_ss_step_matlab) + siso_ss_step_matlab.sys = siso_ss_step_matlab.sys[1, 0] + siso_ss_step_matlab.step_info = siso_ss_step_matlab.step_info[1][0] - @pytest.fixture - def mimo_tf_step_info(self, - siso_tf_kpos, siso_tf_kneg, - siso_tf_step_matlab): Ta = [[siso_tf_kpos, siso_tf_kneg, siso_tf_step_matlab], - [siso_tf_step_matlab, siso_tf_kpos, siso_tf_kneg]] - T = TSys(TransferFunction( + [siso_tf_step_matlab, siso_tf_kpos, siso_tf_kneg]] + mimo_tf_step_info = TSys(TransferFunction( [[Ti.sys.num[0][0] for Ti in Tr] for Tr in Ta], [[Ti.sys.den[0][0] for Ti in Tr] for Tr in Ta])) - T.step_info = [[Ti.step_info for Ti in Tr] for Tr in Ta] + mimo_tf_step_info.step_info = [[Ti.step_info for Ti in Tr] + for Tr in Ta] # enforce enough sample points for all channels (they have different # characteristics) - T.kwargs['step_info'] = {'T_num': 2000} - return T - + mimo_tf_step_info.kwargs['step_info'] = {'T_num': 2000} - @pytest.fixture - def tsystem(self, - request, - siso_ss1, siso_ss2, siso_tf1, siso_tf2, - mimo_ss1, mimo_ss2, mimo_tf2, - siso_dtf0, siso_dtf1, siso_dtf2, - siso_dss1, siso_dss2, - mimo_dss1, mimo_dss2, mimo_dtf1, - pole_cancellation, no_pole_cancellation, siso_tf_type1, - siso_tf_kpos, siso_tf_kneg, - siso_tf_step_matlab, siso_ss_step_matlab, - mimo_ss_step_matlab, mimo_tf_step_info, - siso_tf_asymptotic_from_neg1): - systems = {"siso_ss1": siso_ss1, - "siso_ss2": siso_ss2, - "siso_tf1": siso_tf1, - "siso_tf2": siso_tf2, - "mimo_ss1": mimo_ss1, - "mimo_ss2": mimo_ss2, - "mimo_tf2": mimo_tf2, - "siso_dtf0": siso_dtf0, - "siso_dtf1": siso_dtf1, - "siso_dtf2": siso_dtf2, - "siso_dss1": siso_dss1, - "siso_dss2": siso_dss2, - "mimo_dss1": mimo_dss1, - "mimo_dss2": mimo_dss2, - "mimo_dtf1": mimo_dtf1, - "pole_cancellation": pole_cancellation, - "no_pole_cancellation": no_pole_cancellation, - "siso_tf_type1": siso_tf_type1, - "siso_tf_kpos": siso_tf_kpos, - "siso_tf_kneg": siso_tf_kneg, - "siso_tf_step_matlab": siso_tf_step_matlab, - "siso_ss_step_matlab": siso_ss_step_matlab, - "mimo_ss_step_matlab": mimo_ss_step_matlab, - "mimo_tf_step": mimo_tf_step_info, - "siso_tf_asymptotic_from_neg1": siso_tf_asymptotic_from_neg1, - } - return systems[request.param] + systems = locals() + if isinstance(request.param, str): + return systems[request.param] + else: + return [systems[sys] for sys in request.param] @pytest.mark.parametrize( "kwargs", @@ -402,11 +321,12 @@ def tsystem(self, {'X0': np.array([0, 0])}, {'X0': 0, 'return_x': True}, ]) - def test_step_response_siso(self, siso_ss1, kwargs): + @pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True) + def test_step_response_siso(self, tsystem, kwargs): """Test SISO system step response""" - sys = siso_ss1.sys - t = siso_ss1.t - yref = siso_ss1.ystep + sys = tsystem.sys + t = tsystem.t + yref = tsystem.ystep # SISO call out = step_response(sys, T=t, **kwargs) tout, yout = out[:2] @@ -414,19 +334,21 @@ def test_step_response_siso(self, siso_ss1, kwargs): np.testing.assert_array_almost_equal(tout, t) np.testing.assert_array_almost_equal(yout, yref, decimal=4) - def test_step_response_mimo(self, mimo_ss1): - """Test MIMO system, which contains ``siso_ss1`` twice""" - sys = mimo_ss1.sys - t = mimo_ss1.t - yref = mimo_ss1.ystep + @pytest.mark.parametrize("tsystem", ["mimo_ss1"], indirect=True) + def test_step_response_mimo(self, tsystem): + """Test MIMO system, which contains ``siso_ss1`` twice.""" + sys = tsystem.sys + t = tsystem.t + yref = tsystem.ystep _t, y_00 = step_response(sys, T=t, input=0, output=0) _t, y_11 = step_response(sys, T=t, input=1, output=1) np.testing.assert_array_almost_equal(y_00, yref, decimal=4) np.testing.assert_array_almost_equal(y_11, yref, decimal=4) - def test_step_response_return(self, mimo_ss1): - """Verify continuous and discrete time use same return conventions""" - sysc = mimo_ss1.sys + @pytest.mark.parametrize("tsystem", ["mimo_ss1"], indirect=True) + def test_step_response_return(self, tsystem): + """Verify continuous and discrete time use same return conventions.""" + sysc = tsystem.sys sysd = c2d(sysc, 1) # discrete time system Tvec = np.linspace(0, 10, 11) # make sure to use integer times 0..10 Tc, youtc = step_response(sysc, Tvec, input=0) @@ -434,10 +356,9 @@ def test_step_response_return(self, mimo_ss1): np.testing.assert_array_equal(Tc.shape, Td.shape) np.testing.assert_array_equal(youtc.shape, youtd.shape) - @pytest.mark.parametrize("dt", [0, 1], ids=["continuous", "discrete"]) def test_step_nostates(self, dt): - """Constant system, continuous and discrete time + """Constant system, continuous and discrete time. gh-374 "Bug in step_response()" """ @@ -515,7 +436,7 @@ def test_step_info(self, tsystem, systype, time_2d, yfinal): @pytest.mark.parametrize( "tsystem", ['mimo_ss_step_matlab', - pytest.param('mimo_tf_step', marks=slycotonly)], + pytest.param('mimo_tf_step_info', marks=slycotonly)], indirect=["tsystem"]) def test_step_info_mimo(self, tsystem, systype, yfinal): """Test step info for MIMO systems.""" @@ -551,13 +472,15 @@ def test_step_info_invalid(self): with pytest.raises(ValueError, match="matching time vector"): step_info(np.ones((2, 2, 15))) # no time vector - def test_step_pole_cancellation(self, pole_cancellation, - no_pole_cancellation): + @pytest.mark.parametrize("tsystem", + [("no_pole_cancellation", "pole_cancellation")], + indirect=True) + def test_step_pole_cancellation(self, tsystem): # confirm that pole-zero cancellation doesn't perturb results # https://github.com/python-control/python-control/issues/440 - step_info_no_cancellation = step_info(no_pole_cancellation) - step_info_cancellation = step_info(pole_cancellation) - self.assert_step_info_match(no_pole_cancellation, + step_info_no_cancellation = step_info(tsystem[0].sys) + step_info_cancellation = step_info(tsystem[1].sys) + self.assert_step_info_match(tsystem[0].sys, step_info_no_cancellation, step_info_cancellation) @@ -570,7 +493,7 @@ def test_step_pole_cancellation(self, pole_cancellation, ("siso_dtf0", {})], indirect=["tsystem"]) def test_impulse_response_siso(self, tsystem, kwargs): - """Test impulse response of SISO systems""" + """Test impulse response of SISO systems.""" sys = tsystem.sys t = tsystem.t yref = tsystem.yimpulse @@ -581,12 +504,13 @@ def test_impulse_response_siso(self, tsystem, kwargs): np.testing.assert_array_almost_equal(tout, t) np.testing.assert_array_almost_equal(yout, yref, decimal=4) - def test_impulse_response_mimo(self, mimo_ss2): - """"Test impulse response of MIMO systems""" - sys = mimo_ss2.sys - t = mimo_ss2.t + @pytest.mark.parametrize("tsystem", ["mimo_ss2"], indirect=True) + def test_impulse_response_mimo(self, tsystem): + """"Test impulse response of MIMO systems.""" + sys = tsystem.sys + t = tsystem.t - yref = mimo_ss2.yimpulse + yref = tsystem.yimpulse _t, y_00 = impulse_response(sys, T=t, input=0, output=0) np.testing.assert_array_almost_equal(y_00, yref, decimal=4) _t, y_11 = impulse_response(sys, T=t, input=1, output=1) @@ -599,19 +523,21 @@ def test_impulse_response_mimo(self, mimo_ss2): @pytest.mark.skipif(StrictVersion(sp.__version__) < "1.3", reason="requires SciPy 1.3 or greater") - def test_discrete_time_impulse(self, siso_tf1): + @pytest.mark.parametrize("tsystem", ["siso_tf1"], indirect=True) + def test_discrete_time_impulse(self, tsystem): # discrete time impulse sampled version should match cont time dt = 0.1 t = np.arange(0, 3, dt) - sys = siso_tf1.sys + sys = tsystem.sys sysdt = sys.sample(dt, 'impulse') np.testing.assert_array_almost_equal(impulse_response(sys, t)[1], impulse_response(sysdt, t)[1]) - def test_impulse_response_warnD(self, siso_ss1): + @pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True) + def test_impulse_response_warnD(self, tsystem): """Test warning about direct feedthrough""" with pytest.warns(UserWarning, match="System has direct feedthrough"): - _ = impulse_response(siso_ss1.sys, siso_ss1.t) + _ = impulse_response(tsystem.sys, tsystem.t) @pytest.mark.parametrize( "kwargs", @@ -621,12 +547,13 @@ def test_impulse_response_warnD(self, siso_ss1): {'X0': np.array([[0.5], [1]])}, {'X0': np.array([0.5, 1]), 'return_x': True}, ]) - def test_initial_response(self, siso_ss1, kwargs): + @pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True) + def test_initial_response(self, tsystem, kwargs): """Test initial response of SISO system""" - sys = siso_ss1.sys - t = siso_ss1.t + sys = tsystem.sys + t = tsystem.t x0 = kwargs.get('X0', 0) - yref = siso_ss1.yinitial if np.any(x0) else np.zeros_like(t) + yref = tsystem.yinitial if np.any(x0) else np.zeros_like(t) out = initial_response(sys, T=t, **kwargs) tout, yout = out[:2] @@ -634,12 +561,13 @@ def test_initial_response(self, siso_ss1, kwargs): np.testing.assert_array_almost_equal(tout, t) np.testing.assert_array_almost_equal(yout, yref, decimal=4) - def test_initial_response_mimo(self, mimo_ss1): + @pytest.mark.parametrize("tsystem", ["mimo_ss1"], indirect=True) + def test_initial_response_mimo(self, tsystem): """Test initial response of MIMO system""" - sys = mimo_ss1.sys - t = mimo_ss1.t + sys = tsystem.sys + t = tsystem.t x0 = np.array([[.5], [1.], [.5], [1.]]) - yref = mimo_ss1.yinitial + yref = tsystem.yinitial yref_notrim = np.broadcast_to(yref, (2, len(t))) _t, y_00 = initial_response(sys, T=t, X0=x0, input=0, output=0) @@ -667,16 +595,23 @@ def test_forced_response_step(self, tsystem): [np.zeros((10,), dtype=float), 0] # special algorithm ) - def test_forced_response_initial(self, siso_ss1, u): - """Test forced response of SISO system as intitial response""" - sys = siso_ss1.sys - t = siso_ss1.t - x0 = np.array([[.5], [1.]]) - yref = siso_ss1.yinitial - - tout, yout = forced_response(sys, t, u, X0=x0) - np.testing.assert_array_almost_equal(tout, t) - np.testing.assert_array_almost_equal(yout, yref, decimal=4) + @pytest.mark.parametrize("tsystem", ["siso_ss1", "siso_tf2"], + indirect=True) + def test_forced_response_initial(self, tsystem, u): + """Test forced response of SISO system as intitial response.""" + sys = tsystem.sys + t = tsystem.t + x0 = tsystem.X0 + yref = tsystem.yinitial + + if isinstance(sys, StateSpace): + tout, yout = forced_response(sys, t, u, X0=x0) + np.testing.assert_array_almost_equal(tout, t) + np.testing.assert_array_almost_equal(yout, yref, decimal=4) + else: + with pytest.warns(UserWarning, match="Non-zero initial condition " + "given for transfer function"): + tout, yout = forced_response(sys, t, u, X0=x0) @pytest.mark.parametrize("tsystem, useT", [("mimo_ss1", True), @@ -718,6 +653,66 @@ def test_forced_response_legacy(self): t, y = ct.forced_response(sys, T, U) t, y, x = ct.forced_response(sys, T, U, return_x=True) + @pytest.mark.parametrize( + "tsystem, fr_kwargs, refattr", + [pytest.param("siso_ss1", + {'T': np.linspace(0, 1, 10)}, 'yinitial', + id="ctime no U"), + pytest.param("siso_dss1", + {'T': np.arange(0, 5, 1,)}, 'yinitial', + id="dt=True, no U"), + pytest.param("siso_dtf1", + {'U': np.ones(5,)}, 'ystep', + id="dt=True, no T"), + pytest.param("siso_dtf2", + {'U': np.ones(25,)}, 'ystep', + id="dt=0.2, no T"), + pytest.param("siso_ss2_dtnone", + {'U': np.ones(10,)}, 'ystep', + id="dt=None, no T"), + pytest.param("siso_dtf3", + {'U': np.ones(10,)}, 'ystep', + id="dt with rounding error, no T"), + ], + indirect=["tsystem"]) + def test_forced_response_T_U(self, tsystem, fr_kwargs, refattr): + """Test documented forced_response behavior for parameters T and U.""" + if refattr == 'yinitial': + fr_kwargs['X0'] = tsystem.X0 + t, y = forced_response(tsystem.sys, **fr_kwargs) + np.testing.assert_allclose(t, tsystem.t) + np.testing.assert_allclose(y, getattr(tsystem, refattr), rtol=1e-3) + + @pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True) + def test_forced_response_invalid_c(self, tsystem): + """Test invalid parameters.""" + with pytest.raises(TypeError, + match="StateSpace.*or.*TransferFunction"): + forced_response("not a system") + with pytest.raises(ValueError, match="T.*is mandatory for continuous"): + forced_response(tsystem.sys) + with pytest.raises(ValueError, match="time values must be equally " + "spaced"): + forced_response(tsystem.sys, [0, 0.1, 0.12, 0.4]) + + @pytest.mark.parametrize("tsystem", ["siso_dss2"], indirect=True) + def test_forced_response_invalid_d(self, tsystem): + """Test invalid parameters dtime with sys.dt > 0.""" + with pytest.raises(ValueError, match="can't both be zero"): + forced_response(tsystem.sys) + with pytest.raises(ValueError, match="must have same elements"): + forced_response(tsystem.sys, + T=tsystem.t, U=np.random.randn(1, 12)) + with pytest.raises(ValueError, match="must have same elements"): + forced_response(tsystem.sys, + T=tsystem.t, U=np.random.randn(12)) + with pytest.raises(ValueError, match="must match sampling time"): + forced_response(tsystem.sys, T=tsystem.t*0.9) + with pytest.raises(ValueError, match="must be multiples of " + "sampling time"): + forced_response(tsystem.sys, T=tsystem.t*1.1) + # but this is ok + forced_response(tsystem.sys, T=tsystem.t*2) @pytest.mark.parametrize("u, x0, xtrue", [(np.zeros((10,)), @@ -848,10 +843,11 @@ def test_default_timevector_functions_d(self, fun, dt): @pytest.mark.parametrize("tsystem", ["siso_ss2", # continuous "siso_tf1", - "siso_dss1", # no timebase + "siso_dss1", # unspecified sampling time "siso_dtf1", "siso_dss2", # matching timebase "siso_dtf2", + "siso_ss2_dtnone", # undetermined timebase "mimo_ss2", # MIMO pytest.param("mimo_tf2", marks=slycotonly), "mimo_dss1", @@ -876,9 +872,9 @@ def test_time_vector(self, tsystem, fun, squeeze, matarrayout): kw['T'] = t if fun == forced_response: kw['U'] = np.vstack([np.sin(t) for i in range(sys.ninputs)]) - elif fun == forced_response and isctime(sys): + elif fun == forced_response and isctime(sys, strict=True): pytest.skip("No continuous forced_response without time vector.") - if hasattr(tsystem.sys, "nstates"): + if hasattr(sys, "nstates"): kw['X0'] = np.arange(sys.nstates) + 1 if sys.ninputs > 1 and fun in [step_response, impulse_response]: kw['input'] = 1 @@ -892,6 +888,8 @@ def test_time_vector(self, tsystem, fun, squeeze, matarrayout): if hasattr(tsystem, 't'): # tout should always match t, which has shape (n, ) np.testing.assert_allclose(tout, tsystem.t) + elif fun == forced_response and sys.dt in [None, True]: + np.testing.assert_allclose(np.diff(tout), 1.) if squeeze is False or not sys.issiso(): assert yout.shape[0] == sys.noutputs @@ -899,20 +897,22 @@ def test_time_vector(self, tsystem, fun, squeeze, matarrayout): else: assert yout.shape == tout.shape - if sys.dt > 0 and sys.dt is not True and not np.isclose(sys.dt, 0.5): + if sys.isdtime(strict=True) and sys.dt is not True and not \ + np.isclose(sys.dt, 0.5): kw['T'] = np.arange(0, 5, 0.5) # incompatible timebase with pytest.raises(ValueError): fun(sys, **kw) @pytest.mark.parametrize("squeeze", [None, True, False]) - def test_time_vector_interpolation(self, siso_dtf2, squeeze): - """Test time vector handling in case of interpolation + @pytest.mark.parametrize("tsystem", ["siso_dtf2"], indirect=True) + def test_time_vector_interpolation(self, tsystem, squeeze): + """Test time vector handling in case of interpolation. Interpolation of the input (to match scipy.signal.dlsim) gh-239, gh-295 """ - sys = siso_dtf2.sys + sys = tsystem.sys t = np.arange(0, 10, 1.) u = np.sin(t) x0 = 0 @@ -928,7 +928,8 @@ def test_time_vector_interpolation(self, siso_dtf2, squeeze): assert yout.shape == tout.shape assert np.allclose(tout[1:] - tout[:-1], sys.dt) - def test_discrete_time_steps(self, siso_dtf2): + @pytest.mark.parametrize("tsystem", ["siso_dtf2"], indirect=True) + def test_discrete_time_steps(self, tsystem): """Make sure rounding errors in sample time are handled properly These tests play around with the input time vector to make sure that @@ -936,7 +937,7 @@ def test_discrete_time_steps(self, siso_dtf2): gh-332 """ - sys = siso_dtf2.sys + sys = tsystem.sys # Set up a time range and simulate T = np.arange(0, 100, 0.2) @@ -968,10 +969,11 @@ def test_discrete_time_steps(self, siso_dtf2): with pytest.raises(ValueError): step_response(sys, T) - def test_time_series_data_convention_2D(self, siso_ss1): + @pytest.mark.parametrize("tsystem", ["siso_ss1"], indirect=True) + def test_time_series_data_convention_2D(self, tsystem): """Allow input time as 2D array (output should be 1D)""" tin = np.array(np.linspace(0, 10, 100), ndmin=2) - t, y = step_response(siso_ss1.sys, tin) + t, y = step_response(tsystem.sys, tin) assert isinstance(t, np.ndarray) and not isinstance(t, np.matrix) assert t.ndim == 1 assert y.ndim == 1 # SISO returns "scalar" output diff --git a/control/timeresp.py b/control/timeresp.py index fd2e19c91..989a832cb 100644 --- a/control/timeresp.py +++ b/control/timeresp.py @@ -211,39 +211,45 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, 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 float, optional - Input array giving input at each time `T` (default = 0). + If None, `U` must be given and `len(U)` time steps of sys.dt are + simulated. If sys.dt is None or True (undetermined time step), a time + step of 1.0 is assumed. - If `U` is ``None`` or ``0``, a special algorithm is used. This special - algorithm is faster than the general algorithm, which is used - otherwise. + U : array_like or float, optional + Input array giving input at each time `T`. + If `U` is None or 0, `T` must be given, even for discrete + time systems. In this case, for continuous time systems, a direct + calculation of the matrix exponential is used, which is faster than the + general interpolating algorithm used otherwise. - X0 : array_like or float, optional - Initial condition (default = 0). + X0 : array_like or float, default=0. + Initial condition. - transpose : bool, optional + transpose : bool, default=False If True, transpose all input and output arrays (for backward - compatibility with MATLAB and :func:`scipy.signal.lsim`). Default - value is False. + compatibility with MATLAB and :func:`scipy.signal.lsim`). - interpolate : bool, optional (default=False) + interpolate : bool, default=False If True and system is a discrete time system, the input will be interpolated between the given time steps and the output will be given at system sampling rate. Otherwise, only return the output at the times given in `T`. No effect on continuous - time simulations (default = False). + time simulations. - return_x : bool, optional - If True (default), return the the state vector. Set to False to - return only the time and output vectors. + return_x : bool, default=None + - If False, return only the time and output vectors. + - If True, also return the the state vector. + - If None, determine the returned variables by + config.defaults['forced_response.return_x'], which was True + before version 0.9 and is False since then. squeeze : bool, optional By default, if a system is single-input, single-output (SISO) then the output response is returned as a 1D array (indexed by time). If - squeeze=True, remove single-dimensional entries from the shape of - the output even if the system is not SISO. If squeeze=False, keep + `squeeze` is True, remove single-dimensional entries from the shape of + the output even if the system is not SISO. If `squeeze` is False, keep the output as a 2D array (indexed by the output number and time) - even if the system is SISO. The default value can be set using + even if the system is SISO. The default behavior can be overridden by config.defaults['control.squeeze_time_response']. Returns @@ -252,13 +258,15 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, Time values of the output. yout : array - Response of the system. If the system is SISO and squeeze is not + Response of the system. If the system is SISO and `squeeze` is not True, the array is 1D (indexed by time). If the system is not SISO or - squeeze is False, the array is 2D (indexed by the output number and + `squeeze` is False, the array is 2D (indexed by the output number and time). xout : array - Time evolution of the state vector. Not affected by squeeze. + Time evolution of the state vector. Not affected by `squeeze`. Only + returned if `return_x` is True, or `return_x` is None and + config.defaults['forced_response.return_x'] is True. See Also -------- @@ -277,7 +285,8 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, -------- >>> T, yout, xout = forced_response(sys, T, u, X0) - See :ref:`time-series-convention`. + See :ref:`time-series-convention` and + :ref:`package-configuration-parameters`. """ if not isinstance(sys, (StateSpace, TransferFunction)): @@ -294,10 +303,17 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, "return_x specified for a transfer function system. Internal " "conversion to state space used; results may meaningless.") + # If we are passed a transfer function and X0 is non-zero, warn the user + if isinstance(sys, TransferFunction) and np.any(X0 != 0): + warnings.warn( + "Non-zero initial condition given for transfer function system. " + "Internal conversion to state space used; may not be consistent " + "with given X0.") + sys = _convert_to_statespace(sys) A, B, C, D = np.asarray(sys.A), np.asarray(sys.B), np.asarray(sys.C), \ np.asarray(sys.D) -# d_type = A.dtype + # d_type = A.dtype n_states = A.shape[0] n_inputs = B.shape[1] n_outputs = C.shape[0] @@ -306,65 +322,64 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, if U is not None: U = np.asarray(U) if T is not None: + # T must be array-like T = np.asarray(T) # Set and/or check time vector in discrete time case - if isdtime(sys, strict=True): + if isdtime(sys): if T is None: - if U is None: - raise ValueError('Parameters ``T`` and ``U`` can\'t both be' + if U is None or (U.ndim == 0 and U == 0.): + raise ValueError('Parameters ``T`` and ``U`` can\'t both be ' 'zero for discrete-time simulation') # Set T to equally spaced samples with same length as U if U.ndim == 1: n_steps = U.shape[0] else: n_steps = U.shape[1] - T = np.array(range(n_steps)) * (1 if sys.dt is True else sys.dt) + dt = 1. if sys.dt in [True, None] else sys.dt + T = np.array(range(n_steps)) * dt else: # Make sure the input vector and time vector have same length - # TODO: allow interpolation of the input vector if (U.ndim == 1 and U.shape[0] != T.shape[0]) or \ (U.ndim > 1 and U.shape[1] != T.shape[0]): - ValueError('Pamameter ``T`` must have same elements as' - ' the number of columns in input array ``U``') + raise ValueError('Pamameter ``T`` must have same elements as' + ' the number of columns in input array ``U``') + if U.ndim == 0: + U = np.full((n_inputs, T.shape[0]), U) + else: + if T is None: + raise ValueError('Parameter ``T`` is mandatory for continuous ' + 'time systems.') # Test if T has shape (n,) or (1, n); - # T must be array-like and values must be increasing. - # The length of T determines the length of the input vector. - if T is None: - raise ValueError('Parameter ``T``: must be array-like, and contain ' - '(strictly monotonic) increasing numbers.') T = _check_convert_array(T, [('any',), (1, 'any')], 'Parameter ``T``: ', squeeze=True, transpose=transpose) - dt = T[1] - T[0] - if not np.allclose(T[1:] - T[:-1], dt): - raise ValueError("Parameter ``T``: time values must be " - "equally spaced.") + n_steps = T.shape[0] # number of simulation steps + # equally spaced also implies strictly monotonic increase, + dt = (T[-1] - T[0]) / (n_steps - 1) + if not np.allclose(np.diff(T), dt): + raise ValueError("Parameter ``T``: time values must be equally " + "spaced.") + # create X0 if not given, test if X0 has correct shape X0 = _check_convert_array(X0, [(n_states,), (n_states, 1)], 'Parameter ``X0``: ', squeeze=True) - # If we are passed a transfer function and X0 is non-zero, warn the user - if isinstance(sys, TransferFunction) and np.any(X0 != 0): - warnings.warn( - "Non-zero initial condition given for transfer function system. " - "Internal conversion to state space used; may not be consistent " - "with given X0.") - xout = np.zeros((n_states, n_steps)) xout[:, 0] = X0 yout = np.zeros((n_outputs, n_steps)) # Separate out the discrete and continuous time cases - if isctime(sys): + if isctime(sys, strict=True): # Solve the differential equation, copied from scipy.signal.ltisys. dot = np.dot # Faster and shorter code # Faster algorithm if U is zero - if U is None or (isinstance(U, (int, float)) and U == 0): + # (if not None, it was converted to array above) + if U is None or np.all(U == 0): # Solve using matrix exponential expAdt = sp.linalg.expm(A * dt) for i in range(1, n_steps): @@ -410,7 +425,11 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, else: # Discrete type system => use SciPy signal processing toolbox - if sys.dt is not True: + + # sp.signal.dlsim assumes T[0] == 0 + spT = T - T[0] + + if sys.dt is not True and sys.dt is not None: # Make sure that the time increment is a multiple of sampling time # First make sure that time increment is bigger than sampling time @@ -420,12 +439,21 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, # Now check to make sure it is a multiple (with check against # sys.dt because floating point mod can have small errors - elif not (np.isclose(dt % sys.dt, 0) or - np.isclose(dt % sys.dt, sys.dt)): + if not (np.isclose(dt % sys.dt, 0) or + np.isclose(dt % sys.dt, sys.dt)): raise ValueError("Time steps ``T`` must be multiples of " "sampling time") sys_dt = sys.dt + # sp.signal.dlsim returns not enough samples if + # T[-1] - T[0] < sys_dt * decimation * (n_steps - 1) + # due to rounding errors. + # https://github.com/scipyscipy/blob/v1.6.1/scipy/signal/ltisys.py#L3462 + scipy_out_samples = int(np.floor(spT[-1] / sys_dt)) + 1 + if scipy_out_samples < n_steps: + # parantheses: order of evaluation is important + spT[-1] = spT[-1] * (n_steps / (spT[-1] / sys_dt + 1)) + else: sys_dt = dt # For unspecified sampling time, use time incr @@ -434,7 +462,8 @@ def forced_response(sys, T=None, U=0., X0=0., transpose=False, # Use signal processing toolbox for the discrete time simulation # Transpose the input to match toolbox convention - tout, yout, xout = sp.signal.dlsim(dsys, np.transpose(U), T, X0) + tout, yout, xout = sp.signal.dlsim(dsys, np.transpose(U), spT, X0) + tout = tout + T[0] if not interpolate: # If dt is different from sys.dt, resample the output diff --git a/doc/conventions.rst b/doc/conventions.rst index 4a3d78926..8a4d11992 100644 --- a/doc/conventions.rst +++ b/doc/conventions.rst @@ -83,17 +83,17 @@ The timebase argument can be given when a system is constructed: * dt = 0: continuous time system (default) * dt > 0: discrete time system with sampling period 'dt' * dt = True: discrete time with unspecified sampling period -* dt = None: no timebase specified +* dt = None: no timebase specified Only the :class:`StateSpace`, :class:`TransferFunction`, and :class:`InputOutputSystem` classes allow explicit representation of discrete time systems. -Systems must have compatible timebases in order to be combined. A discrete time -system with unspecified sampling time (`dt = True`) can be combined with a system +Systems must have compatible timebases in order to be combined. A discrete time +system with unspecified sampling time (`dt = True`) can be combined with a system having a specified sampling time; the result will be a discrete time system with the sample time of the latter system. Similarly, a system with timebase `None` can be combined with a system having a specified -timebase; the result will have the timebase of the latter system. For continuous +timebase; the result will have the timebase of the latter system. For continuous time systems, the :func:`sample_system` function or the :meth:`StateSpace.sample` and :meth:`TransferFunction.sample` methods can be used to create a discrete time system from a continuous time system. See :ref:`utility-and-conversions`. The default value of 'dt' can be changed by @@ -179,6 +179,10 @@ can be computed like this:: ft = D * U + +.. currentmodule:: control +.. _package-configuration-parameters: + Package configuration parameters ================================ @@ -207,25 +211,25 @@ on standard configurations. Selected variables that can be configured, along with their default values: * bode.dB (False): Bode plot magnitude plotted in dB (otherwise powers of 10) - + * bode.deg (True): Bode plot phase plotted in degrees (otherwise radians) - + * bode.Hz (False): Bode plot frequency plotted in Hertz (otherwise rad/sec) - + * bode.grid (True): Include grids for magnitude and phase plots - + * freqplot.number_of_samples (None): Number of frequency points in Bode plots - + * freqplot.feature_periphery_decade (1.0): How many decades to include in the frequency range on both sides of features (poles, zeros). - + * statesp.use_numpy_matrix (True): set the return type for state space matrices to `numpy.matrix` (verus numpy.ndarray) - * statesp.default_dt and xferfcn.default_dt (None): set the default value of dt when + * statesp.default_dt and xferfcn.default_dt (None): set the default value of dt when constructing new LTI systems - * statesp.remove_useless_states (True): remove states that have no effect on the + * statesp.remove_useless_states (True): remove states that have no effect on the input-output dynamics of the system Additional parameter variables are documented in individual functions