systprop.m

matlab的FDC工具箱

function [prop,eigA]=systprop(A,B,C,D)
% SYSTPROP - computes the following properties of a linear system:
%
%          1. time constant                             tau
%          2. natural frequency of undamped system      w0
%          3. eigenfrequency of the system              wn
%          4. period                                    P
%          5. damping factor                            zeta
%          6. percentage overshoot                      PO
%          7. peak-time                                 Tpeak
%          8. settling time to a 2% band                Tset
%          9. halve-time                                Thalve
%
% Conclusions are displayed on screen and stored in a user-specified
% diary-file (by default SYSTPROP.LOG). This includes information about
% the stability of the system and possible oscillatory behaviour, and
% a short summary of the meaning of some system properties. 
%
% The default location in which the diary-file will be stored is the 
% FDC data directory; if SYSTPROP is not able to locate that directory, 
% it will default to the current working directory. You can change the 
% destination directory and filename in the save dialog. Creation of 
% the diary-file can be skipped altogether by selecting 'Cancel' or 
% by clicking the Close button in the save dialog. 
%
% [X,Y] = SYSTPROP(A,B,C,D) or [X,Y] = SYSTPROP(A,B,C) returns the 
%    system properties in the matrix X and the eigenvalues of A in the 
%    columnvector Y. In addition, this call returns information about 
%    the controllability and observability of the system.
%
%    The colums of the matrix X correspond with the system properties 
%    1 through 9 in the table above: 
%
%      X = [tau, w0, wn, P, zeta, PO, Tpeak, Tset, Thalve]
%
%    If the elements of the D matrix are all zero, it isn't necessary to
%    enter the D matrix; in that case you may use [X,Y] = SYSTPROP(A,B,C).
%
% [X,Y] = SYSTPROP(A) returns the system properties in the matrix X and 
%    the eigenvalues of A in the columvector Y; see the definitions above.
%
% [X,Y] = SYSTPROP(num,den) evaluates the equivalent transfer function
%    representation, again returning the matrix X and vector Y. It also
%    checks whether there are zeros in the right half of the s-plane, 
%    to identify minimum-phase and non-minimum-phase systems.
%
% X = SYSTPROP(Y) returns the system properties in the matrix X if Y 
%    contains the eigenvalues of a system in a columnvector. If Y is a 
%    scalar, it will be treated as a single eigenvalue, instead of a
%    1x1 A-matrix. To analyze a system consisting of a 1x1 A-matrix, 
%    use the call [X,Y] = SYSTPROP(A,0,0) instead.
%
% Note: this function requires the Control System Toolbox! In addition, 
% be aware that SYSTPROP uses the routine NUM2STR2 for output formatting.
%
% See also DAMP, CTRB, OBSV.


% ---------------------------------------------------------------------
% REFERENCE: J. van de Vegte, 'Feedback Control Systems', Prentice Hall
%            International Editions, London, 2nd edition, 1990.
% ---------------------------------------------------------------------


% Disable Matlab warning messages and clear the command-window, to prevent the 
% screen output from being distorted. This function was originally written for 
% Matlab 3.5, and although it has been updated since then, it may cause inadver-
% tent warning messages due to slightly obsolete coding. 
% ------------------------------------------------------------------------------
warning off
clc

nargs=nargin;

% If number of inputs is zero: return help information only. Else, compute
% the system properties and store the results in a diary-file.
% ------------------------------------------------------------------------
if nargs == 0
    helpwin systprop
