System Identification for FOPDT using Tangent Line method

This is a basic example of identifying First-order plus dead time (FOPDT) using tangent line method, with the real system is already given

Example:

The order of real transfer function is 3, combined with time delay = 2s.

This is just an example. In real, we only have the data describing the inputs and outputs of the system. Our mission is to find its mathematic model (state-space equation, transfer function, differentiate equation…)

A simple way we can identify the system is using identification toolbox available in MATLAB, with the input is table of data from excel. In this article, I used tangent line method to identify it.

The above figure describe three values that we need to get from the real data: process time constant t, process dead time theta, process gain k

Our mission is to find k, t and theta to form the estimated transfer function above. After that, we can use Sum of square errors (SSE) to assess the quality of estimated model compared with real model.

Firstly, we need to find the inflection point and draw a tangent line on this inflection point.

% Calculate the numerical first and second derivatives
dy_real = diff(y_real) ./ diff(t);            % First derivative
d2y_real = diff(dy_real) ./ diff(t(1:end-1)); % Second derivative

% Find the inflection point (where second derivative crosses zero)
inflection_index = find(d2y_real(1:end-1) .* d2y_real(2:end) < 0); % Find sign change
inflection_time = t(inflection_index(1) + 1); % Time at the inflection point
inflection_value = y_real(inflection_index(1) + 1); % Step response value at the inflection point

% Compute the slope (dy/dt) at the inflection point using the first derivative
inflection_slope = dy_real(inflection_index(1)); 

% Plot the inflection point on the step response
hold on;
plot(inflection_time, inflection_value, 'go', 'MarkerSize', 8, 'LineWidth', 4);

With the above code, we can identify the inflection point. Next, we identify the tangent line

%% Define the tangent line using the point-slope form y = m(x - x0) + y0

tangent_time = linspace(t(1), t(end), 100); % Time range for tangent
tangent_line = inflection_slope * (tangent_time - inflection_time) + inflection_value; % Tangent equation

% Plot the tangent line
plot(tangent_time, tangent_line, 'g--', 'LineWidth', 2);

After all, we can identify the process time constant t and process dead time theta easily.

Process dead time is the time when tangent line cross X - axis, means when y = 0

%% Find time delay 
% Find the X-axis intercept (where y = 0 for the tangent line)
x_axis_intersection_time = inflection_time - inflection_value / inflection_slope;
deadtime = x_axis_intersection_time; 
% Plot the inflection point on the step response
hold on;
plot(deadtime, 0, 'ro', 'MarkerSize', 8, 'LineWidth', 4);

To find process gain k, we just need to find the final steady state value of y_real from real model

%% Find k = y_steady_state
kx = y_real(end)
k = round(kx, 1); 

To find process time constant t, firstly we need to find the point when response reach 63.2% of the steady-state response, then align it to X - axis.

Time constant = that point - deadtime

%% Find time constant t 
y_max = round(y_real(end),1)
y_63 = 0.632 * y_max;

% Define the tolerance range
tolerance = 0.005;

% Find the approximate x value with the allowed tolerance
% Initialize variables to store the x' value and the smallest difference
x_63 = NaN;
min_diff = Inf;

for i = 1:length(y_real)
    % Calculate the difference between the current value of y and y'
    diff = abs(y_real(i) - y_63);
    
    % Check if the difference is within the allowed tolerance
    if diff <= tolerance
        % Update x' and the smallest difference
        x_63 = t(i);
        min_diff = diff;
        break; % Exit the loop when the first satisfying point is found
    end
end

% Plot the time constant point on the step response
hold on;
plot(x_63, 0, 'bo', 'MarkerSize', 8, 'LineWidth', 4);

timeconstant = x_63 - deadtime;

You can see in the code, we have to use for ...end function, because in fact the data is not 100% continous, some points are not identified, so we can consider a tolerance about 0.005 as a acceptible error from that x_63 to the real point

Finally, we just write out the estimated model

s = tf('s');  % Define the Laplace variable
G_est = k / (1 + timeconstant * s) * exp(-timedelay * s);  % FOPDT transfer function
% Generate the step response
[y_est, t_est] = step(G_est);

Here is the assesment based on Sum of square errors:

%% Calculate error between real and estimated model

% Resample y_est to match the time vector of y_real
y_est_resampled = interp1(t_est, y_est, t);  % Interpolate y_est at the time points of t
N = length(y_real);
% Calculate error between real and estimated model
error = sum((y_real - y_est_resampled).^2) / N

The final result after running the MATLAB code:

Besides, there are some other methods to identify the system. I will upload it soon.

I thankfully respect my lecturer Dr. Bich Thanh Truong Thi for teaching me this subject.