design optimization of a welded beam -凯发k8网页登录

this example shows how to examine the tradeoff between the strength and cost of a beam. several publications use this example as a test problem for various multiobjective algorithms, including deb et al. [1] and ray and liew [2].

for a video overview of this example, see .

problem description

the following sketch is adapted from ray and liew [2].

this sketch represents a beam welded onto a substrate. the beam supports a load p at a distance l from the substrate. the beam is welded onto the substrate with upper and lower welds, each of length l and thickness h. the beam has a rectangular cross-section, width b, and height t. the material of the beam is steel.

the two objectives are the fabrication cost of the beam and the deflection of the end of the beam under the applied load p. the load p is fixed at 6,000 lbs, and the distance l is fixed at 14 in.

the design variables are:

x(1) = h, the thickness of the welds
x(2) = l, the length of the welds
x(3) = t, the height of the beam
x(4) = b, the width of the beam

the fabrication cost of the beam is proportional to the amount of material in the beam, $(l l) t b$ , plus the amount of material in the welds, $l h^{2}$ . using the proportionality constants from the cited papers, the first objective is

$f_{1} (x) = 1.10471 x_{1}^{2} x_{2} 0.04811 x_{3} x_{4} (14 x_{2}) .$

the deflection of the beam is proportional to p and inversely proportional to $b t^{3}$ . again, using the proportionality constants from the cited papers, the second objective is

$f_{2} (x) = \frac{p}{x_{4} x_{3}^{3}} c$ , where $c = \frac{4 (14)^{3}}{30 \times 1 0^{6}} \approx 3.6587 \times 1 0^{- 4}$ and $p = 6, 000$ .

the problem has several constraints.

the weld thickness cannot exceed the beam width. in symbols, x(1) <= x(4). in toolbox syntax:

aineq = [1,0,0,-1];
bineq = 0;

the shear stress $τ (x)$ on the welds cannot exceed 13,600 psi. to calculate the shear stress, first calculate preliminary expressions:

$τ_{1} = \frac{1}{\sqrt{2} x_{1} x_{2}}$

$r = \sqrt{x_{2}^{2} (x_{1} x_{3})^{2}}$