else

    % Diary-file specification. The while-loop gives the user the possibility to
    % return to filename-specification when the cancel-button of the uiputfile GUI
    % is pressed inadvertantly. The variable usediary determines whether or not
    % a diary-file will be used. By default it is 'on'. This variable can be 
    % changed only by editing systprop.m. 
    % ----------------------------------------------------------------------------
    usediary = 'on';
    ok = 'No';
    while not(strcmp(ok,'Yes'))

        % Let the user specify a suitable name and location for the SYSTPROP diary-file.
        % The default file-name is SYSPROP.LOG; the default location is the FDC data-
        % directory, or, if that directory can't be found, the current working directory.
        % -------------------------------------------------------------------------------
        defaultdir = feval('datadir');
        workingdir = pwd;

        if exist(defaultdir,'dir')
            [filename,pathname] = uiputfile([defaultdir filesep 'systprop.log'], ...
                                            'Select name and location of diary-file:');
        else
            [filename,pathname] = uiputfile([workingdir filesep 'systprop.log'], ...
                                            'Select name and location of diary-file:');
        end

        % If cancel-button was clicked, ask whether or not to omit creating the diary-file.
        % If the user doesn't want a diary-file, change the variable usediary to 'off'.
        % ---------------------------------------------------------------------------------
        if not(ischar(filename) | ischar(pathname))
            ok= questdlg(['No diary-file has been specified. Continue without diary-file?'], ...
                          'WARNING', ...
                          'Yes','No','No');
            if strcmp(ok,'Yes')
                usediary = 'off';
            end
        else
            ok = 'Yes';
        end

    end

    % If required, activate specified diary-file. If file already exists, delete
    % it (a question whether or not to replace an existing file will be created
    % by uiputfile itself, so we don't need to take additional precautions here).
    % ---------------------------------------------------------------------------
    if strcmp(usediary,'on')
        location = [pathname filename];
        if exist(location,'file')
            delete(location);
        end
        diary(location);
    end


    % Now determine which actions to take, depending on the number of input arguments
    % -------------------------------------------------------------------------------
    if nargs >= 3              % Three or four inputs: treat inputs as system matrices
                               % of a linear state-space system [A,B,C,D]

        % If the D matrix is not specified, create D-matrix as a matrix of zeros
        % ----------------------------------------------------------------------
        if nargs == 3;
            [m1,n1] = size(B);
            [m2,n2] = size(C);
            D = zeros(m2,n1);
        end

        % Display system information as specified in SYSTPROP inputs
        % ----------------------------------------------------------
        disp('Linear state-space system with system matrices A, B, C, D:');
        A,B,C,D
        disp(' ');

        % Display error message if the dimensions of the system matrices don't match
        % --------------------------------------------------------------------------
        error(abcdchk(A,B,C,D));

        % Display information about controllability and observability of the system
        % -------------------------------------------------------------------------
        Co = ctrb(A,B);
        Ob = obsv(A,C);
        rankA = rank(A);

        disp(' ');
        disp(' ');
        if rank(Co) == rankA & rank(Ob) == rankA
            disp('System is controllable and observable.');
        elseif rank(Co) == rankA
            disp('System is controllable.');
        elseif rank(Ob) == rankA
            disp('System is observable.');
        else
            disp('System is neither controllable nor observable');
        end;

        % Determine eigenvalues of the system matrix A
        % --------------------------------------------
        eigA = eig(A);

    elseif nargs == 1          % Only one input: treat input as system matrix A or as a vector
                               % with eigenvalues Y. One scalar input will always be treated as
                               % a single eigenvalue instead of a 1x1 matrix A. It is still
                               % possible to input a 1x1 matrix A by defining the matrices B,
                               % C, and D too, e.g.: SYSTPROP(A,0,0,0) is a valid call.
        [m,n] = size(A);

        if n==1                % Treat input A as a vector with eigenvalues.

            % Display system information as specified in SYSTPROP input
            % ---------------------------------------------------------
            Y = A;
            disp('Vector with eigenvalues Y:');
            Y

            % Store vector with eigenvalues in eigA
            % -------------------------------------
            eigA = A;

        elseif m == n          % Treat input A as state matrix.

            % Display system information as specified in SYSTPROP input
            % ---------------------------------------------------------
            disp('System matrix A:');
            A

            % Determine eigenvalues of the system matrix A
            % --------------------------------------------
            eigA = eig(A);

        else
            error('Error:  Incorrect input-type for SYSTPROP');
        end

    elseif nargs == 2          % Two inputs: treat inputs A and B as numerator and
                               % denominator of a transfer function. In this case, A
                               % is a vector with coefficients of the numerator and B
                               % is a vector with coefficients of the denominator.

        % Display system information as specified in SYSTPROP inputs
        % ----------------------------------------------------------
        num = A;
        den = B;
        tf(num,den)

        % Determine poles of the system
        % -----------------------------
        eigA = roots(den);
        eigB = roots(num);
        if max(sign(eigB)) == 1
            disp('Zeros found in right half-plane => non-minimum phase system.')
        elseif max(sign(eigB)) == -1
            disp('All zeros in left half-plane => minimum phase system.')
        end
        disp(' '); 

    else                       % More than four inputs: not correct.
        error('Error: too many inputs for SYSTPROP')
    end


    % Find additional system properties by analyzing the eigenvalues of the 
    % system. First, determine the number of system modes.
    % ---------------------------------------------------------------------
    [m,n] = size(eigA);
    
    if m == 1
       disp('There is one system mode:');
    else