基于最大熵原理的复合形法及其在边坡稳定分析中的应用

摘要

将最大熵原理应用于基本复合形法的优化计算中，旨在提高算法的寻优成功率。具体步骤为：选取复形中熵最大的顶点为寻优点；在复形其余各点的加权中心点与寻优点的连线上，搜索比寻优点改善的新顶点；并替换掉寻优点构成新的复形；如此迭代直至复形中的信息被完全利用。采用基本复合形法和基于最大熵原理的复合形法对两个复杂土坡的最小安全系数进行了搜索，结果表明基于最大熵原理的复合形法的全局搜索能力有了明显改善。

正文

《1 引言》

1 引言

用一种简单、有效的搜索方法寻找土坡最危险滑动面及其对应的安全系数, 是极限平衡法分析土坡稳定的关键。对均质土坡, 经验方法^[1]及传统的搜索算法均可确定最危险的滑动面;但对工程实践中广泛存在的复杂土坡、堆石坝等构筑物而言, 无论是圆弧滑动面还是非圆滑动面, 其安全系数对应的目标函数往往存在多个局部极小值, 利用传统搜索方法如枚举法、模式搜索法、共轭梯度法等搜索到的最危险滑动面往往依赖于初始值或陷入局部极优。陈祖煜等采用单纯形法和DFP优化方法搜索边坡的最小安全系数^[2], 张永生等采用复合形法计算土坡的稳定性^[3]。近来有学者利用遗传算法^[4]、模拟退火算法^[5]、蚁群算法等随机搜索算法来研究边坡的稳定性^[6], 均取得了较好的效果。但遗传算法存在早熟收敛、计算量大等缺陷;模拟退火算法在实际应用过程中不能保证退火充分而陷入局部最优。复合形法在迭代计算中不必计算目标函数的一阶或二阶导数, 是工程设计中较为常用的一种有约束的直接求解方法。但对复杂优化问题, 复合形法在迭代过程中会出现迭代失败现象, 对任意生成的初始复形, 基本复合形法搜索到全局最优解的概率较低。本文利用最大熵原理改进了基本复合形法的寻优思路, 提高了算法的寻优成功率, 因此具有很强的工程应用价值。

《2 边坡最小安全系数搜索的优化模型》

2 边坡最小安全系数搜索的优化模型

给定土坡剖面及土性力学参数后, 任一圆弧滑动面就可由设计变量 $\vec{X} = (X, Y, Ζ)$ $\vec{X} = (X, Y, Ζ)$ 唯一确定, 其对应的抗滑稳定安全系数为S ( $\vec{X}$ $\vec{X}$ ) , 寻找最小安全系数的优化模型表述如下:

$\begin{array}{l} Μ i n i m i z e : S (\vec{X}) \\ S . Τ g (\vec{X}) \leq 0; \vec{X}_{L} \leq \vec{X} \leq \vec{X}_{U} \end{array}$ $\begin{array}{l} Μ i n i m i z e : S (\vec{X}) \\ S . Τ g (\vec{X}) \leq 0; \vec{X}_{L} \leq \vec{X} \leq \vec{X}_{U} \end{array}$

上式中, X为圆弧滑动面圆心的x坐标, Y为圆弧滑动面圆心的y坐标, Z为滑动面在坡脚滑出点的x坐标。g ( $\vec{X}$ $\vec{X}$ ) 代表设计变量需满足的隐式约束条件, 即由设计变量值确定的潜在滑动面要实际可行, $\vec{X}$ $\vec{X}$ _L, $\vec{X}$ $\vec{X}$ _U表示设计变量的取值范围即区间约束, 由 $\vec{X}$ $\vec{X}$ _L, $\vec{X}$ $\vec{X}$ _U包围的空间为搜索域。

为简单起见, 采用瑞典条分法计算每一潜在滑动面的抗滑稳定安全系数^[7], 计算示意图见图1。

