[SciPy-User] Discretisation of continuous state space model almost 20 times slower than manual code

[Clancy Rowley]

One issue is that the first routine below will not work for systems of dimension > 1.  You'd need to use the matrix exponential, for instance with scipy.linalg.expm, instead of math.exp, which does elementwise exponentiation.

> On Oct 25, 2019, at 6:29 AM, Jonas Bülow <jonas.bulow at gmail.com> wrote:
> 
> Hi,
> 
> I'm new to scipy. I implemented a simple state space model where the A and B matrix is specified and C and D is zero. Before I knew about scipy.signal.StateSpace I manually discretised the model using scipy.integrate.quad as shown in the code below. 
> 
> Someone suggested I should use scipy.signal.StateSpace. The code become a lot slower using signal.StateSpace, almost 20 times slower. Any hints on what I'm doing wrong or how I should change my use of signal.StateSpace to get better performance? Or should I just stick to my simplified manual method (variant 1 on the code below)?
> 
> The code:
> 
> ----8<--------
> 
> import math
> import numpy as np
> from scipy import integrate, signal
> 
> import timeit as ti
> 
> def cont2disc_1(A, B, sample_time):
>     Ad = math.exp(A * sample_time)
>     def fq(x):
>         return integrate.quad(lambda t: math.exp(A * t) * B[x], 0, sample_time)[0]
>     Bd = np.array([fq(0), fq(1)])
>     return Ad, Bd
> 
> def cont2disc_2(A, B, sample_time):
>     ds = signal.StateSpace(A, B, 0, [0, 0]).to_discrete(sample_time)
>     Ad = ds.A[0][0]
>     Bd = ds.B[0]
>     return Ad, Bd
> 
> if __name__ == "__main__":
>     A,B,h = -0.0004622051395583765, [0.00046221, 0.01629134], 60
>     print(f'1: {cont2disc_1(A, B, h)}')
>     print(f'2: {cont2disc_2(A, B, h)}')
> 
>     SETUP=f'''
> from __main__ import cont2disc_1, cont2disc_2
> A,B,h = {A}, {B}, {h} 
>     '''
>     t_1 = ti.timeit('cont2disc_1(A, B, h)', SETUP, number=1000)
>     t_2  = ti.timeit('cont2disc_2(A, B, h)', SETUP, number=1000)
>     print(f'Method 1 is {t_2 / t_1} times faster then method 2')
> 
> 
> ----8<--------
> 
>     
> 
> 
> 
> 
>     
> 
> 
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at python.org
> https://mail.python.org/mailman/listinfo/scipy-user

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-user/attachments/20191025/68d092d6/attachment.html>