systprop.m

matlab的FDC工具箱

       disp(['There are ' num2str(m) ' system modes:']);
    end
    disp(' ');

    % Next, determine the system properties for each individual mode
    % --------------------------------------------------------------
    for i = 1:1:m
        ra = real(eigA(i,1));                             % Real parts of the eigenvalues
        ia = imag(eigA(i,1));                             % Imaginary parts of the eigenvalues

        prop(i,9) =  log(0.5)/ra;                         % halve time [s]
        prop(i,1) = -1/ra;                                % time constant [s]
        prop(i,2) =  abs(eigA(i,1));                      % undamped natural frequency [rad/s]
        prop(i,5) = -ra/prop(i,2);                        % damping coefficient [-]
        prop(i,8) =  4*prop(i,1);                         % settling time [s]

        % If the eigenvalues have a non-zero imaginary part, determine the properties of the
        % oscillatory signals. For non-oscillatory behavior, the eigenfrequency equals zero,
        % the period is infinite, as is the peak time (which equals pi/eigenfrequency = pi/0).
        % ------------------------------------------------------------------------------------
        if ia ~= 0
            prop(i,3)= prop(i,2)*sqrt(1-prop(i,5)^2);     % eigenfrequency [rad]
            prop(i,4)= abs(2*pi/prop(i,3));               % period [s]
            prop(i,7)= pi/prop(i,3);                      % peak time [s]
        else
            prop(i,3)= 0;                                 % eigenfrequency [rad]
            prop(i,4)= inf;                               % 'period' [s]
            prop(i,7)= inf;                               % 'peak time' [s]
        end

        % Determine the overshoot percentage
        % ----------------------------------
        if prop(i,5) > 0 & prop(i,5) < 1;                 % damped oscillatory
            prop(i,6) = 100*exp(-pi*prop(i,5)/sqrt(1-prop(i,5)));
        elseif prop(i,5) > 1;                            % non-oscillatory and stable
            prop(i,6) = 0;
        else
            prop(i,6) = inf;                              % unstable
        end 

        % Analyze the results (categorize the characteristic responses for the different
        % eigenvalues into corresponding 'system modes' depending on their eigenvalues).
        %-------------------------------------------------------------------------------
        if ia == 0 & ra < 0
            mode ='stable, non-oscillatory';
        elseif ia == 0 & ra > 0
            mode ='unstable, non-oscillatory';
        elseif ia ~= 0 & ra < 0
            mode ='stable, oscillatory';
        elseif (ia ~= 0 & ra > 0) | (prop(i,5) < 0 & prop(i,5) > -1)
            mode ='unstable, oscillatory';
        elseif prop(i,5) == 1 | ra == 0
            mode ='critically damped, indifferent';
        elseif prop(i,5) < -1
            mode ='unstable, non-oscillatory';
        elseif prop(i,5) > 1
            mode ='non-oscillatory';
        else
            mode ='unknown (!!)';
        end;

        % Display system modes conclusions
        % --------------------------------
        if ia == 0
            disp(['Mode ' num2str(i) ' is ' mode ' with eigenvalue: ' num2str(ra) ]);
        elseif ia < 0
            disp(['Mode ' num2str(i) ' is ' mode ' with eigenvalue: ' num2str(ra), ...
                  ' ' num2str(ia) 'i']);
        else
            disp(['Mode ' num2str(i) ' is ' mode ' with eigenvalue: ' num2str(ra), ...
                  ' +' num2str(ia) 'i']);
        end

    end

    NewPage(usediary); % local function, see below

    % Display results. Note: the routine NUM2STR2 is called instead of the standard Matlab 
    % routine NUM2STR in order to obtain better screen formatting. For 1, 2, or 3 system
    % modes, the page formatting is tuned for optimal readibility. If there are more than
    % 3 system modes, a short-hand notation will be used.
    %-------------------------------------------------------------------------------------
    if m == 1

        disp('                          Mode 1')
        disp('                          ======')
        disp(['Time constant       :   ' num2str2(prop(1,1),8) '   [s]']);
        disp(['Undamped nat. freq. :   ' num2str2(prop(1,2),8) '   [rad/s]']);
        disp(['Eigenfrequency      :   ' num2str2(prop(1,3),8) '   [rad/s]']);
        disp(['Period              :   ' num2str2(prop(1,4),8) '   [s]']);
        disp(['Damping coefficient :   ' num2str2(prop(1,5),8) '   [-]']);
        disp(['Overshoot           :   ' num2str2(prop(1,6),8) '   [%]']);
        disp(['Peak time           :   ' num2str2(prop(1,7),8) '   [s]']);
        disp(['Settling time (2%)  :   ' num2str2(prop(1,8),8) '   [s]']);
        disp(['Halve time          :   ' num2str2(prop(1,9),8) '   [s]']);

    elseif m == 2

        disp('                          Mode 1       Mode 2')
        disp('                          ======       ======')
        disp(['Time constant       :   ' num2str2(prop(1,1),8), ...
              '     ' num2str2(prop(2,1),8) '   [s]']);
        disp(['Natural frequency   :   ' num2str2(prop(1,2),8), ...
              '     ' num2str2(prop(2,2),8) '   [rad/s]']);
        disp(['Eigenfrequency      :   ' num2str2(prop(1,3),8), ...
              '     ' num2str2(prop(2,3),8) '   [rad/s]']);
        disp(['Period              :   ' num2str2(prop(1,4),8), ...
              '     ' num2str2(prop(2,4),8) '   [s]']);
        disp(['Damping coefficient :   ' num2str2(prop(1,5),8), ...
              '     ' num2str2(prop(2,5),8) '   [-]']);
        disp(['Overshoot           :   ' num2str2(prop(1,6),8), ...
              '     ' num2str2(prop(2,6),8) '   [%]']);
        disp(['Peak time           :   ' num2str2(prop(1,7),8), ...
              '     ' num2str2(prop(2,7),8) '   [s]']);
        disp(['Settling time (2%)  :   ' num2str2(prop(1,8),8), ...
              '     ' num2str2(prop(2,8),8) '   [s]']);
        disp(['Halve time          :   ' num2str2(prop(1,9),8), ...
              '     ' num2str2(prop(2,9),8) '   [s]']);

    elseif m == 3

        disp('                          Mode 1       Mode 2       Mode 3')
        disp('                          ======       ======       ======')
        disp(['Time constant       :   ' num2str2(prop(1,1),8), ...
              '     ' num2str2(prop(2,1),8) '     ' num2str2(prop(3,1),8) '   [s]']);
        disp(['Natural frequency   :   ' num2str2(prop(1,2),8), ...
              '     ' num2str2(prop(2,2),8) '     ' num2str2(prop(3,2),8) '   [rad/s]']);
        disp(['Eigen frequency     :   ' num2str2(prop(1,3),8), ...
              '     ' num2str2(prop(2,3),8) '     ' num2str2(prop(3,3),8) '   [rad/s]']);
        disp(['Period              :   ' num2str2(prop(1,4),8), ...
              '     ' num2str2(prop(2,4),8) '     ' num2str2(prop(3,4),8) '   [s]']);
        disp(['Damping coefficient :   ' num2str2(prop(1,5),8), ...
              '     ' num2str2(prop(2,5),8) '     ' num2str2(prop(3,5),8) '   [-]']);
        disp(['Overshoot           :   ' num2str2(prop(1,6),8), ...
              '     ' num2str2(prop(2,6),8) '     ' num2str2(prop(3,6),8) '   [%]']);
        disp(['Peak time           :   ' num2str2(prop(1,7),8), ...
              '     ' num2str2(prop(2,7),8) '     ' num2str2(prop(3,7),8) '   [s]']);
        disp(['Settling time (2%)  :   ' num2str2(prop(1,8),8), ...
              '     ' num2str2(prop(2,8),8) '     ' num2str2(prop(3,8),8) '   [s]']);
        disp(['Halve time          :   ' num2str2(prop(1,9),8), ...
              '     ' num2str2(prop(2,9),8) '     ' num2str2(prop(3,9),8) '   [s]']);

    elseif m > 3

        disp('           tau [s]  w0 [rad/s]  wn [rad/s]       P [s]    zeta [-]');
        for i = 1:1:m
            disp(['Mode ' int2str(i) ':   ' num2str2(prop(i,1),8), ...
                  '    ' num2str2(prop(i,2),8) '    ' num2str2(prop(i,3),8), ...
                  '    ' num2str2(prop(i,4),8) '    ' num2str2(prop(i,5),8)]);
        end
        disp(' ');
        disp('            PO [%]   Tpeak [s]    Tset [s]  Thalve [s]');
        for i=1:1:m
            disp(['Mode ' int2str(i) ':   ' num2str2(prop(i,6),8), ...
                  '    ' num2str2(prop(i,7),8) '    ' num2str2(prop(i,8),8), ...
                  '    ' num2str2(prop(i,9),8)]);
        end
        disp(' ');
        disp('tau   : time constant                      PO    : percentage overshoot');
        disp('w0    : natural freq. of undamped system   Tpeak : peak-time');
        disp('wn    : eigenfrequency of the system       Tset  : settling-time (2%)');
        disp('P     : period                             Thalve: halve-time');
        disp('zeta  : damping factor');
        
    else               % Not 1, 2, 3, or more system modes? Then something must be wrong. 
        error('Invalid number of system modes!');
    end

    disp(' ');

    NewPage(usediary); % local function, see below

    disp('Notes:')
    disp('======')
    disp('  1. The system is stable (i.e. the transient solution will decay to zero)')
    disp('     if all poles are in the left half plane of the s-plane')
    disp('  2. To avoid excessive overshoot and unduly oscillator behaviour, the')
    disp('     damping ratio must be adequate, thus the angle phi (the angle between the')
    disp('     real-axis of the s-plane and the line through the origin and the pole)')
    disp('     must not be too close to zero.')
    disp('  3. The time constant and settling time can be reduced (response speed')
    disp('     increased) by increasing the negative real part of the poles.')
    disp('  4. Undamped frequency = distance from pole to origin; moving the poles out')
    disp('     radially (i.e. with constant damping factor) increases the speed of');
    disp('     response, thus reducing settling time, peak time and rise time, while the')
    disp('     percentage overshoot remains constant.')
    disp('  5. Natural frequency = resonant frequency or damped natural frequency. This')
    disp('     is the frequency of transient oscillations, which equals the imaginary')
    disp('     part of the pole positions.')
    disp('  6. Peak time and rise time are reduced by increasing the imaginary part')
    disp('     of the pole locations.')
    disp(' ');
    disp('Refer to source-code of SYSTPROP.M for more details about the computations.')
    disp('See also: J. van de Vegte, ''Feedback Control Systems'', Prentice Hall')
    disp('International Editions, London, 2nd edition, 1990.')
    
    NewPage(usediary); % local function, see below
    
    disp('Ready.')
    disp(' ')

    if strcmp(usediary,'on')
        diary off
        disp('The results have been stored in the diary-file:');
        disp(location);
    end
end

% Enable Matlab warning messages again
% ------------------------------------
warning on


% ====================
% BEGIN LOCAL FUNCTION
% ====================

    function NewPage(usediary)
    % Small local function, turns off diary function when required (that is: when
    % usediary equals 'on'), waits for user to press a key, clears the screen and
    % turns diary function back on when required. This prevents the 'press a key to
    % continue' messages to be included in the SYSTPROP logfile. 

    disp(' ');
    if strcmp(usediary,'on')
        diary off
    end
    disp('<<< Press a key to continue >>>');
    pause
    clc
    if strcmp(usediary,'on')
        diary on
    end

% ==================
% END LOCAL FUNCTION
% ==================


%-----------------------------------------------------------------------
% This program is based on a tool written by E.A. van der Zwan in 1993.
%
% The Flight Dynamics and Control Toolbox version 1.4.0. 
% (c) Copyright Marc Rauw, 1994-2005. Licensed under the Open Software 
% License version 2.1; see COPYING.TXT and LICENSE.TXT for more details.
% Last revision of this program: July 18, 2005.
%-----------------------------------------------------------------------