Skip to content

Improved default time vector and handling for time response functions step, impulse, and initial #420

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 11 commits into from
Jul 11, 2020
Merged

Improved default time vector and handling for time response functions step, impulse, and initial #420

Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
b39184b
fixed default response time for time response of discrete-time functi…
sawyerbfuller Jun 18, 2020
32b610b
to pass tests, added convenient ability to specify simulation time an…
sawyerbfuller Jun 19, 2020
9a8da65
eliminated deprecation warnings by importing certain functions from n…
sawyerbfuller Jun 19, 2020
db9f197
eliminated deprecation warnings by importing certain functions from n…
sawyerbfuller Jun 19, 2020
6d62910
small fix to pass unit tests
sawyerbfuller Jun 19, 2020
9f02119
adjusted sisotool test so tests pass with new default step response t…
sawyerbfuller Jun 19, 2020
1bc60ba
added functionality to automatically choose dt in timeresp.py based o…
sawyerbfuller Jun 21, 2020
512f822
removed some leftover code and comments
sawyerbfuller Jun 21, 2020
79951a1
explanation in docstrings for how time vector T is auto-computed in t…
sawyerbfuller Jun 28, 2020
5a34fc0
Merge branch 'master' into fix-discrete-step-response
murrayrm Jul 11, 2020
e155a15
fix unit test precision in sisotool_test.py
murrayrm Jul 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
4 changes: 2 additions & 2 deletions control/freqplot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -822,10 +822,10 @@ def default_frequency_range(syslist, Hz=None, number_of_samples=None,

# Set the range to be an order of magnitude beyond any features
if number_of_samples:
omega = sp.logspace(
omega = np.logspace(
lsp_min, lsp_max, num=number_of_samples, endpoint=True)
else:
omega = sp.logspace(lsp_min, lsp_max, endpoint=True)
omega = np.logspace(lsp_min, lsp_max, endpoint=True)
return omega


Expand Down
31 changes: 17 additions & 14 deletions control/matlab/timeresp.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,8 +21,9 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False):
sys: StateSpace, or TransferFunction
LTI system to simulate

T: array-like object, optional
Time vector (argument is autocomputed if not given)
T: array-like or number, optional
Time vector, or simulation time duration if a number (time vector is
autocomputed if not given)

X0: array-like or number, optional
Initial condition (default = 0)
Expand Down Expand Up @@ -59,7 +60,7 @@ def step(sys, T=None, X0=0., input=0, output=None, return_x=False):
from ..timeresp import step_response

T, yout, xout = step_response(sys, T, X0, input, output,
transpose = True, return_x=True)
transpose=True, return_x=True)

if return_x:
return yout, T, xout
Expand All @@ -75,8 +76,9 @@ def stepinfo(sys, T=None, SettlingTimeThreshold=0.02, RiseTimeLimits=(0.1,0.9)):
sys: StateSpace, or TransferFunction
LTI system to simulate

T: array-like object, optional
Time vector (argument is autocomputed if not given)
T: array-like or number, optional
Time vector, or simulation time duration if a number (time vector is
autocomputed if not given)

SettlingTimeThreshold: float value, optional
Defines the error to compute settling time (default = 0.02)
Expand Down Expand Up @@ -127,9 +129,10 @@ def impulse(sys, T=None, X0=0., input=0, output=None, return_x=False):
sys: StateSpace, TransferFunction
LTI system to simulate

T: array-like object, optional
Time vector (argument is autocomputed if not given)

T: array-like or number, optional
Time vector, or simulation time duration if a number (time vector is
autocomputed if not given)

X0: array-like or number, optional
Initial condition (default = 0)

Expand Down Expand Up @@ -182,9 +185,10 @@ def initial(sys, T=None, X0=0., input=None, output=None, return_x=False):
sys: StateSpace, or TransferFunction
LTI system to simulate

T: array-like object, optional
Time vector (argument is autocomputed if not given)

T: array-like or number, optional
Time vector, or simulation time duration if a number (time vector is
autocomputed if not given)

X0: array-like object or number, optional
Initial condition (default = 0)

