## Loren on the Art of MATLABTurn ideas into MATLAB

Note

Loren on the Art of MATLAB has been archived and will not be updated.

# Optimizing a Structural Design

Today you'll see a new demonstration of applying optimization techniques. Today's guest is Takafumi Ohbiraki. This demonstration was part of the contents of the MATLAB EXPO which was held in Tokyo last year (2016).

### Contents

Optimization is a universal topic in both engineering and data science. It is a technology with high needs. Speaking of optimization in engineering, parameter estimation is often performed on a model showing a physical phenomenon.

Generally, when expressing a physical phenomenon by a model, MATLAB has functions of ordinary differential equations and partial differential equations,and performs simulation. For example, since the problem of heat conduction can be represented by parabolic partial differential equations on the second floor, it is possible to simulate with the function pdepe.

In reality, it is often difficult to convert a physical phenomenon into a closed-form mathematical formula. So it is common to utilize a modeling tool to convert the problem into a mathematical model for analysis.

At this time, we will introduce how to optimize the design by using Partial Differential Equation Toolbox and Optimization Toolbox with an example of structural analysis.

#### Design of an umbrella hook

A cool umbrella hook is always seen in Japanese public restrooms. As a designer of an umbrella hook, we will consider how to reduce the cost by reducing the materials used for an umbrella hook. However, in order to keep the functionality of the umbrella hook, it is necessary to for the hook not to bend.

First, we define the design specifications. We consider opening an elliptical cavity in the member. The material (A or B) of the member, the coordinates (X, Y) of the ellipse, the radius of the ellipse (RX, RY), and the rotation angle of the ellipse (theta) are design parameters. The number of parameters is seven.

We will set the downward load on the tightening portion at the upper right and the fixed end at the left end. We perform an optimization with the aim of minimizing the displacement amount of the lower right point. All parameters have the following constraints. #### FEM model (Using Partial Differential Equation (PDE) Toolbox)

We create a function pdesample to actually simulate with parameters using Partial Differential Equation Toolbox (PDE Toolbox).

in = [1 5,2,1.5,1.0,0.3,0];
lh = pdesample
out = lh{1}(in)

lh =
5×1 cell array
@pdesample_input
@pdesample_input_doe_sn
@myobjfun
@myobjfun2
@mynoncon
out =
0.40131 This is the initial condition.

Output is about 0.4013. We do not understand whether this result is good or bad. We will seek a good answer using optimization techniques.

#### Process of Optimization

To actually perform the optimization calculation, if there are too many parameters, it takes time to reach the optimum value. This time, we will utilize a two-step optimization method that is often used in quality engineering. According to the two-step design method, first, make efforts to reduce the number of variables to be optimized by utilizing the experimental design method. In this method, we adopt the Taguchi design. The orthogonal table is uploaded to the File Exchange. In variable selection, we paid attention to the signal to noise ratio (S/N) ratio and the nominal-is-best response. In the second step, we use optimization algorithms to find good parameters. #### Design of Experiments (DOE) (Taguchi design)

We will exploit which parameters should be reduced by utilizing the orthogonal array famous from the Taguchi design. In the orthogonal array of Taguchi design, there are 2 factors, 3 factors, and mixed orthogonal tables. We chose the famous L18 model for mixing both 2 factors and 3 factors.

L18 model : 18 experiments to make the combination of 8 variables efficient; one variable (material) is 2 levels , seven variables are 3 levels. The variable design created by my helper function orthogonalarray is a MATLAB table.

design = orthogonalarray('L18','Bounds',[1 2;2.0 8.0;1.0 3.0;0.5 1.9;0.2 0.9;0.1 0.3;-pi/12 pi/12;0 0],...
'VariableNames',{'material','x','y','rx','ry','dh','theta','xxx'})

