Main Content

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

遺伝的アルゴリズムを使用して混合整数エンジニアリング設計問題を解く

スクリプトを開く

この例では、Global Optimization Toolbox の遺伝的アルゴリズム (ga) ソルバーを使用して、混合整数エンジニアリング設計問題を解決する方法を示します。

この例で示す問題は、段付き片持ち梁の設計に関するものです。特に、梁は規定の端部荷重を支えることができなければなりません。さまざまなエンジニアリング設計の制約の下で梁の体積を最小限に抑える問題を解決します。

この例では、[1]で公開された問題の2つの境界バージョンを解きます。

段付き片持ち梁の設計問題

下の図に示すように、段付き片持ち梁は一端で支えられ、自由端に荷重がかかります。梁は、サポートから一定の距離 $L$ で、指定された荷重 $P$ を支えることができなければなりません。梁の設計者は、各セクションの幅 ($b_i$) と高さ ($h_i$) を変更できます。カンチレバーの各セクションの長さは同じで、$l$ であると仮定します。

ビームの体積

梁の体積$V$は、個々のセクションの体積の合計です。

$$V = l(b_1h_1 + b_2h_2 + b_3h_3 + b_4h_4 + b_5h_5)$$

デザイン上の制約:1 - 曲げ応力

座標の中心が梁の自由端の断面の中心にある単一の片持ち梁を考えます。梁の点$(x, y, z)$における曲げ応力は、次の式で与えられる。

$$\sigma_b = M(x)y/I$$

ここで、$M(x)$ は $x$ における曲げモーメント、$x$ は端部荷重からの距離、$I$ は梁の面積慣性モーメントです。

さて、図に示す段付き片持ち梁では、梁の各セクションの最大モーメントは $PD_i$ です。ここで、$D_i$ は梁の各セクションの端部荷重 $P$ からの最大距離です。したがって、梁の$i$番目のセクション$\sigma_i$の最大応力は次のように表される。

$$\sigma_i = PD_i(h_i/2)/I_i$$

最大応力は梁の端で発生します（$y = h_i/2$）。梁の$i$番目の断面の断面慣性モーメントは次のように表される。

$$I_i = b_ih_i^3/12$$

これを$\sigma_i$の式に代入すると、

$$\sigma_i = 6PD_i/b_ih_i^2$$

カンチレバーの各部分の曲げ応力は、最大許容応力 $\sigma_{max}$ を超えてはなりません。その結果、最終的に5つの曲げ応力制約（カンチレバーの各ステップに1つずつ）を述べることができる。

$$\frac{6Pl}{b_5h_5^2} \leq \sigma_{max}$$

$$\frac{6P(2l)}{b_4h_4^2} \leq \sigma_{max}$$

$$\frac{6P(3l)}{b_3h_3^2} \leq \sigma_{max}$$

$$\frac{6P(4l)}{b_2h_2^2} \leq \sigma_{max}$$

$$\frac{6P(5l)}{b_1h_1^2} \leq \sigma_{max}$$

デザイン上の制約:2 - 端のたわみ

カンチレバーの端のたわみは、カスティリアーノの第2定理を使用して計算できます。

$$\delta = \frac{\partial U}{\partial P}$$

ここで、$\delta$ は梁のたわみ、$U$ は加えられた力によって梁に蓄えられたエネルギー、$P$ です。

片持ち梁に蓄えられたエネルギーは次のように表される。

$$U = \int_0^L \! M^2/2EI \, \mathrm{d} x$$

ここで、$M$ は $x$ における適用された力のモーメントです。

片持ち梁の$M = Px$を考えると、上記の式は次のように書けます。

$$U =
P^2/2E \int_0^l \! [(x+4l)^2/I_1 \,
 + (x+3l)^2/I_2 \,
 + (x+2l)^2/I_3 \,
 + (x+l)^2/I_4 \,
 + x^2/I_5 ]\, \mathrm{d} x$$

