actrim.m

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

   end

   gammatype = input('Use specified manifold pressure or flight-path angle ([m]/f)? ','s');
   if isempty(gammatype)
      gammatype = 'm';
   end
   if gammatype == 'f'
      gamma = input('Give flightpath angle [deg], default = 0: ')*pi/180;
      if isempty(gamma)
         gamma = 0;
      end
   end

   phidot   = 0;
   psidot   = 0;
   thetadot = 0;

   rolltype = 's';   % No rolling, so default setting, which is stability-
                     % axes roll, will be used.

elseif opt == 2                                       % STEADY TURNING FLIGHT
                                                      % ---------------------
   clc
   disp('Steady turning flight.');
   disp('======================');

   turntype = input('Do you want a coordinated or uncoordinated turn ([c]/u)? ','s');
   if isempty(turntype)
      turntype = 'c';
   end

   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 initial heading [deg], default = 0: ')*pi/180;
   if isempty(psi)
      psi = 0;
   end

   gammatype = input('Use specified manifold pressure or flight-path angle ([m]/f)? ','s');
   if isempty(gammatype)
      gammatype = 'm';
   end
   if gammatype == 'f'
      gamma = input('Give flightpath angle [deg], default = 0: ');
      if isempty(gamma)
         gamma = 0;
      end
   end

   phidot = 0;
   thetadot = 0;

   answ = input('Do you want to specify the turnrate (1) or turnradius (2)? ');
   if answ == 1
      psidot = input('Give desired rate of turn [deg/s], default = 0: ')*pi/180;
      if isempty(psidot)
         psidot = 0;
      end
   elseif answ == 2
      R = input('Give desired turn radius (>0) [m], default = straight & level: ');
      if isempty(R) ~= 1
         psidot = V/R;
      else
         psidot = 0;
      end
   end
   clear answ

   rolltype = 's';   % No rolling, so default setting, which is stability-
                     % axes roll, will be used.

elseif opt == 3                                              % STEADY PULL-UP
                                                             % --------------
   clc
   disp('Steady pull-up.');
   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 initial heading [deg], default = 0: ')*pi/180;
   if isempty(psi)
      psi = 0;
   end

   gammatype = 'f' % Use specified flightpath angle and numerically
                   % adjust manifold pressure
   gamma = input('Give initial flightpath angle [deg], default = 0: ')*pi/180;
   if isempty(gamma)
      gamma = 0;
   end

   phidot = 0;
   psidot = 0;

   thetadot = input('Give pull-up rate [deg/s], default = 0: ')*pi/180;
   if isempty(thetadot)
      thetadot = 0;
   end

   rolltype = 's';   % No rolling, so default setting, which is stability-
                     % axes roll, will be used.

elseif opt == 4                                                 % STEADY ROLL
                                                                % -----------
   clc
   disp('Steady roll.');
   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 initial heading [deg], default = 0: ')*pi/180;
   if isempty(psi)
      psi = 0;
   end

   gammatype = input('Use specified manifold pressure or flight-path angle ([m]/f)? ','s');
   if isempty(gammatype)
      gammatype = 'm';
   end
   if gammatype == 'f'
      gamma = input('Give flightpath angle [deg], default = 0: ');
      if isempty(gamma)
         gamma = 0;
      end
   end

   thetadot = 0;
   psidot = 0;

   phidot = input('Give desired roll-rate [deg/s], default = 0: ')*pi/180;
   if isempty(phidot)
      phidot = 0;
   end

   ok = 0;
   while ok~= 1
      if phidot ~= 0
         rolltype = input('Roll in stability or body axes reference frame (b/s)? ','s');
      end
      if rolltype == 'b' | rolltype == 's'
         ok = 1;
      end
   end
   clear ok

else
   % Set helpvariable skip = 1, to ensure that the aircraft configuration
   % does not have to be entered if the user chooses option 5 = QUIT.
   % --------------------------------------------------------------------
   skip = 1;
end

%%%%

if skip ~= 1         % DEFINE CONFIGURATION OF THE AIRPLANE, IF NOT QUITTING
                     % -----------------------------------------------------

   % For the 'Beaver' model, the flap angle and engine speed define the
   % configuration of the aircraft (engine speed is selected by the pilot
   % and maintained by the regulator of the propeller). Other aircraft
   % may have other definitions, like gear in/out, slats, settings of trim
   % surfaces, multiple engines, etc., so change the following lines if
   % required!
   % ---------------------------------------------------------------------
   deltaf = input('Give flap angle [deg], default = 0: ')*pi/180;
   if isempty(deltaf)
      deltaf = 0;
   end

   n = input('Give engine speed [RPM], default = 1800: ');
   if isempty(n)
      n = 1800;
   end

   % For the 'Beaver' aircraft, the engine power is determined by the engine
   % speed, which has been set above, and the manifold pressure, which will
   % be involved in the trim process if the flightpath angle gamma is user-
   % specified, or kept constant if pz is user-specified (defined in string-
   % variable gammatype). If the manifold pressure is to be adjusted by the
   % numerical trim algorithm, it is important to specify a meaningful esti-
   % mation of the manifold pressure, because otherwise the trim process may
   % not converge, or may give unrealistic results if the optimization pro-
   % cess converges to a LOCAL minimum.
   % -----------------------------------------------------------------------
   if gammatype == 'f'
      pz = input('Give estimate of manifold pressure pz ["Hg], default = 20: ');
   else  % gammatype == 'm'
      pz = input('Give manifold pressure pz ["Hg], default = 20: ');
   end
   if isempty(pz)
      pz = 20;
   end

