Repository navigation
lsim is slow #48
Description
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.
Activity
murrayrm commentedon Mar 21, 2015
Since I think we already depend on
scipy.signal, can we just call through to that?cwrowley commentedon Mar 21, 2015
Unfortunately, no, because the LTI objects in
scipy.signalare totally different from theTransferFunctionandStateSpaceobjects incontrol.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 commentedon Mar 21, 2015
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:
cwrowley commentedon Mar 21, 2015
I was just planning to call
scipy.linalg.expm(A * dt)whereAis from the state-space representation of the system anddtis the (constant) timestep. I believe that should work fine even if there is a zero at the origin.blinkminster commentedon Mar 21, 2015
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:
Speed up forced_response (lsim)