MATLAB Examples

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.

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

```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) ```