%   % Hdot is the rate-of-climb, which is a function of the flightpath angle
%   % (change the program if you want to specify Hdot itself in stead of
%   % the flightpath angle gamma).
%   % ----------------------------------------------------------------------
%   Hdot = V*sin(gamma);

   % G is the centripetal acceleration, used for the coordinated turn
   % constraint.
   % ----------------------------------------------------------------
   G = psidot*V/9.80665;

   % If steady, uncoordinated, turning condition must be determined, the
   % roll angle can be specified freely; equilibrium will be obtained for
   % a trimmed-flight sideslip angle, which in this case usually will be
   % quite large.
   % --------------------------------------------------------------------
   phi = [];
   if opt == 2 & turntype == 'u';  % Steady turn, uncoordinated.
      phi = input('Give desired roll angle phi [deg], default = 0: ')*pi/180;
   end

   if isempty(phi)
      phi = 0;
   end

   % All variables which specify the flight condition will be put in the
   % vector ctrim (constants for trim process). The variables which are
   % varied by the minimization routine FMINS will be put into the vector
   % vtrim (variables for trim process). These vectors will have to be
   % redefined for aircraft which use other state or input vectors.
   %
   % The variable phi in ctrim is used for uncoordinated turns only; if a
   % coordinated turn is required, the coordinated turn constraint determines
   % the values of phi and beta. Type HELP ACCONSTR for more info.
   %
   % Definition:
   %
   % if gammatype =='f' (flight-path angle user-specified and manifold
   % pressure numerically adjusted):
   %    vtrim = [alpha beta deltae deltaa deltar pz]'
   % else
   %    vtrim = [alpha beta deltae deltaa deltar gamma]'
   %
   % Note: vtrim contains the first estimation of the non-constant states
   % and inputs. The final values will be determined iteratively.
   % -----------------------------------------------------------------------
   if gammatype == 'f' % gamma in ctrim, pz in vtrim
      ctrim = [V H psi gamma G psidot thetadot phidot deltaf n phi]';
      vtrim = [0 0 0 0 0 pz]';
   else % gamma in vtrim, pz in ctrim
      ctrim = [V H psi pz G psidot thetadot phidot deltaf n phi]';
      vtrim = [0 0 0 0 0 0]';  % <-- initial guess for gamma = 0!
   end

   % Options for the minimization routine FMINS. Type HELP FMINS or
   % HELP FOPTIONS for more info. Note: opts contains the desired values
   % of the options, the actual values will be returned in the vector
   % options when calling the FMINS routine later on.
   % -------------------------------------------------------------------
   options     = foptions([]);  % Use default options for FMINS...
   options(2)  = 1e-10;         % ... except for the error tolerance, which
   options(3)  = 1e-10;         %     is set to a smaller value.

   % The Simulink system BEAVER or an equivalent aircraft model is called by
   % the function call xdot = feval('beaver',t,x,u,'outputs'), followed by
   % xdot = feval('beaver',t,x,u,'derivs') to obtain the time-derivatives of
   % the state vector x. Here t = 0 (system is time-independent). The
   % function call is created here, and stored in the variable modfun.
   % Note: the aircraft model itself will be compiled before applying these
   % function calls!
   % -----------------------------------------------------------------------
   modfun   = ['xdot = feval(''',sysname,''',0,x,u,''outputs'');'];
   modfun   = [modfun 'xdot = feval(''',sysname,''',0,x,u,''derivs'');'];

   % Iteration loop in which the steady-state flight condition is searched.
   % If the variable iterate is set to 0, the iteration will be stopped.
   % If no stable solution has been found before reaching the maximum num-
   % ber of iterations, a warning message will be displayed and the user
   % will be asked if the routine should continue with the iterations, or
   % stop. Usually a stable solution will be found before reaching the
   % maximum number of iterations.
   % ----------------------------------------------------------------------

   iterate = 1;
   while iterate
      clc;
      disp('Searching for stable solution. Wait a moment...');
      disp(' ');
      disp(' ');
      pause(0.5)
      comment
      home
      pause(0.5)

      % First pre-compile the aircraft model.
      % -------------------------------------
      feval(sysname,[],[],[],'compile');