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: http://www.mathworks.com/matlabcentral/fileexchange/3938-flutter-analysis-model

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

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; sim(model) make2plots(sigsOut.plunge.Time,sigsOut.plunge.Data,sigsOut.pitch.Data) flutter = sigsOut.isStable.Data(end) 
flutter = 1 

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

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.

figure ax1 = subplot(2,1,1); x = (mu(1)-3*sigma(1)):sigma(1)/10:(mu(1)+3*sigma(1)); plot(x,pdf(M,x),'r') xlabel('Mach Number'), ylabel('Probability Density') ax2 = subplot(2,1,2); x = (m-3*sqrt(v)):sqrt(v)/10:(m+3*sqrt(v)); plot(x,pdf(A,x),'r') 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); setup s = sim(model,'SrcWorkspace','current'); sigsOut = get(s,'sigsOut'); %% Extract Failure Measure try failed = sigsOut.isStable.Data(end); catch disp('Stopped here because of error') end 

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 = @(x)1-flutterLimitFunction(x,model) 

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

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

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 = @(x)1-flutterLimitFunction(x,model) 

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); figure scatterhist(xM,xA) 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 = []; figure files = dir; matlabpool('open','speedy-cluster',16,... 'FileDependencies',{files(3:end).name}) % use 4 core machine cv = 1; block = 1; tic 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]); end 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 subplot(2,1,1) hist(mcG) 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)]) subplot(2,1,2) plot(N,mpf,'r') hold on plot(N,mpf-1.96*spf,'b-',N,mpf + 1.96*spf,'b-') hold off legend('Probability Estimate','Lower Bound','Upper Bound',... 'Location','Best') title('Monte Carlo Convergence with 95% confidence interval') drawnow % protect against premature convergence due to low sample size if N(end) <= blkSize cv( cv(1:end) <= 0.10) = 1; end end tmc = toc; axLim = axis; cvIndex = find(cv <= 0.10); hold on plot([cvIndex(1),cvIndex(1)],axLim(3:4),'k--') text(cvIndex(1),0.5*mean(axLim(3:4)),... {['CV = 0.10, pf = ',num2str(pf(cvIndex(1)),4),'\rightarrow'],... ['Iterations = ',num2str(cvIndex(1)),'\rightarrow']},... 'HorizontalAlignment','right'); legend('Probability Estimate','Lower Bound','Upper Bound',... 'Location','Best') hold off Pfmc = pf(cvIndex(1)) 
Starting matlabpool using the 'speedy-cluster' configuration ... connected to 16 labs. Pfmc = 0.0258 

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 figure index = 1:length(pf); X = [xM(index)' xA(index)']; U = T(X(logical(mcG),:)); Uf = T(X(~mcG,:)); scatterhist(U(:,1),U(:,2)) axis equal xlabel('Normalized Mach Number'),ylabel('Normalized Altitude') title('Standard Normal Space') hold on plot(Uf(:,1),Uf(:,2),'ro') 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 tic [u,beta,~,info] = patternsearch(objFunction,u0,[],[],[],[],... lb,ub,@(u)nonlinConstraint(u,G,Tinv),... psoptimset('TolFun',1e-3,'Disp','iter',... 'UseParallel','always','CompletePoll','on',... 'PollMethod','GPSPositiveBasisNp1',... 'Vectorized','off','Search',{@searchlhs})) tform = toc; Pform = normcdf(-beta) matlabpool close % plot hold on plot([0,u(1)],[0,u(2)],'g-o') legend('Test Points','Failed Points','FORM Distance') 
objFunction = @(u)norm(u) 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') end [row,col] = size(u); if row>col u = u'; end x = Tinv(u); x(:,1) = x(:,1).*10; % Mach Scaling x(:,2) = x(:,2)./10000; % Altitude Scaling ceq = Gx(x); max 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 = 1.8821 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 = 0.0299 Sending a stop signal to all the labs ... stopped. 

Results

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.

pf(info.funccount) 
ans = 0.0291