# [SciPy-User] Discretisation of continuous state space model almost 20 times slower than manual code Classic List Threaded 4 messages Open this post in threaded view
|

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

 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 mathimport numpy as npfrom scipy import integrate, signalimport timeit as tidef 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)    Bd = np.array([fq(0), fq(1)])    return Ad, Bddef cont2disc_2(A, B, sample_time):    ds = signal.StateSpace(A, B, 0, [0, 0]).to_discrete(sample_time)    Ad = ds.A    Bd = ds.B    return Ad, Bdif __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_2A,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 [hidden email] https://mail.python.org/mailman/listinfo/scipy-user
Open this post in threaded view
|

## Re: Discretisation of continuous state space model almost 20 times slower than manual code

 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 <[hidden email]> 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 mathimport numpy as npfrom scipy import integrate, signalimport timeit as tidef 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)    Bd = np.array([fq(0), fq(1)])    return Ad, Bddef cont2disc_2(A, B, sample_time):    ds = signal.StateSpace(A, B, 0, [0, 0]).to_discrete(sample_time)    Ad = ds.A    Bd = ds.B    return Ad, Bdif __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_2A,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[hidden email]https://mail.python.org/mailman/listinfo/scipy-user_______________________________________________ SciPy-User mailing list [hidden email] https://mail.python.org/mailman/listinfo/scipy-user