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 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<--------
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@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@python.org https://mail.python.org/mailman/listinfo/scipy-user
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 <clancyr@gmail.com> 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 <jonas.bulow@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@python.org https://mail.python.org/mailman/listinfo/scipy-user
_______________________________________________ SciPy-User mailing list SciPy-User@python.org https://mail.python.org/mailman/listinfo/scipy-user
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 <ilhanpolat@gmail.com> 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 <clancyr@gmail.com> 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 <jonas.bulow@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@python.org https://mail.python.org/mailman/listinfo/scipy-user
_______________________________________________ SciPy-User mailing list SciPy-User@python.org https://mail.python.org/mailman/listinfo/scipy-user
_______________________________________________ SciPy-User mailing list SciPy-User@python.org https://mail.python.org/mailman/listinfo/scipy-user
participants (3)
-
Clancy Rowley -
Ilhan Polat -
Jonas Bülow