このページは機械翻訳を使用して翻訳されました。元の英語を参照するには、ここをクリックします。

溶接梁の設計最適化

この例では、梁の強度とコストのトレードオフを調べる方法を示します。Debら[1]やRayとLiew[2]など、いくつかの出版物では、この例をさまざまな多目的アルゴリズムのテスト問題として使用しています。

この例のビデオ概要については、多目的最適化のためのパレートセットを参照してください。

問題の説明

以下のスケッチはRayとLiew[2]から改変したものである。

このスケッチは、基板上に溶接された梁を表しています。ビームは基板から距離 L にある荷重 P を支えます。梁は、それぞれ長さ l と厚さ h の上部溶接部と下部溶接部によって基板に溶接されます。梁の断面は長方形で、幅は b、高さは t です。梁の材質は鋼鉄です。

2 つの目的は、梁の製造コストと、適用された荷重 P の下での梁の端部のたわみです。荷重 P は 6,000 ポンドに固定され、距離 L は 14 インチに固定されています。

設計変数は次のとおりです。

x(1) = h、溶接部の厚さ
x(2) = l、溶接部の長さ
x(3) = t、梁の高さ
x(4) = b、ビームの幅

梁の製造コストは、梁の材料の量 $(l + L) t b$ と溶接部の材料の量 $l h^{2}$ の合計に比例します。引用した論文の比例定数を用いて、最初の目的は

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

ビームのたわみはPに比例し、 $b t^{3}$ に反比例します。再び、引用した論文の比例定数を使用して、2番目の目的は

$F_{2} (x) = \frac{P}{x_{4} x_{3}^{3}} C$ 、ここで $C = \frac{4 (14)^{3}}{30 \times 1 0^{6}} \approx 3.6587 \times 1 0^{- 4}$ と $P = 6, 000$ です。

問題にはいくつかの制約があります。

溶接厚さは梁の幅を超えることはできません。記号では、x(1) <= x(4) です。ツールボックス構文の場合:

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

溶接部のせん断応力 $τ (x)$ は 13,600 psi を超えることはできません。せん断応力を計算するには、まず予備的な式を計算します。

$τ_{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}} .$

要約すると、溶接部のせん断応力には制約 $τ (x)$ <= 13600 があります。

溶接部の通常の応力 $σ (x)$ は 30,000 psi を超えることはできません。通常の応力は $P \frac{6 L}{x_{4} x_{3}^{2}} \leq 30 \times 1 0^{3}$ です。
垂直方向の座屈荷重容量は、適用荷重 6,000 ポンドを超える必要があります。ヤング率 $E = 30 \times 1 0^{6}$ psi と $G = 12 \times 1 0^{6}$ psi の値を使用すると、座屈荷重制約は $\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$ になります。数値的には、これは不等式 $64, 746.022 (1 - 0.0282346 x_{3}) x_{3} x_{4}^{3} \geq 6000$ になります。
変数の境界は0.125 <= x(1) <= 5、0.1 <= x(2) <= 10、0.1 <= x(3) <= 10、0.125 <= x(4) <= 5です。ツールボックス構文の場合:

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

目的関数は、この例の最後にある関数 objval(x) に現れます。非線形制約は、この例の最後の関数 nonlcon(x) に表示されます。

多目的問題の定式化と`paretosearch`解決

この問題を最適化するには、いくつかの方法があります。

最大たわみを設定し、最大たわみ制約を満たす設計に対して単一目的の最小製造コストを見つけます。
最大製造コストを設定し、製造コストの制約を満たす設計に対して単一目的の最小たわみを見つけます。
2 つの目的間のトレードオフを視覚化しながら、多目的問題を解決します。

問題についてのより多くの情報を提供する多目的アプローチを使用するには、目的関数と非線形制約関数を設定します。

fun = @objval;
nlcon = @nonlcon;

'psplotparetof' プロット関数と paretosearch を使用して問題を解決します。診断表示情報の量を減らすには、Display オプションを '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);

Figure paretosearch contains an axes object. The axes object with title Pareto Front, xlabel Objective 1, ylabel Objective 2 contains an object of type scatter.

disp("Total Function Count: " + psoutput1.funccount);

Total Function Count: 1467

より滑らかなパレート面を実現するには、より多くのポイントを使用してみてください。

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

Figure paretosearch contains an axes object. The axes object with title Pareto Front, xlabel Objective 1, ylabel Objective 2 contains an object of type scatter.

disp("Total Function Count: " + psoutput2.funccount);

Total Function Count: 4697

このソリューションはより滑らかな曲線のように見えます。60 個ではなく 160 個のパレートポイントを使用する場合、ソルバーは 3 倍以上の関数評価を実行します。

`gamultiobj` ソリューション

ソルバーが違いを生むかどうかを確認するには、問題に gamultiobj ソルバーを試してください。前のソリューションと同じオプションを設定します。gamultiobj ソルバーは、そのソリューションの半分未満を最良のパレートフロント上に保持するため、以前の 2 倍のポイントを使用します。

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

Figure Genetic Algorithm contains an axes object. The axes object with title Pareto Front, xlabel Objective 1, ylabel Objective 2 contains an object of type scatter.

disp("Total Function Count: " + gaoutput1.funccount);

Total Function Count: 33281

gamultiobj は数万回の関数評価を必要としますが、paretosearch は数千回しか必要としません。gamultiobj は、特に目的 2 の場合、目的関数値の範囲がわずかに広くなります。

解の比較

gamultiobj ソリューションは paretosearch ソリューションとは異なるように見えますが、プロットされたスケールが異なるため、区別するのは困難です。同じスケールを使用して、2 つの解を同じプロットにプロットします。

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

