actrim.m

MATLAB在飞行动力学和控制中应用的工具

% The FDC toolbox. Trim routine ACTRIM.
% =====================================
% ACTRIM is a Matlab program, which computes steady-state trimmed
% flight conditions for nonlinear aircraft models in SIMULINK,
% which must have the same structure as the system BEAVER from
% the toolbox FDC 1.3. This trim routine is based upon an algo-
% rithm from ref.[1].
%
% The routine can be applied to all aircraft models which use
% the same definitions of the input and state vectors as the
% system BEAVER, being:
%
% x = [V alpha beta p q r psi theta phi xe ye H]',
% u = [deltae deltaa deltar deltaf n pz,
%                                   ,uw vw ww uwdot vwdot wwdot]'
%
% where: V      = airspeed [m/s]
%        alpha  = angle of attack [rad]
%        beta   = sideslip angle [rad]
%        p      = roll rate [rad/s]
%        q      = pitch rate [rad/s]
%        r      = yaw rate [rad/s]
%        psi    = yaw angle [rad]
%        theta  = pitch angle [rad]
%        phi    = roll angle [rad]
%        xe     = X-coordinate in Earth-fixed reference frame [m]
%        ye     = Y-coordinate in Earth-fixed reference frame [m]
%        H      = altitude above sea-level [m]
%
% and:   deltae = elevator deflection [rad]
%        deltaa = ailerons deflection [rad]
%        deltar = rudder deflection [rad]
%        deltaf = flap setting [rad]
%        n      = engine speed [RPM]
%        pz     = manifold pressure ["Hg]
%        uw     = wind & turbulence speed along Xb-axis [m/s]
%        vw     = wind & turbulence speed along Yb-axis [m/s]
%        ww     = wind & turbulence speed along Zb-axis [m/s]
%        uwdot  = d(ut)/dt [m/s^2]
%        vwdot  = d(vt)/dt [m/s^2]
%        wwdot  = d(vt)/dt [m/s^2]
%
% Although the wind and turbulence velocities and their time-
% derivatives are not involved in the trim process, these varia-
% bles cannot be ignored, because they are part of the input-
% vector to the nonlinear aircraft model. The trim program
% ACTRIM sets these variables to zero during trimming.
%
% It is necessary to edit ACTRIM if another definition of the in-
% put or state vector is used. To facilitate this, many comment
% lines have been included in the source file ACTRIM.M.
%
% Subroutines: ACCONSTR.M, evaluates flightpath constraints
%              ACCOST.M, evaluates cost function
%              COMMENT.M, evaluates... (well, you better try it yourself)
%              DATADIR.M, determines the default directory where FDC
%                  stores its datafiles
%              LOADER.M is called to load parameter matrices and vectors
%                  for the nonlinear aircraft model (type HELP LOADER for
%                  more info).
%
% Results from ACTRIM.
% ====================
% ACTRIM stores its results in the following variables:
%
%    xinco   : state vector in trimmed condition (see the defini-
%              tion, used in system BEAVER, given above),
%    uaero0  : inputvector to aerodynamic model, which is defined
%              by:  uaero = [deltae deltaa deltar deltaf]' in the
%              system BEAVER,
%    uprop0  : inputvector to engine model, which is defined by:
%              uprop = [n pz]' in the system BEAVER,
%    xdot0   : time-derivative of state vector in trimmed-flight
%              condition (Vdot, alphadot, betadot, pdot, qdot,
%              and rdot should be zero in a perfect steady-state
%              trim),
%    trimdef : text matrix, which contains the exact definition
%              of the trimmed flight condition.
%
% Note: in the system BEAVER, the total inputvector equals:
%
%     u = [ua' ut' uwind']',
%     where: uwind = [uw vw ww uwdot vwdot wwdot]'.
%
%
% Steady-state value of output vector
% ===================================
% ACTRIM does not return the initial value of the OUTPUTVECTOR
% from the system that is trimmed, because the definition of y
% is not standardized in the FDC package. If you do need the
% outputs for the trimmed-flight condition, run ACTRIM first,
% and type:
%
%     y0 = sfun(0, xinco, [uaero0; uprop0; 0;0;0;0;0;0], 3)
%
% where 'sfun' is the name of the system, e.g. sfun = beaver. The
% six zeros correspond with the wind velocities and their time-
% derivatives, which for trimming purposes are set zero (ACTRIM
% needs to be edited if you need trimmed-flight conditions for
% non-zero wind velocities!); the vector [uaero0;uprop0;0;0;0;0;0;0]
% is the total inputvector to the Simulink system sfun which con-
% tains the nonlinear aircraft model.
%
% Type HELP SFUNC for more info about this function-call.
%
% References
% ==========
% [1]   Stevens, B.L., Lewis, F.L: 'Aircraft Control and Simula-
%       tion'. John Wiley & Sons Inc., 1992.


% List of variables for ACTRIM
% ============================
% alpha       angle of attack [rad]
% answ        general: string variable used for storing user answers to yes/no
%             questions
% answ1       same as answ for nested questions
% beta        sideslip angle [rad]
% datadir     not a variable but a Matlab function, determines the default
%             data-directory for the FDC toolbox
% defdir      default directory where the FDC toolbox stores and searches for
%             datafiles
% deltaa      ailerons deflection [rad]
% deltae      elevator angle [rad]
% delttaf     flap angle [rad]
% deltar      rudder angle [rad]
% dirname     string variable with name of data-directory, entered by the user
% filename    string variable with name of datafile, entered by the user
% H           altitude [m]
% ii          counter
% jj          counter
% line#       string vectors, used to fill the text matrix 'trimdef' with infor-
%             mation about the linearized model
% loadcmmnd   string variable in which the load-command is defined (used if the
%             user wants to load an operating point from file, where the names
%             of the file, file-extension, and data-directory are obtained from
%             the variables filename and dirname)
% n           engine speed [RPM]
% ok          in general: flag used for while-loops in which the user must enter
%             certain parameters, etc.; while-loop is not quitted until this
%             variable has got the value 1 (for instance if the user has acknow-
%             ledged that the user-input is correct)
% ok1,ok2,... other flags, used for nested while-loops
% opt         variable used in user-menu's, which determines which menu option
%             has been selected by the user (exact meaning explained in source
%             text)
% p           roll-rate [rad/s]
% phi         roll-angle [rad]
% proceed     flag which is set to 0 if the user wants to quit the manual
%             operating-point definition; used to 'break' the corresponding
%             while-loop
% psi         yaw-angle [rad]
% pz          manifold pressure ["Hg]
% q           pitch-rate [rad/s]
% r           yaw-rate [rad/s]
% savecmmnd   string variable which specifies which results to save in which
%             file
% skip        flag which is set to 1 if user wants to quit the program; in
%             that case some parts of the program are skipped
% sysname     string variable, contains the name of the system to be
%             linearized
% t           time (t = clock, type HELP CLOCK for more info)
% t1          hours, extracted from t
% t2          minutes, extracted from t
% t3          seconds, extracted from t
% theta       pitch-angle [rad]
% trimdef     text matrix with information about the trimmed-flight operating
%             point (exists only if the user retrieves the operating point
%             from a file, created by ACTRIM, or if ACTRIM itself is called
%             during the operating point definition for the linearization
%             routine)
% turntype    variable for selecting coordinated or uncoordinated turns
%             (turntype == 'c' or turntype == 'v', respectively)
% uaero0      vector with initial values of aerodynamic control inputs,
%             defines operating point for the linearization routine
% udef        vector with numbers of the input variables which the user wants
%             to extract from the total of six control inputs and six wind &
%             turbulence inputs from the linearized model
% uinco       initial input vector, consists of uaero0, uprop0, and a zero-
%             value wind & turbulence vector
% uprop0      vector with initial values of engine control inputs, defines
%             operating point for the linearization routine
% V           airspeed [m/s]
% xdef        vector with numbers of the state variables which the user wants
%             to extract from the total of twelve state variables from the
%             linearized model
% xe          X-coordinate [m]
% xfix        gain variable which is used to specify states that should be
%             artificially fixed to their initial values (xfix is either equal
%             to 1 or equal to a vector of length twelve with elements that
%             equal 0 or 1, see the help-file XFIX.HLP or type HELP FIXSTATE
%             at the Matlab command line)
% xinco       initial value of the state vector, defines operating point for
%             linearization routine
% ye          Y-coordinate [m]
%-----------------------------------------------------------------------------

% Initialization commands
% -----------------------
format short e;
options = [];
turntype = 'c';

% Display header; welcome to ACTRIM!
% ----------------------------------
clc
disp('The FDC toolbox - ACTRIM');
disp('========================');
disp(' ');
disp('This program searches determines a steady-state trimmed-flight condition');
disp('for a non-linear aircraft model in Simulink.');
disp(' ');
disp(' ');

% Parameters for aircraft model:
%
% AM, EM, GM1, and GM2: matrices, containing parameters for the aircraft
%                       model. See also the source-code of MODBUILD.M
%                       (contained in the directory AIRCRAFT)
%
% The other model parameters will either be explained when they are used,
% or are self-explaining (at least if you know something about aircraft
% stability and control).
%
% Check if AM, EM, GM1, and GM2 have been defined in the Matlab workspace.
% If not, run LOADER to load them from file.
% -------------------------------------------------------------------------
if exist('AM')==0 | exist('EM')==0 | exist('GM1')==0 | exist('GM2')==0
   h = warndlg('First, the model parameters need to be retrieved from file (e.g. AIRCRAFT.DAT). Click ''OK'' to continue.');
   uiwait(h);
   loader
end

% Set xinco (= initial value of the state vector). Note: the value of xinco
% which is defined here doesn't really matter, because this variable won't
% be used during the trim process. The variable xinco must be defined in the
% Matlab Workspace, however, because Simulink needs it in order to define
% the system internally. The variable xinco is called by the Integrator
% block in the heart of the aircraft model (level 3: 'Aircraft Equations of
% Motion').
%
% In the aircraft model, some variables are made dimensionless by dividing
% by the dynamic pressure, 1/2*rho*V^2. In order to avoid 'Division By Zero'
% errors, the first element of xinco (=V) has been equalled to 45 [m/s],
% which is the mean velocity of the region in which the aerodynamic model
% of the 'Beaver' is valid. It could have been set to any value unequal to
% zero, however.
% -------------------------------------------------------------------------
xinco = [45 0 0 0 0 0 0 0 0 0 0 0]';

% The system BEAVER uses a gain block for arbitrarily setting time-
% derivatives of the state variables to zero. This may be useful for
% some purposes, for instance to eliminate longitudinal-lateral cross-
% coupling, or to fix the airspeed to its initial value if an altitude
% controller is to be evaluated, but no autothrottle is available yet.
% The gain value is xfix, which is set to 1 if all states may vary, i.e.,
% if the complete six degree-of-freedom model is used without restrictions.
% When states have to be fixed, the routine FIXSTATE.M can be called
% (type HELP FIXSTATE for more info).
%
% Here, all state variables are allowed to vary, so xfix will be set to
% one. Replace this line by:
%
%   fixstate;
%
% if you want to set some time-derivatives to zero, even when trimming
% the aircraft model.
% -------------------------------------------------------------------------
xfix = 1;
%%%% fixstate;

% The name of the aircraft model to be evaluated will be stored in the
% stringvariable sysname. The aircraft model must use the same definitions
% of the state and input vectors as the system BEAVER in order to function
% properly. Edit this program if your model uses other definitions!
% ------------------------------------------------------------------------
disp('Give name of system with aircraft model (8 characters max.)');
disp('default = beaver');
sysname = input('> ','s');
if isempty(sysname)
   sysname = 'beaver';
end

% Display menu in which the user can choose between a number of trimmed-
% flight conditions and quitting.
% ----------------------------------------------------------------------
clc
opt = menu('Select type of steady-state flight:',...
      'Steady wings-level flight','Steady turning flight',...
      'Steady pull-up','Steady roll','Quit');

skip  = 0; % Do not skip iteration block, unless option 'quit' is used
           % (then skip will be set to 1).

% Define flight condition, depending upon trim option chosen
% ----------------------------------------------------------
if opt == 1                                       % STEADY WINGS-LEVEL FLIGHT
                                                  % -------------------------
   clc
   disp('Steady wings-level flight.');
   disp('==========================');

   V = input('Give desired airspeed [m/s], default = 45: ');
   if isempty(V)
      V = 45;
   end

   H = input('Give (initial) altitude [m], default = 0: ');
   if isempty(H)
      H = 0;
   end

   psi  = input('Give heading [deg], default = 0: ')*pi/180;
   if isempty(psi)
      psi = 0;