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

It's the expm, exp difference as Clancy mentions. However you can instead
use python-control or harold packages for specific control purposes if you
can't make scipy ones work for you.

Admittedly in scipy.signal some LTI functions might indeed be the naive
implementations and can suffer from performance issues.

On Fri, Oct 25, 2019 at 2:51 PM Clancy Rowley <[email protected]> wrote:

> 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 <[email protected]> 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
> [email protected]
> https://mail.python.org/mailman/listinfo/scipy-user
>
>
> _______________________________________________
> SciPy-User mailing list
> [email protected]
> https://mail.python.org/mailman/listinfo/scipy-user
>

_______________________________________________
SciPy-User mailing list
[email protected]
https://mail.python.org/mailman/listinfo/scipy-user