# Second Order Linear Systems¶

## Standard form¶

The standard form for a second order linear system is given by

$$\tau^2\frac{d^2y}{dt^2} + 2\zeta\tau\frac{dy}{dt} + y = K u(t)$$

where the symbols are defined in the following table.

Quantity Symbol$\$
Input $u(t)$
Response $y(t)$
Steady State Gain $K$
Time Constant $\tau$
Damping Factor $\zeta$
Natural Frequency $\omega_n = \frac{1}{\tau}$

The standard form assumes that a zero input (i.e, $u(t) = 0$) results in a zero response ($y(t) = 0$). In practice, the nominal or quiescent value of $y$ or $u$ may different from zero. In that case we would write

$$\tau^2\frac{d^2y}{dt^2} + 2\zeta\tau\frac{dy}{dt} + y - y_{ref} = K\left(u(t) - u_{ref}\right)$$

where $u_{ref}$ and $y_{ref}$ represent constant reference values.

## Step Response¶

The step response of a system corresponds to a system that is initially at steady with $u = u_{ref}$ with a steady response $y = y_{ref}$. At time $t=0$ the input is incremented by a constant value U, i.e. $u = u_{ref} + U$ for $t \geq 0$. The subsequent response $y(t) - y_{ref}$ is the step response.

Second order linear systems have rather elegant analytica solutions expressed using exponential and trignometric functions. There are three distinct cases that depend on the value of the damping factor $\zeta$.

### Overdamped ($\zeta > 1$)¶

In the case of a large damping factor $\zeta > 1$, the step response is given by

$$y(t) = y_{ref} + KU\left(1 - \frac{\tau_1e^{-t/\tau_1} - \tau_2e^{-t/\tau_2}}{\tau_1 - \tau_2}\right)$$

where $\tau_1$ and $\tau_2$ are found by factor the polynomial

$$\tau^2s^2 + 2\zeta\tau s + 1 = (\tau_1s + 1)(\tau_2s + 1)$$

which for $\zeta \geq 1$ are given by

\begin{align*} \tau_1 & = \frac{\tau}{\zeta - \sqrt{\zeta^2 -1}} \\ \tau_2 & = \frac{\tau}{\zeta + \sqrt{\zeta^2 -1}} \end{align*}

### Critically Damped ($\zeta = 1$)¶

$$y(t) = y_{ref} + KU\left(1 - \left(1 + \frac{t}{\tau}\right)e^{-t/\tau}\right)$$

### Underdamped ($0 \leq \zeta < 1$)

$$y(t) = y_{ref} + KU\left[1 - e^{-\zeta t/\tau}\left(\cos\left(\omega t\right) + \frac{\zeta}{\sqrt{1-\zeta^2}}\,\sin\left(\omega t\right) \right)\right]$$

where

$$\omega = \frac{\sqrt{1-\zeta^2}}{\tau}$$

## Interactive Simulation¶

\begin{align*} \frac{dy}{dt} & = v \\ \frac{dv}{dt} & = -\frac{1}{\tau^2}(y-y_{ref}) - \frac{2\zeta}{\tau}v + K\left(u(t)-u_{ref}\right) \end{align*}
In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from ipywidgets import interact

def simulation(yref=0, U=1, K=1, tau=1, zeta=0.2):

def deriv(X,t):
y,v = X
ydot = v
vdot = -(y-yref)/tau/tau - 2*zeta*v/tau + K*U/tau/tau
return[ydot,vdot]

# simulation
t = np.linspace(0,20*tau,1000)
y = odeint(deriv, [yref,0], t)[:,0]

# plot steady state line and bounds
plt.figure(figsize=(12,6))

# plot solution
plt.plot(t,y,lw=3)
plt.title('Step Response of a Second Order System')
plt.xlabel('Time')
plt.ylabel('y')

# plot limits
plt.ylim(plt.ylim()[0],1.1*plt.ylim()[1])
plt.xlim(t[0],t[-1])
dy = np.diff(plt.ylim())

# arrow props
ap1 = dict(arrowstyle="->")
ap2 = dict(arrowstyle="<->")

if zeta < 1:
#overshoot
os = np.exp(-np.pi*zeta/np.sqrt(1-zeta**2))

# time to first peak
tp = np.pi*tau/np.sqrt(1-zeta**2)
yp = (1+os)*K*U + yref