$S (\vec{X}) = \frac{\sum_{i = 1}^{Ν} (c_{i} l_{i} + γ_{i} b_{i} h_{i} \cos α_{i} \tan ϕ_{i}}{\sum_{i = 1}^{Ν} γ_{i} b_{i} h_{i} \sin α_{i}} (1)$ $S (\vec{X}) = \frac{\sum_{i = 1}^{Ν} (c_{i} l_{i} + γ_{i} b_{i} h_{i} \cos α_{i} \tan ϕ_{i}}{\sum_{i = 1}^{Ν} γ_{i} b_{i} h_{i} \sin α_{i}} (1)$

《图1》

图1 滑动面安全系数计算模型示意

Fig.1 Illustration of computation model for factor of safety of slip surface

式 (1) 及图1中N为条分数目;c_i为滑动面上土层的粘聚力;l_i为第i土条滑动面的弧长;γ_i为土的容重;b_i, h_i为土条的宽度和高度;α_i为第i土条滑动面法线与竖直线的夹角;ϕ_i为滑动面上土层的内摩擦角。

《3 基本复合形法》

3 基本复合形法

复合形法是在单纯形法基础上发展起来的。该方法是在n维受约束的设计空间内由k (n+1≤k≤2n) 个顶点构成多面体 (复形) , 然后对复形顶点的函数值逐一进行比较, 不断丢掉函数值最劣的顶点, 代之以满足约束条件且函数值有所改善的新顶点, 如此重复, 逐步逼近最优点为止, 见图2。

若目标函数最坏点b为 $\vec{X}_{b} = (X_{b 1}, X_{b 2}, \dots, X_{b n})$ $\vec{X}_{b} = (X_{b 1}, X_{b 2}, \dots, X_{b n})$ , 除最坏点b外其余顶点的中心点o为 $\vec{X}^{o} = (X_{1}^{o}, X_{2}^{o}, \dots, X_{n}^{o})$ $\vec{X}^{o} = (X_{1}^{o}, X_{2}^{o}, \dots, X_{n}^{o})$ , 改善点r由下式寻找

$\vec{X}^{r} = \vec{X}^{o} + α (\vec{X}^{o} - \vec{X}_{b}) (2)$ $\vec{X}^{r} = \vec{X}^{o} + α (\vec{X}^{o} - \vec{X}_{b}) (2)$

式中α为映射系数, 开始寻找时通常给定初值α_ini=constan t (>1) , 然后不断减半 (α=α/2) 收缩至搜索到的r点比b点的目标函数值有所改善为止。若当α≤ξ时仍然没有找到改善点r, 则称“关于最坏点b映射失败现象”, ξ为一很小的正数。此时就选取目标函数值次大的顶点为最坏点b, 重新进行搜索。若o点不可行, 就不能利用式 (2) 来寻找改善点r, 一般处理方法为在以最好点 (目标函数值最小的顶点) g和中心点o为界的超立方体中随机产生新的初始复形重新计算。

《图2》

图2 基本复合形法流程示意

Fig.2 The illustrator of flowchart of basic complex method

因此, 基本复合形法实质上是一种群体搜索算法, 多个顶点构成的初始复形类似于遗传算法中的初始群体, 最坏点b相当于遗传算法中欲变异的个体, 其余各点综合为几何中心点o与遗传算法中交叉算子生成两个子体雷同, 只不过前者利用了多个顶点的信息, 而后者仅仅利用了两个父体的信息而已。遗传算法中变异算子对未知的解空间进行开发, 交叉算子对已知的解空间进行更细致的探索。在基本复合形法中, 首先规定了一个最坏点b为欲改善的点, 然后将剩余顶点综合为一点o, 在bo连线上搜索新点以代替b点, 可见基本复合形法通过直线搜索将对解空间的探索以及对解空间的开发有机结合起来, 初步具有了群体进化的特征。笔者以为, 既然是对解空间进行探索、开发, 就应该最大限度地利用当前复形的信息, 为此本文提出了基于最大熵原理的复合形法 (CMBOMEP) 。

《4 基于最大熵原理的复合形法》

4 基于最大熵原理的复合形法

《4.1 复形顶点熵的度量》

4.1 复形顶点熵的度量

由信息理论可知^[8], 信息量蕴含于不确定性之中。而不确定性在概率论中是用随机变量或随机事件来描述的。设随机变量为y, 其可能的取值有m个, 即y₁, y₂, …, y_m, p_i为y取为y_i的概率, 则随机变量y的信息熵为

$Η (y) = Η (p_{1}, p_{2}, \dots, p_{m}) = - d \sum_{i = 1}^{m} p_{i} l o g a^{p_{i}} (3)$ $Η (y) = Η (p_{1}, p_{2}, \dots, p_{m}) = - d \sum_{i = 1}^{m} p_{i} l o g a^{p_{i}} (3)$

它度量了随机变量y的不确定程度, 信息熵可理解为物质系统内部状态的丰富度或复杂程度。式中d为一常数, 具体值取决于熵的单位, 选取不同的对数底, 熵就有不同的单位, 本文取e为底, 此时d=1, 熵的单位为奈特。

复形中每个顶点包含的信息量是不同的。给定设计变量 $\vec{X} = (x_{1}, x_{2}, \dots, x_{n})$ $\vec{X} = (x_{1}, x_{2}, \dots, x_{n})$ , 复形中k个顶点为 $\vec{X}$ $\vec{X}$ ₁, $\vec{X}$ $\vec{X}$ ₂, …, $\vec{X}$ $\vec{X}$ _k, 其中 $\vec{X}_{i} = (X_{i 1}, X_{i 2}, \dots, X_{i n}), i = 1, 2, \dots, k$ $\vec{X}_{i} = (X_{i 1}, X_{i 2}, \dots, X_{i n}), i = 1, 2, \dots, k$ 。将x_i, i=1, 2, …, n, 视为随机变量, 其可能的取值有X_1i, X_2i, …, X_kik个, p_ji代表x_i取X_ji的概率, 则第i个顶点 $\vec{X}$ $\vec{X}$ _i的信息熵为

$Η (\vec{X}_{i}) = - \sum_{m = 1}^{n} p_{i m} l n p_{i m} (4)$ $Η (\vec{X}_{i}) = - \sum_{m = 1}^{n} p_{i m} l n p_{i m} (4)$

式中 $p_{j i} = \frac{b o t t o m_{i}}{d i s_{j i}}$ $p_{j i} = \frac{b o t t o m_{i}}{d i s_{j i}}$ ; $b o t t o m_{i} = \min (d i s_{j i}) ‚ d i s_{j i} = a b s (X_{j i} - e q_{i}), e q_{i} = \frac{1}{k} \sum_{j = 1}^{k} X_{j i}, j = 1, 2, \dots, k$ $b o t t o m_{i} = \min (d i s_{j i}) ‚ d i s_{j i} = a b s (X_{j i} - e q_{i}), e q_{i} = \frac{1}{k} \sum_{j = 1}^{k} X_{j i}, j = 1, 2, \dots, k$ 。

基本复合形法利用下式计算几何中心点o,

$\vec{X}^{o} = \frac{1}{k - 1} \sum_{\begin{array}{l} i = 1 \\ i \neq b \end{array}}^{k} \vec{X}_{i} (5)$ $\vec{X}^{o} = \frac{1}{k - 1} \sum_{\begin{array}{l} i = 1 \\ i \neq b \end{array}}^{k} \vec{X}_{i} (5)$

可见, 几何中心点o平均利用每个顶点的信息, 这不符合最大限度利用信息的原则。对优化目标贡献大的顶点要多加利用, 相反贡献小的顶点对o点的作用就小, 根据各个顶点对优化目标的贡献大小赋予每个顶点不同的权重, 由下式来计算中心点o:

$\vec{X}^{o} = \sum_{\begin{array}{l} i = 1 \\ i \neq h \end{array}}^{k} w_{i} \vec{X}_{i} (6)$ $\vec{X}^{o} = \sum_{\begin{array}{l} i = 1 \\ i \neq h \end{array}}^{k} w_{i} \vec{X}_{i} (6)$

式中:w_i为第i个顶点的权重, 由以下公式给出

$w_{i} = \frac{\frac{1.0}{S (\vec{X}_{i})} * Η (\vec{X}_{i})}{\sum_{\frac{j = 1}{j \neq h}}^{k} \frac{1.0}{S (\vec{X}_{j})} * Η (\vec{X}_{j})} ‚ i = 1, 2, \dots, k; i \neq h (7)$ $w_{i} = \frac{\frac{1.0}{S (\vec{X}_{i})} * Η (\vec{X}_{i})}{\sum_{\frac{j = 1}{j \neq h}}^{k} \frac{1.0}{S (\vec{X}_{j})} * Η (\vec{X}_{j})} ‚ i = 1, 2, \dots, k; i \neq h (7)$

式中S ( $\vec{X}$ $\vec{X}$ _i) 为第i个顶点的目标函数值, h为当前复形中熵最大的顶点, 利用式 (7) 计算得到的中心点o是加权中心点, 比几何中心点更加有效地利用了当前复形的信息。

《4.2 动态全域映射收缩算子》

4.2 动态全域映射收缩算子

基本复合形法在bo连线上寻找改善点时, 映射系数始终为正, 这就将寻找范围限制在图3中的实线部分, 而忽略了位于虚线上的可行点, 导致搜索不完全, 可称之为局域映射收缩算子, 为此本文提出了动态全域映射收缩算子。

《图3》

图3 全域映射收缩算子示意

Fig.3 The illustrator of full reflection and retraction operator

用下述程序计算映射系数的最大、最小值α_max, α_min。

for i=1 to n (n为设计变量的个数)

$\begin{array}{l} α_{l_{i}} = (l_{i} - X_{i}^{o}) / (X_{i}^{o} - X_{b i}) \\ α_{u_{i}} = (u_{i} - X_{i}^{o}) / (X_{i}^{o} - X_{b i}) \\ α_{i, 1} = \min (α_{l_{i}}, α_{u_{i}}) \\ α_{i, 2} = \max (α_{l_{i}}, α_{u_{i}}) \\ n e x t i \\ α_{\min} = \max {α_{i, 1}}_{i = 1, 2 \dots, n} \\ α_{\max} = \min {α_{i, 2}}_{i = 1, 2, \dots, n} \end{array}$ $\begin{array}{l} α_{l_{i}} = (l_{i} - X_{i}^{o}) / (X_{i}^{o} - X_{b i}) \\ α_{u_{i}} = (u_{i} - X_{i}^{o}) / (X_{i}^{o} - X_{b i}) \\ α_{i, 1} = \min (α_{l_{i}}, α_{u_{i}}) \\ α_{i, 2} = \max (α_{l_{i}}, α_{u_{i}}) \\ n e x t i \\ α_{\min} = \max {α_{i, 1}}_{i = 1, 2 \dots, n} \\ α_{\max} = \min {α_{i, 2}}_{i = 1, 2, \dots, n} \end{array}$

程序中l_i, u_i为第i设计变量的下、上限。搜索改善点r时, 首先赋α_ini=α_max, 在bo连线的实线部分寻找改善点r, 若当α→0时没有找到, 则赋α_ini=α_min在bo连线的虚线部分再次寻找改善点r, 这种映射收缩算子称为动态全域映射收缩算子。若仍然采取基本复合形法的最坏点替代策略, 则算法称之为动态全域映射收缩复合形法。

《4.3 基于最大熵原理的复合形法步骤》

4.3 基于最大熵原理的复合形法步骤

1) 确定设计变量个数n=3以及设计变量的上下限 $\vec{X}_{U} = (u_{1}, u_{2}, \dots, u_{n}) ‚ \vec{X}_{L} = (l_{1}, l_{2}, \dots, l_{n})$ $\vec{X}_{U} = (u_{1}, u_{2}, \dots, u_{n}) ‚ \vec{X}_{L} = (l_{1}, l_{2}, \dots, l_{n})$ , 在搜索域内随机产生k=6个顶点 $\vec{X}$ $\vec{X}$ ₁, $\vec{X}$ $\vec{X}$ ₂, …, $\vec{X}$ $\vec{X}$ ₆构成初始复形, 并计算各顶点的目标函数值S ( $\vec{X}$ $\vec{X}$ _i) , order=1。

2) 计算复形中H值第order大顶点 (该顶点的H值在所有顶点中按降序排第order位, 下同) $\vec{X}$ $\vec{X}$ _H。