$τ_{2} = \frac{(l x_{2} / 2) r}{\sqrt{2} x_{1} x_{3} (x_{2}^{2} / 3 (x_{1} x_{3})^{2}}$

$τ (x) = p \sqrt{τ_{1}^{2} τ_{2}^{2} \frac{2 τ_{1} τ_{2} x_{2}}{r}} .$

in summary, the shear stress on the welds has the constraint $τ (x)$ <= 13600.

the normal stress $σ (x)$ on the welds cannot exceed 30,000 psi. the normal stress is $p \frac{6 l}{x_{4} x_{3}^{2}} \leq 30 \times 1 0^{3}$ .
the buckling load capacity in the vertical direction must exceed the applied load of 6,000 lbs. using the values of young's modulus $e = 30 \times 1 0^{6}$ psi and $g = 12 \times 1 0^{6}$ psi, the buckling load constraint is $\frac{4.013 e x_{3} x_{4}^{3}}{6 l^{2}} (1 - \frac{x_{3}}{2 l} \sqrt{\frac{e}{4 g}}) \geq 6000$ . numerically, this becomes the inequality $64, 746.022 (1 - 0.0282346 x_{3}) x_{3} x_{4}^{3} \geq 6000$ .
the bounds on the variables are 0.125 <=x(1) <= 5, 0.1 <= x(2) <= 10, 0.1 <= x(3) <= 10, and 0.125 <= x(4) <= 5. in toolbox syntax:

lb = [0.125,0.1,0.1,0.125];
ub = [5,10,10,5];

the objective functions appear at the end of this example in the function objval(x). the nonlinear constraints appear at the end of this example in the function nonlcon(x).

multiobjective problem formulation and `paretosearch` solution

you can optimize this problem in several ways:

set a maximum deflection, and find a single-objective minimal fabrication cost over designs that satisfy the maximum deflection constraint.
set a maximum fabrication cost, and find a single-objective minimal deflection over designs that satisfy the fabrication cost constraint.
solve a multiobjective problem, visualizing the tradeoff between the two objectives.

to use the multiobjective approach, which gives more information about the problem, set the objective function and nonlinear constraint function.

fun = @objval;
nlcon = @nonlcon;

solve the problem using paretosearch with the 'psplotparetof' plot function. to reduce the amount of diagnostic display information, set the display option to 'off'.

opts_ps = optimoptions('paretosearch','display','off','plotfcn','psplotparetof');
rng default % for reproducibility
[x_ps1,fval_ps1,~,psoutput1] = paretosearch(fun,4,aineq,bineq,[],[],lb,ub,nlcon,opts_ps);

disp("total function count: "   psoutput1.funccount);

total function count: 1870

for a smoother pareto front, try using more points.

npts = 160; % the default is 60
opts_ps.paretosetsize = npts;
[x_ps2,fval_ps2,~,psoutput2] = paretosearch(fun,4,aineq,bineq,[],[],lb,ub,nlcon,opts_ps);

disp("total function count: "   psoutput2.funccount);

total function count: 6254

this solution looks like a smoother curve, but it has a smaller extent of objective 2. the solver takes over three times as many function evaluations when using 160 pareto points instead of 60.

`gamultiobj` solution

to see if the solver makes a difference, try the gamultiobj solver on the problem. set equivalent options as in the previous solution. because the gamultiobj solver keeps fewer than half of its solutions on the best pareto front, use two times as many points as before.

opts_ga = optimoptions('gamultiobj','display','off','plotfcn','gaplotpareto','populationsize',2*npts);
[x_ga1,fval_ga1,~,gaoutput1] = gamultiobj(fun,4,aineq,bineq,[],[],lb,ub,nlcon,opts_ga);

disp("total function count: "   gaoutput1.funccount);

total function count: 38401

gamultiobj takes tens of thousands of function evaluations, whereas paretosearch takes only thousands.

compare solutions

the gamultiobj solution seems to differ from the paretosearch solution, although it is difficult to tell because the plotted scales differ. plot the two solutions on the same plot, using the same scale.

fps2 = sortrows(fval_ps2,1,'ascend');
figure
hold on
plot(fps2(:,1),fps2(:,2),'r-')
fga = sortrows(fval_ga1,1,'ascend');
plot(fga(:,1),fga(:,2),'b--')
xlim([0,40])
ylim([0,1e-2])
legend('paretosearch','gamultiobj')
xlabel 'cost'
ylabel 'deflection'
hold off

the gamultiobj solution is better in the rightmost portion of the plot, whereas the paretosearch solution is better in the leftmost portion. paretosearch uses many fewer function evaluations to obtain its solution.

typically, when the problem has no nonlinear constraints, paretosearch is at least as accurate as gamultiobj. however, the resulting pareto sets can have somewhat different ranges. in this case, the presence of a nonlinear constraint causes the paretosearch solution to be less accurate over part of the range.

one of the main advantages of paretosearch is that it usually takes many fewer function evaluations.

start from single-objective solutions

to help the solvers find better solutions, start them from points that are the solutions to minimizing the individual objective functions. the pickindex function returns a single objective from the objval function. use fmincon to find single-objective optima. then use those solutions as initial points for a multiobjective search.

x0 = zeros(2,4);
x0f = (lb   ub)/2;
opts_fmc = optimoptions('fmincon','display','off','maxfunctionevaluations',1e4);
x0(1,:) = fmincon(@(x)pickindex(x,1),x0f,aineq,bineq,[],[],lb,ub,@nonlcon,opts_fmc);
x0(2,:) = fmincon(@(x)pickindex(x,2),x0f,aineq,bineq,[],[],lb,ub,@nonlcon,opts_fmc);

examine the single-objective optima.

objval(x0(1,:))

ans = 1×2
    2.3810    0.0158

objval(x0(2,:))

ans = 1×2
   76.7188    0.0004

the minimum cost is 2.381, which gives a deflection of 0.158. the minimum deflection is 0.0004, which has a cost of 76.7253. the plotted curves are quite steep near the ends of their ranges, meaning you get much less deflection if you take a cost a bit above its minimum, or much less cost if you take a deflection a bit above its minimum.

start paretosearch from the single-objective solutions. because you will plot the solutions later on the same plot, remove the paretosearch plot function.

opts_ps.initialpoints = x0;
opts_ps.plotfcn = [];
[x_psx0,fval_ps1x0,~,psoutput1x0] = paretosearch(fun,4,aineq,bineq,[],[],lb,ub,nlcon,opts_ps);
disp("total function count: "   psoutput1x0.funccount);

total function count: 4839

start ga from the same initial points, and remove its plot function.

opts_ga.initialpopulationmatrix = x0;
opts_ga.plotfcn = [];
[~,fval_ga,~,gaoutput] = gamultiobj(fun,4,aineq,bineq,[],[],lb,ub,nlcon,opts_ga);
disp("total function count: "   gaoutput.funccount);

total function count: 37441

plot the solutions on the same axes.

fps = sortrows(fval_ps1x0,1,'ascend');
figure
hold on
plot(fps(:,1),fps(:,2),'r-')
fga = sortrows(fval_ga,1,'ascend');
plot(fga(:,1),fga(:,2),'b--')
xlim([0,40])
ylim([0,1e-2])
legend('paretosearch','gamultiobj')
xlabel 'cost'
ylabel 'deflection'
hold off

by starting from the single-objective solutions, the gamultiobj solution is slightly better than the paretosearch solution throughout the plotted range. however, gamultiobj takes almost ten times as many function evaluations to reach its solution.

hybrid function

gamultiobj can call the hybrid function fgoalattain automatically to attempt to reach a more accurate solution. see whether the hybrid function improves the solution.

opts_ga.hybridfcn = 'fgoalattain';
[xgah,fval_gah,~,gaoutputh] = gamultiobj(fun,4,aineq,bineq,[],[],lb,ub,nlcon,opts_ga);
disp("total function count: "   gaoutputh.funccount);

total function count: 57478

fgah = sortrows(fval_gah,1,'ascend');
figure
hold on
plot(fps(:,1),fps(:,2),'r-')
plot(fga(:,1),fga(:,2),'b--')
plot(fgah(:,1),fgah(:,2),'g-')
xlim([0,40])
ylim([0,1e-2])
legend('paretosearch','gamultiobj','gamultiobj/fgoalattain')
xlabel 'cost'
ylabel 'deflection'
hold off

the hybrid function provides a slight improvement on the gamultiobj solution, mainly in the leftmost part of the plot.

run `fgoalattain` manually from `paretosearch` solution points

although paretosearch has no built-in hybrid function, you can improve the paretosearch solution by running fgoalattain from the paretosearch solution points. create a goal and weights for fgoalattain by using the same setup for fgoalattain as described in gamultiobj hybrid function.

fmax = max(fval_ps1x0);
nobj = numel(fmax);
fmin = min(fval_ps1x0);
w = sum((fmax - fval_ps1x0)./(1   fmax - fmin),2);
p = w.*((fmax - fval_ps1x0)./(1   fmax - fmin));
xnew = zeros(size(x_psx0));
nsol = size(xnew,1);
fvalnew = zeros(nsol,nobj);
opts_fg = optimoptions('fgoalattain','display','off');
nfv = 0;
for ii = 1:nsol
    [xnew(ii,:),fvalnew(ii,:),~,~,output] = fgoalattain(fun,x_psx0(ii,:),fval_ps1x0(ii,:),p(ii,:),...
        aineq,bineq,[],[],lb,ub,nlcon,opts_fg);
    nfv = nfv   output.funccount;
end
disp("fgoalattain function count: "   nfv)

fgoalattain function count: 14049

fnew = sortrows(fvalnew,1,'ascend');
figure
hold on
plot(fps(:,1),fps(:,2),'r-')
plot(fga(:,1),fga(:,2),'b--')
plot(fgah(:,1),fgah(:,2),'g-')
plot(fnew(:,1),fnew(:,2),'k.-')
xlim([0,40])
ylim([0,1e-2])
legend('paretosearch','gamultiobj','gamultiobj/fgoalattain','paretosearch/fgoalattain')
xlabel 'cost'
ylabel 'deflection'

the combination of paretosearch and fgoalattain creates the most accurate pareto front. zoom in to see.

xlim([3.64 13.69])
ylim([0.00121 0.00442])
hold off

even with the extra fgoalattain computations, the total function count for the combination is less than half of the function count for the gamultiobj solution alone.

fprintf("total function count for gamultiobj alone is %d.\n"   ...
    "for paretosearch and fgoalattain together it is %d.\n",...
    gaoutput.funccount,nfv   psoutput1x0.funccount)

total function count for gamultiobj alone is 37441.
for paretosearch and fgoalattain together it is 18888.

find good parameters from plot

the plotted points show the best values in function space. to determine which parameters achieve these function values, find the size of the beam and size of the weld in order to get a particular cost/deflection point. for example, the plot of paretosearch followed by fgoalattain shows points with a cost of about 6 and a deflection of about 3.5e–3. determine the sizes of the beam and weld that achieve these points.

whichgood = find(fvalnew(:,1) <= 6 & fvalnew(:,2) <= 3.5e-3);
goodpoints = table(xnew(whichgood,:),fvalnew(whichgood,:),'variablenames',{'parameters' 'objectives'})

goodpoints=4×2 table
                   parameters                       objectives     
    ________________________________________    ___________________
    0.63457     1.5187         10    0.67262    5.6974    0.0032637
    0.61635     1.5708         10    0.63165     5.391    0.0034753
    0.63228     1.5251         10     0.6674    5.6584    0.0032892
    0.65077     1.4751         10    0.70999     5.976    0.0030919

four sets of parameters achieve a cost of less than 6 and a deflection of less than 3.5e–3:

weld thickness slightly over 0.6
weld length about 1.5
beam height 10 (the upper bound)
beam width between 0.63 and 0.71

objective and nonlinear constraints

function [cineq,ceq] = nonlcon(x)
sigma = 5.04e5 ./ (x(:,3).^2 .* x(:,4));
p_c = 64746.022*(1 - 0.028236*x(:,3)).*x(:,3).*x(:,4).^3;
tp = 6e3./sqrt(2)./(x(:,1).*x(:,2));
tpp = 6e3./sqrt(2) .* (14 0.5*x(:,2)).*sqrt(0.25*(x(:,2).^2   (x(:,1)   x(:,3)).^2)) ./ (x(:,1).*x(:,2).*(x(:,2).^2 / 12   0.25*(x(:,1)   x(:,3)).^2));
tau = sqrt(tp.^2   tpp.^2   (x(:,2).*tp.*tpp)./sqrt(0.25*(x(:,2).^2   (x(:,1)   x(:,3)).^2)));
cineq = [tau - 13600,sigma - 3e4,6e3 - p_c];
ceq = [];
end
function f = objval(x)
f1 = 1.10471*x(:,1).^2.*x(:,2)   0.04811*x(:,3).*x(:,4).*(14.0 x(:,2));
f2 = 2.1952./(x(:,3).^3 .* x(:,4));
f = [f1,f2];
end
function z = pickindex(x,k)
    z = objval(x); % evaluate both objectives
    z = z(k); % return objective k
end

references

[1] deb, kalyanmoy, j. sundar, udaya bhaskara rao n, and shamik chaudhuri. reference point based multi-objective optimization using evolutionary algorithms. international journal of computational intelligence research, vol. 2, no. 3, 2006, pp. 273–286. available at

[2] ray, t., and k. m. liew. a swarm metaphor for multiobjective design optimization. engineering optimization 34, 2002, pp.141–153.

design optimization of a welded beam -凯发k8网页登录

problem description

multiobjective problem formulation and paretosearch solution

gamultiobj solution