ここで、$I_n$ はカンチレバーの $n$ 番目の部分の面積慣性モーメントです。積分を評価すると、$U$ の次の式が得られます。

$$U = (P^2/2)(l^3/3E)(61/I_1 + 37/I_2 + 19/I_3 + 7/I_4 + 1/I_5)$$

カスティリアーノの定理を適用すると、梁の端部たわみは次のように表される。

$$\delta = Pl^3/3E(61/I_1 + 37/I_2 + 19/I_3 + 7/I_4 + 1/I_5)$$

ここで、カンチレバーの端部たわみ $\delta$ は、最大許容たわみ $\delta_{max}$ よりも小さくなる必要があり、次の制約が得られます。

$$Pl^3/3E(61/I_1 + 37/I_2 + 19/I_3 + 7/I_4 + 1/I_5) \leq
\delta_{max}$$

デザイン上の制約:3 - アスペクト比

カンチレバーの各ステップのアスペクト比は、最大許容アスペクト比 $a_{max}$ を超えてはなりません。つまり、

$h_i/b_i \leq a_{max}$ ($i = 1, ..., 5$ 用)

最適化問題を述べる

これで、規定された制約を前提として、段付き片持ち梁の最適なパラメータを見つけるという問題を定義できるようになりました。

$x_1 = b_1$、$x_2 = h_1$、$x_3 = b_2$、$x_4 = h_2$、$x_5 = b_3$、$x_6 = h_3$、$x_7 = b_4$、$x_8 = h_4$、$x_9 = b_5$、$x_{10} = h_5$ とします。

最小化:

$$V = l(x_1x_2 + x_3x_4 + x_5x_6 + x_7x_8 + x_9x_{10})$$

以下を条件とします:

$$\frac{6Pl}{x_9x_{10}^2} \leq \sigma_{max}$$

$$\frac{6P(2l)}{x_7x_8^2} \leq \sigma_{max}$$

$$\frac{6P(3l)}{x_5x_6^2} \leq \sigma_{max}$$

$$\frac{6P(4l)}{x_3x_4^2} \leq \sigma_{max}$$

$$\frac{6P(5l)}{x_1x_2^2} \leq \sigma_{max}$$

$$\frac{Pl^3}{E}(\frac{244}{x_1x_2^3} + \frac{148}{x_3x_4^3} +
\frac{76}{x_5x_6^3} + \frac{28}{x_7x_8^3} +
\frac{4}{x_9x_{10}^3}) \leq \delta_{max}$$

$x_2/x_1 \leq 20$、$x_4/x_3 \leq 20$、$x_6/x_5 \leq 20$、$x_8/x_7 \leq 20$、$x_{10}/x_9 \leq 20$

ビームの最初のステップは、最も近いセンチメートル単位でのみ機械加工できます。つまり、$x_1$ と $x_2$ は整数である必要があります。残りの変数は連続的です。変数の境界は以下の通りです。

$$1 \leq x_1 \leq 5$$

$$30 \leq x_2 \leq 65$$

$$2.4 \leq x_3, x_5 \leq 3.1$$

$$45 \leq x_4, x_6 \leq 60$$

$$1 \leq x_7, x_9 \leq 5$$

$30 \leq x_8, x_{10} \leq 65$

この問題の設計パラメータ

この例で解決する問題では、梁がサポートしなければならない端部荷重は $P = 50000 N$ です。

梁の長さと最大端部たわみは次のとおりです。

  • ビーム全長、$L = 500 cm$

  • 梁の個々のセクション、$l = 100 cm$

  • 最大ビーム端たわみ、$\delta_{max} = 2.7 cm$

梁の各ステップにおける最大許容応力、$\sigma_{max} = 14000 N/cm^2$

梁の各段のヤング率、$E = 2\times10^{7} N/cm^2$

混合整数最適化問題を解く

