การประมาณ Transfer function ของระบบที่เราไม่ทราบสมการพลศาสตร์ (System Indentification) ด้วยสมการอนุพันธ์อันดับสอง

clc
    clear

    % Declaration of Omega_n Value
    wn = 1.4925;   %***** Our system **** %

    % Zeta Values for Each Variation
    z1 = 1;
    z2 = 2;
    z3 = 0.1708;   %***** Our system **** %
    z4 = 0.5;
    z5 = 0;

    % Numerator and Denominator Values for Each Variation
    num1 = wn^2;
    den1 = [1 2*z1*wn wn^2];

    num2 = wn^2;
    den2 = [1 2*z2*wn wn^2];

    num3 = wn^2;
    den3 = [1 2*z3*wn wn^2];

    num4 = wn^2;
    den4 = [1 2*z4*wn wn^2];

    num5 = wn^2;
    den5 = [1 2*z5*wn wn^2];

    figure;

    % When Zeta = 1
    sys_tf1 = tf(num1, den1);

    % Pole Zero Map
    subplot(2,5,1)
    pzmap(sys_tf1);
    title(['Pole Zero Map (Zeta = ', num2str(z1), ')'])

    % Step Response
    subplot(2,5,6)
    step(sys_tf1);
    title(['Step Response (Zeta = ', num2str(z1), ')'])

    % When Zeta = 2
    sys_tf2 = tf(num2, den2);

    % Pole Zero Map
    subplot(2,5,2)
    pzmap(sys_tf2);
    title(['Pole Zero Map (Zeta = ', num2str(z2), ')'])

    % Step Response
    subplot(2,5,7)
    step(sys_tf2);
    title(['Step Response (Zeta = ', num2str(z2), ')'])

    % When Zeta = 0.1
    sys_tf3 = tf(num3, den3);

    % Pole Zero Map
    subplot(2,5,3)
    pzmap(sys_tf3);
    title(['Pole Zero Map (Zeta = ', num2str(z3), ')'])

    % Step Response
    subplot(2,5,8)
    step(sys_tf3);
    title(['Step Response (Zeta = ', num2str(z3), ')'])

    % When Zeta = 0.5
    sys_tf4 = tf(num4, den4);

    % Pole Zero Map
    subplot(2,5,4)
    pzmap(sys_tf4);
    title(['Pole Zero Map (Zeta = ', num2str(z4), ')'])

    % Step Response
    subplot(2,5,9)
    step(sys_tf4);
    title(['Step Response (Zeta = ', num2str(z4), ')'])

    % When Zeta = 0
    sys_tf5 = tf(num5, den5);

    % Pole Zero Map
    subplot(2,5,5)
    pzmap(sys_tf5);
    title(['Pole Zero Map (Zeta = ', num2str(z5), ')'])

    % Step Response
    subplot(2,5,10)
    step(sys_tf5,5);
    title(['Step Response (Zeta = ', num2str(z5), ')'])

    % Displaying damping mode information and time domain specifications
    for z = [z1, z2, z3, z4, z5]
        if z == 0
            disp(['Mode of damping for Zeta = ', num2str(z), ': Undamped']);
        elseif 0 < z && z < 1
            % Displaying damping mode information
            disp(['Mode of damping for Zeta = ', num2str(z), ': Underdamped']);
            
            % Calculating and displaying time domain specifications
            % Rise Time
            tr = pi / (wn * sqrt(1 - z^2));
            disp(['Rise Time: ', num2str(tr)]);
            
            % Settling Time
            ts = 4 / (z * wn);
            disp(['Settling Time: ', num2str(ts)]);
            
            % Natural Frequency of Oscillation
            fn = wn * sqrt(1 - z^2);
            disp(['Natural Frequency of Oscillation: ', num2str(fn)]);
        elseif z == 1
            disp(['Mode of damping for Zeta = ', num2str(z), ': Critically Damped']);
        elseif z > 1
            disp(['Mode of damping for Zeta = ', num2str(z), ': Overdamped']);
        end
        disp(' '); % Empty line break
    end

figure(2);
bode(sys_tf1,sys_tf2,sys_tf3,sys_tf4,sys_tf5)
title('Bode Plot')
grid on;

figure(3)
step(sys_tf3);
title(['Step Response (Zeta = ', num2str(z3), ')'])
grid on;