openloop.m

matlab的FDC工具箱

% The FDC toolbox. Force/moment coefficient plotting routine PLOTCOEFF.
% =====================================================================
% OPENLOOP is a Matlab macro that determines nonlinear and linear open-loop
% responses of the Beaver for six specific test-inputs. The results were 
% compared with ref.[1] to verify proper implementation of the Beaver 
% dynamic model in Simulink and to verify the results from the linearization
% process by the ACLIN tool. 
%
% The test inputs are:
%
% - block-shaped elevator, aileron, and rudder deflections, maintaining
%   a deflection of 3-degrees during 2 seconds
% - a ramp-shaped flap extension from zero to 3 degrees in 3 seconds
% - a ramp-shaped rise in engine RPM, from plus zero to plus 200 RPM in
%   4 seconds
% - a ramp-shaped rise in manifold pressure, from plus zero to plus 2
%   inch Hg in 2 seconds
%
% These test inputs are identical to the default test inputs used by the
% systems oloop1.mdl and oloop3.mdl. Note: oloop1.mdl and oloop3.mdl 
% provide interactive access to the nonlinear and linearized Beaver model,
% whereas this macro uses batch processing to achieve similar results 
% OPENLOOP uses the Simulink models oloop1a.mdl and oloop3a.mdl, which 
% correspond to oloop1.mdl and oloop3.mdl, except in this case the test 
% inputs are defined in this Matlab script i.s.o. the Simulink models.
%
% All results are based on an initial steady flight-condition at an 
% altitude of 6000 feet (approx. 1850 m) and an airspeed of 45 m/s.
% Nonlinear results shown as a solid line, linear results are shown as
% a dotted line. 
%
% REFERENCE:
% [1] Baarspul, M.: Lecture Notes on Flight Simulation Techniques. 
%     Report LR-596, Delft University of Technology, Faculty of 
%     Aerospace Engineering, Delft, 1989.
% -----------------------------------------------------------------------

% Variables:
% ==========
% Delta_deltaa    block-shaped change in aileron deflection, equal to u_block
% Delta_deltae    block-shaped change in elevator deflection, equal to u_block
% Delta_deltaf    ramp-shaped change in flap deflection (3 degrees in 3 seconds)
% Delta_deltar    block-shaped change in rudder deflection, equal to u_block
% Delta_n         ramp-shaped change in engine RPM (200 RPM in 4 seconds)
% Delta_pz        ramp-shaped change in manifold pressure (2 inch Hg in 2 sec.)

% h               handle to waitbar figure
% h1 ... h6       handles to the six results plots
% h_fig           height of results plots
% hh              handle to current result plot (equals h1, h2, ..., or h6)
% ii              counter
% OldAxes...      used to temporarily store default axes properties
% screen_size     used to store screen dimensions
% time            time-axis
% u_block         general block-signal, equal to 3 degrees (converted to rad) 
%                  during 2 seconds
% u0              vector with zeros, used to fill up u1 ... u6 input matrices
% u1 ... u6       matrices containing the six individual control input 
%                  trajectories
% x_fig           X-coordinate of lower left corner of results plots
% y_fig           Y-coordinate of lower left corner of results plots
% w_fig           width of results plots
% x1 ... x6       matrices with nonlinear state trajectories for the six test-cases
%                  (not used here, as these results are duplicated in y1 ... y6)
% x1lin ... x6lin matrices with linear state trajectories for the six test-cases
%                  (not used here, as these results are duplicated in y1lin ... y6lin)
% y1 ... y6       matrices with output trajectories (including states) for the 
%                  six test-cases
% yy              current output matrix of nonlinear results (equals y1, y2, ..., 
%                  or y6)
% yylin           current output matrix of linear results (equals y1in, y2lin, ..., 
%                  or y6lin)


% Load model parameters
% ---------------------
try
   load('aircraft.dat','-mat');
   load('cr4560.lin','-mat');
catch
   error('Failed to load the model parameters from file.');
end
          

% Define time-axis (using steps of 0.1 seconds results in 1201 data-points)
% -------------------------------------------------------------------------
time = 0:0.1:120;

% Define block-shaped signal for elevator, aileron, and rudder deflections
% (deflection of 3 degrees, which lasts 2 seconds; ten data points per second,
% leaves 1201-20 = 1181 data points after the step signal)
% ----------------------------------------------------------------------------
u_block = 3*pi/180 * [ones(1,20) zeros(1,1181)];
Delta_deltae = u_block;
Delta_deltaa = u_block;
Delta_deltar = u_block;

% Define ramp-shaped signal for flap extension (reaching a deflection of 
% 3 degrees in 3 seconds; ten data points per second, plus one additional
% point, leaves 1201-30-1 = 1170 data points after the ramp signal)
% ----------------------------------------------------------------------
Delta_deltaf = pi/180 * [[0:0.1:3] 3*ones(1,1170)];