ここで、最適化問題を述べるで説明した問題を解決します。

適合度関数と制約関数を定義する

MATLAB® ファイル cantileverVolume.m および cantileverConstraints.m を調べて、適合度関数と制約関数がどのように実装されているかを確認します。

線形制約に関する注意:ga に線形制約を指定する場合、通常は A、b、Aeq、および beq 入力を介して指定します。この場合、非線形制約関数を使用して指定しました。これは、この例の後半で、変数の一部が離散的になるためです。問題に離散変数がある場合、非線形制約関数で線形制約を指定する方がはるかに簡単になります。代替案としては、線形制約行列を変換された変数空間で動作するように変更することですが、これは簡単ではなく、おそらく不可能です。また、混合整数 ga ソルバーでは、線形制約は、指定方法に関係なく、非線形制約とは異なる扱いを受けません。

境界を設定する

下限制約 (lb) と上限制約 (ub) を含むベクトルを作成します。

lb = [1 30 2.4 45 2.4 45 1 30 1 30];
ub = [5 65 3.1 60 3.1 60 5 65 5 65];

オプションを設定する

より正確なソリューションを得るために、PopulationSize および MaxGenerations オプションをデフォルト値から増やし、EliteCount および FunctionTolerance オプションを減らします。これらの設定により、ga はより大きな集団を使用し (PopulationSize が増加)、設計空間の検索を増やし (EliteCount が減少)、最適なメンバーがほとんど変化しなくなるまで (FunctionTolerance が小さくなる) 処理を続けます。また、ga の進行に応じてペナルティ関数の値を監視するためのプロット関数も指定します。

混合整数問題を解くときに使用できる ga オプションのセットは制限されていることに注意してください。詳細については、Global Optimization Toolbox ユーザーズ ガイドを参照してください。

opts = optimoptions(@ga, ...
                    'PopulationSize', 150, ...
                    'MaxGenerations', 200, ...
                    'EliteCount', 10, ...
                    'FunctionTolerance', 1e-8, ...
                    'PlotFcn', @gaplotbestf);

問題を解決するには ga に電話してください

ga を呼び出して、問題を解きます。問題文では、$x_1$ と $x_2$ は整数変数です。これを指定するには、非線形制約入力の後、オプション入力の前にインデックス ベクトル [1 2] から ga を渡します。再現性を確保するために、ここで乱数ジェネレーターのシードと設定も行います。

rng(0, 'twister');
[xbest, fbest, exitflag] = ga(@cantileverVolume, 10, [], [], [], [], ...
    lb, ub, @cantileverConstraints, [1 2], opts);

ga stopped because it exceeded options.MaxGenerations.

結果を分析する

問題に整数制約がある場合、ga はそれを内部的に再定式化します。特に、問題内の適応度関数は、制約を処理するペナルティ関数に置き換えられます。実行可能な集団メンバーの場合、ペナルティ関数は適応度関数と同じです。

ga から返されたソリューションを以下に示します。サポートに最も近いセクションは、幅 ($x_1$) と高さ ($x_2$) が整数値に制限されており、この制約は GA によって尊重されていることに注意してください。

display(xbest);

xbest =

  Columns 1 through 7

    3.0000   60.0000    2.8504   57.0057    2.6114   50.6243    2.2132

  Columns 8 through 10

   44.2349    1.7543   35.0595

ga にビームの最適なボリュームを返すように要求することもできます。

fprintf('\nCost function returned by ga = %g\n', fbest);

Cost function returned by ga = 63408.9

離散非整数変数制約を追加する

エンジニアは、カンチレバーの 2 段目と 3 段目の幅と高さは標準セットから選択できることを知らされます。このセクションでは、この制約を最適化問題に追加する方法を示します。この制約を追加すると、この問題は[1]で解決された問題と同一になることに注意してください。