3) 由式 (6) 计算除 $\vec{X}$ $\vec{X}$ _H外其余顶点的加权中心点 $\vec{X}$ $\vec{X}$ ^o, 若 $\vec{X}$ $\vec{X}$ ^o可行, 则进行动态全域映射收缩算子的计算;若找到改善点r, 则利用r替换掉b构成新复形, 赋order=1转2) , 若没有找到改善点r, 则order=order+1转2) ;若 $\vec{X}$ $\vec{X}$ _o不可行则order=order+1转2) 。

计算终止准则有两个, 其一为order>6;其二为 $| S (\vec{X}_{g}) - S (\vec{X}_{b}) < ε, ε = 0.001 ‚ g$ $| S (\vec{X}_{g}) - S (\vec{X}_{b}) < ε, ε = 0.001 ‚ g$ 代表目标函数值最小点, b代表目标函数值最大点。程序计算过程中满足任一终止准则即停止迭代。

《5 算例》

5 算例

给出两个复杂土坡的剖面及土质力学参数分别如图4、图5及表1、表2所示。

对S1土坡, 给出两组搜索域:

$\begin{array}{l} \vec{X}_{U 1} = (90, 80, 45), \vec{X}_{L 1} = (10, 10, 25) ‚ \\ \vec{X}_{U 2} = (50, 60, 45), \vec{X}_{L 2} = (10, 0, 25) 。 \end{array}$ $\begin{array}{l} \vec{X}_{U 1} = (90, 80, 45), \vec{X}_{L 1} = (10, 10, 25) ‚ \\ \vec{X}_{U 2} = (50, 60, 45), \vec{X}_{L 2} = (10, 0, 25) 。 \end{array}$