Expand Down Expand Up @@ -245,9 +249,8 @@ def lsim(sys, U=0., T=None, X0=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.

T: array-like
Time steps at which the input is defined, numbers must be (strictly
monotonic) increasing.
T: array-like, optional for discrete LTI `sys`
Time steps at which the input is defined; values must be evenly spaced.

X0: array-like or number, optional
Initial condition (default = 0).
Expand Down
2 changes: 1 addition & 1 deletion control/rlocus.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy import array, poly1d, row_stack, zeros_like, real, imag
from numpy import array, poly1d, row_stack, zeros_like, real, imag
import scipy.signal # signal processing toolbox
import pylab # plotting routines
from .xferfcn import _convert_to_transfer_function
Expand Down
2 changes: 1 addition & 1 deletion control/sisotool.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@ def _SisotoolUpdate(sys,fig,K,bode_plot_params,tvect=None):
# Generate the step response and plot it
sys_closed = (K*sys).feedback(1)
if tvect is None:
tvect, yout = step_response(sys_closed)
tvect, yout = step_response(sys_closed, T_num=100)
else:
tvect, yout = step_response(sys_closed,tvect)
ax_step.plot(tvect, yout)
Expand Down
4 changes: 2 additions & 2 deletions control/tests/sisotool_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ def test_sisotool(self):
assert_array_almost_equal(ax_rlocus.lines[2].get_data(),initial_point_2)

# Check the step response before moving the point
step_response_original = np.array([ 0., 0.02233651, 0.13118374, 0.33078542, 0.5907113, 0.87041549, 1.13038536, 1.33851053, 1.47374666, 1.52757114])
step_response_original = np.array([0. , 0.0217, 0.1281, 0.3237, 0.5797, 0.8566, 1.116 , 1.3261, 1.4659, 1.526 ])

@murrayrm murrayrm Jun 28, 2020

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do these values need to be changed? The system is the same and so shouldn't the step response be the same?

@sawyerbfuller sawyerbfuller Jun 28, 2020

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The function in scipy that used to be used to generate the time vector, scipy.signal.ltisys._default_response_times(A, n) uses numpy.linspace without the keyword arg endpoint=False so it returns n+1 data points instead of n. The new function does, so it generates only n data points, which is consistent with how other functions in python-control seem to do it. (thus the slightly changed values). Let me know if you have a preference to go back to n+1 datapoints. My thinking is, if you ask for n datapoints, that ought to be how many you get.

assert_array_almost_equal(ax_step.lines[0].get_data()[1][:10],step_response_original)

bode_plot_params = {
Expand Down Expand Up @@ -59,7 +59,7 @@ def test_sisotool(self):
assert_array_almost_equal(ax_mag.lines[0].get_data()[1][10:20],bode_mag_moved)

# Check if the step response has changed
step_response_moved = np.array([[ 0., 0.02458187, 0.16529784 , 0.46602716 , 0.91012035 , 1.43364313, 1.93996334 , 2.3190105 , 2.47041552 , 2.32724853] ])
step_response_moved = np.array([0. , 0.0239, 0.161 , 0.4547, 0.8903, 1.407 , 1.9121, 2.2989, 2.4686, 2.353 ])

@murrayrm murrayrm Jun 28, 2020

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same question as above. Just want to make sure I understand the need for the change.

assert_array_almost_equal(ax_step.lines[0].get_data()[1][:10],step_response_moved)


Expand Down
71 changes: 69 additions & 2 deletions control/tests/timeresp_test.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
import unittest
import numpy as np
from control.timeresp import *
from control.timeresp import _ideal_tfinal_and_dt, _default_time_vector
from control.statesp import *
from control.xferfcn import TransferFunction, _convert_to_transfer_function
from control.dtime import c2d
Expand Down Expand Up @@ -94,6 +95,7 @@ def test_step_response(self):
np.testing.assert_array_equal(Tc.shape, Td.shape)
np.testing.assert_array_equal(youtc.shape, youtd.shape)


# Recreate issue #374 ("Bug in step_response()")
def test_step_nostates(self):
# Continuous time, constant system
Expand Down Expand Up @@ -346,10 +348,75 @@ 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)
t2, y2 = step_response(sys2, input=0)
t1, y1 = step_response(sys1, input=0, T_num=100)
t2, y2 = step_response(sys2, input=0, 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
p = 0.5
np.testing.assert_array_almost_equal(
_ideal_tfinal_and_dt(TransferFunction(1, [1, .5]))[0],
(7/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
np.testing.assert_array_almost_equal(
_ideal_tfinal_and_dt(TransferFunction(1, [1, .5, 0]))[0],
(7/p))
# confirm a TF with a natural frequency of wn rad/s gets a
# dt of 1/(7.0*wn)
wn = 10
np.testing.assert_array_almost_equal(
_ideal_tfinal_and_dt(TransferFunction(1, [1, 0, wn**2]))[1],
1/(7.0*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))
# 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],
.1)
np.testing.assert_array_almost_equal(
np.diff(initial_response(TransferFunction(1, [1, 2*zeta*wn, wn**2]).sample(.1))[0][0:2]),
.1)
np.testing.assert_array_almost_equal(
_ideal_tfinal_and_dt(TransferFunction(1, [1, 2*zeta*wn, wn**2]))[1],
1/(7.0*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 impose final time
np.testing.assert_array_almost_equal(
100,
step_response(sys, 100)[0][-1],
decimal=.5)
np.testing.assert_array_almost_equal(
100,
step_response(sysdt, 100)[0][-1],
decimal=.5)
np.testing.assert_array_almost_equal(
100,
impulse_response(sys, 100)[0][-1],
decimal=.5)
np.testing.assert_array_almost_equal(
100,
initial_response(sys, 100)[0][-1],
decimal=.5)


def test_time_vector(self):
"Unit test: https://github.com/python-control/python-control/issues/239"
# Discrete time simulations with specified time vectors
Expand Down
Loading