% Define ramp-shaped signal for engine RPM increase (reaching an RPM increase 
% of 200 RPM in 4 seconds; ten data points per second, plus one additional 
% point, leaves 1201-40-1 = 1160 data points after the ramp signal)
% ---------------------------------------------------------------------------
Delta_n = 50*[[0:0.1:4] 4*ones(1,1160)];

% Define ramp-shaped signal for engine manifold pressure increase (reaching 
% an additional 2 inch Hg in 2 seconds; ten points per second, plus one 
% additional point, leaves 1201-20-1 = 1180 data points after the ramp signal)
% ----------------------------------------------------------------------------
Delta_pz = [[0:0.1:2] 2*ones(1,1180)];

% Define a vector of zeros with the same length as the vector defined above
% -------------------------------------------------------------------------
u0 = zeros(1,length(time));

% Define input matrices corresponding with the individual control input 
% trajectories defined above
% ---------------------------------------------------------------------
u1 = [time; Delta_deltae; u0; u0; u0; u0; u0]';
u2 = [time; u0; Delta_deltaa; u0; u0; u0; u0]';
u3 = [time; u0; u0; Delta_deltar; u0; u0; u0]';
u4 = [time; u0; u0; u0; Delta_deltaf; u0; u0]';
u5 = [time; u0; u0; u0; u0; Delta_n; u0]';
u6 = [time; u0; u0; u0; u0; u0; Delta_pz]';

% Run simulations to find open-loop aircraft responses to the individual
% control inputs, and show a waitbar to indicate progress
% ----------------------------------------------------------------------
h = waitbar(0,'Please wait...');
set(h,'Name','Computing open-loop responses');
waitbar(0,h);
[time,x1,y1]=sim('oloop1a',time,[],u1);
[tdummy,x1lin,y1lin]=sim('oloop3a',time,[],u1);
waitbar(17/100,h);
[time,x2,y2]=sim('oloop1a',time,[],u2);
[time,x2lin,y2lin]=sim('oloop3a',time,[],u2);
waitbar(33/100,h);
[time,x3,y3]=sim('oloop1a',time,[],u3);
[time,x3lin,y3lin]=sim('oloop3a',time,[],u3);
waitbar(50/100,h);
[time,x4,y4]=sim('oloop1a',time,[],u4);
[time,x4lin,y4lin]=sim('oloop3a',time,[],u4);
waitbar(67/100,h);
[time,x5,y5]=sim('oloop1a',time,[],u5);
[time,x5lin,y5lin]=sim('oloop3a',time,[],u5);
waitbar(83/100,h);
[time,x6,y6]=sim('oloop1a',time,[],u6);
[time,x6lin,y6lin]=sim('oloop3a',time,[],u6);
waitbar(100/100,h);
pause(0.5); close(h);

% Determine suitable figure dimensions
% ------------------------------------
screen_size = screensize;
x_fig = screen_size(1)+50;
y_fig = screen_size(2)+150;
w_fig = screen_size(3)-250;
h_fig = screen_size(4)-250;

% Temporarily store default axes properties (will be changed for 
% plotting, but restored later)
% --------------------------------------------------------------
OldAxesFontSize       = get(0,'DefaultAxesFontSize');
OldAxesPosition       = get(0,'DefaultAxesPosition');
OldAxesColorOrder     = get(0,'DefaultAxesColor');
OldAxesLineStyleOrder = get(0,'DefaultAxesLineStyleOrder');

% Set fontsize for axes
% ---------------------
set(0,'DefaultAxesFontSize',8);
set(0,'DefaultAxesColorOrder',[0 0 0]);
set(0,'DefaultAxesLineStyleOrder',['-|:|--']);

% Initiate six figures to plot the aircraft responses to the six test-signals
% ---------------------------------------------------------------------------
h1=figure('Name','Responses to a block-shaped elevator deflection', ...
     'NumberTitle','off', ...
     'Position',[x_fig,y_fig,w_fig,h_fig]);
h2=figure('Name','Responses to a block-shaped aileron deflection', ...
     'NumberTitle','off', ...
     'Position',[x_fig+20,y_fig-20,w_fig,h_fig]);
h3=figure('Name','Responses to a block-shaped rudder deflection', ...
     'NumberTitle','off', ...
     'Position',[x_fig+40,y_fig-40,w_fig,h_fig]);
h4=figure('Name','Responses to a ramp-shaped flap extension', ...
     'NumberTitle','off', ...
     'Position',[x_fig+60,y_fig-60,w_fig,h_fig]);
h5=figure('Name','Responses to a ramp-shaped increase in engine RPM', ...
     'NumberTitle','off', ...
     'Position',[x_fig+80,y_fig-80,w_fig,h_fig]);
