lsim is slow

[cwrowley]

The lsim command can be much, much slower than scipy.signal.lsim. For instance, consider the following example (run in IPython for %timeit):

import scipy.signal as sig
import control
import numpy as np
t = np.arange(0,1000,0.5)
u = np.sin(2.*t)
sys1 = sig.lti(1.,[1,0,1])
sys2 = control.tf(1.,[1,0,1])
%timeit sig.lsim(sys1,u,t)
%timeit control.lsim(sys2,u,t)

For me, scipy.signal.lsim gives 17.4 ms per loop, and control.lsim gives 7.64 s per loop.

The current implementation of control.lsim calls a general-purpose ODE solver at each step. I think it is better to assume a linear system with constant-spaced timesteps in t (as done in scipy.signal.lsim and in Matlab's lsim), and then one can evaluate the matrix exponential once and apply it at each step. If non-uniform timesteps are desired (which is very rare, at least in my experience), then we could have an alternative routine control.lsim2 that calls the general-purpose ODE solver--that is what is done in scipy.signal.

[murrayrm]

Since I think we already depend on scipy.signal, can we just call through to that?

[cwrowley]

Unfortunately, no, because the LTI objects in scipy.signal are totally different from the TransferFunction and StateSpace objects in control.

This raises another issue, that ideally it would be nice if these objects were all compatible with one another, since a lot of stuff is already included in scipy.signal. However, those are pretty different, so that would require a lot of work.

[blinkminster]

Defaulting to the general-purpose solver was also motivated by how scipy's
lsim doesn't work for systems with a pole at the origin, which is
relatively common in control systems. So any function that tries to use
scipy.signal.lsim or whatever algorithm that uses should also fall back to
the general-purpose solver if there is a zero at the origin.

Sawyer

---

Sawyer B. Fuller, Ph.D.
Harvard School of Engineering and Applied Science
Wyss Institute for Biologically Inspired Engineering
http://people.seas.harvard.edu/~minster/

On Fri, Mar 20, 2015 at 9:18 PM, Richard Murray notifications@github.com
wrote:

Since I think we already depend on scipy.signal, can we just call through
to that?

Reply to this email directly or view it on GitHub
#48 (comment)
.

[cwrowley]

I was just planning to call scipy.linalg.expm(A * dt) where A is from the state-space representation of the system and dt is the (constant) timestep. I believe that should work fine even if there is a zero at the origin.

[blinkminster]

I've never looked into why scipy's lsim fails with a pole at the origin,
but if what you're talking about is faster than a standard ode solver, is
correct, and works for a zero pole, then that sounds great to me.

Sawyer

---

Sawyer B. Fuller, Ph.D.
Harvard School of Engineering and Applied Science
Wyss Institute for Biologically Inspired Engineering
http://people.seas.harvard.edu/~minster/

On Fri, Mar 20, 2015 at 10:31 PM, Clancy Rowley notifications@github.com
wrote:

I was just planning to call scipy.linalg.expm(A * dt) where A is from the
state-space representation of the system and dt is the (constant)
timestep. I believe that should work fine even if there is a zero at the
origin.

Reply to this email directly or view it on GitHub
#48 (comment)
.