MATLAB Examples

Failure Boundary Identification

This file describes the development of a failure boundary identication algorithm shown in "Using Statistics and Optimization to Support Design Activities" Webinar, July 21, 2009.

(C) Copyright 2009 The MathWorks, Inc.


Use Peaks as the Test Function

Peaks has two valleys which we'll use as our failure regions. Let's denote a failure when z <= -1.

z =  3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ... 
   - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ... 
   - 1/3*exp(-(x+1).^2 - y.^2) 

Sample Peaks with a Uniform Grid

Let's start with a uniform grid sampling of the peaks function.

[x,y,z] = peaks(10);
x = x(:); y = y(:); z = z(:);
fail = z < -1; % specify failure
peaks, hold on
grid on, hold off
z =  3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ... 
   - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ... 
   - 1/3*exp(-(x+1).^2 - y.^2) 

Use Clustering to Identify Number of Differnt Regions

Clustering is a technique that segments the data based upon grouping the data with similar characteristics (a distance measure of how similar or differnent the data is).

maxRegions = 3;
zones = clusterdata([x,y,fail],maxRegions);
zoneNumbers = unique(zones)
axis square equal
zoneNumbers =


We have successfully segmented the failure regions from non-failure regions. Now we need to develop a method to sample near the boundary of the failure region to continue to resolve the boundary.

Use Voronoi Diagram to Find Next Points

Voronoi diagrams add cells around each point that separates each point from the others. For this problem, the cells look like rectangles, but thet can be any shape/sized polygon that best separates the points.

We'll use the cell polygons to find points to use that are close to our regions of interest and use them to update our sampling of the peaks function to resolve the failure boundary.

DT = DelaunayTri(x,y);
[V,R] = voronoiDiagram(DT);
hold on

Use Convex Hull to Identify First Failure Region

We need a way to identify the failure regions boundaries. Convex Hulls will encompase the failure points in a polygon. We can then use the polygon to identify the Voronoi cells we can use to update our sampling of the peaks function.

failZones = unique( fail .* zones);
failZones = failZones(failZones > 0)
out = find(~(zones == failZones(1)));
DT.X(out,:) = [];

% Identify convex hull
ch = convexHull(DT);
failZones =


Show the next points to update.

% Make note of any points that are in the convex hull or on it
in = inpolygon(x,y,DT.X(ch,1),DT.X(ch,2));

% Use the points in the convex hull to find the voronoi cells
vc = R(in);

% Use the voronoi cell points as the next points to evaluate
idxNewPts = [];
for j = 1:length(vc)
    idxNewPts = [idxNewPts, vc{j}];
idxNewPts = unique(idxNewPts);
newXY = V(idxNewPts,:);

% Remove any points in/on the convex hull
in = inpolygon(newXY(:,1),newXY(:,2),DT.X(ch,1),DT.X(ch,2));
newXY = newXY(~in,:);

% Show plot of the points we want to update

Note that Convex Hull was used to find the boundary. This assumes the failure region will be convex in shape. freeBoundary could be used instead of Convex Hull if you expect a non-convex shape to the failure region.

Encapsulate Above Process in Function boundarySearch

Repeat the process above to find the boundary for the other failure region, and continue to add points to the boundary until satisfied with the boundary resolution (in this case, a change in perimeter of less than 0.1 is used)

help boundarySearch
 BOUNDARYSEARCH searches a 2-D space to identify the failure region(s)
   BOUNDARYSEARCH(X, LIMITFUNCTION) searches over a 2D space to identify
   regions where the limit function returns (True) 1 if failure is present.
   X0 (nx2) is initial search starting points of (x,y) pairs. LIMITFUNCTION
   is a handle to the limit state function that that takes X as an input.
   Optional input are:
     'showplot'     show plot update of algorithm (voronoi cells and
                    next points to update)
     'ptol'         Perimeter change tolerance, followed by numeric value.
      [xb,P]   = boundarySearch(X0, limitFunction)
      [xb,P,A] = boundarySearch(X0, limitFunction,'showplot','ptol',val)
     [x,y,z] = peaks(10);
     X0 = [x(:) y(:)];
     ub = [3 3];
     lb = [-3 -3];
     [xb P a] = boundarySearch(X0,lb,ub,@peaksLimitFunction,'showplot','ptol',0.1)
   See also clusterdata, DelaunayTri, voronoi, peaksLimitFunction

Peaks Limit Function

type peaksLimitFunction
function failed = peaksLimitFunction(X)

%   Copyright 2009 The MathWorks, Inc.

% Calcualte Peaks Limit Function for z <= -1
z = peaks(X(:,1),X(:,2));
failed = z < -1;
failed = logical(failed);

Run the algoirthm.

[x,y,z] = peaks(10);
X0 = [x(:) y(:)];
ub = [3 3];
lb = [-3 -3];
[xb P a] = boundarySearch(X0,lb,ub,@peaksLimitFunction,'showplot','ptol',0.1)
xb = 

    [34x2 double]    [35x2 double]

P =

    6.1405    4.4354

a =

    2.8461    1.5107

Note that the boundarySearch function takes in a function handle (@) symbol as an arguement. This allows us to use the boundarySearch for an arbitrary function. For example, we'll use this function on a Simulink® Model in statOptimDesign.

Show Boundary on Peaks

Plot the result on top of the peaks function to see how well the algorithm did.

hold off
hold on

col = 'rb';
for i = 1:length(P)
z =  3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ... 
   - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ... 
   - 1/3*exp(-(x+1).^2 - y.^2)