使用遗传算法解决混合整数工程设计问题

此示例展示如何使用 Global Optimization Toolbox 中的遗传算法(ga)求解器解决混合整数工程设计问题。

本例所说明的问题涉及阶梯悬臂梁的设计。具体来说，梁必须能够承载规定的端部载荷。我们将解决一个问题，即在各种工程设计约束下最小化梁体积。

在这个例子中，我们将解决 [1] 中发表的两个有界版本的问题。

阶梯式悬臂梁设计问题

阶梯悬臂梁的一端受到支撑，并在自由端施加载荷，如下图所示。梁必须能够支撑距支撑固定距离 $P$ 的给定负载 $L$ 。梁的设计者可以改变每个部分的宽度（ $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}$ 。因此，我们最终可以状态五个弯曲应力约束（悬臂的每个台阶一个）

$\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 - 末端挠度

悬臂的末端挠度可以利用卡斯蒂利亚诺第二定理计算，该定理指出

$\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$ 。

梁长度和最大端部挠度为：

梁总长度，
梁的单个部分，
最大梁端挠度， $\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

添加离散非整数变量约束

现在工程师们被告知，悬臂的第二步和第三步只能具有从标准集合中选择的宽度和高度。在本节中，我们展示如何将此约束添加到优化问题中。请注意，增加这个约束后，这个问题与[1]中解决的问题相同。

首先，我们状态将添加到上述优化中的额外约束

梁的第二和第三步的宽度必须从以下集合中选择：- [2.4, 2.6, 2.8, 3.1] 厘米
横梁第二和第三台阶的高度必须从以下集合中选择：- [45, 50, 55, 60] 厘米

为了解决这个问题，我们需要能够将变量 $x_3$ 、 $x_4$ 、 $x_5$ 和 $x_6$ 指定为离散变量。要指定分量 $x_j$ 从集合 $S = {v_1,\ldots,v_k}$ 中取离散值，请使用 $x_j$ 优化取值范围从 1 到 $k$ 的整数变量，并使用 $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., and G. N. Vanderplaats."Survey of Discrete Variable Optimization for Structural Design."Journal of Structural Engineering 121 (3), 1995, pp. 301–306.

另请参阅

主题

混合整数 ga 优化