MATLAB Examples

First Order Reliability Method using a Simulink Flutter Model

This examples illustrates how to perform a FORM analysis on a discrete (0 or 1) failure response. In the example we'll compare a traditional Monte Carlo method with FORM. This example is was created to illustrate how to use apply probabalistic techniques o a Simulink Model (originally developed as a demo for an aerospace customer).


Flutter Model

This Simulink model is adapted from Cole Stephen's MATLAB Central submission:

We're using it here to illustrate how to perform a probability of failure analysis for:

  1. Simulink model
  2. Discrete valued failure response

FORM is traditionally a gradient based linear approximation of failure probability and in this example we extend it to a nonsmooth example using solvers from Global Optimization Toolbox.

model  = 'Flutter';
setup, open(model)

Run Model

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. Flutter occurs at this condition.

Mach = 0.77;
Altitude = 20000;
flutter = sigsOut.isStable.Data(end)
flutter =


Failure Boundary

This model was previously run to map the flutter boundary and is reproduced here to provide context. In this example, we'd like to evaluate the amount of time flutter will occur (without active suppression) for an aircraft at 0.77 Mach Number and 17,000 ft (think of this as a typical cruising altitude for a Boeing 737 on short routes, e.g. Boston to Washington D.C.).

See also the demo on mapping the failure boundary:

open flutterBoundary.fig

Random Variables

We'd like to estimate the probability of flutter (failure) for nominal operating condition of 0.77 Mach at 17,000 ft. Define Mach and Altitude using probability distribution objects and show distributions graphically. We'll also assume that Mach and Altitude are independent.

Mach Number is normally distributed (assume)

mu(1) = 0.77;
sigma(1) = 0.01;
M = ProbDistUnivParam('normal',[mu(1) sigma(1)])
M = 

normal distribution

    mu = 0.77
    sigma = 0.01

Altitude in units {ft} and is lognormally distributed (assume)

m     = 17000;  % mean
v     = 500^2; % variance

Note that the inputs for the lognormal disribution are the mean and standard deviation of a normal distribution see the documentation with (doc lognormal) for more information. So we need to transfrom m,v to equivalent normal values

mu(2)    = log(m^2/sqrt(v+m^2)) % need to transfrom m,v to mu,sigma
sigma(2) = sqrt(log(v/m^2+1))
A = ProbDistUnivParam('lognormal',[mu(2) sigma(2)])
mu =

    0.7700    9.7405

sigma =

    0.0100    0.0294

A = 

lognormal distribution

    mu = 9.74054
    sigma = 0.0294054

Plot probability density functions for each variable.

ax1 = subplot(2,1,1);
x = (mu(1)-3*sigma(1)):sigma(1)/10:(mu(1)+3*sigma(1));
xlabel('Mach Number'), ylabel('Probability Density')
ax2 = subplot(2,1,2);
x = (m-3*sqrt(v)):sqrt(v)/10:(m+3*sqrt(v));
xlabel('Altitude'), ylabel('Probability Density')

Limit State Function

Define the limit function (failure function) as function that accepts inputs M and A (as a matrix [M A]).

type flutterLimitFunction
function failed = flutterLimitFunction(X,model)
% Calcualte Flutter Limit Function

% Copyright 2010 The MathWorks, Inc.
%% Scale X's
% Good practice to scale variables of different magnitudes when using
% optimization, need to scale back to model units.
X(1) = X(1)/10;    % Mach Scaling
X(2) = X(2)*10000; % Altitude Scaling

%% Run Simulink Model
Mach = X(1);
Altitude = X(2);
s = sim(model,'SrcWorkspace','current');
sigsOut = get(s,'sigsOut');

%% Extract Failure Measure
    failed = sigsOut.isStable.Data(end);
    disp('Stopped here because of error')

A 1 indicats flutter and a 0 indicates no flutter from the simulink model. It is standard practice to define the limit state function as G(...) <= 0, so we'll switch the sign to be consistent with convention. So flutter occurs when G = 0. It does not occur when G = 1.

G = @(x) 1-flutterLimitFunction(x,model)
G = 