Figure contains an axes object. The axes object with xlabel Cost, ylabel Deflection contains 2 objects of type line. These objects represent paretosearch, gamultiobj.

paretosearch の解はプロットの左端の部分ではより優れており、残りの部分では解は区別がつかないように見えます。paretosearch は、解を得るために使用する関数評価の回数がはるかに少なくなります。

通常、問題に非線形制約がない場合、paretosearch は少なくとも gamultiobj と同程度に正確です。ただし、結果として得られるパレート集合の範囲は多少異なる場合があります。

paretosearch の主な利点の 1 つは、通常、関数評価の回数が大幅に少なくなることです。

単一目的のソリューションから始める

ソルバーがより良いソリューションを見つけられるように、個々の目的関数を最小化するソリューションであるポイントからソルバーを始めます。pickindex 関数は、objval 関数から単一の目的を返します。fmincon を使用して単一目的の最適値を見つけます。次に、それらのソリューションを多目的検索の初期ポイントとして使用します。

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

単一目的最適値を調べます。

objval(x0(1,:))

ans = 1×2

    2.3810    0.0158

objval(x0(2,:))

ans = 1×2

   76.7188    0.0004

最小コストは 2.381 で、偏向は 0.158 になります。最小たわみは 0.0004 で、コストは 76.7188 です。プロットされた曲線は範囲の端近くでかなり急勾配になっています。つまり、コストを最小値より少しだけ高くすると偏向が大幅に少なくなり、偏向を最小値より少しだけ高くするとコストが大幅に少なくなります。

単一目的のソリューションから paretosearch を開始します。後で同じプロットにソリューションをプロットするため、paretosearch プロット関数を削除します。

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: 5007

同じ初期点から ga を開始し、そのプロット関数を削除します。

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: 33281

解を同じ軸上にプロットします。

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

Figure contains an axes object. The axes object with xlabel Cost, ylabel Deflection contains 2 objects of type line. These objects represent paretosearch, gamultiobj.

単一目的のソリューションから開始すると、gamultiobj ソリューションはプロットされた範囲全体で paretosearch ソリューションに類似します。ただし、gamultiobj では、解に到達するまでにほぼ 10 倍の関数評価が必要になります。

ハイブリッド機能

gamultiobj は、ハイブリッド関数 fgoalattain を自動的に呼び出して、より正確なソリューションに到達しようとします。ハイブリッド関数によってソリューションが改善されるかどうかを確認します。

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: 44923

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

Figure contains an axes object. The axes object with xlabel Cost, ylabel Deflection contains 3 objects of type line. These objects represent paretosearch, gamultiobj, gamultiobj/fgoalattain.

ハイブリッド関数は、主にプロットの左端の部分で、gamultiobj ソリューションにわずかな改善をもたらします。

`paretosearch` ソリューションポイントから `fgoalattain` を手動で実行する

paretosearch にはハイブリッド関数が組み込まれていませんが、paretosearch ソリューションポイントから fgoalattain を実行することで paretosearch ソリューションを改善できます。gamultiobj ハイブリッド機能で説明したのと同じ fgoalattain の設定を使用して、fgoalattain の目標と重みを作成します。

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: 10438

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'

Figure contains an axes object. The axes object with xlabel Cost, ylabel Deflection contains 4 objects of type line. These objects represent paretosearch, gamultiobj, gamultiobj/fgoalattain, paretosearch/fgoalattain.

paretosearch と fgoalattain の組み合わせにより、最も正確なパレート面が作成されます。拡大してご覧ください。

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

追加の fgoalattain 計算があっても、組み合わせの合計関数数は、gamultiobj ソリューションのみの関数数の半分未満になります。

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 33281.
For paretosearch and fgoalattain together it is 15445.

プロットから適切なパラメータを見つける

プロットされた点は関数空間における最良の値を示します。どのパラメータがこれらの関数値を達成するかを判断するには、特定のコスト/たわみポイントを取得するために梁のサイズと溶接のサイズを見つけます。たとえば、paretosearch の後に fgoalattain が続くプロットは、コストが約 6 で偏向が約 3.5e–3 のポイントを示しています。これらのポイントを達成する梁と溶接のサイズを決定します。

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

goodpoints=7×2 table
                   Parameters                       Objectives     
    ________________________________________    ___________________

    0.63052       1.53         10     0.6634    5.6286     0.003309
    0.63907     1.5063         10     0.6829    5.7741    0.0032145
    0.64311     1.4954         10    0.69221    5.8435    0.0031713
    0.61497     1.5749         10     0.6286    5.3681    0.0034922
    0.62035     1.5591         10    0.64056    5.4577     0.003427
    0.63655     1.5132         10    0.67714    5.7311    0.0032419
    0.64834     1.4814         10    0.70433    5.9338    0.0031167

いくつかのパラメータセットでは、コストが 6 未満、偏向が 3.5e–3 未満になります。

溶接厚さは0.6をわずかに超える
溶接長さ約1.5
梁の高さ 10 (上限)
ビーム幅0.63～0.71

客観的制約と非線形制約

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

参照

[1] Deb、Kalyanmoy、J. Sundar、Udaya Bhaskara Rao N、Shamik Chaudhuri。進化的アルゴリズムを使用した参照ポイントベースの多目的最適化。国際計算知能研究ジャーナル、第2巻第3号、2006年、273〜286頁。https://www.softcomputing.net/ijcir/vol2-issu3-paper4.pdf から入手可能

[2] Ray, T.、およびK. M. Liew.多目的設計最適化のための群集メタファー。エンジニアリング最適化34、2002年、pp.141-153。