actrim.m

matlab的FDC工具箱

   %
   % 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

   % 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'');'];

   % First pre-compile the aircraft model (temporarily suppress warnings
   % to prevent harmless, but inconvenient messages when the model is 
   % already compiled).
   %--------------------------------------------------------------------
   warning off
   feval(sysname,[],[],[],'compile');
   clear ans % (the above feval statement somehow generates an 'ans' variable)
   warning on

   % Call minimization routine FMINSEARCH. A scalar cost function is contained
   % in the M-file ACCOST, which also calls the constraints for coordinated
   % turns and rate-of-climb of the aircraft. FMINSEARCH returns the trim values
   % of alpha, beta, deltae, deltaa, deltar, and pz in the vector vtrimmed.
   % The convergence tolerance has been set to 1e-10 (default 1e-4), the 
   % maximum number of function evaluations has been set to 5000 (default 1200),
   % and the maximum number of iterations has been set to 3000 (default 1200)
   % ----------------------------------------------------------------------
   clc;
   disp('Searching for stable solution. Wait a moment...');
   disp(' ');
   
   options = optimset('Display','off','TolX',1e-10,'MaxFunEvals',5000,'MaxIter',3000);
   [vtrimmed,fval,exitflag,output] = fminsearch('accost',vtrim,options,...
      ctrim,rolltype,turntype,gammatype,modfun);
   
   % Display error message when maximum number of iterations is 
   % reached before finding converged solution
   % ----------------------------------------------------------
   if exitflag == 0
      warning('Maximum number of iterations was exceeded!');
      disp(['- number of function evaluations: ' num2str(output.funcCount)]);
      disp(['- number of iterations: ' num2str(output.iterations)]);
   else
      disp('Converged to a trimmed solution...');
   end

   disp(' ');
   disp('<<< Press a key to get results >>>');
   pause

   % Call subroutine ACCONSTR, which contains the flight-path constraints
   % once more to obtain the final values of the inputvector and states.
   % --------------------------------------------------------------------
   [x,u] = acconstr(vtrimmed,ctrim,rolltype,turntype,gammatype);
   
   % Call a/c model once more, to obtain the time-derivatives of the states
   % in trimmed-flight condition. By examining xdot, the user can decide if
   % the trimmed-flight condition is accurate enough. If not, either the
   % error tolerance of the minimization routine needs to be tightened, or
   % the cost function in ACCOST.M needs to be changed.
   % -----------------------------------------------------------------------
   eval(modfun);
   
   % Release compiled aircraft model now that all results are known
   % --------------------------------------------------------------
   feval(sysname,[],[],[],'term');
   
   % Display the final results. x, u, and xdot contain the inputvector,
   % state vector, and time-derivative of the state vector, respectively.
   % --------------------------------------------------------------------
   clc
   disp('State vector (trimmed): ');
   x
   disp(' ');
   disp(' ');
   disp('<<< Press a key >>>');
   pause

   clc
   disp('Input vector (trimmed): ');
   u
   disp(' ');
   disp(' ');
   disp('<<< Press a key >>>');
   pause

   clc
   disp('Time derivative of state vector (trimmed): ');
   xdot
   disp(' ');
   disp(' ');
   disp('<<< Press a key >>>');
   pause

   % The integrator block in the aircraft model uses the Matlab variable
   % xinco to store the initial value of x. The models from the library
   % EXAMPLES use the variables uaero0 and uprop0 to store the initial values
   % of the inputs to the aerodynamic and engine models, respectively.
   % Usually, trimmed-flight conditions are used as initial conditions.
   % Therefore, the trimmed-flight conditions will be stored in these'
   % variables, to ensure compatibility with the Simulink systems. The
   % trimmed-flight value of xdot = dx/dt will be stored in the variable
   % xdot0, which indicates that this value of xdot is obtained for the
   % INITIAL flight condition, defined by xinco, uaero0, and uprop0.
   %-------------------------------------------------------------------
   xinco  = x;          % Initial value of state vector
   xdot0  = xdot;       % Initial value of time-derivative of state vector
   uaero0 = u(1:4);     % Initial value of input vector for aerodynamic model
   uprop0 = u(5:6);     % Initial value of input vector for engine
                        %                            forces and moments model
                        
   % Now, a string-matrix will be created, which contains a description of
   % the trimmed-flight condition. Note: the function NUM2STR2 is just a
   % sligtly changed version of NUM2STR. Type HELP NUM2STR2 at the Matlab
   % command line for more info.
   % ---------------------------------------------------------------------
   line1  = ['Trimmed-flight condition for S-function ' sysname ':'];
   line2 = '';

   if opt == 1
      line3 = ['Steady-state wings-level flight'];
   elseif opt == 2
      line3 = ['Steady-state turning flight'];
   elseif opt == 3
      line3 = ['Steady pull-up'];
   elseif opt == 4
      if rolltype == 'b'                            % Body-axis roll
         line3 = ['Steady roll along the body-axis'];
      else                                          % Stability-axis roll (default)
         line3 = ['Steady roll along the stability-axis'];
      end
   end

   line4  = '';
   line5  = ['User-specified definition of aircraft and flight condition:'];
   line6  = '';
   line7  = ['Flap angle:       deltaf    =  ' num2str2(deltaf,10) ' [deg]  '];
   line8  = ['Engine speed:     n         =  ' num2str2(n,10)      ' [RPM]  '];
   line9  = ['Airspeed:         V         =  ' num2str2(V,10)      ' [m/s]  '];
   line10 = ['Altitude:         H         =  ' num2str2(H,10)      ' [m]    '];
   line11 = ['Yaw-angle:        psi       =  ' num2str2(psi,10)    ' [deg]  '];

   if gammatype == 'f'
      line12 = ['Flightpath angle: gamma     =  ' num2str2(gamma,10)  ' [deg]  '];
   else
      line12 = ['Manifold pressure: pz       =  ' num2str2(pz,10)     ' ["Hg]  '];
   end
   line13 = ['Yaw-rate:         psi dot   =  ' num2str2(psidot,10)   ' [deg/s]'];
   line14 = ['Pitch-rate:       theta dot =  ' num2str2(thetadot,10) ' [deg/s]'];
   line15 = ['Roll-rate:        phi dot   =  ' num2str2(phidot,10)   ' [deg/s]'];

   if opt == 2 & turntype == 'u' % uncoordinated turn
      line16 = ['Uncoordinated turn,'];
      line17 = ['Roll-angle:       phi       =  ' num2str2(phi,10) ' [deg]'];
   elseif opt == 2 & turntype == 'c' % coordinated turn
      line16 = ['Coordinated turn'];
      line17 = '';
   else % no turning flight
      line16 = '';
      line17 = '';
   end

   % Some general remarks, need to be changed manually by the user!
   % --------------------------------------------------------------
   line18 = 'Definitions of state and input vectors as in system BEAVER:';
   line19 = 'x = [V alpha beta p q r psi theta phi xe ye H]''';
   line20 = 'u = [deltae deltaa deltar deltaf n pz uw vw ww uwdot vwdot wwdot]''';
   line21 = '';

   % Add info about the variables in the Matlab Workspace.
   % -----------------------------------------------------
   line22 = 'Variables: xdot=initial state, uaero0=initial inputs to aeromodel';
   line23 = 'uprop0=initial input to engine forces/moments model, xdot0=dx/dt(0)';
   line24 = '';

   % Add time and date to text matrix.
   % ---------------------------------
   t  = clock;
   t1 = num2str(t(4));
   t2 = num2str(t(5));
   t3 = num2str(t(6));

   line25 = ['date: ', date, '   time: ', t1, ':', t2, ':', t3];
   
   % Enhance explanatory lines to 80 characters and collect the results in the 
   % string matrix 'trimdef'. Clear the individual line variables.
   % -------------------------------------------------------------------------
   trimdef = [];
   for ii = 1:25
      linestr = ['line' num2str(ii)];
      eval([linestr '= [' linestr ' blanks(79-length(' linestr '))];']);
      eval(['trimdef = [trimdef; ' linestr '];']);
      eval(['clear ' linestr ]);
   end
   
   % Results can now be saved to a file, which will be filled with the
   % trim condition xinco, uaero0, uprop0, xdot0, and the explanation matrix
   % trimdef.
   %
   % Note: files with steady-state trimmed flight condition have
   % extension .TRI!
   % -----------------------------------------------------------------
   clc
   answ = input('Save trimmed condition to file (y/[n])? ','s');
   if isempty(answ)
      answ = 'n';
   end
   if answ == 'y'

      % Destination directory for storing the results is the default FDC 
      % data-directory, which is obtained from the function DATADIR.M. 
      % If that directory does not exist, the current directory will be 
      % used instead.
      % ----------------------------------------------------------------
      defdir = datadir;
      currentdir = chdir;

      % Go to default directory if that directory exists (if not, start
      % save-routine from the current directory).
      % ---------------------------------------------------------------
      try
        chdir(defdir);
      catch
        chdir(currentdir);
      end

      % Obtain filename and path
      % ------------------------
      [filename,dirname] = uiputfile('*.tri','Save trimmed flight condition');

      % Save results
      % ------------
      save([dirname,filename],'xinco','uaero0','uprop0','xdot0','trimdef');
     
      % Display file location
      % ---------------------
      clc;
      disp('Results saved to the file: ');
      disp(' ');
      disp(['     ',dirname,filename]);
      disp(' ');
      disp('<<< Press a key >>>');
      pause
   end

   % Display user-information.
   % -------------------------
   clc
   disp('The results have been stored in the following variables:');
   disp(' ');
   disp('xinco = [V alpha beta p q r psi theta phi xe ye H]'' = state vector');
   disp('xdot0 = dx/dt(0)');
   disp('uaero0= [deltae deltaa deltar deltaf]'' = initial input vector for');
   disp('                                          aerodynamic model');
   disp('uprop0= [n pz]'' = initial input vector for engine model');
   disp(' ');
   disp('The text-matrix ''trimdef'' contains more info about the trimmed');
   disp('flightcondition.');
   disp(' ');
   
   % Repeat warning message if the maximum number of iterations was 
   % reached before finding a converged solution
   % --------------------------------------------------------------
   if exitflag == 0
   disp(' ');
      disp('WARNING!');
      disp('  The solution may not be accurate, as the maximum number of iterations'); 
      disp('  was reached before finding a converged solution. You can increase the');
      disp('  number of iterations by changing the MaxIter parameter of the simulation');
      disp('  options. This involves editing actrim.m!');
   end

end                % END OF TRIM-LOOP (SKIPPED IF OPTION 5 = QUIT IS SELECTED)


% Clear unneccessary variables
% ----------------------------
clear h ctrim vtrim vtrimmed modfun output options exitflag fval
clear opt V H psi gamma phidot psidot thetadot rolltype turntype gammatype
clear deltaf n pz Hdot G phi sysname x xdot u exitflag t t1 t2 t3 ii linestr
clear ok answ dirname filename savecmmnd defdir currentdir skip

disp(' ');
disp('Ready.');
disp(' ');


%-----------------------------------------------------------------------
% 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 file: September 7, 2005.  
%-----------------------------------------------------------------------