A second order model is given by the following equation:

$$ \begin{equation} \frac{Y(s)}{U(s)}=\frac{K_{ss}\,c}{s^2+bs+c}e^{-s\theta} \end{equation} $$

where $b$, and $c$ are the parameters of the system. For a general second-order system, $c$ corresponds to the system's natural frequency, while b corresponds to its damping ratio. We can predict that we will have a very large b which corresponds to its damping ratio. Since we will use an optimization technique available in MATLAB for parameter guessing, it is essential to ensure that the model’s parameters do not have redundancies. That is why we select a second-order model equation in its simplest form.

To find the values for $b$, and $c$, we will use an optimization technique available in MATLAB. Again, we need to collect the step response data, similar to what we did for the FOPTD modeling. Re-using the same data is also acceptable.

To implement optimization technique for parameter identification, we must know how to program $\frac{Y(s)}{U(s)}=\frac{K_{ss}c}{s^2+bs+c}e^{-s\theta}$ in MATLAB. Therefore, we apply the following discretization process.

First, take its time domain alternative by performing an inverse Laplace operation:

$$ \begin{align} \begin{split} cK_{ss}u(t-\theta)&=\ddot{y}(t)+b\dot{y}(t)+cy(t) \\

\ddot{y}(t)&=cK_{ss}u(t-\theta)-b\dot{y}(t)-cy(t) \end{split} \end{align} $$

Next, take:

$$ \begin{equation} \dot{y}(t) =\frac{y(k)-y(k-1)}{\Delta t} \end{equation} $$

Also take:

$$ \begin{equation} \begin{split} \ddot{y}(t)&=\frac{\dot{y}(k)-\dot{y}(k-1)}{\Delta t} \\

\ddot{y}(t)&=\frac{\frac{y(k)-y(k-1)}{\Delta t}-\frac{y(k-1)-y(k-2)}{\Delta t}}{\Delta t}\\

\ddot{y}(t)&=\frac{y(k)-2y(k-1)+y(k-2)}{\Delta t^2} \\

t &=k \Delta t \end{split} \end{equation} $$

Now that we have re-expressed $\dot{x}$ and $\ddot{x}$, we can then rewrite Eq. 2 by substituting Eq. 3 and Eq. 4 to Eq. 2. The result is as follows.

$$ \begin{equation} \begin{split} \ddot{y}(t)&=cK_{ss}u(t-\theta)-b\dot{y}(t)-cy(t) \\

\frac{y(k)-2y(k-1)+y(k-2)}{\Delta t^2}&=cK_{ss}u(t-\theta)-b\,\frac{y(k)-y(k-1)}{\Delta t}-cy(k) \\

\end{split} \end{equation} $$

Next, we rearrange Eq. 5 such that only $y(k)$ sits in the left side of the equation. The result is as follows.

$$ \begin{equation} \begin{split}

y(k)&=\frac{(2+b\Delta t) y(k-1)-y(k-2)+cK_{ss} \Delta t^2 u\bigg(k-\frac{\theta}{\Delta t}\bigg)}{\bigg[1+b\Delta t+c\Delta t^2\bigg]}

\end{split} \end{equation} $$

We can now program Eq. 6 in MATLAB, as follows.

sopdt.m

function y = sopdt(u, u0, dt, b, c, theta, Kss, y0, y0dot)
% This simulates a second-order system plus time delay
% The system model is:
% Y(s)/U(s) = Kss / (s^2 + b*s + c) * exp(-theta*s)
%
% y is the discretized system's output (vector)
%
% u is the discretized system's input (vector)
%
% u0 is the initial value of u, u0 will be used  as the
%    input u when 0 <= t <= theta
%
% y0 is the initial output (scalar)
% y0dot is the first time derivative of the initial output (scalar)
%
% dt, b, c, theta, and Kss are all scalar parameters
%
% See the explanation at: 
%   <https://crystal-rain-1f4.notion.site/Second-Order-Plus-Time-Delay-SOPTD-7fa1d65c898540cb9b0e7fecc96c322e>

% Apply the delay
y = zeros(1, length(u));
y(1:2) = [y0 y0+y0dot*dt];
u_d = zeros(1, length(u));
n = floor(theta/dt);

% Simulate the system
for k = 1 : length(u)
    if k > n
        u_d(k) = u(k-n);        
    else
        u_d(k) = u0;        
    end

    if k > 3
        A = 1 + b * dt + c * dt^2;
        y(k) = (c*Kss*dt^2*u_d(k) + b*dt*y(k-1) + 2*y(k-1) - y(k-2)) / A;
    end
end

end

The codes above is a MATLAB function. For sanity check, use the following codes.