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
- 1.1-Define Bounds
- 1.2-Run Model Over Bounds
- 1.3-Design Exploration Via Design of Experiments
- 1.4-Scale to Engineering Units
- 1.5-Run Simulink at DOE Points
- 1.6-Find Failure Boundary Using boundarySearch
- 2-Sensitivity Analysis
- 2.1-Define Allowable Bounds
- 2.2-Run Simulink with Given Dataset and Bounds (Initial Run)
- 2.3-Define Our Factors and Responses of Interest
- 2.4-Sensitivity Analysis by Gradient (Local Method)
- 2.5-Sensitivity Analysis by Response Surfaces (more Global Method)
- 2.6-Scale DOE to Engineering Units
- 2.7-Run the Model
- 2.8-Visualize the Results
- 2.9-Fit a Quadratic Model and Determine Which Parameters are Signficant
- 3-Find Design Parameters that Minimize Failure Region Area
- 3.1-Define Allowable Bounds
- 3.2-Run Model Over Bounds and Intial Point
- 3.3-Initialize Design Space with Design of Experiments
- 3.2-Define Limit Function and Objective Function
- 3.3-Optimize
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)