まず、上記の最適化に追加される追加の制約を述べる。

  • 梁の2段目と3段目の幅は、次のセットから選択する必要があります:- [2.4、2.6、2.8、3.1] cm

  • 梁の2段目と3段目の高さは、次のセットから選択する必要があります:- [45、50、55、60] cm

この問題を解決するには、変数 $x_3$、$x_4$、$x_5$、$x_6$ を離散変数として指定できる必要があります。コンポーネント $x_j$ がセット $S = {v_1,\ldots,v_k}$ から離散値を取るように指定するには、1 から $k$ までの値を取る整数変数 $x_j$ を使用して最適化し、離散値として $S(x_j)$ を使用します。範囲 (1 から $k$) を指定するには、下限として 1、上限として $k$ を設定します。

そこで、まず離散変数の境界を変換します。各セットには4つのメンバーがあり、離散変数を範囲[1, 4]の整数にマッピングします。したがって、これらの変数を整数にマッピングするには、各変数の下限を 1 に、上限を 4 に設定します。

lb = [1 30 1 1 1 1 1 30 1 30];
ub = [5 65 4 4 4 4 5 65 5 65];

ga ソルバーが呼び出されると、$x_3$、$x_4$、$x_5$、$x_6$ の変換された (整数) バージョンがフィットネス関数と制約関数に渡されるようになりました。これらの関数を正しく評価するには、$x_3$、$x_4$、$x_5$、$x_6$ をこれらの関数内の指定された離散セットのメンバーに変換する必要があります。これがどのように行われるかを確認するには、MATLAB ファイル cantileverVolumeWithDisc.m、cantileverConstraintsWithDisc.m、および cantileverMapVariables.m を調べてください。

これで、ga を呼び出して離散変数の問題を解決できます。この場合、$x_1, ..., x_6$ は整数です。つまり、インデックス ベクトル 1:6 を ga に渡して整数変数を定義します。

rng(0, 'twister');
[xbestDisc, fbestDisc, exitflagDisc] = ga(@cantileverVolumeWithDisc, ...
    10, [], [], [], [], lb, ub, @cantileverConstraintsWithDisc, 1:6, opts);

ga stopped because it exceeded options.MaxGenerations.

結果を分析する

xbestDisc(3:6) は ga から整数として (つまり、変換された状態で) 返されます。エンジニアリング単位で値を取得するには、変換を逆にする必要があります。

xbestDisc = cantileverMapVariables(xbestDisc);
display(xbestDisc);

xbestDisc =

  Columns 1 through 7

    3.0000   60.0000    3.1000   55.0000    2.6000   50.0000    2.2430

  Columns 8 through 10

   44.8603    1.8279   36.5593

前と同様に、ga から返されるソリューションは、$x_1$ と $x_2$ が整数であるという制約を尊重します。また、$x_3$、$x_5$ は [2.4, 2.6, 2.8, 3.1] cm のセットから選択され、$x_4$、$x_6$ は [45, 50, 55, 60] cm のセットから選択されていることもわかります。

変数 x(3)、x(4)、x(5)、x(6) に追加の制約を追加したことを思い出してください。予想どおり、これらの変数に追加の離散制約がある場合、最適ソリューションの最小ボリュームは高くなります。さらに、[1]で報告された解は$64558 cm^3$の最小体積を持ち、[1]で報告されたものとほぼ同じ解が見つかることにも注意してください。

fprintf('\nCost function returned by ga = %g\n', fbestDisc);

Cost function returned by ga = 64795

まとめ

この例では、遺伝的アルゴリズム ソルバー ga を使用して、整数制約を持つ制約付き非線形最適化問題を解決する方法を示します。この例では、問題の定式化に離散変数が含まれる問題を処理する方法も示します。

参照

[1] Thanedar, P. B.、およびG. N. Vanderplaats。「構造設計のための離散変数最適化の調査」構造工学ジャーナル121（3）、1995年、pp.301-306。

関連するトピック