Test at the mean values (scaled for limit function inputs)

X0 = [mean(M)*10,mean(A)/10000];
ans =


Close Simulink model to run in background (faster)

close_system(model,false)   % close model
close all                   % close figures

Use the model without animation

model = 'Flutter_noAnimation';
G = @(x) 1-flutterLimitFunction(x,model)
G = 


Monte Carlo Simulation

We'll first run a Monte Carlo simulation to provide us with a reasonable estimate of the failure probability. Then we'll compare it to FORM. This simulation will also limit the number of simulation points by establising a convergence criteria based upon the coefficient of varation (in essence, stop the simulation once we achieve 90% confidence in our estimate of the failure probability).

Generate random numbers from the distribution objects for simulation and show a scatterplot with marginal histograms

nSamples = 32000;
xM = random(M,1,nSamples);
xA = random(A,1,nSamples);

xlabel('Mach Number'),ylabel('Altitude')

Run over all points and calculate the failure probability with confidence intervals. The simulation will stop when the coefficient of varation reaches 0.1, or 90% confidence.

blkSize = 160;
mcGBlk = zeros(blkSize,1);
mcG = [];
files = dir;
    'FileDependencies',{files(3:end).name})  % use 4 core machine
cv = 1;
block = 1;
while cv(end) > 0.10
    % Run a block in parallel
    parfor n = 1:blkSize
        mcGBlk(n) = G([xM(block-1+n)*10,xA(block-1+n)/10000]);
    block = block+blkSize;

    mcG = [mcG; mcGBlk];

    N   = 1:length(mcG);
    pf  = cumsum(mcG <=0)./N';
    mpf = cumsum(pf)./N';
    spf = sqrt(mpf.*(1-mpf)./N');
    cv  = sqrt((1-mpf)./(N'.*mpf)); % coefficient of varation

    % Convergence plot
    xlabel('Limit Function G(M,A)'), ylabel('Frequency')
    title(['Probability of Failure = ',num2str(pf(end),4),' Iteration = ',...
        num2str(length(mcG)),' CV = ',num2str(cv(end),4)])

    hold on
    plot(N,mpf-1.96*spf,'b-',N,mpf + 1.96*spf,'b-')
    hold off
    legend('Probability Estimate','Lower Bound','Upper Bound',...
    title('Monte Carlo Convergence with 95% confidence interval')
    % protect against premature convergence due to low sample size
    if N(end) <= blkSize
        cv( cv(1:end) <= 0.10) = 1;
tmc = toc;
axLim = axis;
cvIndex = find(cv <= 0.10);
hold on
    {['CV = 0.10, pf = ',num2str(pf(cvIndex(1)),4),'\rightarrow'],...
    ['Iterations = ',num2str(cvIndex(1)),'\rightarrow']},...
    legend('Probability Estimate','Lower Bound','Upper Bound',...
hold off
Pfmc = pf(cvIndex(1))
Starting matlabpool using the 'speedy-cluster' configuration ... connected to 16 labs.

Pfmc =


First Order Reliability Method (FORM)

FORM takes the input space, X, and trasfroms it to a standard normal space, then solves for the point in the failure region that is closest to the origin.

Transform to Standard Normal Space for Independent Random Variables X --> U. Since the marginals (the input random variables) are independent (not correlated), the transfomation is a simple mapping using the CDF of the marginal distribution. Note that for the case when the marginals are correlated, we'd need to do a more invovled transfromation, using a Nataf or Rosenblat transformation.

% Tranform X to U
T = @(x) [ norminv( normcdf(x(:,1),mu(1),sigma(1)), 0,1 ),...
           norminv( logncdf(x(:,2),mu(2),sigma(2)), 0,1 ) ];

% Inverse Transform from U to X
Tinv = @(u)[ norminv( normcdf(u(:,1),0,1), mu(1),sigma(1) ),...
             logninv( normcdf(u(:,2),0,1), mu(2),sigma(2) ) ];

% Plot the mapping
index = 1:length(pf);
X = [xM(index)' xA(index)'];
U = T(X(logical(mcG),:));
Uf = T(X(~mcG,:));
axis equal
xlabel('Normalized Mach Number'),ylabel('Normalized Altitude')
title('Standard Normal Space')
hold on
legend('Test Points','Failed Test Points')

There are a few failure points that show in outside the failure region. These failures are a result of the numeric stiffness at these locations and the solver limits on time-steps that were used. I didn't look into changing the solver type/options for this problem.

First Order Reliability Problem Formulation

FORM expressed as an optimization problem is:

Objective: Minimize norm(U)

Constraint: G(U) = 0

and we'll solve this using patternsearch since the limit state function is discrete (0 or 1).

Essentially what FORM is doing is minimizing the distance to the failure region from the origin and using that distance as an approximation for the probability of failure.

objFunction    = @(u) norm(u)
type nonlinConstraint

% solve with patternsearch
u0 = [0,0];
lb = [-5 -5]; ub = [5 5];  % keep from searching too far out from origin
                           % and to prevent a simulation crash
[u,beta,~,info] = patternsearch(objFunction,u0,[],[],[],[],...
tform = toc;
Pform = normcdf(-beta)
matlabpool close
% plot
hold on
legend('Test Points','Failed Points','FORM Distance')
objFunction = 


function [c,ceq] = nonlinConstraint(u,Gx,Tinv)
%% Nonlinear Constratint Function for FORM
% inputs:
%   u  - vector of inputs in normal space
%   Gx - Limit State Function in standard space
%   Tinv - Inverse tranfrom (U --> X)

% Copyright 2010 The MathWorks, Inc.
c = []; % no nonlinear inequality constraints

% input to Tinv must a row vector to eval appropriately
if ~isvector(u)
    error('Input must be a vector')
[row,col] = size(u);
if row>col
    u = u';
x = Tinv(u);
x(:,1) = x(:,1).*10;    % Mach Scaling
x(:,2) = x(:,2)./10000; % Altitude Scaling
ceq = Gx(x);

  Iter   f-count      f(x)      constraint   MeshSize      Method
    0     1            0            1            1    
    1    33       3.0738            0     0.009772   Update multipliers
    2   101      1.88482            0     0.000955   Update multipliers
    3   175      1.88281            0   9.333e-005   Update multipliers
    4   258       1.8821            0    9.12e-006   Update multipliers
    5   344       1.8821            0   8.913e-007   Update multipliers
Optimization terminated: mesh size less than options.TolMesh
 and constraint violation is less than options.TolCon.

u =

   -1.3661    1.2947

beta =


info = 

         function: @(u)norm(u)
      problemtype: 'nonlinearconstr'
       pollmethod: 'gpspositivebasisnp1'
     searchmethod: @searchlhs
       iterations: 5
        funccount: 344
         meshsize: 8.9125e-007
    maxconstraint: 0
          message: [1x115 char]

Pform =


Sending a stop signal to all the labs ... stopped.


Monte Carlo and FORM produce similar estimates of the probablity of failure for this problem. FORM is a linear method which isn't as accurate (in general), but can achieve a reasonable approximation, with fewer function calls. This problem had a discontinous failure surface (0 or 1) which required more function evaluations than would be needed in a continous case.

fprintf('            16 MATLAB Worker Case\n');
fprintf('            Monte Carlo     FORM \n');
fprintf('Pf          %2.4f          %2.4f \n',Pfmc,Pform);
fprintf('Func Evals  %d           %d\n',cvIndex(1),info.funccount);
fprintf('Time(s)     %4.1f          %4.1f \n',tmc,tform);
            16 MATLAB Worker Case
            Monte Carlo     FORM 
Pf          0.0258          0.0299 
Func Evals  3532           344
Time(s)     244.6          86.7 

The trade-off is between accuracy in the estimate and the number of function evaluations required. In many cases, Monte Carlo simulation can be too costly or used only to verify the final design. FORM is cheaper and good for ranking designs early on in the design process and is generally conservative in the failure estimate. The estimate that would be attained from the Monte Carlo simulation with the same number of iterations is shown below.

ans =