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())

sim_utility_state-space1_simple1

 

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())

sim_utility_state-space2_spring1

 

Leave a Reply