h6=figure('Name','Responses to a ramp-shaped increase in manifold pressure', ...
     'NumberTitle','off', ...
     'Position',[x_fig+100,y_fig-100,w_fig,h_fig]);

% Plot all relevant output trajectories and the relevant test-sigal for the
% six test-cases, using the figures defined above
% -------------------------------------------------------------------------
for ii = 1:6
  
   % Determine ii-th output matrix and figure handle
   % -----------------------------------------------
   yy    = eval(['y' num2str(ii)]);
   yylin = eval(['y'num2str(ii) 'lin']);
   hh    = eval(['h' num2str(ii)]);

   % Open ii-th figure and plot the relevant outputs
   % -----------------------------------------------
   figure(hh);
   
   set(0,'DefaultAxesPosition',[0.08 0.06 0.93 0.88]);

   subplot(6,2,1);  plot(time,yy(:,1),time,yylin(:,1));
   ylabel('V [m s^{-1}]','FontAngle','italic');
   subplot(6,2,3);  plot(time,yy(:,2)*180/pi,time,yylin(:,2)*180/pi);
   ylabel('\alpha [deg]','FontAngle','italic');
   subplot(6,2,5);  plot(time,yy(:,3)*180/pi,time,yylin(:,3)*180/pi);
   ylabel('\beta [deg]','FontAngle','italic');
   subplot(6,2,7);  plot(time,yy(:,4)*180/pi,time,yylin(:,4)*180/pi);
   ylabel('p [deg s^{-1}]','FontAngle','italic');
   subplot(6,2,9);  plot(time,yy(:,5)*180/pi,time,yylin(:,5)*180/pi);
   ylabel('q [deg s^{-1}]','FontAngle','italic');
   subplot(6,2,11); plot(time,yy(:,6)*180/pi,time,yylin(:,6)*180/pi);
   xlabel('t [s]','FontAngle','italic'); ylabel('r [deg s^{-1}]','FontAngle','italic');

   set(0,'DefaultAxesPosition',[0.02 0.06 0.93 0.88]);

   subplot(6,2,2);  plot(time,yy(:,7)*180/pi,time,yylin(:,7)*180/pi);
   ylabel('\psi [deg]','FontAngle','italic');
   subplot(6,2,4);  plot(time,yy(:,8)*180/pi,time,yylin(:,8)*180/pi);
   ylabel('\theta [deg]','FontAngle','italic');
   subplot(6,2,6);  plot(time,yy(:,9)*180/pi,time,yylin(:,9)*180/pi);
   ylabel('\phi [deg]','FontAngle','italic');
   subplot(6,2,8);  plot(time,yy(:,12),time,yylin(:,12));
   ylabel('H [m]','FontAngle','italic');
   subplot(6,2,10); plot(time,yy(:,13),time,yylin(:,13));
   ylabel('dH/dt [ms^{-1}]','FontAngle','italic');
end   

% Include the relevant test-inputs for the six test-cases
% -------------------------------------------------------
figure(h1); subplot(6,2,12); plot(time,(Delta_deltae+uaero0(1))*180/pi);
   ylabel('\delta_e [deg]','FontAngle','italic');
figure(h2); subplot(6,2,12); plot(time,(Delta_deltaa+uaero0(2))*180/pi);
   ylabel('\delta_a [deg]','FontAngle','italic');
figure(h3); subplot(6,2,12); plot(time,(Delta_deltar+uaero0(3))*180/pi);
   ylabel('\delta_r [deg]','FontAngle','italic');
figure(h4); subplot(6,2,12); plot(time,(Delta_deltaf+uaero0(4))*180/pi);
   ylabel('\delta_f [deg]','FontAngle','italic');
figure(h5); subplot(6,2,12); plot(time,(Delta_n+uprop0(1)+.01)); % fudge-factor for better plot-scaling 
   ylabel('n [RPM]','FontAngle','italic');
figure(h6); subplot(6,2,12); plot(time,(Delta_pz+uprop0(2)));
   ylabel('p_z ["Hg]','FontAngle','italic');
   
% Display a message box with a legend (this does not distort the plots,
% contrary to using the 'legend' command)
% ---------------------------------------------------------------------
newMsgBox({'Solid lines: results from nonlinear simulations', ...
           'Dotted lines: results from linear simulations', ...
           ' ', ...
           'Relevant test inputs have been shown in the lower-left', ...
           'subplots, and mentioned in the figure titles.'},'Plot legend');

% Restore old default axes properties
% -----------------------------------
set(0,'DefaultAxesFontSize',OldAxesFontSize);
set(0,'DefaultAxesPosition',OldAxesPosition);
set(0,'DefaultAxesColorOrder',OldAxesColorOrder);
set(0,'DefaultAxesLineStyleOrder',OldAxesLineStyleOrder);

%-----------------------------------------------------------------------
% 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: May 19, 2005.
%-----------------------------------------------------------------------