Setup

Let us begin by setting up our simulation time and a sample input.

% Initialize a time vector for our investigation

t_start = 0;

t_stop = 50;

t_step = 0.1;

t = 0:0.01:50;

% Initialize a step input for our investigation

u_mag = 1;

u_start = 5;

u = t;

u(u<=u_start) = 0;

u(u>u_start) = u_mag;

% Plot our input

f = figure(); f.Position = [0 0 500 300];

plot(t, u); title('sample input'); xlabel('time (s)'); ylabel('magnitude'); ylim([0 u_mag*1.1]);

Sample Second Order System:

To begin our analysis we first generate a random second order system.

For practical purposes, this system will approximate transfer functions of the form:

% Random system to approximate (second order; no zero)

k_mystery = rand;

w_mystery = rand;

z_mystery = rand;

a_mystery = rand;

sys_mystery = tf(k_mystery*[1 a_mystery], [1 2*z_mystery*w_mystery w_mystery^2]);

% Save the input frequencies magnitude & frequency response of the mystery system for later use

[mag_mystery, phase_mystery, freq_mystery] = bode(sys_mystery);

% Display the frequency response of our mystery system

f = figure(); f.Position = [0 0 800 500];

bode(sys_mystery); title('Mystery System Response');

Optimization

We can manually alter k, w, z, & a values until the system response matches the mystery response.

Note: Nesting multiple for-loops is grossly inefficient computionally. It takes approximately 5-10 minues to run this script with 20 values per parameter.

k = 0.05:0.05:1;         % range of potential values for gain
w = 0.05:0.05:1;         % range of potential values for natural frequency
z = 0.05:0.05:1;         % range of potential values for damping ratio
a = 0.05:0.05:1;         % range of potential values for the zero
 
error_all      = [];
error_combined = [];
 
% loop through every potential combination of k, w, & z
for i = 1:length(k)
    for ii = 1:length(w)
        for iii = 1:length(z)
            for iiii = 1:length(a)
                sys_candidate = tf(k(i)*[1 a(iiii)], [1 2*z(iii)*w(ii) w(ii)^2]);
                [mag_candidate, phase_candidate, freq_candidate] = bode(sys_candidate, freq_mystery);
                    
                error_mag      = sqrt(mean((mag_candidate-mag_mystery).^2));        % root mean square error
                error_phase    = sqrt(mean((phase_candidate-phase_mystery).^2));    % root mean square error
                error_all      = [error_all;      error_mag error_phase, k(i), w(ii), z(iii), a(iiii)];
                error_combined = [error_combined; error_mag+error_phase, k(i), w(ii), z(iii), a(iiii)];
            end
        end
    end
end
 
[min_error, location] = min(error_combined(:,1), [], 'ComparisonMethod','abs');
results = [error_combined(location, 2:end); k_mystery w_mystery z_mystery a_mystery];
 
error       = error_combined(location, 1);
k_estimated = error_combined(location, 2);
w_estimated = error_combined(location, 3);
z_estimated = error_combined(location, 4);
a_estimated = error_combined(location, 5);
fprintf('k_est = %g \nw_est = %g \nz_est = %g \na_est = %g\nerror = %g', ...
    k_estimated, w_estimated, z_estimated, a_estimated, error);
k_est = 0.2 
w_est = 0.5 
z_est = 0.45 
a_est = 0.7
error = 0.752665

Second Order with Zero (frequency domain)

Housekeeping