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

I did try python-control as well. The state space discretization in
pyton-control is just a thin layer on top of scipy.signal.

Besides performance difference in expm and exp there is also a
noticeable performance difference in math.exp and numpy.exp.

To summarize, there seems to be a simple way to optimize the 1D state
space case  giving a performance boost of about 15-20x by implementing
it using the trivial math operations.

Some performance measurements:

In[35]: math.exp(1)
Out[35]: 2.718281828459045
In [36]: np.exp(1)
Out[36]: 2.718281828459045
In [39]: linalg.expm([[1]])
Out[39]: array([[2.71828183]])
In [44]: m = timeit.timeit('math.exp(1)', setup='import math', number=1000)
In [46]: n = timeit.timeit('numpy.exp(1)', setup='import numpy', number=1000)
In [47]: l = timeit.timeit('linalg.expm([[1]])', setup='from scipy
import linalg', number=1000)
In [48]: l / n
Out[48]: 5.557622207273845
In [49]: n / m
Out[49]: 12.922172140014915
In [50]: l / m
Out[50]: 71.81655085156227

On Sat, Oct 26, 2019 at 12:22 PM Ilhan Polat <[email protected]> wrote:
>
> 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
_______________________________________________
SciPy-User mailing list
[email protected]
https://mail.python.org/mailman/listinfo/scipy-user