Re: Discretisation of continuous state space model almost 20 times slower than manual code
Clancy Rowley <[email protected]> Fri, 25 Oct 2019 08:50:55 -0400
| Newsgroups | gmane.comp.python.scientific.user |
|---|---|
| Message-ID | <[email protected]> |
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