对S2土坡, 给出两组搜索域:

$\begin{array}{l} \vec{X}_{U 1} = (100 ‚ 80 ‚ 58) ‚ \vec{X}_{L 1} = (0, 10, 28) ‚ \\ \vec{X}_{U 2} = (50 ‚ 50 ‚ 58) ‚ \vec{X}_{L 2} = (10, 0, 28) 。 \end{array}$ $\begin{array}{l} \vec{X}_{U 1} = (100 ‚ 80 ‚ 58) ‚ \vec{X}_{L 1} = (0, 10, 28) ‚ \\ \vec{X}_{U 2} = (50 ‚ 50 ‚ 58) ‚ \vec{X}_{L 2} = (10, 0, 28) 。 \end{array}$

对每组搜索域由程序随机生成50组初始复形 (每组复形的顶点数为6) , 这样共有200 组初始复形产生。限于篇幅, 本文不给出200组初始复形顶点的具体位置, 仅简要介绍程序的伪代码。

《图4》

图4 S1土坡剖面图

Fig.4 The cross-section of first slope

《图5》

图5 S2土坡剖面图

Fig.5 The cross-section of second slope

for I=1, 300 (50组初始复形顶点总数)

loop X=l₁+rand* (u₁-l₁)

Y=l₂+rand* (u₂-l₂)