plt.text(tp,yp+0.02*dy,"Overshoot\n b/a = {0:0.2f}".format(os), ha='center')
plt.annotate('',xy=(tp,K*U+yref),xytext=(tp,yp),arrowprops=ap2)
plt.text(tp,(K*U+yref+yp)/2,' b')
plt.annotate('',xy=(tp,yref),xytext=(tp,K*U+yref),arrowprops=ap2)
plt.text(tp,K*U/2+yref,' a')
plt.annotate("Time to first\n peak = {0:.2f}".format(tp),
xy=(tp,yref), xytext=(1.2*tp,0.2*K*U+yref),arrowprops=ap1)

# rise time
tr = t[np.where(np.diff(np.sign(y-yref-K*U))*np.sign(K*U)>0)[0][0]]
if tr < plt.xlim()[1]:
plt.plot([tr,tr],[0.3*K*U+yref,K*U+yref],'r:')
plt.annotate('',xy=(plt.xlim()[0],0.4*K*U+yref),xytext=(tr,0.4*K*U+yref),
arrowprops=ap2)
plt.text(plt.xlim()[0]+tr/2,0.42*K*U+yref+0.02*dy,
'Rise Time\n = {0:.2f}'.format(tr),ha='center')

# period
P = 2*np.pi*tau/np.sqrt(1-zeta**2)
if tr + P < plt.xlim()[1]:
plt.plot([tr,tr],[0.3*K*U+yref,K*U+yref],'r:')
plt.plot([tr+P,tr+P],[0.3*K*U+yref,K*U+yref],'r:')
plt.annotate('',xy=(tr,0.4*K*U+yref),xytext=(tr+P,0.4*K*U+yref),arrowprops=ap2)
plt.text(tr+P/2,0.42*K*U+yref+0.02*dy,'Period = {0:.2f}'.format(P), ha='center')

# second peak
if tp + P < plt.xlim()[1]:
plt.annotate('',xy=(tp+P,K*U+yref),xytext=(tp+P,K*U*(1+os**3)+yref),
arrowprops=ap2)
plt.text(tp+P,K*U*(1+os**3/2)+yref,' c')
plt.text(tp+P,K*U*(1+os**3)+yref+0.02*dy,
'Decay Ratio\n c/b = {0:.2f}'.format(os**2),va='bottom',ha='center')

# settling time
ts = -np.log(0.05)*np.sqrt(1-zeta**2)*tau/zeta
if ts < plt.xlim()[1]:
plt.fill_between(t[t>ts],0.95*K*U+yref,1.05*K*U+yref,alpha=0.4,color='y')
plt.text(ts,1.05*K*U+yref+0.02*dy,
'Settling Time\n = {0:.2f}'.format(ts),ha='center')

plt.plot(plt.xlim(),[yref,yref],'k--')
plt.plot(plt.xlim(),[K*U+yref,K*U+yref],'k--')

interact(simulation, yref = (-10,10,0.1), U=(0.01,5,0.01),
K = (-5,5,0.01), zeta=(0.01,3,0.01), tau = (0.1,5.0,0.01));


## Performance Indicators for Underdamped Systems¶

For an underdamped second order system, the desired performance metrics are given by the following by formulas in the following table.

Quantity Symbol Expression/Value
Rise Time $t_r$ Time to first SS crossing
Time to first peak $t_p$ $\frac{\pi\tau}{\sqrt{1-\zeta^2}}$
Overshoot OS $\exp\left(-\frac{\pi\zeta}{\sqrt{1-\zeta^2}}\right)$
Decay Ratio DR $\exp\left(-\frac{2\pi\zeta}{\sqrt{1-\zeta^2}}\right)$
Period $\frac{2\pi\tau}{\sqrt{1-\zeta^2}}$
Setting Time $t_s$ Time to +/- 5% of SS
In [2]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

zeta = np.linspace(0,0.999,100)
os = np.exp(-np.pi*zeta/np.sqrt(1-zeta**2))
dr = np.exp(-2*np.pi*zeta/np.sqrt(1-zeta**2))
pd = np.sqrt(1-zeta**2)

plt.figure(figsize=(8,8))
plt.plot(zeta,os)
plt.plot(zeta,dr)
plt.plot(zeta,pd)
plt.axis('square')
plt.xlim(0,1)
plt.ylim(0,1)
plt.title('Performance Characteristics of Underdamped Second Order Systems')
plt.xlabel('$\zeta$')
plt.ylabel('Performance Characteristic')
plt.text(0.35,0.4,'Overshoot')
plt.text(0.05,.2,'Decay Ratio')
plt.text(0.7,0.8,'Natural Period / Period')
plt.gca().set_xticks(np.arange(0,1,0.1),minor=True)
plt.gca().set_yticks(np.arange(0,1,0.1),minor=True)
plt.grid(b=True,which='major')
plt.grid(b=True,which='minor')

In [ ]: