State-Space Methods Part 2: Harmonic Oscillation

A simple harmonic oscillator is a standard introductory physics problem that we’ve all seen in high school. It consists of a block with mass sitting on a frictionless surface. This mass is attached to a spring, which is then attached to a wall. If the block is moved away from the Equilibrium point and then released, its trajectory will follow a sinusoidal oscillation known as Simple Harmonic Motion.

The differential equations of this mass-spring system are:
$$mv’=-kx$$
or
$$m\frac{d^2u}{d^2x}+ku=0$$
with \(m\) being the mass and \(k\) is the spring constant.

When the spring constant is 1 we get the following equations for the undamped oscillator:
$$\begin{matrix}
x’=v \\
v=-1x-0v
\end{matrix}$$

Which can be represented in the matrix form:
$$\begin{bmatrix}
x’ \\ v’
\end{bmatrix}=
\begin{bmatrix}
0 & 1 \\ -1 & 0
\end{bmatrix}
\begin{bmatrix}
x \\ v
\end{bmatrix}$$

When we look at the eigenvectors of the state matrix for this system we can see that they don’t have a real component. The imaginary eigenvectors are:
$$\lambda_{1,2} = 0 \pm 1i$$
This means that the poles of the system lie on the vertical axis in the complex s-domain. The time response is therefore a constant amplitude oscillation.

Spring Mass Damper

The spring-mass-damper builds on this fundamental introductory problem and turns it into something much more useful that has wide-ranging applications. A spring-mass-damper system is a 2nd-order system with an energy-dissipating or damping component. This is important because most physical systems have some sort of energy loss. This includes motors, structures, gyroscopes, tuning forks, etc. This system is identical to the simple harmonic oscillation except for adding a damping term. This may represent friction, inertia, drag, or resistance, or any other energy-dissipating force or component.

A mass is attached to a spring and damper, the spring and damper are attached to a fixed reference point, and there is a force applied to the mass, just like in the previous example. The total force acting on the mass is the combined force of the Spring, Damper, and external force
$$F=F_{spring}+F_{damper}+F_{applied}$$
where the spring force is
$$F_{spring}=-k_sx$$
with \(k_s\) being the spring constant. The damper force can be described with a different constant:
$$F_{damper}=-k_d \dot{x}$$
where \(k_d\) is the damping constant. When combined with Newton’s second law \(F=ma\) then the system dynamic can be described by
$$m \ddot{x}=F_{applied} – k_d\dot{x}-k_sx$$
which can be rearranged to
$$F_{applied}=m\ddot{x}+k_d\dot{x}+k_sx$$
Equations of harmonic oscillators always take the form
$$\ddot{x}+2\zeta\omega\dot{x}+\omega^2x=0$$
This gives us a second-order linear ordinary differential equation, which has the general solution of
$$x=f(t,x(0),\dot{x}(0),F(t))$$
where \(t\) is the time, \(x(0)\) is the initial position, \(\dot{x}(0)\) is the initial velocity, and \(F(t)\) is the external applied force. We can then use this to make a state-space model where our states are position and velocity, with the input being force. The state variables are:
$$x_1=x$$
$$x_2=\dot{x}$$
We then differentiate these to states by using the equations above.
$$\dot{x}_1=\dot{x}$$
$$\dot{x_2}=\ddot{x}=-\frac{k_d}{m}\dot{x}-\frac{k_s}{m}x+\frac{F}{m}$$
We can then put these equations into the following state space form :
$$\begin{bmatrix}
\dot{x}_1 \\
\dot{x}_2
\end{bmatrix} =
\begin{bmatrix}
0 & 1 \\
-\frac{k_s}{m} & -\frac{k_d}{m}
\end{bmatrix}
\begin{bmatrix}
x_1 \\ x_2
\end{bmatrix}+
\begin{bmatrix}
0 \\ \frac{1}{m}
\end{bmatrix}F$$
$$y=
\begin{bmatrix}
1 & 0
\end{bmatrix}x$$

Now let’s look at two implementations of this system. I developed a python class that can be used for any state-space implementation and a Matlab script for the same problem so that we can cross-check the answers.

Python Model Implementation

This project used a state-space model of a spring-mass damper system using the forward Euler integration method to model the time domain response of the system. The python model used in this example is available on GitHub here. I will be going through snippets of the code in this article, but I encourage you to view it in its entirety. The first part that we will look at is shown below

cuda = torch.device('cuda')
cpu = torch.device('cpu')

class StateSpace():
    def __init__(self, A, B, C, D, x, t):
        self.A = A
        self.B = B
        self.C = C
        self.D = D
        self.x = x
        self.t = t
        
    def step(self, u):
        x_dot = torch.mm(self.A,self.x)+torch.mm(self.B,u)
        self.x = torch.mm(x_dot,self.t)+self.x
        return torch.mm(self.C,self.x)+torch.mm(self.D,u)

In this class I have the basic state-space matrices, the state vector, and the time. I also defined a step function that implements the forward Euler discretization scheme (you can read more about that here). I used the common Pytorch library due to its ability to run matrix operations efficiently, but I won’t go into too much detail as that is beyond the scope of this example and we are running all of our computations on the CPU. The torch.mm() function is the Pytorch implementation of matrix multiplication. It can multiply vectors with matrices, and matrices with matrices. Now let’s look at setting up our system.

 ks = 2
kd = 2
m = 4
t = [[0.001]]
v0 = 0
x0 = 1
x = [[x0],[v0]]
A = [[0,1],[-ks/m,-kd/m]]
B = [[0],[1/m]]
C = [[1,0]]
D = [[0]]

These lines define the initial conditions of the state vector x, and the matrices for our state-space model. The mass is 4, and the spring and damping constants are both 2. Next we can build the model using our class constructor that we made previously.

# build the model
SS = StateSpace(A, B, C, D, x, t)

This simulation does not have any input forces, just the initial conditions. We can then run through our step method for 25 seconds to see the output of our system.

u = [[0]] #input force
u = torch.tensor(u,dtype=float,device=cpu)
out = []

start=time()
# run simulation
for i in range(25000):
    out.append(SS.step(u))
end=time()

The unit for the x-axis is ms, but we can see out expected damped second-order response with its overshoot and settling to zero.

Matlab Model Implementation

Matlab is absolutely wonderful and makes it insanely easy to build and analyze systems in state-space form. The code is here if you would like to try it out on your own. We will simulate the same system parameters and compare the results.

ks = 2; % spring constant
kd = 2; % damping constant
m = 4; % mass
t = 0.001; % timestep
v0 = 0; % initial velocity
x0 = 1; % initial position

A = [[0 1];[-ks/m -kd/m]] 
B = [[0];[1/m]]
C = [1 0]
D = [0]

model = ss(A,B,C,D)

initial(model,[1;0])

This is much more compact as we don’t have to write any of integration or simulation steps ourselves. We use the state space method ss() to build a state-space system from the A, B, C, and D matrices. Then we can use the initial() method to simulate the model with the same initial conditions as our python system that we listed above. The results are shown in the figure below. We can see that the responses for both of our two implementations match.

Up Next

Now this was a very basic example, but it lays the groundwork for may other items that we need to be able to build practical simulations. Next we will look at multi-dimensional spring-mass damper systems which we can use to model, cloth and structures. Stay tuned!