design =
18×8 table
material    x    y    rx      ry     dh      theta     xxx
________    _    _    ___    ____    ___    _______    ___
1           2    1    0.5     0.2    0.1    -0.2618    0
1           2    2    1.2    0.55    0.2          0    0
1           2    3    1.9     0.9    0.3     0.2618    0
1           5    1    0.5    0.55    0.2     0.2618    0
1           5    2    1.2     0.9    0.3    -0.2618    0
1           5    3    1.9     0.2    0.1          0    0
1           8    1    1.2     0.2    0.3          0    0
1           8    2    1.9    0.55    0.1     0.2618    0
1           8    3    0.5     0.9    0.2    -0.2618    0
2           2    1    1.9     0.9    0.2          0    0
2           2    2    0.5     0.2    0.3     0.2618    0
2           2    3    1.2    0.55    0.1    -0.2618    0
2           5    1    1.2     0.9    0.1     0.2618    0
2           5    2    1.9     0.2    0.2    -0.2618    0
2           5    3    0.5    0.55    0.3          0    0
2           8    1    1.9    0.55    0.3    -0.2618    0
2           8    2    0.5     0.9    0.1          0    0
2           8    3    1.2     0.2    0.2     0.2618    0


#### Calculation with setting parameters for L18

At experimental points, we calculate the outputs (strain) and S/N ratios.

output = zeros(18,1);
sn = zeros(18,1);

for i = 1:18
output(i,1) = lh{1}(table2array(design(i,1:8)));
sn(i,1) = lh{2}(table2array(design(i,1:8)));
end

results = [design table(output,sn)]

results =
18×10 table
material    x    y    rx      ry     dh      theta     xxx    output       sn
________    _    _    ___    ____    ___    _______    ___    _______    ______
1           2    1    0.5     0.2    0.1    -0.2618    0        0.346    108.16
1           2    2    1.2    0.55    0.2          0    0      0.36042    88.027
1           2    3    1.9     0.9    0.3     0.2618    0      0.42694    57.447
1           5    1    0.5    0.55    0.2     0.2618    0      0.35443    88.063
1           5    2    1.2     0.9    0.3    -0.2618    0      0.37875    78.837
1           5    3    1.9     0.2    0.1          0    0      0.36595    84.527
1           8    1    1.2     0.2    0.3          0    0      0.33836     121.1
1           8    2    1.9    0.55    0.1     0.2618    0      0.38466    59.585
1           8    3    0.5     0.9    0.2    -0.2618    0      0.34978    96.122
2           2    1    1.9     0.9    0.2          0    0       1.8659    32.802
2           2    2    0.5     0.2    0.3     0.2618    0      0.92894    132.08
2           2    3    1.2    0.55    0.1    -0.2618    0      0.97881    90.175
2           5    1    1.2     0.9    0.1     0.2618    0       1.2566    41.031
2           5    2    1.9     0.2    0.2    -0.2618    0       1.0356     78.99
2           5    3    0.5    0.55    0.3          0    0      0.94161       113
2           8    1    1.9    0.55    0.3    -0.2618    0        1.019     59.27
2           8    2    0.5     0.9    0.1          0    0      0.95126    106.67
2           8    3    1.2     0.2    0.2     0.2618    0      0.95404    88.904 #### Factor Effect Chart

The factor effect diagram is a graph expressing the effect of each factor, with the horizontal axis representing the level of the factor and the vertical axis representing the S / N ratio and the characteristic value. From this graph, find the condition with high S / N ratio and select the one with better characteristics.

factoreffectchart is my helper function.

#### S/N ratio

According to this chart, the S/N ratio is very good at both the 1st level of the 4th and 5th parameters (or variables).

The 4th variable is the x radius and the 5th variables is the y radius.

d = orthogonalarray('L18');
Osn = factoreffectchart(d,results.sn) % S/N

Osn =
86.873       82.546          NaN
84.781       80.741       88.608
75.07       90.698       88.362
107.35       84.679       62.103
102.29       83.019       68.818
81.69       78.818       93.621
85.258        91.02       77.851
80.213       83.509       90.407 #### Output

To minimize the strain, the 1st level of the 1st variable is good.

The 1st variable is material.

