Python: State-Space model test with numpy

This is state space model test with numpy. Sorry if it includes mistakes.
The gist for the code is python:state-space model test with numpy.

You can use a module for space state model with quantecon.LinearStateSpace in QuantEcon.

Simple state-space model:

Equations:

Differential equation：

$y'' = ay' + by$

Turn into 1st order differential equation：

$x_{2}' = ax_{2} + bx_{1},\; x_{1}=y,\; x_{2}=y', \; x_{2}=x_{1}'$

State:

$\binom{x_{1}'}{x_{2}'}= \begin{pmatrix}0&1 \\ a&b \end{pmatrix}\binom{x_{1}}{x_{2}}$

Output:

$y = \binom{0}{1}\binom{x_{1}}{x_{2}}$

Code:

# Differential equation：　y'' = ay' + by
# Turn into 1st order differential equation：　x2' = a*x2 + b*x1; x1=y, x2=y', x2=x1'
# State：　[x1', x2'] = [[0,1], [a,b]] * [x1, x2]
# Output：　y = [1,0]*[x1, x2]

class SpaceState(object):
    def __init__(self):
        self.a = 0.2
        self.b = 1.3
        
        self.x1 = 1.3
        self.x2 = 0.6

    def calc(self):
        vec = np.array([self.x1, self.x2])
        mat = np.array([[0,1], [self.a, self.b]]) #x2 goes to x1. x2=x2*a + x1*b.
        
        y = []
        for i in range(10):
            vec = mat.dot(vec)
            y.append(vec[0]) # Output x1.
        return y

ss1 = SpaceState()
plt.plot(ss1.calc())

# Change x2.
ss2 = SpaceState()
ss2.x2 = 2.3
plt.plot(ss2.calc())

Spring model:

Equations:

Differential equation:

$m{y}''\left ( t\right ) = -c{y}'\left ( t \right ) -ky\left ( t \right ) + u\left ( t \right )$

Turn into 1st order differential equation：

$x_{2}' = (-c/m)x_{2} + (-k/m)x_{1},\; x_{2}=x_{1}'$

State:

$\binom{x_{1}'}{x_{2}'} = \begin{pmatrix}0 & 1\\-k/m & -c/m\end{pmatrix}\binom{x_{1}}{x_{2}} + u\binom{0}{1/m}$

Output:

$y = \binom{1}{0}\binom{x_{1}}{x_{2}}$

# DIfferential equation：　my’’(t) =  -cy’(t) -ky(t) + u(t)
# Turn into 1st order differential equation：　x2' = -c/m*x2 + -k/m*x1; x2=x1'
# State：　[x1’,x2’] = [[0,1], [-k/m,-c/m]]*[x1,x2] + [0,1/m]*u
# Output：　y = [1,0]*[x1,x2]

class Spring(object):
    def __init__(self):
        self.time = 10
        
        self.k = 0.2
        self.c = 1.3
        self.m = 3
        
        self.x1 = 1.3
        self.x2 = 0.6
        self.u = 2.3
    
    def calc(self):
        vec = np.array([self.x1, self.x2])
        mat = np.array([[0,1], [-self.k/self.m, -self.c/self.m]])
        u_mat = np.array([0, 1/self.m])
        
        y = []
        for i in range(self.time):
            vec = mat.dot(vec) + u_mat.dot(self.u) # State.
            y.append(vec[0]) # Output x1.
        return y

spring1 = Spring()
spring1.m = 1.3
plt.plot(spring1.calc())

# Change value.
spring2 = Spring()
spring2.x2 = 2.1
plt.plot(spring2.calc())

Simple state-space model:

Equations:

Code:

Spring model:

Equations:

Share it!

Related posts:

Leave a Reply Cancel reply