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

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

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

Jonas Bülow
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
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: 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 <[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 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
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

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

Ilhan Polat
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 <[hidden email]> 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 <[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 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
[hidden email]
https://mail.python.org/mailman/listinfo/scipy-user

_______________________________________________
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
Reply | Threaded
Open this post in threaded view
|

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

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