Code covered by the BSD License

# Automated Failure Boundary Mapping

### Stuart Kozola (view profile)

18 Aug 2009 (Updated )

Demo files from July 21, 2009 webinar

Failure Boundary Identification

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

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

```[x,y,z] = peaks(10);
x = x(:); y = y(:); z = z(:);
fail = z < -1; % specify failure
peaks, hold on
plot3(x,y,z,'ro','MarkerFaceColor','r')
view(2)
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)
scatter(x,y,50,zones,'filled')
axis square equal
```
```zoneNumbers =

1
2
3

```

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
voronoi(x,y)
```

## 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);
plot(DT.X(ch,1),DT.X(ch,2),'r')
```
```failZones =

1
3

```

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}];
end
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
plot(newXY(:,1),newXY(:,2),'r+')
```

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.

SYNTAX
[xb,P]   = boundarySearch(X0, limitFunction)
[xb,P,A] = boundarySearch(X0, limitFunction,'showplot','ptol',val)

Example
[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)

```

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
peaks
hold on
view(2);

col = 'rb';
for i = 1:length(P)
plot(xb{i}(:,1),xb{i}(:,2),col(i),'LineWidth',2)
end
```
```
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)

```