Simulating a mass attached to a spring with control theory

Control system theory is very useful when it comes to simulate physical systems and their behaviour.
Suppose you have a mass attached to a spring on an horizontal flat surface. The force diagram associated with this physical system can be sketched as follows:



By analysing the force diagram, you can actually see that on the x axis the only forces acting on the mass are the elastic force (Fe), the friction (Fr) due to the air resistance and the external force applied to the mass (u). On the y axis the normal force and the force of gravity balance out, therefore there is no motion.

Setting up the physical system

The elastic force can be assumed to be proportional to the spring elongation as such:
$$\hat F_e = -k x(t) \hat u_x $$
where k is the spring constant (N/m). I am assuming the spring is ideal. This assumption can be a good approximation if the spring is not streched too far.
The force due to air friction could be assumed to be proportional to the mass velocity
$$\hat F_r = - \nu \dot x(t) \hat u_x $$
in fact, we could even assume that the spring is in water, so that the friction due to the fluid is not neglectable. This model works very good as a first approximation of a such a system in a fluid, be it air, water or another fluid.
The external applied force will be denoted as $u(t)$.
If we apply Newton's second law, these are the equations that we obtain on $x$ and $y$:
$$N - mg = 0$$
$$-k x(t) - \nu \dot x(t) + u(t) = m \ddot x(t) $$
Now, equations such as this one can be solved in the time domain or in the frequency domain immediately, however, using a control system approach could help us understanding better what is going on.

State space representation

For starters, we need to put the system in a state space representation as follows
$$  x_1(t) = x(t)         \\
    x_2(t) = \dot x(t)    \\
    \begin{cases}
    \dot x_1(t) = x_2(t)   \\
    \dot x_2(t) = \frac{-k}{m}x_1 - \frac{-\nu}{m}x_2 + \frac{1}{m}u(t) \\
    y(t) = x_1(t) \\
    \end{cases}
$$
Or even better, using matrices
$$\left\{\begin{matrix}
\begin{pmatrix}
\dot x_1\\
\dot x_2
\end{pmatrix}
=
\begin{pmatrix}
0 & 1 \\
\frac{-k}{m} & \frac{-\nu}{m}
\end{pmatrix}
\begin{pmatrix}
x_1\\
x_2
\end{pmatrix}
+
\begin{pmatrix}
0\\
1
\end{pmatrix}
u(t)=AX+BU
\\
\\
y(t) = \begin{pmatrix}
1 & 0
\end{pmatrix}
\begin{pmatrix}
x_1\\
x_2
\end{pmatrix}
=CX+DU
\end{matrix}\right.
$$
After uploading the model in Matlab we need to check whether the system is stable or unstable. A good guess is that the system should be stable since there is no reason for it to be unstable, when you pull the spring, it should oscillate and then, after a certain period of time, it should get back to the equilibrium state.
For linear systems such as this one, asymptotical stability can be checked by making sure that all the eigenvalues of the matrix $A$ have negative real part. I will not dig deep into this property however you can easily look it up on Wikipedia.
Asymptotically stable systems are very easy to deal with. In case the system is not stable there are methods to make it stable.

We can easily note that the system is indeed stable.

Frequency domain analysis

By calculating the transfer function of the system, $G(s)$ we will be able to study the system in the frequency domain.



The Bode diagram shows a resonance peak at about 0.071 Hz (remember $w = 2 \pi f$) and a sharp change of phase of 180 degrees at the natural frequency. This is what we expected since the system has a couple of complex conjungate poles.
What does that resonance peak means? Well, if we were to set $u(t)=sin(w_c t)$ then after a certain period of time, because of the frequency response theorem, we expect the system to behave exactly like this
$$y_{asymptotical}(t) = \left | G(jw_c) \right |sin(w_ct+\angle G(jw_c))$$
that is, $y$ oscillates with the same frequency as $u$ and (not necessarily) a different phase. In our case the phase of $y$ is about 90 degrees off the phase of $u$ when $w_c = 0.441 rad/s$. We can simulate such $u$ and $y$ with Matlab if you do not take my word for this ;-)

Furthermore, we know that the amplitude of the wave at the exit should be $10^{1 \frac{28.9Db}{20}}$ that is about 27.86.



And this is exactly what we see. To think about it, with relatively little effort since $u$ has amplitude 1, you can make the system move with an amplitude 28 times bigger. 
Finally, we can see that if the frequency is too high, the amplitude of the sine wave at the exit will be very low. This is somewhat familiar, try to move a mass attached to a spring with very high frequency. You will see that the spring will move very little.



Note how much smaller is the amplitude. It is very close to zero after about 150 seconds.
Now we can simulate how the system is going to react when $u$ is an impulse or a step:




Note that the mass oscillates: this is exactly what we expected.
When $u(t)$ is an impulse the system oscillates and then it gets back to its original equilibrium state since there is no more external force applied. However, when $u(t)=1$ that is, $u(t)$ is a unit step, the system slowly reaches a new constant value which is exactly $G(0)=5$ the gain of the transfer function.

Putting the system in a closed loop

In order to spice things up, let's put the system in a closed loop configuration with a proportional controller $k$ and see what happens when $k$ varies.



By sketching the root locus for $k > 0$ we actually see that there is no value of $k>0$ for which the system goes unstable, that is, the closed loop system poles are all in the left half complex plane. By the way, whatever $k$ the system will always oscillate since the poles are always complex conjugates.



Note that by increasing $k$ you are actually moving the poles further away from the real axis, this basically means that as $k$ gets bigger, the frequency of the oscillations gets higher. We can simulate
the closed loop system behaviour to the unit step with $k=50$




As you can see the frequency of the oscillations has got higher. In fact, the following relationship holds:
$$T = \frac{2 \pi}{w \sqrt{1-\xi ^2}}$$
where $\xi$ is the dumping coefficient.
This application of control theory was very easy and simple, however it is easy to see how this could be extended to any linear physical system. For instance, it is very straight forward to extend this to a circuit with capacitors and inductors.

I hope this was useful and interesting.