Z=l₃+rand* (u₂-l₃)

表1 S1土坡的土质力学参数

Table 1 Soil properties of first slope

《表1》

层次	土类	γ/kN·m^-3	C/kPa	ϕ/ (°)
1	素填土	19.22	0.0	35.0
2	粉土	19.6	6.0	38.5
3	淤泥质粘土	18.5	6.2	0.0
4	粉土	19.8	5.0	36.0

表2 S2土坡的土质力学参数

Table 2 Soil properties of the second slope

《表2》

层次	土类	γ/kN·m^-3	C/kPa	ϕ/ (°)
1	素填土	19.22	0.0	35.0
2	粉土	19.6	8.0	32.5
3	淤泥质粘土	18.5	10.5	0.0
4	粉土	19.8	10.0	30.0

IF (X, Y, Z) 确定的潜在滑动面可行 then

保存该点

else

go to loop

end if

next I

若记基本复合形法 (α_ini=1.3) 为F1, 基于最大熵原理的复合形法为F3, 动态全域映射收缩复合形法为F2, 分别用F1, F2, F3对每组搜索域下的50组初始复形进行了计算, 以得到的最小安全系数为横轴, 搜索到某安全系数的初始复形数为纵轴, 绘成的土坡复形浓度分布如图6所示。

由图6 中a, b, c, d图可以看出, F1计算下, 每种搜索域生成的50组初始复形只有20～30组寻找到最优解, 局部最优解2.2附近集中了10多个初始复形, 若定义寻优成功率=集中在最优解附近的复形数/50, 则基本复合形法的寻优成功率介于40%～60%之间。F2计算下, 每个搜索域下的50组初始复形有30～35组集中在最优解附近, 局部最优解2.2附近也集中了五六个初始复形, 其寻优成功率达到60%～70%, 这有力地说明了基本复合形法因为搜索不完全而导致的较低寻优成功率。而采用F3计算下, 每个搜索域的50组初始复形均有40～45组集中在最优解附近, 局部最优解2.2附近没有集中初始复形, 只有少数的几个初始复形集中在局部最优解1.7～1.9附近。其寻优成功率达到80%～90%, 这说明利用信息量最大的顶点比利用目标函数值最大的顶点来搜索改进点更能有效地开发解空间, 收到较好效果。

《图6》

图6 土坡复形浓度分布

Fig.6 The distribution of complex density for soil slopes

《6 结语》

6 结语

从更深层次上对基本复合形法的寻优机理进行了剖析, 结合最大熵原理确定了复形中信息量最大的顶点, 在信息量最大的顶点以及其余顶点的加权中心点的连线上开发未知的解空间, 当复形的全部信息被利用后算法结束。同时针对基本复合形法搜索不完全的缺陷提出了动态全域映射收缩算子来对已知的解空间进行更为穷尽的搜索, 改进后算法的寻优成功率提高了40%。

展示更多