We succeed at reducing the parameters from 7 variables to 4.

O = factoreffectchart(d,results.output) % Output

O =
0.36726       1.1035          NaN
0.81783      0.72215      0.66618
0.86336      0.67327      0.66952
0.64534      0.71116      0.84967
0.66149      0.67315      0.87153
0.71387      0.82003      0.67226
0.68465      0.80391       0.7176
0.81182      0.71343       0.6809 #### Optimization using patternsearch solver

We will optimize 4 parameters using the patternsearch solver.

First of all, we are going to optimize the location of the hole that minimizes strain.

The objective function is the code from 192 to 197 lines in "pdesample.m".

dbtype pdesample 193:197

close all

193   function out = myobjfun(in)
194
195   out = pdesample_input([1 in(1) in(2) 0.5 0.2 in(3) in(4)]);
196
197   end

opt = optimoptions('patternsearch');
opt.MeshTolerance = 0.001;
opt.FunctionTolerance = 0.00001;
opt.MaxIterations = 500;
opt.TimeLimit = 10*60;
opt.UseParallel = true;
opt.Display = 'iter';
[x,obj] = patternsearch(lh{3}, [7,3,0.2,0],[],[],[],[],[2.0,0.3,0.1,-pi/4],[9.5,4.5,0.3,pi/4],lh{5},opt)

                                      Max
Iter   Func-count       f(x)      Constraint   MeshSize      Method
0         1     0.338279            0            1
1         6      0.33769            0     0.009772   Update multipliers
2       329     0.334631            0     0.000955   Update multipliers
Optimization terminated: mesh size less than options.MeshTolerance
and constraint violation is less than options.ConstraintTolerance.
x =
9.3125       2.9688      0.14141         0.25
obj =
0.33463 #### Optimization using ga solver

We will also optimize 4 parameters using the ga> solver.

gaoptions = gaoptimset(@ga);
gaoptions.CrossoverFcn = @crossoverintermediate;
gaoptions.CrossoverFraction = 0.5;
gaoptions.TolFun = 0.00001;
gaoptions.PopulationSize = 20; % population size

gaoptions.InitialPopulation = [7,3,0.2,0]; % initial population
gaoptions.Generations = 500; %
gaoptions.TimeLimit = 10*60;
gaoptions.UseParallel = true;
gaoptions.Display = 'iter';

[x,fval]  = ga(lh{3},4,[],[],[],[],[2.0,0.3,0.1,-pi/4],[9.5,4.5,0.3,pi/4],lh{5},gaoptions)

                              Best       Max        Stall
Generation  Func-count        f(x)     Constraint  Generations
1           1078       0.33505            0      0
2           1598      0.334569    0.0008296      0
Optimization terminated: average change in the fitness value less than options.FunctionTolerance
and constraint violation is less than options.ConstraintTolerance.
x =
9.3345       3.1148      0.19511      0.36912
fval =
0.33457 The values of the output by both the patternsearch and ga solvers are reduced from the inital conditions. Both results give especially similiar results for the x,y positions. We conclude that both results are good. However, the area of hole is the smallest possible because of considering the S/N ratio. In order to lighten the weight, we will enlarge the area of hole.

#### Multi-objective optimization using the gamultiobj solver

We would like to optimize both the size and location of the hole in order to minmize the strain and lighten the weight.

In the previous case, radius x and y are the minimum values. In order to increase the size of the hole, we will use a multi-objective optimization gamultiobj solver.

The relationship between the value of strain and the area of hole is inverse proportional and there is a trade-off.

The objective function is the code from 200 to 205 lines in "pdesample.m".

dbtype pdesample 200:205

200   function out = myobjfun2(in)
201
202   out(1) = pdesample_input([1 in 0.2 0*pi/12]);
203   out(2) = -pi*in(3)*in(4);
204
205   end

