Code covered by the BSD License  

Highlights from
Automated Failure Boundary Mapping

image thumbnail

Automated Failure Boundary Mapping

by

 

18 Aug 2009 (Updated )

Demo files from July 21, 2009 webinar

Flutter Model Demo

Flutter Model Demo

From: "Using Statistics and Optimization to Support Design Activities" Webinar, July 21, 2009.

This demo shows how to use the boundary mapping algorithm developed in boundaryMapping and encapsulated in the function boundarySearch on a Simulink® Model of wing flutter. This demo also highlights two approaches to sensitivity analysis and using optimization to minimize the flutter failure region.

NOTE: This demo uses Simulink® for the modeling environment, but the techniques shown here can be applied to any model that can be called from a MATLAB® function.

The Flutter Analysis Model was originally developed by Cole Stephens and can be found on MATLAB® Central.

(C) Copyright 2009 The MathWorks, Inc.

Contents

1-Failure Boundary Mapping

The goal of this analysis is to map the failure boundary (Flutter Region) using boundarySearch. First let's define the model and parameters that we want to map the failure region over (Mach and Altitude).

setup
% Simulink model and blocks
model  = 'Flutter';
blocks = {'Mach', 'Altitude'};

% Create DataSet Array extracting parameters from blocks (Statistics
% Toolbox(TM) required here).
ds = getParamNames(model,blocks)
ds = 

               Mach    Altitude
    Initial    0.77    30000   

1.1-Define Bounds

Define the range of values that we are interested in.

NOTE: The Simulink Model has been verified for Mach Number of 0.77 and Altitude of 30,000 ft. The aerodynamic model defined in the Simulink model may not be appropriate for other Mach Number and Altitude ranges. We are assuming for the purpose of this example that they are. If you plan to use this Simulink model, verification of the aerodynamic model should be performed. The goal of this analysis is to show how to perform a failure boundary mapping.

ds.Mach('lowerBound',1)= 0.65;
ds.Mach('upperBound',1) = 0.85;
ds.Altitude('lowerBound',1) = 10000;
ds.Altitude('upperBound',1) = 50000
ds = 

                  Mach    Altitude
    Initial       0.77    30000   
    lowerBound    0.65    10000   
    upperBound    0.85    50000   

1.2-Run Model Over Bounds

Open Simulink Model and run it over the bounds and initial setting. Flutter occurs when the plunge (up and down motion) and pitch (twist of airfoil) frequencies align. The response recorded from the model is a true (1) or false (0) indication that flutter occurs (unstable). Also recorded in the dataset array are the plunge and pitch frequencies, as well as the time series data.

open(model)
ds = runModel(ds)
ds = 

                  Mach    Altitude    unstable    time                    pitch               
    Initial       0.77    30000       1           [ 3417x1 double]        [ 3417x1 double]    
    lowerBound    0.65    10000       0           [10370x1 double]        [10370x1 double]    
    upperBound    0.85    50000       0           [10010x1 double]        [10010x1 double]    


                  plunge                  pitchFreq    plungeFreq
    Initial       [ 3417x1 double]        3.5494        3.6718   
    lowerBound    [10370x1 double]        3.2593       0.31644   
    upperBound    [10010x1 double]        4.3985        4.3985   

Let's take a look at the responses for the initial setting. It is clear that flutter is occuring (unstable oscillations as the unstable flag inidicated in ds.unstable).

make2plots(ds.time{'Initial'},ds.pitch{'Initial'},ds.plunge{'Initial'})

1.3-Design Exploration Via Design of Experiments

Let's seed our failure boundary search with a set of points that span the range of interests. We'll use Design of Experiements for this, but other approaches could be used. DOE will return a nondimensional test matrix.

factors = {'Mach','Altitude'};
doe = ccdesign(length(factors),'center',1,'type','faced')

close all %figures
doe =

    -1    -1
    -1     1
     1    -1
     1     1
    -1     0
     1     0
     0    -1
     0     1
     0     0

1.4-Scale to Engineering Units

Scale the DOE to work with the units in the Simulink Model

doeRuns    = length(doe);
doeRuns    = strrep( strcat({'DOE'},num2str((1:doeRuns)','%d')), ' ','');
doeRange   = [min(doe(:,1)) max(doe(:,1))];

for i = 1:length(factors)
    range = ds.(factors{i})({'lowerBound','upperBound'});
    ds.(factors{i})(doeRuns) = interp1(doeRange, range, doe(:,i));
end
ds(:,factors)

plot(ds.Mach(doeRuns),ds.Altitude(doeRuns),'bo','MarkerFaceColor','b')
axis tight
xlabel('Mach Number')
ylabel('Altitude (feet)')
ans = 

                  Mach    Altitude
    Initial       0.77    30000   
    lowerBound    0.65    10000   
    upperBound    0.85    50000   
    DOE1          0.65    10000   
    DOE2          0.65    50000   
    DOE3          0.85    10000   
    DOE4          0.85    50000   
    DOE5          0.65    30000   
    DOE6          0.85    30000   
    DOE7          0.75    10000   
    DOE8          0.75    50000   
    DOE9          0.75    30000   

1.5-Run Simulink at DOE Points

ds(doeRuns,:) = runModel(ds(doeRuns,:))
ds = 

                  Mach    Altitude    unstable    time                    pitch               
    Initial       0.77    30000       1           [ 3417x1 double]        [ 3417x1 double]    
    lowerBound    0.65    10000       0           [10370x1 double]        [10370x1 double]    
    upperBound    0.85    50000       0           [10010x1 double]        [10010x1 double]    
    DOE1          0.65    10000       0           [10370x1 double]        [10370x1 double]    
    DOE2          0.65    50000       0           [10010x1 double]        [10010x1 double]    
    DOE3          0.85    10000       0           [10098x1 double]        [10098x1 double]    
    DOE4          0.85    50000       0           [10010x1 double]        [10010x1 double]    
    DOE5          0.65    30000       1           [ 3105x1 double]        [ 3105x1 double]    
    DOE6          0.85    30000       1           [ 3875x1 double]        [ 3875x1 double]    
    DOE7          0.75    10000       0           [10141x1 double]        [10141x1 double]    
    DOE8          0.75    50000       0           [10010x1 double]        [10010x1 double]    
    DOE9          0.75    30000       1           [ 3368x1 double]        [ 3368x1 double]    


                  plunge                  pitchFreq    plungeFreq
    Initial       [ 3417x1 double]        3.5494        3.6718   
    lowerBound    [10370x1 double]        3.2593       0.31644   
    upperBound    [10010x1 double]        4.3985        4.3985   
    DOE1          [10370x1 double]        3.2593       0.31644   
    DOE2          [10010x1 double]          3.36         4.765   
    DOE3          [10098x1 double]        3.2046        18.334   
    DOE4          [10010x1 double]        4.3985        4.3985   
    DOE5          [ 3105x1 double]        3.9176        3.9176   
    DOE6          [ 3875x1 double]        3.4259        3.4259   
    DOE7          [10141x1 double]         3.342        4.7655   
    DOE8          [10010x1 double]        4.5818        4.5818   
    DOE9          [ 3368x1 double]        3.6719        3.6719   

1.6-Find Failure Boundary Using boundarySearch

We need to pass in values of Mach Number and Altitude that have similar magnitudes, otherwise the triangulation methods used in boundarySearch will run into numeric (scaling) problems.

X0 = double(ds(doeRuns,factors));
lb = double(ds('lowerBound',factors));
ub = double(ds('upperBound',factors));

X0(:,1) = X0(:,1).*10;    % Mach Scaling
X0(:,2) = X0(:,2)./10000; % Altitude Scaling
lb = lb.*[10 1/10000];
ub = ub.*[10 1/10000];

close_system(model,false)   % close model (will be run in the background)
close all                   % close figures

% Define the failure function (often called the limit or limit state
% function) and map the bounary.
limitFunction = @(x)flutterLimitFunction(x,factors);
xb = boundarySearch(X0,lb,ub,limitFunction,'showplot','ptol',0.15);

We can see that the algorithm found the boundary, and the cluster of points tested are predominately located near the boundary region. The length of xb indicates the number of points used to define the boundary with the change in perimeter of <= 0.15 (ptol argument).

length(xb)
xb
ans =

     2


xb = 

    [14x2 double]    [5x2 double]

2-Sensitivity Analysis

Let's take a look at the sensitivity of the the flutter boundary relative to the design parameters we can control (in the BACT Wing & PAPA Mount block).

% Define model and parameters of interest
clear all, close all, close_system(find_system,false)
% Simulink model and blocks
model  = 'Flutter';
blocks = 'BACT Wing & PAPA Mount';

% Create DataSet Array extracting parameters from blocks (and masks)
ds = getParamNames(model,blocks)

% Open block to show parameters
open_system([model,'/',blocks])
ds = 

               m       iyy    Kp      zetaP      Kt      zetaT     xcg       
    Initial    6.08    2.8    2686    0.35782    3000    0.1833    -0.0023275

2.1-Define Allowable Bounds

We'll set the bounds to +/- 10% of the initial value

ds('upperBound',:) = datasetfun(@(x) 1.1*x, ds('Initial',:), 'DatasetOutput', true);
ds('lowerBound',:) = datasetfun(@(x) 0.1*x, ds('Initial',:), 'DatasetOutput', true);

% Swap any bounds that might have been negative
swapIndx = double(ds('upperBound',:)) < double(ds('lowerBound',:));
if sum(swapIndx) ~= 0
    tmp = ds('upperBound',swapIndx);
    ds('upperBound',swapIndx) = ds('lowerBound',swapIndx);
    ds('lowerBound',swapIndx) = tmp;
end

2.2-Run Simulink with Given Dataset and Bounds (Initial Run)

close_system(model,false) % Close model before running
ds = runModel(ds)
ds = 

                  m        iyy     Kp        zetaP       Kt      zetaT      xcg            unstable
    Initial        6.08     2.8      2686     0.35782    3000     0.1833     -0.0023275    1       
    upperBound    6.688    3.08    2954.6      0.3936    3300    0.20163    -0.00023275    1       
    lowerBound    0.608    0.28     268.6    0.035782     300    0.01833     -0.0025603    0       


                  time                    pitch                   plunge                  pitchFreq
    Initial       [ 3417x1 double]        [ 3417x1 double]        [ 3417x1 double]        3.5494   
    upperBound    [ 3336x1 double]        [ 3336x1 double]        [ 3336x1 double]         3.672   
    lowerBound    [10523x1 double]        [10523x1 double]        [10523x1 double]        3.0505   


                  plungeFreq
    Initial       3.6718    
    upperBound     3.672    
    lowerBound    35.611    

2.3-Define Our Factors and Responses of Interest

Since flutter occurs when the plunge and pitch occur at the same frequency, we'll perform the sensitivity analysis around the pitch and plunge frequency. Alternatively we could perform the sensitivity analysis on the failure (flutter) region area, but this would require several Simulink Runs.

factors = {'m','iyy','Kp','zetaP','Kt','zetaT','xcg'};
response= {'pitchFreq','plungeFreq'};

2.4-Sensitivity Analysis by Gradient (Local Method)

Gradient estimation is often used to estimate the sensitivity of a given response to the input factors. This is a local measure of sensitivity since it does not account for how the response may change when we move away from the current point. We'll use central finite differences to estimate the gradient (uses Optimization Toolbox™) around the current factor settings. We'll also set the change in the parameters for the finite difference approximation to be 1% (default values are on the order of 1e-6).

[gradf,ds] = saGradientMethod(ds,factors,response,0.1);
factors
gradf = gradf'
factors = 

    'm'    'iyy'    'Kp'    'zetaP'    'Kt'    'zetaT'    'xcg'


gradf =

   -0.2017   -0.2194    0.0002         0    0.0006    0.0001         0

Note that GRADF is the gradient (or slope) at the initial point and is valid only near the initial point.

ds('Initial',:)
ans = 

               m       iyy    Kp      zetaP      Kt      zetaT     xcg           unstable
    Initial    6.08    2.8    2686    0.35782    3000    0.1833    -0.0023275    1       


               time                   pitch                  plunge                 pitchFreq
    Initial    [3417x1 double]        [3417x1 double]        [3417x1 double]        3.5494   


               plungeFreq
    Initial    3.6718    

2.5-Sensitivity Analysis by Response Surfaces (more Global Method)

We can extend the gradient based method by developing a model over the bounds of interest that approximates the response surface of the true model. Design of Experiments and response surface modeling (quadratic modeling) is a common approach and we'll demonstrate that here.

doe = ccdesign(length(factors),'center',1,'type','faced')
doe =

    -1    -1    -1    -1    -1    -1     1
    -1    -1    -1    -1    -1     1    -1
    -1    -1    -1    -1     1    -1    -1
    -1    -1    -1    -1     1     1     1
    -1    -1    -1     1    -1    -1    -1
    -1    -1    -1     1    -1     1     1
    -1    -1    -1     1     1    -1     1
    -1    -1    -1     1     1     1    -1
    -1    -1     1    -1    -1    -1    -1
    -1    -1     1    -1    -1     1     1
    -1    -1     1    -1     1    -1     1
    -1    -1     1    -1     1     1    -1
    -1    -1     1     1    -1    -1     1
    -1    -1     1     1    -1     1    -1
    -1    -1     1     1     1    -1    -1
    -1    -1     1     1     1     1     1
    -1     1    -1    -1    -1    -1    -1
    -1     1    -1    -1    -1     1     1
    -1     1    -1    -1     1    -1     1
    -1     1    -1    -1     1     1    -1
    -1     1    -1     1    -1    -1     1
    -1     1    -1     1    -1     1    -1
    -1     1    -1     1     1    -1    -1
    -1     1    -1     1     1     1     1
    -1     1     1    -1    -1    -1     1
    -1     1     1    -1    -1     1    -1
    -1     1     1    -1     1    -1    -1
    -1     1     1    -1     1     1     1
    -1     1     1     1    -1    -1    -1
    -1     1     1     1    -1     1     1
    -1     1     1     1     1    -1     1
    -1     1     1     1     1     1    -1
     1    -1    -1    -1    -1    -1    -1
     1    -1    -1    -1    -1     1     1
     1    -1    -1    -1     1    -1     1
     1    -1    -1    -1     1     1    -1
     1    -1    -1     1    -1    -1     1
     1    -1    -1     1    -1     1    -1
     1    -1    -1     1     1    -1    -1
     1    -1    -1     1     1     1     1
     1    -1     1    -1    -1    -1     1
     1    -1     1    -1    -1     1    -1
     1    -1     1    -1     1    -1    -1
     1    -1     1    -1     1     1     1
     1    -1     1     1    -1    -1    -1
     1    -1     1     1    -1     1     1
     1    -1     1     1     1    -1     1
     1    -1     1     1     1     1    -1
     1     1    -1    -1    -1    -1     1
     1     1    -1    -1    -1     1    -1
     1     1    -1    -1     1    -1    -1
     1     1    -1    -1     1     1     1
     1     1    -1     1    -1    -1    -1
     1     1    -1     1    -1     1     1
     1     1    -1     1     1    -1     1
     1     1    -1     1     1     1    -1
     1     1     1    -1    -1    -1    -1
     1     1     1    -1    -1     1     1
     1     1     1    -1     1    -1     1
     1     1     1    -1     1     1    -1
     1     1     1     1    -1    -1     1
     1     1     1     1    -1     1    -1
     1     1     1     1     1    -1    -1
     1     1     1     1     1     1     1
    -1     0     0     0     0     0     0
     1     0     0     0     0     0     0
     0    -1     0     0     0     0     0
     0     1     0     0     0     0     0
     0     0    -1     0     0     0     0
     0     0     1     0     0     0     0
     0     0     0    -1     0     0     0
     0     0     0     1     0     0     0
     0     0     0     0    -1     0     0
     0     0     0     0     1     0     0
     0     0     0     0     0    -1     0
     0     0     0     0     0     1     0
     0     0     0     0     0     0    -1
     0     0     0     0     0     0     1
     0     0     0     0     0     0     0

2.6-Scale DOE to Engineering Units

doeRuns    = length(doe);
doeRuns    = strrep( strcat({'DOE'},num2str((1:doeRuns)','%d')), ' ','');
doeRange   = [min(doe(:,1)) max(doe(:,1))];
for i = 1:length(factors)
    range = ds.(factors{i})({'lowerBound','upperBound'});
    ds.(factors{i})(doeRuns) = interp1(doeRange, range, doe(:,i));
end
ds(:,factors)
ans = 

                  m        iyy     Kp        zetaP       Kt      zetaT      xcg        
    Initial        6.08     2.8      2686     0.35782    3000     0.1833     -0.0023275
    upperBound    6.688    3.08    2954.6      0.3936    3300    0.20163    -0.00023275
    lowerBound    0.608    0.28     268.6    0.035782     300    0.01833     -0.0025603
    deltaX11      5.472     2.8      2686     0.35782    3000     0.1833     -0.0023275
    deltaX12      6.688     2.8      2686     0.35782    3000     0.1833     -0.0023275
    deltaX21       6.08    2.52      2686     0.35782    3000     0.1833     -0.0023275
    deltaX22       6.08    3.08      2686     0.35782    3000     0.1833     -0.0023275
    deltaX31       6.08     2.8    2417.4     0.35782    3000     0.1833     -0.0023275
    deltaX32       6.08     2.8    2954.6     0.35782    3000     0.1833     -0.0023275
    deltaX41       6.08     2.8      2686     0.32204    3000     0.1833     -0.0023275
    deltaX42       6.08     2.8      2686      0.3936    3000     0.1833     -0.0023275
    deltaX51       6.08     2.8      2686     0.35782    2700     0.1833     -0.0023275
    deltaX52       6.08     2.8      2686     0.35782    3300     0.1833     -0.0023275
    deltaX61       6.08     2.8      2686     0.35782    3000    0.16497     -0.0023275
    deltaX62       6.08     2.8      2686     0.35782    3000    0.20163     -0.0023275
    deltaX71       6.08     2.8      2686     0.35782    3000     0.1833     -0.0023336
    deltaX72       6.08     2.8      2686     0.35782    3000     0.1833     -0.0023214
    DOE1          0.608    0.28     268.6    0.035782     300    0.01833    -0.00023275
    DOE2          0.608    0.28     268.6    0.035782     300    0.20163     -0.0025603
    DOE3          0.608    0.28     268.6    0.035782    3300    0.01833     -0.0025603
    DOE4          0.608    0.28     268.6    0.035782    3300    0.20163    -0.00023275
    DOE5          0.608    0.28     268.6      0.3936     300    0.01833     -0.0025603
    DOE6          0.608    0.28     268.6      0.3936     300    0.20163    -0.00023275
    DOE7          0.608    0.28     268.6      0.3936    3300    0.01833    -0.00023275
    DOE8          0.608    0.28     268.6      0.3936    3300    0.20163     -0.0025603
    DOE9          0.608    0.28    2954.6    0.035782     300    0.01833     -0.0025603
    DOE10         0.608    0.28    2954.6    0.035782     300    0.20163    -0.00023275
    DOE11         0.608    0.28    2954.6    0.035782    3300    0.01833    -0.00023275
    DOE12         0.608    0.28    2954.6    0.035782    3300    0.20163     -0.0025603
    DOE13         0.608    0.28    2954.6      0.3936     300    0.01833    -0.00023275
    DOE14         0.608    0.28    2954.6      0.3936     300    0.20163     -0.0025603
    DOE15         0.608    0.28    2954.6      0.3936    3300    0.01833     -0.0025603
    DOE16         0.608    0.28    2954.6      0.3936    3300    0.20163    -0.00023275
    DOE17         0.608    3.08     268.6    0.035782     300    0.01833     -0.0025603
    DOE18         0.608    3.08     268.6    0.035782     300    0.20163    -0.00023275
    DOE19         0.608    3.08     268.6    0.035782    3300    0.01833    -0.00023275
    DOE20         0.608    3.08     268.6    0.035782    3300    0.20163     -0.0025603
    DOE21         0.608    3.08     268.6      0.3936     300    0.01833    -0.00023275
    DOE22         0.608    3.08     268.6      0.3936     300    0.20163     -0.0025603
    DOE23         0.608    3.08     268.6      0.3936    3300    0.01833     -0.0025603
    DOE24         0.608    3.08     268.6      0.3936    3300    0.20163    -0.00023275
    DOE25         0.608    3.08    2954.6    0.035782     300    0.01833    -0.00023275
    DOE26         0.608    3.08    2954.6    0.035782     300    0.20163     -0.0025603
    DOE27         0.608    3.08    2954.6    0.035782    3300    0.01833     -0.0025603
    DOE28         0.608    3.08    2954.6    0.035782    3300    0.20163    -0.00023275
    DOE29         0.608    3.08    2954.6      0.3936     300    0.01833     -0.0025603
    DOE30         0.608    3.08    2954.6      0.3936     300    0.20163    -0.00023275
    DOE31         0.608    3.08    2954.6      0.3936    3300    0.01833    -0.00023275
    DOE32         0.608    3.08    2954.6      0.3936    3300    0.20163     -0.0025603
    DOE33         6.688    0.28     268.6    0.035782     300    0.01833     -0.0025603
    DOE34         6.688    0.28     268.6    0.035782     300    0.20163    -0.00023275
    DOE35         6.688    0.28     268.6    0.035782    3300    0.01833    -0.00023275
    DOE36         6.688    0.28     268.6    0.035782    3300    0.20163     -0.0025603
    DOE37         6.688    0.28     268.6      0.3936     300    0.01833    -0.00023275
    DOE38         6.688    0.28     268.6      0.3936     300    0.20163     -0.0025603
    DOE39         6.688    0.28     268.6      0.3936    3300    0.01833     -0.0025603
    DOE40         6.688    0.28     268.6      0.3936    3300    0.20163    -0.00023275
    DOE41         6.688    0.28    2954.6    0.035782     300    0.01833    -0.00023275
    DOE42         6.688    0.28    2954.6    0.035782     300    0.20163     -0.0025603
    DOE43         6.688    0.28    2954.6    0.035782    3300    0.01833     -0.0025603
    DOE44         6.688    0.28    2954.6    0.035782    3300    0.20163    -0.00023275
    DOE45         6.688    0.28    2954.6      0.3936     300    0.01833     -0.0025603
    DOE46         6.688    0.28    2954.6      0.3936     300    0.20163    -0.00023275
    DOE47         6.688    0.28    2954.6      0.3936    3300    0.01833    -0.00023275
    DOE48         6.688    0.28    2954.6      0.3936    3300    0.20163     -0.0025603
    DOE49         6.688    3.08     268.6    0.035782     300    0.01833    -0.00023275
    DOE50         6.688    3.08     268.6    0.035782     300    0.20163     -0.0025603
    DOE51         6.688    3.08     268.6    0.035782    3300    0.01833     -0.0025603
    DOE52         6.688    3.08     268.6    0.035782    3300    0.20163    -0.00023275
    DOE53         6.688    3.08     268.6      0.3936     300    0.01833     -0.0025603
    DOE54         6.688    3.08     268.6      0.3936     300    0.20163    -0.00023275
    DOE55         6.688    3.08     268.6      0.3936    3300    0.01833    -0.00023275
    DOE56         6.688    3.08     268.6      0.3936    3300    0.20163     -0.0025603
    DOE57         6.688    3.08    2954.6    0.035782     300    0.01833     -0.0025603
    DOE58         6.688    3.08    2954.6    0.035782     300    0.20163    -0.00023275
    DOE59         6.688    3.08    2954.6    0.035782    3300    0.01833    -0.00023275
    DOE60         6.688    3.08    2954.6    0.035782    3300    0.20163     -0.0025603
    DOE61         6.688    3.08    2954.6      0.3936     300    0.01833    -0.00023275
    DOE62         6.688    3.08    2954.6      0.3936     300    0.20163     -0.0025603
    DOE63         6.688    3.08    2954.6      0.3936    3300    0.01833     -0.0025603
    DOE64         6.688    3.08    2954.6      0.3936    3300    0.20163    -0.00023275
    DOE65         0.608    1.68    1611.6     0.21469    1800    0.10998     -0.0013965
    DOE66         6.688    1.68    1611.6     0.21469    1800    0.10998     -0.0013965
    DOE67         3.648    0.28    1611.6     0.21469    1800    0.10998     -0.0013965
    DOE68         3.648    3.08    1611.6     0.21469    1800    0.10998     -0.0013965
    DOE69         3.648    1.68     268.6     0.21469    1800    0.10998     -0.0013965
    DOE70         3.648    1.68    2954.6     0.21469    1800    0.10998     -0.0013965
    DOE71         3.648    1.68    1611.6    0.035782    1800    0.10998     -0.0013965
    DOE72         3.648    1.68    1611.6      0.3936    1800    0.10998     -0.0013965
    DOE73         3.648    1.68    1611.6     0.21469     300    0.10998     -0.0013965
    DOE74         3.648    1.68    1611.6     0.21469    3300    0.10998     -0.0013965
    DOE75         3.648    1.68    1611.6     0.21469    1800    0.01833     -0.0013965
    DOE76         3.648    1.68    1611.6     0.21469    1800    0.20163     -0.0013965
    DOE77         3.648    1.68    1611.6     0.21469    1800    0.10998     -0.0025603
    DOE78         3.648    1.68    1611.6     0.21469    1800    0.10998    -0.00023275
    DOE79         3.648    1.68    1611.6     0.21469    1800    0.10998     -0.0013965

2.7-Run the Model

ds(doeRuns,:) = runModel(ds(doeRuns,:))
ds = 

                  m        iyy     Kp        zetaP       Kt      zetaT      xcg            unstable
    Initial        6.08     2.8      2686     0.35782    3000     0.1833     -0.0023275    1       
    upperBound    6.688    3.08    2954.6      0.3936    3300    0.20163    -0.00023275    1       
    lowerBound    0.608    0.28     268.6    0.035782     300    0.01833     -0.0025603    0       
    deltaX11      5.472     2.8      2686     0.35782    3000     0.1833     -0.0023275    1       
    deltaX12      6.688     2.8      2686     0.35782    3000     0.1833     -0.0023275    1       
    deltaX21       6.08    2.52      2686     0.35782    3000     0.1833     -0.0023275    1       
    deltaX22       6.08    3.08      2686     0.35782    3000     0.1833     -0.0023275    1       
    deltaX31       6.08     2.8    2417.4     0.35782    3000     0.1833     -0.0023275    1       
    deltaX32       6.08     2.8    2954.6     0.35782    3000     0.1833     -0.0023275    1       
    deltaX41       6.08     2.8      2686     0.32204    3000     0.1833     -0.0023275    1       
    deltaX42       6.08     2.8      2686      0.3936    3000     0.1833     -0.0023275    1       
    deltaX51       6.08     2.8      2686     0.35782    2700     0.1833     -0.0023275    1       
    deltaX52       6.08     2.8      2686     0.35782    3300     0.1833     -0.0023275    1       
    deltaX61       6.08     2.8      2686     0.35782    3000    0.16497     -0.0023275    1       
    deltaX62       6.08     2.8      2686     0.35782    3000    0.20163     -0.0023275    1       
    deltaX71       6.08     2.8      2686     0.35782    3000     0.1833     -0.0023336    1       
    deltaX72       6.08     2.8      2686     0.35782    3000     0.1833     -0.0023214    1       
    DOE1          0.608    0.28     268.6    0.035782     300    0.01833    -0.00023275    0       
    DOE2          0.608    0.28     268.6    0.035782     300    0.20163     -0.0025603    0       
    DOE3          0.608    0.28     268.6    0.035782    3300    0.01833     -0.0025603    1       
    DOE4          0.608    0.28     268.6    0.035782    3300    0.20163    -0.00023275    0       
    DOE5          0.608    0.28     268.6      0.3936     300    0.01833     -0.0025603    0       
    DOE6          0.608    0.28     268.6      0.3936     300    0.20163    -0.00023275    0       
    DOE7          0.608    0.28     268.6      0.3936    3300    0.01833    -0.00023275    1       
    DOE8          0.608    0.28     268.6      0.3936    3300    0.20163     -0.0025603    0       
    DOE9          0.608    0.28    2954.6    0.035782     300    0.01833     -0.0025603    0       
    DOE10         0.608    0.28    2954.6    0.035782     300    0.20163    -0.00023275    0       
    DOE11         0.608    0.28    2954.6    0.035782    3300    0.01833    -0.00023275    0       
    DOE12         0.608    0.28    2954.6    0.035782    3300    0.20163     -0.0025603    1       
    DOE13         0.608    0.28    2954.6      0.3936     300    0.01833    -0.00023275    0       
    DOE14         0.608    0.28    2954.6      0.3936     300    0.20163     -0.0025603    0       
    DOE15         0.608    0.28    2954.6      0.3936    3300    0.01833     -0.0025603    1       
    DOE16         0.608    0.28    2954.6      0.3936    3300    0.20163    -0.00023275    1       
    DOE17         0.608    3.08     268.6    0.035782     300    0.01833     -0.0025603    0       
    DOE18         0.608    3.08     268.6    0.035782     300    0.20163    -0.00023275    0       
    DOE19         0.608    3.08     268.6    0.035782    3300    0.01833    -0.00023275    0       
    DOE20         0.608    3.08     268.6    0.035782    3300    0.20163     -0.0025603    0       
    DOE21         0.608    3.08     268.6      0.3936     300    0.01833    -0.00023275    0       
    DOE22         0.608    3.08     268.6      0.3936     300    0.20163     -0.0025603    0       
    DOE23         0.608    3.08     268.6      0.3936    3300    0.01833     -0.0025603    0       
    DOE24         0.608    3.08     268.6      0.3936    3300    0.20163    -0.00023275    0       
    DOE25         0.608    3.08    2954.6    0.035782     300    0.01833    -0.00023275    0       
    DOE26         0.608    3.08    2954.6    0.035782     300    0.20163     -0.0025603    0       
    DOE27         0.608    3.08    2954.6    0.035782    3300    0.01833     -0.0025603    0       
    DOE28         0.608    3.08    2954.6    0.035782    3300    0.20163    -0.00023275    0       
    DOE29         0.608    3.08    2954.6      0.3936     300    0.01833     -0.0025603    0       
    DOE30         0.608    3.08    2954.6      0.3936     300    0.20163    -0.00023275    0       
    DOE31         0.608    3.08    2954.6      0.3936    3300    0.01833    -0.00023275    0       
    DOE32         0.608    3.08    2954.6      0.3936    3300    0.20163     -0.0025603    0       
    DOE33         6.688    0.28     268.6    0.035782     300    0.01833     -0.0025603    0       
    DOE34         6.688    0.28     268.6    0.035782     300    0.20163    -0.00023275    0       
    DOE35         6.688    0.28     268.6    0.035782    3300    0.01833    -0.00023275    0       
    DOE36         6.688    0.28     268.6    0.035782    3300    0.20163     -0.0025603    0       
    DOE37         6.688    0.28     268.6      0.3936     300    0.01833    -0.00023275    0       
    DOE38         6.688    0.28     268.6      0.3936     300    0.20163     -0.0025603    0       
    DOE39         6.688    0.28     268.6      0.3936    3300    0.01833     -0.0025603    0       
    DOE40         6.688    0.28     268.6      0.3936    3300    0.20163    -0.00023275    0       
    DOE41         6.688    0.28    2954.6    0.035782     300    0.01833    -0.00023275    0       
    DOE42         6.688    0.28    2954.6    0.035782     300    0.20163     -0.0025603    0       
    DOE43         6.688    0.28    2954.6    0.035782    3300    0.01833     -0.0025603    0       
    DOE44         6.688    0.28    2954.6    0.035782    3300    0.20163    -0.00023275    0       
    DOE45         6.688    0.28    2954.6      0.3936     300    0.01833     -0.0025603    0       
    DOE46         6.688    0.28    2954.6      0.3936     300    0.20163    -0.00023275    0       
    DOE47         6.688    0.28    2954.6      0.3936    3300    0.01833    -0.00023275    0       
    DOE48         6.688    0.28    2954.6      0.3936    3300    0.20163     -0.0025603    0       
    DOE49         6.688    3.08     268.6    0.035782     300    0.01833    -0.00023275    0       
    DOE50         6.688    3.08     268.6    0.035782     300    0.20163     -0.0025603    0       
    DOE51         6.688    3.08     268.6    0.035782    3300    0.01833     -0.0025603    1       
    DOE52         6.688    3.08     268.6    0.035782    3300    0.20163    -0.00023275    1       
    DOE53         6.688    3.08     268.6      0.3936     300    0.01833     -0.0025603    0       
    DOE54         6.688    3.08     268.6      0.3936     300    0.20163    -0.00023275    0       
    DOE55         6.688    3.08     268.6      0.3936    3300    0.01833    -0.00023275    1       
    DOE56         6.688    3.08     268.6      0.3936    3300    0.20163     -0.0025603    1       
    DOE57         6.688    3.08    2954.6    0.035782     300    0.01833     -0.0025603    0       
    DOE58         6.688    3.08    2954.6    0.035782     300    0.20163    -0.00023275    0       
    DOE59         6.688    3.08    2954.6    0.035782    3300    0.01833    -0.00023275    1       
    DOE60         6.688    3.08    2954.6    0.035782    3300    0.20163     -0.0025603    1       
    DOE61         6.688    3.08    2954.6      0.3936     300    0.01833    -0.00023275    0       
    DOE62         6.688    3.08    2954.6      0.3936     300    0.20163     -0.0025603    0       
    DOE63         6.688    3.08    2954.6      0.3936    3300    0.01833     -0.0025603    1       
    DOE64         6.688    3.08    2954.6      0.3936    3300    0.20163    -0.00023275    1       
    DOE65         0.608    1.68    1611.6     0.21469    1800    0.10998     -0.0013965    0       
    DOE66         6.688    1.68    1611.6     0.21469    1800    0.10998     -0.0013965    0       
    DOE67         3.648    0.28    1611.6     0.21469    1800    0.10998     -0.0013965    0       
    DOE68         3.648    3.08    1611.6     0.21469    1800    0.10998     -0.0013965    0       
    DOE69         3.648    1.68     268.6     0.21469    1800    0.10998     -0.0013965    0       
    DOE70         3.648    1.68    2954.6     0.21469    1800    0.10998     -0.0013965    0       
    DOE71         3.648    1.68    1611.6    0.035782    1800    0.10998     -0.0013965    0       
    DOE72         3.648    1.68    1611.6      0.3936    1800    0.10998     -0.0013965    0       
    DOE73         3.648    1.68    1611.6     0.21469     300    0.10998     -0.0013965    0       
    DOE74         3.648    1.68    1611.6     0.21469    3300    0.10998     -0.0013965    1       
    DOE75         3.648    1.68    1611.6     0.21469    1800    0.01833     -0.0013965    0       
    DOE76         3.648    1.68    1611.6     0.21469    1800    0.20163     -0.0013965    0       
    DOE77         3.648    1.68    1611.6     0.21469    1800    0.10998     -0.0025603    0       
    DOE78         3.648    1.68    1611.6     0.21469    1800    0.10998    -0.00023275    0       
    DOE79         3.648    1.68    1611.6     0.21469    1800    0.10998     -0.0013965    0       


                  time                    pitch                   plunge                  pitchFreq
    Initial       [ 3417x1 double]        [ 3417x1 double]        [ 3417x1 double]         3.5494  
    upperBound    [ 3336x1 double]        [ 3336x1 double]        [ 3336x1 double]          3.672  
    lowerBound    [10523x1 double]        [10523x1 double]        [10523x1 double]         3.0505  
    deltaX11      [ 3318x1 double]        [ 3318x1 double]        [ 3318x1 double]         3.6721  
    deltaX12      [ 3499x1 double]        [ 3499x1 double]        [ 3499x1 double]         3.5492  
    deltaX21      [ 3320x1 double]        [ 3320x1 double]        [ 3320x1 double]         3.6721  
    deltaX22      [ 3496x1 double]        [ 3496x1 double]        [ 3496x1 double]         3.5492  
    deltaX31      [ 3491x1 double]        [ 3491x1 double]        [ 3491x1 double]         3.5492  
    deltaX32      [ 3610x1 double]        [ 3610x1 double]        [ 3610x1 double]         3.6713  
    deltaX41      [ 3417x1 double]        [ 3417x1 double]        [ 3417x1 double]         3.5494  
    deltaX42      [ 3417x1 double]        [ 3417x1 double]        [ 3417x1 double]         3.5494  
    deltaX51      [ 3875x1 double]        [ 3875x1 double]        [ 3875x1 double]         3.4259  
    deltaX52      [ 3224x1 double]        [ 3224x1 double]        [ 3224x1 double]         3.7948  
    deltaX61      [ 3417x1 double]        [ 3417x1 double]        [ 3417x1 double]         3.5494  
    deltaX62      [ 3416x1 double]        [ 3416x1 double]        [ 3416x1 double]         3.5494  
    deltaX71      [ 3417x1 double]        [ 3417x1 double]        [ 3417x1 double]         3.5494  
    deltaX72      [ 3417x1 double]        [ 3417x1 double]        [ 3417x1 double]         3.5494  
    DOE1          [10568x1 double]        [10568x1 double]        [10568x1 double]          3.515  
    DOE2          [10444x1 double]        [10444x1 double]        [10444x1 double]         2.2946  
    DOE3          [ 1158x1 double]        [ 1158x1 double]        [ 1158x1 double]         11.605  
    DOE4          [10014x1 double]        [10014x1 double]        [10014x1 double]         11.581  
    DOE5          [10450x1 double]        [10450x1 double]        [10450x1 double]         3.2844  
    DOE6          [10461x1 double]        [10461x1 double]        [10461x1 double]         2.5218  
    DOE7          [ 1158x1 double]        [ 1158x1 double]        [ 1158x1 double]         11.605  
    DOE8          [10014x1 double]        [10014x1 double]        [10014x1 double]         11.581  
    DOE9          [10626x1 double]        [10626x1 double]        [10626x1 double]         12.127  
    DOE10         [10590x1 double]        [10590x1 double]        [10590x1 double]         11.052  
    DOE11         [10507x1 double]        [10507x1 double]        [10507x1 double]         10.645  
    DOE12         [  965x1 double]        [  965x1 double]        [  965x1 double]         12.869  
    DOE13         [10473x1 double]        [10473x1 double]        [10473x1 double]         11.665  
    DOE14         [10537x1 double]        [10537x1 double]        [10537x1 double]         10.578  
    DOE15         [  965x1 double]        [  965x1 double]        [  965x1 double]         12.869  
    DOE16         [  965x1 double]        [  965x1 double]        [  965x1 double]         12.869  
    DOE17         [10052x1 double]        [10052x1 double]        [10052x1 double]         3.2207  
    DOE18         [10072x1 double]        [10072x1 double]        [10072x1 double]         3.2886  
    DOE19         [10040x1 double]        [10040x1 double]        [10040x1 double]         3.3394  
    DOE20         [10048x1 double]        [10048x1 double]        [10048x1 double]          3.434  
    DOE21         [10058x1 double]        [10058x1 double]        [10058x1 double]          3.284  
    DOE22         [10103x1 double]        [10103x1 double]        [10103x1 double]          3.237  
    DOE23         [10049x1 double]        [10049x1 double]        [10049x1 double]         3.3117  
    DOE24         [10058x1 double]        [10058x1 double]        [10058x1 double]         3.4375  
    DOE25         [10055x1 double]        [10055x1 double]        [10055x1 double]         11.352  
    DOE26         [10063x1 double]        [10063x1 double]        [10063x1 double]         11.208  
    DOE27         [10010x1 double]        [10010x1 double]        [10010x1 double]         3.4821  
    DOE28         [10010x1 double]        [10010x1 double]        [10010x1 double]         3.4821  
    DOE29         [10043x1 double]        [10043x1 double]        [10043x1 double]         11.278  
    DOE30         [10070x1 double]        [10070x1 double]        [10070x1 double]         11.277  
    DOE31         [10010x1 double]        [10010x1 double]        [10010x1 double]         3.4821  
    DOE32         [10010x1 double]        [10010x1 double]        [10010x1 double]         3.4821  
    DOE33         [10232x1 double]        [10232x1 double]        [10232x1 double]         1.0303  
    DOE34         [10459x1 double]        [10459x1 double]        [10459x1 double]         1.0213  
    DOE35         [10013x1 double]        [10013x1 double]        [10013x1 double]        0.97773  
    DOE36         [10013x1 double]        [10013x1 double]        [10013x1 double]        0.97773  
    DOE37         [10410x1 double]        [10410x1 double]        [10410x1 double]         1.0165  
    DOE38         [10490x1 double]        [10490x1 double]        [10490x1 double]         0.9603  
    DOE39         [10013x1 double]        [10013x1 double]        [10013x1 double]        0.97773  
    DOE40         [10013x1 double]        [10013x1 double]        [10013x1 double]        0.97773  
    DOE41         [10399x1 double]        [10399x1 double]        [10399x1 double]         3.3636  
    DOE42         [10570x1 double]        [10570x1 double]        [10570x1 double]         3.3222  
    DOE43         [10013x1 double]        [10013x1 double]        [10013x1 double]          3.361  
    DOE44         [10013x1 double]        [10013x1 double]        [10013x1 double]          3.361  
    DOE45         [10344x1 double]        [10344x1 double]        [10344x1 double]         3.4089  
    DOE46         [10488x1 double]        [10488x1 double]        [10488x1 double]         3.4244  
    DOE47         [10013x1 double]        [10013x1 double]        [10013x1 double]          3.361  
    DOE48         [10013x1 double]        [10013x1 double]        [10013x1 double]          3.361  
    DOE49         [10020x1 double]        [10020x1 double]        [10020x1 double]         1.0396  
    DOE50         [10037x1 double]        [10037x1 double]        [10037x1 double]        0.94945  
    DOE51         [ 4972x1 double]        [ 4972x1 double]        [ 4972x1 double]         3.4853  
    DOE52         [ 8140x1 double]        [ 8140x1 double]        [ 8140x1 double]         3.4829  
    DOE53         [10058x1 double]        [10058x1 double]        [10058x1 double]         1.0128  
    DOE54         [10041x1 double]        [10041x1 double]        [10041x1 double]        0.94983  
    DOE55         [ 4973x1 double]        [ 4973x1 double]        [ 4973x1 double]         3.4853  
    DOE56         [ 8137x1 double]        [ 8137x1 double]        [ 8137x1 double]         3.4829  
    DOE57         [10041x1 double]        [10041x1 double]        [10041x1 double]         3.2784  
    DOE58         [10047x1 double]        [10047x1 double]        [10047x1 double]         1.7475  
    DOE59         [ 3336x1 double]        [ 3336x1 double]        [ 3336x1 double]          3.672  
    DOE60         [ 3334x1 double]        [ 3334x1 double]        [ 3334x1 double]          3.672  
    DOE61         [10032x1 double]        [10032x1 double]        [10032x1 double]          3.398  
    DOE62         [10036x1 double]        [10036x1 double]        [10036x1 double]         3.3687  
    DOE63         [ 3334x1 double]        [ 3334x1 double]        [ 3334x1 double]          3.672  
    DOE64         [ 3336x1 double]        [ 3336x1 double]        [ 3336x1 double]          3.672  
    DOE65         [10372x1 double]        [10372x1 double]        [10372x1 double]         0.1266  
    DOE66         [10255x1 double]        [10255x1 double]        [10255x1 double]         2.5034  
    DOE67         [10648x1 double]        [10648x1 double]        [10648x1 double]         3.1517  
    DOE68         [10340x1 double]        [10340x1 double]        [10340x1 double]         3.3445  
    DOE69         [10061x1 double]        [10061x1 double]        [10061x1 double]         1.3508  
    DOE70         [10487x1 double]        [10487x1 double]        [10487x1 double]         4.5761  
    DOE71         [10333x1 double]        [10333x1 double]        [10333x1 double]         3.3738  
    DOE72         [10200x1 double]        [10200x1 double]        [10200x1 double]         3.3615  
    DOE73         [10078x1 double]        [10078x1 double]        [10078x1 double]          3.352  
    DOE74         [ 2163x1 double]        [ 2163x1 double]        [ 2163x1 double]         4.7806  
    DOE75         [10164x1 double]        [10164x1 double]        [10164x1 double]         3.3496  
    DOE76         [10152x1 double]        [10152x1 double]        [10152x1 double]         3.3457  
    DOE77         [10226x1 double]        [10226x1 double]        [10226x1 double]         3.3701  
    DOE78         [10164x1 double]        [10164x1 double]        [10164x1 double]         3.3806  
    DOE79         [10230x1 double]        [10230x1 double]        [10230x1 double]         3.3402  


                  plungeFreq
    Initial         3.6718  
    upperBound       3.672  
    lowerBound      35.611  
    deltaX11        3.7945  
    deltaX12        3.5492  
    deltaX21        3.6721  
    deltaX22        3.5492  
    deltaX31        3.5492  
    deltaX32        3.6713  
    deltaX41        3.6718  
    deltaX42        3.6718  
    deltaX51        3.4259  
    deltaX52        3.7948  
    deltaX61        3.6718  
    deltaX62        3.6718  
    deltaX71        3.6718  
    deltaX72        3.6718  
    DOE1            25.831  
    DOE2             19.95  
    DOE3            11.605  
    DOE4            11.581  
    DOE5            39.892  
    DOE6            25.218  
    DOE7            11.605  
    DOE8            11.581  
    DOE9             56.29  
    DOE10           37.421  
    DOE11           54.441  
    DOE12           12.869  
    DOE13           34.579  
    DOE14           55.239  
    DOE15           12.869  
    DOE16           12.869  
    DOE17           3.2207  
    DOE18         0.030734  
    DOE19           15.625  
    DOE20           15.606  
    DOE21           10.497  
    DOE22            3.237  
    DOE23           15.639  
    DOE24           15.561  
    DOE25           16.108  
    DOE26           8.9971  
    DOE27           3.4821  
    DOE28           3.4821  
    DOE29           16.334  
    DOE30           16.224  
    DOE31           3.4821  
    DOE32           3.4821  
    DOE33           41.401  
    DOE34           37.883  
    DOE35           11.458  
    DOE36           11.458  
    DOE37             71.6  
    DOE38           55.377  
    DOE39           11.458  
    DOE40           11.458  
    DOE41           51.977  
    DOE42           36.608  
    DOE43           11.458  
    DOE44           11.458  
    DOE45           7.1335  
    DOE46           36.356  
    DOE47           11.458  
    DOE48           11.458  
    DOE49           11.313  
    DOE50            11.73  
    DOE51           3.4853  
    DOE52           3.4829  
    DOE53           15.284  
    DOE54           6.6488  
    DOE55           3.4853  
    DOE56           3.4829  
    DOE57           11.674  
    DOE58           8.6149  
    DOE59            3.672  
    DOE60            3.672  
    DOE61           16.071  
    DOE62           20.304  
    DOE63            3.672  
    DOE64            3.672  
    DOE65           0.1266  
    DOE66           5.6953  
    DOE67           13.972  
    DOE68          0.66259  
    DOE69           5.6796  
    DOE70            1.504  
    DOE71           5.8963  
    DOE72           5.7892  
    DOE73           7.6881  
    DOE74           4.7806  
    DOE75          0.34117  
    DOE76           17.472  
    DOE77         0.093613  
    DOE78           16.996  
    DOE79         0.031216  

2.8-Visualize the Results

The main effects plot shows how the average response (ptich frequency in this case) changes with the design factors and indicates the response is nonlinear (kink in the curves).

The interaction plot shows how each factor changes relative to settings in other factors. If the plot is parallel, it indicates that there are no interaction effects (i.e. factor a * factor b). If the plots cross, or converge/diverge it indicates an interaction is present. From the plot we can see that interactions appear to be important. However, from this plot it is hard to tell wich ones are significant.

X = double(ds(doeRuns,factors));
Y = double(ds(doeRuns,response));
figure, maineffectsplot(Y(:,1),doe,'varnames',factors)
figure, interactionplot(Y(:,1),doe,'varnames',factors)

2.9-Fit a Quadratic Model and Determine Which Parameters are Signficant

We fit a quadratic model to data using stepwise regression. Stepwise regression automatically adds terms to the model and arrives at a model that best captures our data with the least number. It adds a term to the model if it is significant (with 95% confidence).

The pareto plot shows the parameters ordered by factors the have the most influence on the response, pitch frequency and the amount of uncertainty captured in the black line by adding in each factor. We can see that the terms that are important and account for up to 85% of the variation observed. The interaction terms that are important become apparent as well.

D = x2fx(doe,'quadratic');
[b,se,pval,inmodel] = stepwisefit(D,Y(:,1));

% generate labels
linearTerms      = factors;
squaredTerms     = strcat(linearTerms,'^2');
interactionTerms = {};
for i = 1:length(linearTerms)-1
    interactionTerms = [interactionTerms strcat(linearTerms(i),'*',linearTerms(i+1:end))];
end
fullQuadraticModel   = [{'const.'} linearTerms, interactionTerms, squaredTerms];

% Show Paretoplot
pareto(abs(b(inmodel)),fullQuadraticModel(inmodel))

% list variables that are important
termNames = fullQuadraticModel(inmodel)
terms = b(inmodel)'
Initial columns included:  none
Step 1, added column 2, p=1.12763e-008
Step 2, added column 4, p=1.74993e-005
Step 3, added column 9, p=0.000223949
Step 4, added column 21, p=0.00025782
Step 5, added column 3, p=0.00140525
Step 6, added column 17, p=0.000672302
Step 7, added column 34, p=0.000730152
Step 8, added column 10, p=0.00391999
Final columns included:  2 3 4 9 10 17 21 34 
    'Coeff'          'Std.Err.'    'Status'    'P'          
    [          0]    [  1.7537]    'Out'       [          1]
    [    -2.3729]    [  0.2143]    'In'        [4.8951e-017]
    [    -0.8698]    [  0.2143]    'In'        [1.2690e-004]
    [     1.5154]    [  0.2143]    'In'        [9.2825e-010]
    [     0.0461]    [  0.2158]    'Out'       [     0.8315]
    [     0.3334]    [  0.2121]    'Out'       [     0.1205]
    [    -0.0558]    [  0.2158]    'Out'       [     0.7966]
    [    -0.0428]    [  0.2158]    'Out'       [     0.8433]
    [     1.1958]    [  0.2176]    'In'        [5.9629e-007]
    [    -0.6494]    [  0.2176]    'In'        [     0.0039]
    [     0.0081]    [  0.2192]    'Out'       [     0.9706]
    [     0.0749]    [  0.2190]    'Out'       [     0.7333]
    [9.6672e-004]    [  0.2192]    'Out'       [     0.9965]
    [     0.0015]    [  0.2192]    'Out'       [     0.9946]
    [    -0.2285]    [  0.2175]    'Out'       [     0.2971]
    [     0.0053]    [  0.2192]    'Out'       [     0.9808]
    [    -0.8782]    [  0.2176]    'In'        [1.3744e-004]
    [     0.0068]    [  0.2192]    'Out'       [     0.9752]
    [     0.0073]    [  0.2192]    'Out'       [     0.9735]
    [     0.0514]    [  0.2191]    'Out'       [     0.8153]
    [    -1.0883]    [  0.2176]    'In'        [4.0694e-006]
    [     0.0054]    [  0.2192]    'Out'       [     0.9805]
    [    -0.0655]    [  0.2191]    'Out'       [     0.7659]
    [     0.0210]    [  0.2192]    'Out'       [     0.9239]
    [    -0.0028]    [  0.2192]    'Out'       [     0.9899]
    [     0.0618]    [  0.2191]    'Out'       [     0.7787]
    [     0.1323]    [  0.2186]    'Out'       [     0.5472]
    [    -0.0242]    [  0.2192]    'Out'       [     0.9123]
    [     0.0245]    [  0.2192]    'Out'       [     0.9112]
    [    -0.4319]    [  0.9188]    'Out'       [     0.6398]
    [     0.6326]    [  0.9171]    'Out'       [     0.4926]
    [     0.4759]    [  0.9185]    'Out'       [     0.6060]
    [     0.6985]    [  0.9164]    'Out'       [     0.4486]
    [     1.9670]    [  0.5283]    'In'        [3.9461e-004]
    [     0.6874]    [  0.9165]    'Out'       [     0.4558]
    [     0.7027]    [  0.9164]    'Out'       [     0.4458]


termNames = 

    'm'    'iyy'    'Kp'    'm*iyy'    'm*Kp'    'iyy*Kt'    'Kp*Kt'    'Kt^2'


terms =

   -2.3729   -0.8698    1.5154    1.1958   -0.6494   -0.8782   -1.0883    1.9670

3-Find Design Parameters that Minimize Failure Region Area

Let's now identify the design parameters, starting from the initial design paraemeteres, that will minimize the failure region. This section extends the anlaysis we've done in 1 and 2.

clear all, close all, close_system(find_system,false)
model  = 'Flutter';
blocks = {'Mach', 'Altitude', 'BACT Wing & PAPA Mount';};

% Create DataSet Array extracting parameters from blocks
ds = getParamNames(model,blocks)
ds = 

               Mach    Altitude    m       iyy    Kp      zetaP      Kt      zetaT     xcg       
    Initial    0.77    30000       6.08    2.8    2686    0.35782    3000    0.1833    -0.0023275

3.1-Define Allowable Bounds

We'll set the bounds to +/- 10% of the initial value

ds('upperBound',:) = datasetfun(@(x) 1.5*x, ds('Initial',:), 'DatasetOutput', true);
ds('lowerBound',:) = datasetfun(@(x) 0.5*x, ds('Initial',:), 'DatasetOutput', true);

% Swap any bounds that might have been negative
swapIndx = double(ds('upperBound',:)) < double(ds('lowerBound',:));
if sum(swapIndx) ~= 0
    tmp = ds('upperBound',swapIndx);
    ds('upperBound',swapIndx) = ds('lowerBound',swapIndx);
    ds('lowerBound',swapIndx) = tmp;
end

% Define Bounds for mach/alt
ds.Mach('lowerBound',1)= 0.65;
ds.Mach('upperBound',1) = 0.85;
ds.Altitude('lowerBound',1) = 10000;
ds.Altitude('upperBound',1) = 50000
ds = 

                  Mach    Altitude    m       iyy    Kp      zetaP      Kt      zetaT   
    Initial       0.77    30000       6.08    2.8    2686    0.35782    3000      0.1833
    upperBound    0.85    50000       9.12    4.2    4029    0.53673    4500     0.27495
    lowerBound    0.65    10000       3.04    1.4    1343    0.17891    1500    0.091652


                  xcg       
    Initial       -0.0023275
    upperBound    -0.0011637
    lowerBound    -0.0034912

3.2-Run Model Over Bounds and Intial Point

ds = runModel(ds)
ds = 

                  Mach    Altitude    m       iyy    Kp      zetaP      Kt      zetaT   
    Initial       0.77    30000       6.08    2.8    2686    0.35782    3000      0.1833
    upperBound    0.85    50000       9.12    4.2    4029    0.53673    4500     0.27495
    lowerBound    0.65    10000       3.04    1.4    1343    0.17891    1500    0.091652


                  xcg           unstable    time                    pitch               
    Initial       -0.0023275    1           [ 3417x1 double]        [ 3417x1 double]    
    upperBound    -0.0011637    0           [10011x1 double]        [10011x1 double]    
    lowerBound    -0.0034912    0           [10169x1 double]        [10169x1 double]    


                  plunge                  pitchFreq    plungeFreq
    Initial       [ 3417x1 double]        3.5494       3.6718    
    upperBound    [10011x1 double]        3.3603       4.6739    
    lowerBound    [10169x1 double]        3.3202       5.8647    

3.3-Initialize Design Space with Design of Experiments

Set an initial response for boundarySearch to start from.

factors = {'Mach','Altitude'};
doe = ccdesign(length(factors),'center',1,'type','faced')
doe =

    -1    -1
    -1     1
     1    -1
     1     1
    -1     0
     1     0
     0    -1
     0     1
     0     0

%Scale to Engineering Units
doeRuns    = length(doe);
doeRuns    = strrep( strcat({'DOE'},num2str((1:doeRuns)','%d')), ' ','');
doeRange   = [min(doe(:,1)) max(doe(:,1))];
for i = 1:length(factors)
    range = ds.(factors{i})({'lowerBound','upperBound'});
    ds.(factors{i})(doeRuns) = interp1(doeRange, range, doe(:,i));
end

3.2-Define Limit Function and Objective Function

We want to minimize flutter boundary area. We need to run boundarySearch for each optimization iteration and measure it's area.

% Define design variables and bounds
designVariables = {'m','iyy','Kp','zetaP','Kt','zetaT','xcg'};
for i = 1:length(designVariables)
    ds.(designVariables{i})(doeRuns) = ds('Initial',designVariables{i});
end
X0 = double(ds(doeRuns,factors));
lb = double(ds('lowerBound',factors));
ub = double(ds('upperBound',factors));

% Scale factors to be similar in magnitude for |boundarySearch|
X0(:,1) = X0(:,1).*10;    % Mach Scaling
X0(:,2) = X0(:,2)./10000; % Altitude Scaling
lb = lb.*[10 1/10000];
ub = ub.*[10 1/10000];

% Define initial conditions (starting point) and bounds
Y0 = double(ds('Initial',designVariables));
lbY = double(ds('lowerBound',designVariables));
ubY = double(ds('upperBound',designVariables));

% Define the limit function
limitFunction = @(X,Y)flutterLimitFunction(X,factors,designVariables,Y);

% And the objective function (minimize area)
A = @(Y) boundaryArea(Y,X0,lb,ub,limitFunction)

% Initial Area
A(Y0)
A = 

    @(Y)boundaryArea(Y,X0,lb,ub,limitFunction)


ans =

    4.7424

3.3-Optimize

Changes in the flutter boundary are expected to be discontinuos (jumps or noisy) and nonlinear. So we'll use the pattern search solver from Genetic Algorithm and Direct Search Toolbox (TM) to find the minimum.

*Note: This example limits the iterations to 3, to show the solver progress. The solver will continue to improve the solution (but it takes a long time). Remove the MaxIter Option to run longer.

options = optimset('Disp','iter','MaxIter',3);
tic
[bestParameters FlutterArea] = patternsearch(A,Y0,[],[],[],[],...
                                             lbY,ubY,[],options)
toc

Iter     f-count          f(x)      MeshSize     Method
    0        1        4.74243             1      
    1        1        4.74243           0.5     Refine Mesh
    2        2        3.07774             1     Successful Poll
    3        3        3.07774           0.5     Refine Mesh
    4        7        1.53391             1     Successful Poll
Maximum number of iterations exceeded: increase options.MaxIter.

bestParameters =

  1.0e+003 *

    0.0081    0.0028    2.6860    0.0004    4.0240    0.0002   -0.0000


FlutterArea =

    1.5339

Elapsed time is 1759.629238 seconds.

Close all models

close_system(find_system,false)

Contact us