mulgaoptions = gaoptimset(@gamultiobj);
mulgaoptions.CrossoverFraction = 0.5;
mulgaoptions.TolFun = 0.001;
mulgaoptions.PopulationSize = 100;
mulgaoptions.Generations = 3000;
mulgaoptions.InitialPopulation = [7 2.5 0.5 0.5];
mulgaoptions.UseParallel = false;
mulgaoptions.TimeLimit = 20*60;
[x,fval,exitflag,output]  = gamultiobj(lh{4},4,[1 0 1 0;0 1 0 1;-1 0 1 0; 0 -1 0 1],[9.9;4.9;-0.2;-0.2],[],[],[2.0,0.3,0.2,0.2],[9.6,4.5,3,2.6],mulgaoptions);

Optimization terminated: time limit exceeded. Here is the biggest hole.

lh{4}(x(1,:))

ans =
0.33632     -0.62854 In this next plot, we show the pareto frontier which is popular for design engineering. There you can see the trade-off between area and strain for making design decisions.

scatter(fval(:,1),50.4 + fval(:,2),'m')
set(gca,'FontSize',16,'FontName','Meiryo UI')
xlabel('strain');ylabel('area');title('pareto frontier') Pareto frontier results.

T1 = array2table([fval(:,1) 50.4+fval(:,2) x],'VariableNames',{'Strain','Area','x','y','rx','ry'})

T1 =
35×6 table
Strain      Area       x          y         rx         ry
_______    ______    ______    _______    _______    _______
0.33632    49.771    8.4674    0.83566    0.54447    0.36746
2.9137    29.829    6.1874     2.4294     2.9804      2.197
0.8458    36.385    6.1587     2.0938     2.7957     1.5957
3.9636    29.123    6.1812      2.593     2.9923     2.2634
1.553    32.151     5.883     2.3126     2.9526     1.9674
0.83451    37.023    6.7695     2.1714     2.4951     1.7065
0.42155    44.452     7.107     2.1224     1.3364     1.4167
1.2357    33.732    6.0206     2.5458      2.879     1.8429
0.93663    35.477    6.1315     2.3105     2.9743     1.5971
0.75125    37.875    6.7413     2.0409     2.5471     1.5653
1.0758     34.69     6.385     2.4479     2.8222     1.7719
3.077    29.714    6.1778     2.4877     2.9777     2.2113
2.8688    29.875     6.022     2.4611     2.9817     2.1911
0.39498    45.353    7.3136     1.8261     1.4593     1.1008
1.1526    34.152    6.2626     2.1237     2.9141     1.7747
0.35257     47.25    8.2497     1.2095     1.0588    0.94702
3.4396    29.299    6.3585     2.4648          3     2.2389
1.4782    32.434    6.1279     2.4721     2.9814     1.9181
2.1484    30.683    5.9299     2.4064     2.9995     2.0924
0.50403    41.862    7.0168     2.0862     1.8298     1.4853
2.3041    30.529    5.9573     2.4324     2.9834     2.1201
1.6972    31.718    5.9601     2.4295     2.9809      1.995
3.2141    29.467    6.1923     2.4677          3      2.221
0.34045    48.795    8.1318     1.0284     0.8714    0.58627
1.0626    35.115    6.6304     2.3413     2.6315     1.8489
1.3309    33.168     6.023     2.5945     2.9801     1.8405
0.38738    45.534    7.4165     1.5702     1.5472     1.0012
1.4303    32.719    5.9897      2.279     2.8842     1.9514
0.55291    40.816    5.8008     1.8469      2.672     1.1418
1.8278    31.425    5.9551      2.491     2.9899     2.0201
0.3506    48.091    8.7273     1.5772    0.62704      1.172
1.9699     31.13    6.0712     2.4664       2.98     2.0583
3.7843    29.133    6.0624     2.4919     2.9961     2.2595
2.7137    30.065     6.191     2.4831     2.9761     2.1749
3.6181    29.252    5.7451     2.5008     2.9962     2.2467


#### Conclusion

We showed how to use FEM tools and Optimization Toolbox to perform structural design. What can you shape with these tools? Let us know your thoughts here.

Published with MATLAB® R2017a