第4章 确定最小安全系数的最优化方法
4,1 概述
4,1,1 应用最优化方法确定最小安全系数边坡稳定分析应包含下面两个步骤
1) 对滑坡体内某一滑裂面按第2或第3章介绍的方法确定其抗滑稳定安全系数
2) 在所有可能的滑裂面中重复上述步骤找出相应最小安全系数的临界滑裂面在讨论了计算一个滑裂面的安全系数的方法后本章介绍边坡稳定分析极限平衡法的第二步就是要寻找最小安全系数的方法如果滑裂面曲线为y(x)那么这个问题具体化为寻找下列泛函的极值
)(yFF = (4.1)
岩土工程中的边坡的几何形状各异材料通常是非均质性纯解析的变分原理很难进行极值计算用最优化方法通过数值方法求解是一个比较现实可行的途径我国学者较早开展应用数值规划方法求解安全系数极值的问题闫中华,1983张天宝,
1981周文通,1984 20世纪80年代初期孙君实(1983)和Nguyen (1985)分别提出了使用复形法和单形法搜索任意形状和圆弧滑裂面的最小安全系数的方法80年代中期有更多的学者发表了有关研究工作(Celestino and Duncan,1981; Li,and White,1987) Chen &
Shao(1988)采用单形法和牛顿法进行任意形状滑裂面搜索的研究成果这篇论文所附的几道例题后来在国外多篇论文中被引用作为检验新的优化方法的考题借此书的机会作者将这些例题和原数据文件在本章介绍国内外有关的研究成果和STAB软件十几年在全国推广的实践证明应用计算机自动搜索临界滑裂面是可行的同时作者也发现相对于三维或二维斜分条的极限平衡法垂直条分法的极小值搜索问题比较简单采用任何一种优化计算方法配合本章介绍的随机搜索法即可快速地找到临界滑裂面因此目前在STAB程序中向用户提供的只是单形法这个功能以后作者还将把本章中介绍的各种方法补充进程序供用户选用
4,1,2 最优化方法最优化方法是近代数学规划中十分活跃的一个领域目前已有许多十分成熟的计算方法这些计算方法总的来看可以分为以下三大类
1,枚举法枚举法的基本思想是根据一定的模式比较不同自变量的目标函数经过筛选最
88 土质边坡稳定分析?原理
方法
程序终找到最小值这是最原始简单的方法如图4.1中示任一圆弧可用其圆心坐标(x
0
,y
0
)和半径r确定其相应的安全系数F可表达为
),,(
00 s
DyxfF = (4.2)
式中D
s
为滑弧深度即圆弧最低点的坐标可知
0
yrD
s
+= (4.3)
显然这是三个自由度的问题采用枚举法不断地改变x
0
,y
0
和D
s
的数值逐一比较相应的安全系数最终找到最小的安全系数在具体操作中先固定一个D
s
然后在圆心可能的位置中布置一个网格见第12章图12.1设网格的中心坐标为x
c
和y
c
在左右方向各布置了N
x
格在上下方向各布置了N
y
格则共计有(2N
x
+1)×(2N
y
+1)个网格点分别以该网格点为圆心以D
s
为滑弧深度计算相应安全系数找出最小的安全系数然后改变一个D
s
值重复相同的步骤在这一过程中有可能出现圆弧和边坡不相交的情况则应令程序自动抛弃该圆弧同样D
s
也是以一个中心值起算在其上下各布置N
d
层这样总计计算(2N
x
+1)×(2N
y
+1)×(2N
d
+1)个圆弧在枚举法的基础上用0.618法或其他方法提高搜索最小安全系数工作的效率这方面的工作可参阅有关文献图 4,1 圆弧滑裂面
2,数值分析方法随着计算机的发展数值分析方法逐步形成一门完整的学科统称为最优化方法(Method
of optimization)这一方法又可分为两大类第一类称模式搜索法(Pattern search method)其基本思想是根据一定的模式比较不同自变量的目标函数确定最优的搜索方向最终找到最小值另一类称牛顿法它要通过解析手段寻找使目标函数F对自变量z
i
的偏导数为零的极值点(?F/?z
i
=
0,i
=
1,2,...,n)同时从理论上讲还需要满足由二阶导数形成的Hessian矩阵正定这个达到极小值的充分条件此类方法中以导数为研究的主要对象因此也称为以导数为基础的方法(Gradient?based?method)一般认为当自由度较多时直接搜索法效率较低此时需要考虑牛顿法体系的分析方法第4章 确定最小安全系数的最优化方法 89
本节介绍直接搜索法中的单纯形法和Powell法以及牛顿法中负梯度法和DFP法
(Davidon?Fletcher?Powell method)由于这些方法的原理在众多的文献及教科书中都有所介绍这里不拟全面阐述其原理而着重讨论在边坡稳定分析领域需解决的特殊问题
3,非数值分析方法第三类方法是在近期计算机发展基础上形成的称为非数值方法这一方法被广泛应用于管理科学计算机科学分子物理学及超大规模集成电路设计中用于解决组合优化问题非数值分析利用计算机具有容量大计算速度快的优点通过大量随机采样来找到目标函数的最优值近期涌现了诸如模拟退火遗传算法神经网络和蚂蚁算法等均属此类在边坡稳定分析中上述算法均有人尝试本章4.7 节将对此类方法作一简要介绍
4,2 任意形状滑裂面的模拟和目标函数的确立
4,2,1 任意形状滑裂面的模拟最优化问题的提法是对于一个具有n个自变量的向量 z
=
(z
1
,z
2
,...,z
n
)确定其目标函数F的最小安全系数F
m
相应的自变量为z
m
为此需要对式(4.1)中的曲线y(x)用若干参数来模拟也就是说需要将任意形状滑裂面y(x)用z来近似表达将滑裂面曲线用m个点A
1
,A
2
,…,A
m
离散见图4.2也就是将此m个点用直线或光滑的曲线连起来以近似模拟此曲线此m个点的坐标用z
i
(i
=
1,2,,m)
表示
=
i
i
i
y
x
z (4.4)
图 4,2 任意形状滑裂面一旦这种连接的模式确定安全系数F即可表达成此m个点的坐标x
1
,y
1
,x
2
,y
2,
...,x
n
,y
n
的函数
()
mm
yxyxyxFF,,,,,
2211
L= (4.5)
90 土质边坡稳定分析?原理
方法
程序在进行最优化搜索过程中A
1
,A
2
,…,A
m
将移到临界滑裂面的位置
m21
B,...,B,B ′′见图
4.2此处m
=
6其中端点A
1
,A
m
原来在边坡线上有可能移到边坡线外或内如图4.2中的
m1
B,B ′′为此需要通过一定的处理方式分别找到它们和边坡线的交点B
1
和B
m
仍以
m21
B,...,B,B作为研究的对象对均匀的土质边坡通常希望滑裂面比较光滑此时可采用三次或更高次的样条函数连接这些点第11章将详细介绍使用样条函数构筑光滑曲线滑裂面的方法当然也可以采用直线和光滑曲线的组合构筑滑裂面例如图4.2,A
3
,A
4
,A
5
,
A
6
用曲线相联A
2
,A
3
用直线相联应用样条函数构筑光滑滑裂面对于减少自由度提高数值计算效率具有重要意义
[例4.1] 折线和光滑滑裂面比较如图4.3所示例如果用四个点模拟滑裂面那么折线滑裂面和曲线滑裂面相应的最小安全系数分别为1.489和1.364差别颇大换句话说如果使用折线模式要得到与曲线模式相同的精度就要更多的离散点这意味着更多的自由度而数值分析的机时和收敛难度是随自由度按指数增加的值得注意的是在笔者所知的所有这方面的文献中还没有其他研究者使用过曲线模拟任意形状滑裂面在稳定分析中对土质比较均匀的边坡常采用圆弧滑裂面此时自变量可以是圆弧的圆心的x,y坐标和半径r由于自变量仅三个最优化方法搜索最小安全系数收敛性能很好因此将不作为重点在本章中讨论图 4,3 说明用不同模式模拟滑裂面的差别滑裂面1?直线模式滑裂面2?曲线模式
4,1,2 目标函数的确定当滑裂面z的离散模型确定后安全系数便是m个控制点A
1
,A
2
,...,A
m
的函数在优化计算过程中这m个点中有n个点各沿某一设定方向 β
i
向临界滑裂面移动或者不规定方向任其自由移动其余m?n个点由于问题本身的要求可以固定滑裂面上任意一点的z
i
可以用相应于一个初始滑裂面的相对坐标来代表见图4.2
=
o
i
o
io
i
y
x
z (4.6)
+=?+=
i
i
i
o
it
o
ii
d
β
β
sin
cos
zzzz (4.7)
式中i
=
1,2,…,n,d
i
为第i个点沿β
i
移动的距离第4章 确定最小安全系数的最优化方法 91
对于固定的点其自由度为零如图4.2中的A
2
,A
3
当某点需沿一个软弱夹层移动时则其自由度为1如对该点的移动方向无特殊要求如图4.2中的其它点该点的自由度为
2问题的自由度为各点自由度的总和搜索最小安全系数的问题具体化为求下列函数的最小值问题
),...,,,,,,(
2121 nn
dddFF βββL= (4.8)
如前所述滑裂面的一部分可能是光滑曲线另一部分受软弱夹层的控制为直线STAB
程序是这样规定这一定义格式的设滑裂面被m个点A
1
,A
2
,……,A
m
分为m?1段从上交点向下交点编号为1,2,…,m?1在STAB程序中用LNO代表此m?1段中为直线段的线段总数在数组LOO(I),I
=
1,2,…,LNO中存入这些直线段的编号其它线段则默认为曲线当LNO为零时滑裂面为一光滑曲线当LNO为一大于m?1的数时则程序默认为全部为直线段这两种情况都无须填写LOO(I)
4,3 模式搜索法
4,3,1 单形法朱伯芳等,1984
1,建立初始单形对某一初始向量z
0
按下面模式构筑n个向量z
i
(i
=
1,2,..,n)组成n+1个顶点
T
m
n
T
m
T
m
pzqzqz
qzpzqz
qzqzpz
],,,[
],,,[
],,,[
00
2
0
1
00
2
0
1
2
00
2
0
1
1
+++=
+++=
+++=
LL
LL
LL
LL
z
z
z
(4.9)
其中
a
n
nn
p
2
1)1(?++
= (4.10)
a
n
n
q
2
1)1(?+
=
(4.11)
式中a为初始步长可根据实际情况进行调试
2,确定搜索方向计算单纯形n+1个顶点的目标函数值比较其大小从中找出目标函数最大点z
H
和次大点z
G
和最小点z
L
然后按下式确定下一步的搜索方向
∑=
+
v
H
v
i
n
i
v
n
n
zzz
1
2
1
(4.12)
式中上标v表示迭代次数
92 土质边坡稳定分析?原理
方法
程序显然从
ν
H
z到
ν
2+n
z目标函数是下降的因此以
ν
H
z和
ν
2+n
z两点的联线作为搜索方向
3,优化计算
(1) 反射沿搜索方向
ν
H
z和
ν
2+n
z点向前再走一步步长为)(
2
v
H
v
n
zz?
+
α到达
ν
3+n
z点
)(
223
v
H
v
n
v
n
v
n
zzzz?+=
+++
α (4.13)
ν
3+n
z称为反射点α>0称为反射系数计算反射点的目标函数并根据它的大小来决定下一步的走法
(2) 扩张如果反射点的函数值小于z
L
点的函数值即
)()(
v
L
FF zz
v
3n
<
+
(4.14)
则表明反射后情况有所改善沿搜索方向
v
n
v
H 2+
z,z还可以试探一下是否可以走得更远一些即是否可以扩张到
v
n 4+
z点
)(
2324
v
n
v
n
v
n
v
n
v
++++
+= zzzz (4.15)
式中v>1为扩张系数如果)()(
4
v
L
v
n
FF zz <
+
则以
v
n 4+
z点替换原来的最坏点
H
z构成新的单纯形转入第4
步进行收敛判断如果)()(
4
v
L
v
n
VV zz ≥
+
则以
v
n 3+
z点替换
v
H
z构成新的单纯形也转入第(4)步进行收敛判断
(3) 收缩如果反射点的函数值大于次坏点的函数值即
)()(
3
v
G
v
n
VV zz >
+
(4.16)
表明反射点走得太远需要缩回去按下式计算
)(
225
v
n
v
H
v
n
v
n +++
+= zzzz β (4.17)
式中0
<
β
<
1是收缩系数用
v
n 5+
z代替
v
H
z构成新的单纯形并转入第(4)步
(4) 缩小边长如果反射点函数值大于最坏点函数值即)()(
3
v
H
v
n
zVzV >
+
则缩小单纯形的边长以最好点
L
z为顶点其它各顶点向
L
z移近一半距离即按下式计算
ni
v
L
v
i
v
i
,...,2,1,0 ),(5.0 =?+= zzzz
v
L
(4.18)
得到新的单纯形转入第(4)步重复计算
4,收敛判断按照一定的方式通过反射扩充和收缩使单形不断更新逼近极值点收敛准则为
[]ε≤?∑
+
+
2
2
0
)()(
1
1
v
n
v
i
n
i
vv
n
zz (4.19)
第4章 确定最小安全系数的最优化方法 93
式中ε为要求的计算精度如式(4.19)满足则结束计算并以
v
L
z作为极小点否则置v
=
v+1转至第3步重复计算据已有经验可取α=1; 0.4
≤
β
≤
0.6; 2.0
≤
γ
≤
3.0经Nelder和Mead论证为使单纯形适应函数的性态及便于收敛α不宜比1大很多而α
<
1的计算次数比α=1要多故折衷取α
=
1另外认为β的数值变化对搜索效率的影响比γ要大他们推荐采用
6.04.0 ≤≤β (4.20)
0.38.2 ≤≤γ (4.21)
[例4.2] 说明单形法计算过程例为了形象直观地了解单形法在搜索最小安全系数时的工作状况我们考察图4.4所示的一个有两个自由度的例子滑裂面它由ABC组成计算时令C点固定不动A,B两点沿水平线移动则该滑裂面的安全系数F由A点的x坐标x
1
和B点的x坐标x
2
决定图4.5
示F相应x
1
,x
2
的等值线图根据枚举法可以发现在x
1
=
92.0,x
2
=
143.0时安全系数获得最小值1.257相应临界滑裂面如图4.4中标5的那个滑裂面图 4,4 具有两个自由度的算例
1,2,3,4?初始滑裂面5?临界滑裂面如果使用单形法则按式(4.9)初始生成的单形即三个滑裂面如图4.5在左下角三角形所示最大次大和最小分别为A B C三点D为按式(4.12)所得的新的顶点AD代表了两点联线的方向D代表了第一次反射和扩张后达到的点第一次迭代后形成了B
C D构成的新的单形开始新的一轮迭代依次循环最终达到安全系数的极值E搜索过程如图4.5中折线1所示最终收敛到F
m
=
1.257相应的z
m
=
(92.00,143.00)
T
94 土质边坡稳定分析?原理
方法
程序图 4,5 对图4.4示例用单形法计算过程图 4,6 具有软弱夹层的算例
[例4.3] 使用单形法计算最小安全系数图4.6具有分层剖面和一个倾斜的软弱夹层图中示4个不同的初始滑裂面点A和
B可以不受限制地自由移动故自由度为2同时此两点有意未设在边坡面上以说明程序在搜索过程中将进行判断和调整自动找到滑面与坡面的交点C和B点位于夹层上因此自由度为l该两点在搜索过程中将沿一定方向β移动β为该夹层的倾角图4.7中滑裂面3为采用单纯形法获得的临界滑裂面相应安全系数为1.025
第4章 确定最小安全系数的最优化方法 95
图 4,7 对图4.6示例用最优化方法计算过程临界滑裂面1?改良后的DFP法F
m
=1.025 2?负梯度法F
m
=1.025
3?单形法F
m
=1.025 4?DFP法F
m
=1.025
4,1,2 Powell法
1,基本原理有关此法的详细原理可参阅文献朱伯芳等1984其基本算法如下
(1) 给定初始点z
0
并令
Ti
)0,...,1,...,0,0(=S (4.22)
式(4.22)右端的1处于第i位i
=
1,2,…,n
(2) 相应迭代步i
=
1进行一次一维搜索第一次搜索方向p
i
=
S
i
使极小=+
)(
1 iii
V pz α (4.23)
并令
iiii
pzz α+=
1
(4.24)
然后令i
=
i+1重复上述过程直至i
=
n为止
(3) 取p
i
=
p
i+1
,i
=
1,2,…,n?1并令
0
zzp?=
nn
(4.25)
通过这一步增加
0
zz?
n
丢弃原来的坐标方向
1
p
(4) 沿z
n
z
0
作一次一维搜索求α
n
使V(z
n
+α
n
p
n
)
=极小并令
nnn
pzz α+=
0
(4.26)
然后转向第(2)步
(5) 重复进行上述运算n次即沿n个互相共轭的方向都进行了一维搜索然后作如下判断
εα <
nn
p (4.27)
96 土质边坡稳定分析?原理
方法
程序如上述判据满足则结束计算此时z
0
即为所求极小点否则转向第(1)步
2,应用Powell法时所做的改进在采用Powell法进行计算时第一步需要确定n个独立的搜索方向在大多数情况下都采用与n个自变量z
i
坐标轴方向一致的单位向量即由式(4.22)定义的S
i
这种处理方式并不适于边坡稳定分析参见图4.8(a)所示一个有两个自由度例子按图 4,8 在边坡稳定分析中使用Powell法作的改进
(a) 未作改进(b) 改进后
Powell法第一步固定B点不同让A
点沿水平方向移动到A′使安全系数达到极小值显然能够移动的范围是很有限的第二步固定A′不动B在其移动方向到达B′再次得一个极小值其移动范围更为有限这样的搜索效率很低当自由度较多时事实上无法实现有效的搜索设想如果在A点移动时让B点也向前移动一定值第一次一维搜索就可以使初值z
0
向极值z
m
迈进一大步如图4.8 (b)
因此在构筑S
i
时我们的思路是当某一控制点在移动方向搜索最小值时要求其邻近的点也按一定的比例移动使滑裂面整体而不是局部在优化计算中不断地调整以迅速地达到极值点为此按下式构筑n个向量( ),...,3,2,1 ni
i
=S
Τ
= ),...,,(
21
i
n
iii
eeeS (4.28)
按下式确定( ),...,3,2,1 nie
i
j
=
>
≤
=
时当时当
ji
xx
xx
ji
xx
xx
e
in
jn
ij
i
j
ij
i
j
)cos(
)cos(
1
1
ββ
ββ
(4.29)
当i
=
1时总是用(4.29)的下式i
=
n时总是用式(4.29)的上式
i
j
e的物理意义是当第i个变量沿其搜索方向移动一个单位值时邻近的点j都将在其移动方向获得一个增量
i
j
e这个增量值取决于第j点距端点与j同侧的那个的水平距离在i
=
0和i
=
n时等于零i
=
j时等于1其余点按线性分布原则确定内插参见图4.9
第4章 确定最小安全系数的最优化方法 97
图 4,9 Powell法搜索方向
i
j
e的改进
(a) 1<i<n情况(b) i=1情况(c) i=n情况
[例4.4] 说明对Powell法所做的改进例本例和[例4.2]具有相同的计算条件采用折线滑裂面参见图4.4如果初始滑裂面设为1图4.5中线2和3分别代表未经改进和改进后的计算过程可见改进后搜索效率大大提高采用这种改进的方法要求各控制点都按同一方向朝临界滑裂面逼近例如图4.2所示例可以将初始滑裂面布置在离估计的临界滑裂面稍远的地方但A,B,C点都采用水平向左的移动方向即式(4.29)中的β均为零通常可以很快地找到极值
4,4 牛顿法
4,4,1 负梯度法负梯度法的基本思想是对一个初始滑裂面寻找一个使安全系数减少速率最大的方向数学上就是
T
n
z
F
z
F
z
F
,,,
21
L这个向量在这个方向上进行一维搜索找到这一方向安全系数的低谷点完成了这第一次迭代后再在这个新的起点即上述低谷区重复这样的运算直到收敛至极值点
[例4.5] 对[例4.2]使用负梯度法计算最小安全系数在使用负梯度法时分别以(84.0,160.0),(70.0,145.0),(112.0,150.0)作为起点相应滑裂面如图4.4中1,2,3示使用负梯度法的搜索路径如图4.10中1,2,3三条折线所示可见每一次搜索均是沿着等值线的法线方向即下降速率最大的方向进行的对这样一个比较简单的算例负梯度法相应不同初始输入的滑裂面都可以成功地找到极值
98 土质边坡稳定分析?原理
方法
程序图 4,10 对[例4.2]使用负梯度法和DFP法计算最小安全系数
4,1,2 DFP法
1,基本原理
DFP法为Davidon?Flefchen?Powell法的简称其原理和详细分析方法参见文献朱伯芳等1984
数学上的极值点是目标函数对各自变量的一次导数为零的点也就是说在极值点
),,,(
21 n
z
F
z
F
z
F
= LG均为零同样要求由二阶导数构成的海色(Hessian)矩阵正定海色矩阵由下式定义
=
2
2
2
2
2
2
2
1
2
21
2
2
1
2
.,,
.,,
n
n
n
z
F
zz
F
z
F
zz
F
zz
F
z
F
称对MO
H (4.30)
第4章 确定最小安全系数的最优化方法 99
DFP法试图通过一系列迭代步骤使假定的海色矩阵的逆阵A
0
逐渐过渡到其真值每一次迭代的方向均按下式决定
vvv
GAS?= (4.31)
在迭代步 v +1通过下式确定A
ν+1
vvvv
DCAA?+=
+1
(4.32)
其中 v
=
0,1,2,…
v
T
v
T
vvv
Yz
zz
C
= (4.33)
v
vT
v
T
v
v
v
v
v
YAY
YAYA
D
)(
= (4.34)
vv
v
GGY?=
+1
(4.35)
vv
v
zzz?=?
+1
(4.36)
DFP法即是基于这一基本思路的一种效率较高的方法DFP法的基本思路是通过迭代使A
0
,A
1
,…,A
n
逼近极值点处Hessian阵的逆阵A通常取A
0
为单位矩阵
2,对DFP法的改进本书作者在20世纪80年代开始研究使用DFP法计算安全系数的极值问题发现在使用这一方法时还需要针对边坡问题的特点解决量纲问题和导数的准确计算问题下面分述此两个问题
(1) 关于量纲问题如前所述DFP法的基本思路是通过迭代使A
0
,A
1
,…,A
n
逼近极值点处Hessian阵的逆阵在边坡稳定分析中?
2
F/?z
i
z
j
的量纲通常在10
2
～10
5
的范围内H
1
的各元素数值的量纲大约在10
2
～10
5
范围内常规的DFP法采用A
0
为单位矩阵那么要从A
0
过渡到10
2
~10
5
的量级就变得十分困难为此改进的方法取
p
101
0
×=A (4.37)
式中p为对H
1
量纲的平均值的一个估值可按下式确定
1
1
122
10
|/|
10
+
=
<
<
∑
p
n
i
ip
n
zF
(4.38)
[例4.6] 对[例4.2]使用改进后的DFP法计算最小安全系数同样我们通过[例4.2]计算以图4.4中的第4个滑裂面为初始滑裂面如果不采用本节引入的改进则根据式(4.32)第一次迭代获得搜索方向为
T1
)7.1,3.19(=S以图4.10中
BC
1
代表沿着该方向显然无法到达极值点D如果看一下计算这个方向的细节
100 土质边坡稳定分析?原理
方法
程序
=
+
=
+=
24.3 217.7-
217.7- 2303.8
0.408 0.492-
0.492- 5917.0
23.7 218.2-
218.2- 4.2003
1 0
0 1
0001
DCAA
可以看到(A
0
D
0
)与C不是同一量级的量(A
0
D
0
)对A
1
的贡献非常小如果使用本节提出的改进那么根据式(4.38)可得到p
=
3,A
1
计算细节为
=
+
=
+=
615.7 273.8
273.8 2411.7
408.0 492.0-
492.0- 5917
23.7 218.2-
218.2- 4.2003
1000.0 0.0
0.0 1000.0
0001
DCAA
可以看到使用了改进的方法后(A
0
D
0
)和C
0
具有相同的量级据此算得的
,)45.3,07.3(
1 T
=S为图4.10中BC所示方向这就保证了第二次迭代快速地向极值点逼进并顺利地达到极值点
(2) 关于导数的计算问题牛顿法体系的最优化方法均是建立在寻找目标函数的导数的基础上的准确地确定?F/?z
i
的数值是个十分重要的问题我们遇到的困难是
(a)?F/?z
i
不能通过一个公式直接算得因此只能通过差分法计算?F/?z
i
的近似值?F/?z
i
(b) 由于F本身只能通过迭代求解获得因此F的有效位很少在差分计算时步长?z
i
也就不能取得很小这就使算得的导数值的误差较大另一方向当趋近极值点时导数值本身却趋近于零这两个因素迭加在一起就使计算所处得的导数在极值点附近相对误差急剧增加导致数值计算无法顺利地收敛到极值点大部分数值分析问题都不可能通过解析方法算得导数值在使用牛顿法进行非线性分析时解决在差分运算中导数计算精度问题是一个具有普遍意义的问题因此本节结合边坡稳定分析的特点讨论使用中心差分时步长的合理选定问题同时提出了一个计算导数的半解析法使用这个方法可以获得较精确的导数值
(c) 用差分法确定导数的步骤首先简单回顾Stewart (1967)讨论的有关差分计算中步长确定的原则对于某一函数)(zφφ =用向前差分计算dφ / dz在z
=
0处值的公式是
δ
φδφφφ )0()(
d
d?
=
≈
zz
(4.39)
中心差分计算公式是第4章 确定最小安全系数的最优化方法 101
δ
δφδφφ
2
)()(
d
d
≈
z
(4.40)
式中δ为选用的步长使用式(4.39)计算dφ/dz其误差来自两个方面
(a) 截断误差?φ的准确值为
...
2
1
=
)0()(
2
++
=?
ρδγδ
φδφφ
(4.41)
式中γ,ρ分别为?φ/?z,?
2
φ/?z
2
在z
=
0处的值因此用式(4.39)计算γ实际上忽略了ρδ
2
/2
由此引起的相对截断误差η
t
近似为
γ
ρδ
η
2
2
1
=
t
(4.42)
(b) 舍入误差任何一个数值在计算机中都是一个有限字长的数字近似代表的如果δ
选取得太小计算机将无法识别φ(δ)和φ(0)之间的区别由于φ的有限精度造成的按式(4.39)
计算的γ的误差称舍入误差舍入误差η
c
按下式估算
η
φ
φ
η
≤
)0(
2
c
(4.43)
式中η为φ(0)的相对误差如果φ(z)的数值是通过某一公式对自变量z进行若干算术运算得到的那么η值首先取决于z本身的相对误差在计算机中z是被一个在区间[ )1(),1(
mm
zz εε +? ]的某个数值近似代表的其中ε
m
称为机器误差它的大小取决于计算机的种类和精度级别(单双精度)
可以按图4.11所示框图编一个简单的程序来测试ε
m
的数值参见Dennis and Shanabl,1983
年对于16位的计算机可以知道ε
m
=
5.96×10
8
图 4,11 确定计算机精度的计算框图
η值还取决于在进行φ(z)的算术运算时积累的误差它使φ(z)的相对误差η比z的相对误差ε
m
大几倍或几十倍如果φ(z)的数值不能用显式计算需要通过迭代求得则η取决于迭代精度例如稳定分析的安全系数F它是通过迭代求解式(2.31)和式(2.32)获得的在单精度运算时迭代收敛标准是10
5
由于安全系数的数值在1左右故可大致认定F的相对误差
η<10
4
Steward提出的确定步长δ的步骤如下
102 土质边坡稳定分析?原理
方法
程序
(a) 要求式(4.42)的和式(4.43)右端相等到即要求舍入误差和截断误差平衡据此计算
δ值
(b) 用式(4.39)计算γ的值并按式(4.42)估算η
t
如果η
t
超过了允许值ε(例如ε = 5%)
则要用中心差分即使用中心差分式(4.40)计算γ
(c) 如果使用中心差分则δ通过解下式求得
0)0(2
2
1
2
=?+φηδγεεδp (4.44)
图4.12示dφ/dz数值其中tanα
0
为准确值tanα
1
为中心差分计算值tanα
2
为向前差分计算值tanα
3
为向后差分计算值式(4.44)是令式(4.43)右端等于ε获得的其中?φ按式(4.41)确定在稳定分析中由于F的精度低故总是用中心差分因此我们不进行
Steward提出的第一步计算在用式(4.44)
计算δ时取ε
=
5%,η
=
10
4
ρ,γ的值近似采用在DFP法计算时上一次迭代已经算得的值假定φ在(?δ,δ)区间呈抛物线分布这样γ值仍按式(4.39)计算ρ的计算公式为图 4,12 向前差分可能获得相反符号的导数示意图
2
2
)0(2)()(
δ
φδφδφ
ρ
+
= (4.45)
采用中心差分的另一个理由是由于F的精度低δ的数值相对较大当快要接近极值时向前差分可能得到与真实值符号截然相反的导数值见图4.12中心差分可在很大程度上避开这一危险
3,半解析法计算导数的步骤通用条分法的特点是在计算F值时各项导数不通过差分而是通过解析公式获得的利用这一特点我们可以将计算?F/?z
i
所包含的差分运算减少到最低量设想由于z
i
的一个增长量?z
i
使式(2.31)和(2.32)中的F和λ得各有一增量?F和?λ则
()0,,=?+?+?+=
iinn
zzFFGG λλ (4.46)
0),,( =?+?+?+=
iinn
zzFFMM λλ (4.47)
全微分式(4.46)和式(4.47)可得第4章 确定最小安全系数的最优化方法 103
i
n
i
n
i
n
z
G
z
G
z
F
F
G
=
+
λ
λ
(4.48)
i
n
ii
n
zz
M
z
F
F
M
Μ?
=
+
λ
λ
(4.49)
由此可得
F
ΜGM
F
G
z
GΜ
z
ΜG
z
F
nnnn
i
nn
i
nn
i
=
λλ
λλ
(4.50)
式(4.50)给出了一个不直接计算F值确定?F/?z
i
的公式式(4.50)右边的?M
n
/?λ,?M
n
/?F,
G
n
/?λ,?G
n
/?F都是在计算一个滑裂面安全系数时已经算得的而且是通过精确的计算公式确定的即第2章中式(2.37)至式(2.40)?G
n
/?z
i
,?M
n
/?z
i
仍然需要通过差分求得但是由于计算G
n
和M
n
的公式(2.31)和式(2.32)不需迭代故其导数的精度远比对F差分精度高用中心差分计算?G
n
/?z
i
和?M
n
/?z
i
时δ的选用原则是要求按式(4.43)确定的舍入误差相对值η
c
小于允许值ε G
n
和M
n
是通过式积分算得的在两个相邻的控制点1,+ii之间被分成了
1+i
N个土条因此要求在进行差分运算时每个土条高度h的相对误差都小于规定值当第i个控制点在其移动方向z有一增量δ 时其(i,i+1),(i?1,i)两段的每个土条的高度
h都将获得一个增量h?其中h获得增量最小的那个土条相应值为δ
/N
m
这里
),max(
1+
=
iim
NNN,
1
,
+ii
NN分别为(i?1,i)和(i,i+1)段的土条数分别以η
c
=
ε,?φ
=
δ
/
N
m
,φ
=
z
m
代入式(4.43)可得一确定δ的近似公式
η
ε
δ
mm
Nz
2≥ (4.51)
式中η为z
m
的相对误差可取η
=
ε
m
选用5%~10% z
m
为x
i
,y
i
中的大值下面列举一个例题以验证本节建议的计算导数的方法的可行性
[例4.7] 测试计算导数精度例题图4.13示一垂直均匀边坡假设滑裂面为通过坡趾A与水平线倾角为α的直线那么沿该滑裂面滑动的安全系数F仅为B点的x坐标值的函数并且安全系数F和dF/dx均可通过理论推导获得:
αγ
αφ
2sin
4
cottan
h
c
F += (4.52)
x
h
c
dx
dF
/)2cot4cot(tan
γ
ααφ+= (4.53)
104 土质边坡稳定分析?原理
方法
程序图 4,13 测试计算导数精度例题当F获得极值时
s/21tan +=α (4.54)
φγ tan
4
h
c
s = (4.55)
对图4.13示边坡取h
=
25m,φ
=
35
°
,c
=
49kPa,γ
=
1.8g/cm
3
分别采用差分和半解析法计算dF/dx并与理论值对照列于表4.1
此例当x
=
12.283时获得安全系数极小值F
m
=
0.9056表4.1所列计算结果表明对于本例半解析解和差分法均可以获得和理论值完全一致的最小安全系数和临界滑裂面表 4,1 理论值和近似法结果对照
x 25 8 5 12.283
理论值式4.52) 1.1446 0.09896 1.2956 0.9056
F
程序计算1.1447 0.9896 1.2956 0.9055
理论值式4.53 0.02801 -0.0499 -0.1854 ≈0
δ 0.019 0.011 0.011 0.011
半解析法
dF/dx 0.02801 -0.0499 -0.1853 0.713×10
-4
δ 0.035 0.078 0.014 0.378
x
F
d
d
数值差分法
dF/dx 0.02806 -0.0499 -0.1854 0.367×10
-4
[例4.8] 对[例4.3]使用DFP法计算成果用经过改进的DFP法计算[例4.3]成果示于图4.14临界滑裂面1至4相应图4.6
中的四个不同的初始滑裂面可见基本上收敛在相同的位置F
m
在1.009~1.025之间
4,5 确定整体极值的随机搜索方法
4,5,1 概述数学中的最优化理论研究的目标函数在其自变量空间是光滑连续的同时要求在极值点附近目标函数是自变量的二次函数而岩土工程中的边坡一般为非均质的因此在实际应用时经常遇到多极值问题在使用数值分析方法进行边坡稳定的最大化分析时笔者发现几乎所有的方法都存在陷入局部极值陷井无法获得整体极值的问题第4章 确定最小安全系数的最优化方法 105
图 4,14 对[例4.3]使用负梯度法和DFP法计算成果
[例4.9] 说明确定整体极值困难一例如图4.15所示已知边坡的整体极值为1.364对于选定的那个初始滑裂面1,F
o
=
1.593
单形法和DFP法计算均终止在线2相应安全系数F
=
1.452其相应临界滑裂面和已知的实际位置相差甚远经验表明在下面的两种情况下数值计算的收敛问题变为极为严峻
(1) 当问题包括较多的自由度时可以设想此时目标函数的自变量空间存在许多局部极值点使搜索整体极值的任务变得十分困难
(2) 当边坡断面较复杂时此时也会出现许多安全系数的局部极值图 4,15 说明确定整体极值困难一例
1?初始滑裂面; 2?应用单形法获得的临界滑裂面如图4.15所示计算终止的那个滑裂面的一部分实际上和两种不同土层的交界线重合很可能是一个局部极值点这道题如果改成均质土坡同样的初始滑裂面位置则计算所得的临界滑裂面可很好地收敛到整体极值点的位置由此可见为了使最优化方法真正成为一个在边坡稳定分析中实用性强的工具还需作进一步的努力有效地解决确定整体极值问题在这方面一个十分简单又十分有效的思路是设法确定一个靠近整体极值的初值显然初值离整体极值愈近丢失整体极值的可能性越小不少教科书中都指出确定一个好的初值的重要性同时指出这个任务可以由随机搜索圆满完成(Fox,1971; Shoup,et,al,1987; Box,1965)
随机搜索的基本思想是应用随机数在研究域内构筑一系列自变量比较其相应的目标函数寻找最小的目标函数这样在进行常规的最优化计算以前对自变量空间进行一次均匀的高密度的扫瞄把整体极值的大体位置确定下来再将这个最小目标函数作为初值进行最优化计算
106 土质边坡稳定分析?原理
方法
程序根据Monte Carlo法的原理当随机搜索的次数趋于无穷大时目标函数的最小值就是整体极值因此随机搜索本身也可以独立地成为一种计算最小安全系数的方法已经有一些文献报导这方面的研究成果(Boutrup,1980; Siegel,1981)在笔者发表了有关使用随机搜索和常规最优化方法相结合的方法求解安全系数极值的论文后Greco (1997)针对笔者的论点作过讨论他认为在一定模式指导下纯随机搜索方法一样可以精确地确定临界滑裂面但是这种纯粹的随机方法通常需要几千几万次的随机搜索才能达到足够的精度而本节所研究的把随机方法和确定性方法结合起来的途径在数值计算效率和收敛性能两个方面均具有明显的优势
4,1,2 最小安全系数的随机搜索对于某一边坡根据问题特点确定一个搜索区域如图4.16所示其轴线用z
o
表示其宽度为D
i
半带宽为d
i
=
D
i
/2这个搜索区域左右边界分别用两个滑裂面z
1
,z
r
来代表即
=
i
i
i
o
i
l
i
d
β
β
sin
cos
zz (4.56)
+=
i
i
i
o
i
r
i
d
β
β
sin
cos
zz (4.57)
式中D
=
(D
1
,D
2
,,D
n
)
T
称为搜索区宽度图 4,16 生成随机滑裂面示意图在搜索区内任意一个滑裂面可用下式表示
+=
i
i
ii
o
ii
dr
β
β
sin
cos
)5.0(zz (4.58)
式中
t
n
rrr ),...,,(
21
=r为该滑裂面和各控制点相对于轴线距离的系数
n
rrr,...,,
21
为伪随机数其值均在(0,1)之间随机搜索的步骤如下第4章 确定最小安全系数的最优化方法 107
1) 计算相应于滑裂面z
1
的安全系数F
0
2) 使用计算机伪随机发生器产生n个随机数
n
rrr,...,,
21
应用式(4.58)确定一个滑裂面z计算其相应的安全系数F
1
3) 比较F
1
和F
0
如果F
1
小于F
0
则F
0
和z
0
用F
1
和z
1
更新否则直接转入步骤(4)
4) 重复步骤2和3直到比较的次数足够大获得的最小安全系数足够小作为最优化法的初值足够好为止由于计算机产生的伪随机数具有很好的均匀性可以认为搜索区域内的滑裂面空间的每个部分都机会均等地被扫描了一遍搜索次数越多扫描密度越高成果越佳应用上述步骤在求解某一边坡的具体问题时需要根据经验确定一个搜索区域和一个搜索次数现在需要解决的问题是如何确定一个合适的随机搜索次数N
设想所寻找的临界滑裂面存在着一个条带称为保证区间任何落在该区间内的滑裂面如果被用作是最优化方法的初值都可以保证计算收敛那么现在的问题就变成按上述的随机方法生成滑裂面要生成多少个即可保证至少有一个落在这个保证区间中我们进一步假定该保证区间的宽度d
i
和搜索区间宽度D
i
存在一个固定和关系
n
n
D
d
D
d
D
d
m,..
2
2
1
1
=== (4.59)
式中m称为置信系数(0
<
m
<
1)
那么当随机数发生器为滑裂面第一个控制点按式(4.58)生成一个相对坐标时该点落在保证区间中的概率为m根据概率论理论n个控制点全部落在保证区间中的概率为m
n
由二项定理可知在随机生成的N个滑裂面中有r个落在保证区间中有N?r个未落在保证区间中的概率是
rNnrnr
N
r
N
mmCp
= )1()()(
(4.60)
那么在N次中至少有一次落在保证区间中的概率称成功率是
Nn
Nnn
N
m
mmCP
)(11
)1()(1
00
=
=
(4.61)
举例在图4.15所示的例中假如用3个控制点并且规定了搜索区域如图4.16中条带部分所示如果假定m
=
0.4那么如果进行了50次随机搜索其成功率是
%96)5.01(1
253
==P
从另一角度说如果取m
=
0.3为了使成功率达到96%我们需要作122次随机搜索
%96)3.01(1
1223
==P
应用上述步骤在求解某一边坡的具体问题是需要根据经验确定一个搜索区域和一个相应此搜索区的置信系数m搜索区域划得越大说明对临界滑裂的位置推测的把握越小m值也越小这样根据式(4.61)所算出的需随机搜索次数也就越多一般来说根据经验是不难划定搜索区域并确定m值的也许反复试算调整几次才能找到较好的结果
108 土质边坡稳定分析?原理
方法
程序使用计算机这也不是一件很难的事因此在实际应用时也可以根据经验来确定搜索次数STAB程序默认的功能则自动设定一个搜索次数
4,1,3 提高随机搜索效率的途径在将这一套随机搜索的方法运用到实际解题中去以前为了提高随机搜索效率还需要解决两个问题
1,滑裂面的几何合理性问题在随机生成的滑裂面中应该删去那些几何形状明显不合理的滑裂面这些形状不合理的滑裂面的存在不仅耗费许多无谓的计算时间而且在计算这种滑裂面的安全系数时还会带来数值分析收敛的麻烦根据边坡稳定分析实际特点对滑裂面和几何形状作出如下限制
(1) 相邻两个控制点的x坐标值应满足:
ii
xx >
+1
(4.62)
(2) 相邻两个控制点所连的直线与x轴的夹角α
i
应控制在一个合理的范围内如
oo
8045 ≤≤?
i
α (4.63)
2,应用简化方法提高随机搜索的效率在自由度较多时为使搜索达到一定保证率需要的随机搜索次数很多因此有必要找到一个途径将随机搜索使用的时间尽量减少一个很有效的方法就是在随机搜索时采用简化方法计算安全系数因为随机搜索的目的是确定一个初始滑裂面过高的安全系数精度是没有必要的在第3.6节中详细介绍的简化法就是用于随机搜索的简化方法简化法1用于光滑形状的滑裂面简化法2用于折线形状滑裂面在一般情况下使用简化法可以使随机搜索的时间减少2/3至1/2
[例4.10] 采用随机搜索重新分析[例4.9]
回到本节一开始提到的那个搜索整体极值失败的[例4.9]如图4.15所示仍使用四个点并用光滑曲线相连来代表滑裂面搜索区间如图4.17所示采用简化法1进行45次随机搜索图 4,17 对[例4.9]采用随机搜索求解示意第4章 确定最小安全系数的最优化方法 109
图4.18示其中的一些典型的随机生成的滑裂面由于使用简化法很快即可在此45
个滑裂面中找到一个相应最小安全系数的滑裂面如图4.17线1所示以此滑裂面作为初始值任何一种最优化方法均很轻松地收敛到相应整体极值的临界滑裂面2上相应安全系数为1.364 Morgenstern?Price法图 4,18 对[例4.9]采用随机生成的滑裂面
4,6 应用实题本节介绍采用上述最优化方法进行边坡稳定分析的工程实例
110 土质边坡稳定分析?原理
方法
程序
[例4.11] 紫坪铺库区左岸堆积体在第2章[例2.1]中讨论了紫坪铺库区左岸堆积体沿夹泥层ABCEFGHIJ的稳定安全系数在进行详细核算时需分析各种可能的滑裂面的安全系数因此要对这一折线型的滑裂面分段核算如图4.19所示初始输入一个滑裂面LGHM相应F
0
=
1.476令G沿GH
方向移动优化计算获得的临界滑裂面为L′G′H′M′相应F
m
=
1.121这一计算确定了以GH
中的一部分为底滑面临界破坏模式下一步还要核算其它可能的破坏模式第13章13.7
节详细介绍了这一核算过程STAB程序为折线型的底滑面提供了分段进行优化计算的工具但不可能一次自动地就把整体极值找出来图 4,19 紫坪铺库区左岸堆积体例
[例4.12] 小浪底水库大坝沿软弱夹层滑动分析高167m的小浪底水库大坝地基中存在缓倾角软弱夹层需要复核大坝下游坝坡在8
度地震情况下沿此软弱夹层滑动的最小安全系数用五个点来模拟滑裂面图4.20中AB
BC均为曲线CD段为软弱夹层故为直线DE段较短采用直线首先划定一个搜索区域如阴影部分所示采用简化方法2作223次随机搜索作为研究作者也曾用严格的通用条分法计算这223个随机滑裂面发现使用简化法2和使用严格法发现的最小安全系数的滑裂面都是第193个但使用时间减少一半根据随机搜索获得的初值最后获得的上下游边坡临界面滑裂面如图4.20中ABCDE所示分析此临界滑裂面图中的点位置是软弱夹层上部的终止点在这样的位置滑裂面中软弱夹层的长度最大因而安全系数最小故计算成果和常识判断是一致的图4.20
中E点的位置肯定了这样一个见识滑裂面通常从坝趾处逸出此处的抗滑力最小表4.2对比了以下三种情况的数值计算效率和精度采用纯粹的随机搜索不采用随机搜索单纯应用最优化方法本节介绍的将两者结合起来的方法第4章 确定最小安全系数的最优化方法 111
可以看到单纯依靠随机搜索当达到2000次后安全系数改善得很慢机时达11170
秒所得到的最小安全系数为1.253与精确值1.199相差颇远说明Boutrup (1980)和Siegel
(1981)建议的纯随机搜索的方法不甚有效单纯使用最优化方法虽然得到较小的最小安全系数1.234但无论从所用的机时和还是从得到的安全系数的精度来衡量都不如经过随机搜索确定初值的方法有效需要说明的是上述讨论系1988年的成果耗费机时系相应
IBM?PC?XT微机图 4,20 随机搜索在小浪底水库边坡稳定分析中的应用表 4,2 三种情况下数值计算效率与精度对比状况单纯随机搜索次数单纯最优化方法两者结合的方法
223 605 1336 2356 0 223
机时秒699 3431 6437 11170 2940 2661
最小安全系数F
m
1.293 1.270 1.261 1.253 1.234 1.199
注 两者结合的方法先进行223次随机搜索然后再以获得的最小安全系数为初值使用单形法确定F
m
耗费机
时系相应IBM?PC?XT微机
[例4.13] 天生桥二级闸首滑坡反演分析该滑坡在水电工程中为一个值得吸取教训的实例发生于1985年12月24日48人丧生由于边坡中存在一层饱和软粘土见图4.21中第5层当边坡开挖将此层粘土揭露出后即发生了滑面滑面下部沿堆积土和基岩的接触面各土层的岩土力学参数如表4.3
所示初始滑裂面如图4.22中A
1
B
1
C
1
D
1
E
1
F
1
G
1
H
1
所示最终获得的临界滑裂面A
2
B
2
…G
2
H
2
临界滑裂面的位置在上部与实测的偏离较大下部则与实测的吻合也是沿基岩接触面下滑
112 土质边坡稳定分析?原理
方法
程序图 4,21 天生桥二级闸首滑坡工程地质图表 4,3 各层土的物理力学特性指标图4.21
抗剪强度指标土层特性密度ρ (g/cm
3
)
φ (°) c (kPa)
1新填土1.85 21.8 19.6
2旧填土1.85 21.8 19.6
3第四系人工堆积1.85 21.8 0.0
4中细砂1.85 20.8 29.4
5砂质淤泥1.81 10.2 34.3
6卵砾石及砂1.90 24.2 0.0
7灰岩2.40 45.0 39.2
第4章 确定最小安全系数的最优化方法 113
图 4,22 天生桥二级闸首滑坡计算
4,7 非数值方法
4,7,1 模拟退火法模拟退火法是模拟金属退火过程的一种最优化计算方法唐立山等1994在金属退火时温度徐徐降低在每一个温度阶段系统的能量都要达到平衡状态使其达到最小值假设在某一随机过程中获得了一个相应某一随机变量的安全系数为F(z
i+1
)如果它与已往已经得到的最优解F(z
i
)的差值?F
=
F(z
i+1
)-F(z
i
)小于零那么自然要用F(z
i+1
)来代替
F(z
i
)但是如果F大于等于零则传统的随机搜索方法就要抛弃这一选择在模拟退火法中根据Monte Carlo法的Metropolis准则还要作一次抽签试验也就是要让计算机产生一个在(0,1)之间的随机数r
i
将其与exp[F/T]比较其中T就是这一阶段的温度在这里T是一个具有安全系数F特征的量如果r
i
<
]/exp[ TF那么仍然要用z
i+1
来代替z
i
这一处理有利于跳过局部极限形成的陷井但是这一做法终究抛弃了一个比现有值更优的自变量因此是带有一定的风险的我们用承担一定的风险作为代价使计算跳出局部极限形成的陷井可以看出在这一决策过程中以下做法是明智的
(1) 在退火过程中温度是逐渐缓慢减少的初期T较大用z
i+1
代替z
i
的可能性也较大因而允许自变量在较大的范围内变动以利于跳过陷阱确定一个较好的方向越到计算接近终止这一替代的可能性越小
(2) 只有在?F的绝对值足够少时即F(z
i+1
)与F(z
i
)十分接近时F(z
i+1
)这个表面上欠优的安全系数才有可能入选
(3) 随机数r
i
足够的小时F(z
i+1
)才有可能入选某种意义上说这意味我们把一部分成功的希望寄托在运气上
114 土质边坡稳定分析?原理
方法
程序这里关键的问题是T的数值在开始的时候取得很大其变小的过程是十分缓慢因此即使这一选择错了以后仍然有机会纠正过来在实际应用这一过程时还采用了一种记忆功能如果在某一温度计算步结束时实践证明把搜索方向从F(z
i
)改为F(z
i+1
)
并没有带来更好的结果我们还是要将原来被抛弃的F(z
i
)取回来进入下一温度计算步根据上述原则模拟退火法的计算步骤如下
(1) 设置初始温度T
0
这一温度应足够高以防止计算落入局部极值的陷井同时过高的初始温度可能降低搜索的效率建议先进行一定数量比如100次的随机搜索找到安全系数的最大值和最小值F
max
和F
min
将初始温度设为
minmax
FFT
o
= (4.64)
将优化计算过程分成m段第k温度段的温度按下式确定
k
Tk
TT )(
0
α?= (4.65)
式中T
k
为第k步迭代时的温度α
T
为温度冷却的阻尼系数α
r
应在[0,1]区间内通过试算确定
(2) 对于一个具有n个自由度的自变量z产生n个随机数v
=
(v
1
,v
2
,…,v
n
)每个ν
i
都在(?1,1)区间并将其换成方向导数u
=
(u
1
,u
2
,…,u
n
)
2/1
1
2
)(
∑
=
=
n
i
i
i
i
v
v
u (4.66)
为了减少自变量空间的粗糙度我们还引入一个系数
)....,,2,1(
)}(,),(,),(max{
)(
11
ni
zzzzzz
zz
nnii
ii
i
=
=
+?+?+
+
κ (4.67)
式中z
和z
+
分别代替第i个自变量在探索区域的最大值和最小值这样相应某一迭代步的自变量搜索为
1
0
)(
=?
k
rk
rr α (4.68)
式中?r
k
为第k步的带宽同样要随着迭代步的增加逐步减少这样自变量的增量为
nidurkz
iikii
,,2,1 L==? (4.69)
第4章 确定最小安全系数的最优化方法 115
[例4.14] 采用模拟退火法重新分析[例4.9]
对于[例4.9]第4.5.3节曾以随机搜索找到整体极值现使用模拟退火法进行分析采用同样的初始滑裂面同样用4个点联成一条光滑的曲线位于坡趾的点在优化过程中不动其它3个点沿水平方向移动因此共有3个自由度计算中采用α
T
=
0.9,α
r
=
0.98首先通过100次随机搜索生成初始温度T
0
随后进行m
k
=
300次退火模拟相应?r
o
=
0.27,T
k
=
1.006和1.001
从图4.23中可见在安全系数从初值1.427开始经过21676次迭代趋于极小值1.364
在向最终值过渡的过程中有时F值反而增加但最终顺利地收敛到整体极值图 4,23 采用模拟退火法分析[例4.9]安全系数收敛到整体极值过程
4,1,2 遗传算法肖专文等(1998)和Goh (1999)分别用遗传算法计算圆弧和任意形状临界滑裂面这里只对基本遗传算法作一概括的介绍有关遗传算法更为详细的内容请参见有关专著遗传算法基于自然界优胜劣汰的生物进化机制是一种自适应搜索技术在特定的竟争环境中依据这种规则适应性强的个体生存下来的可能性大Holland (1975)在其著作中介绍了GA的这一基本原理近来遗传算法作为一种很有应用潜力的优化方法倍受关注原因是它能有效地搜索一个复杂问题的整个定义域全局搜索遗传算法是一个迭代操作过程该迭代操作针对多个待求问题的候选解进行的这些候选解目标函数组成一个固定规模的群体候选解是模拟生物体的染色体将待求问题编码而形成的最初候选解的群体是随机生成的每一个染色体代表给定优化问题的一个可能的解每一个候选解由一串基因组成通常表示为一由0和1组成的二进制串组成染色体的每一个基因二进制子串代表一个待优化的参数整个二进制串称为染色体或个体由代表各个变量的基因子串集合而成使用目标函数可计算一个染色体对应的目标函数值安全系数进而可以
116 土质边坡稳定分析?原理
方法
程序确定与每一个染色体相对应的适应值安全系数的函数染色体通过迭代而进化在每一个迭代步骤中由多个染色体组成的群体称为一代下一代群体中的每一个染色体是由上一代群体中的两个染色体相互结合交叉操作形成也可以通过直接改变上一代群体中的某个染色体变异操作而形成从父代和子代中选择某些适应值大的染色体而淘汰其它适应值小的染色体选择操作可以形成新一代的群体而且通过这种优胜劣汰机制每一代群体中染色体个体的个数保持不变染色体所对应的适应值决定了其生存与产生后代的能力适应值最大最小安全系数最小的染色体最有可能被选择并用于产生下一代染色体这一迭代过程延缓下去直到寻找到最优解或者达到某一个迭代终止的标准图4.24所示的流程图给出了遗传算法的主要操作步骤为了实现遗传算法的一个优化循环需要选择交叉和变异三个基本操作下面对这些基本操作作一简单的介绍
1,遗传算法的基本步骤
(1) 选择操作选择操作可以确定选用当代群体中哪些个体生成下一代群体中新的个体选择操作的基本思想是将适应值超出平均水平的染色体以概率的方式挑选出来比如在轮盘选择策略(Holland,1975)中一个父代个体被选择的概率与其对应的适应值成正比
(2) 交叉操作交叉操作用于确定父代的染色体怎样结合形成子代的染色体交叉操作针对于两个经选择操作选定的父代染色体这两个父代染色体根据给定的交叉概率相互交换基因值二进制数单点交叉是最基本的交叉操作在这一操作中随机选取二进制串上的交叉点然后两个染色体在该点相互交换基因值在双点交叉中随机选取两个交叉点然后两个染色体相互交换该两点间的染色体部分图4.25 (a)和(b)分别给出了单点交叉操作和双点交叉操作的例子图 4,24 遗传算法的流程图图 4,25 交叉操作
(a) 单点交叉(b) 双点交叉第4章 确定最小安全系数的最优化方法 117
(3) 变异操作变异操作是按给定的概率随机改变某个染色体中一个二进制位上的数值例如在图4.26中由0变为1
变异操作的主要目的是提供新的基因以防止迭代落入一个局部峰值区域
(4) 参数初始化群体规模交叉概率变异概率的选取依不同的问题而异通常需要做一些试验确定图 4,26 变异操作
2,与传统方法的比较遗传算法的实现是对一个群体的多方位操作该群体是由许多潜在解组成的而不是仅对单一的解进行操作这一并行性质决定了遗传算法全局优化的特点遗传算法使用适应值而不使用适应函数的梯度或更高次的微分及其它信息所以适用于解决各种不同的优化问题具有鲁棒性强的特点遗传算法属于非确定性算法但是它不同于纯随机的算法遗传算法将随机搜索与方向性搜索结合在一起因此可以平衡优化搜索的两个相互矛盾的策略搜寻最优解与全局搜索爬山法是搜寻最优解但忽视全局搜索的典型例子而随机搜索是全局搜索的典型方法但它忽视了可能存在最优解的局部空间正如其它优化算法一样遗传算法也有其局限性主要缺点是要化费大量的计算时间但随着计算机计算速度的迅速提高这一缺点将被克服
3,遗传算法的实现
Goh (1999)的算例使用楔体分析的程序GWEDGEM该程序使用多楔体方法确定边坡的安全系数GWEDGEM的当时版本将楔体的数量限制在6块以内该程序引入遗传算法搜索临界滑裂面由于遗传算法通常用于搜索最大值故将目标函数转化为适应函数
2
min
max
)(
1
F
f = (4.70)
使用二进制编码表示参数变量二进制串的长度取决于所需的精度例如对于一个四位的二进制编码其精度大约为搜索区间的1/2
4
(1/16)当然根据所需的精度不同的变量可以使用不同长度的二进制串表示在GAWEDGE中每一个变量采用9位长度的二进制串表示这样对于一个含6个楔体的问题共17个变量染色体的总长度为
153位
[例4.15] ACAD考题3(a)
Goh在其论文中介绍的3个考题恰为将在第11章第11.4.2节中ACAD程序中的Exl(e)
Ex3(a)和Ex4现择Ex3(a)为例介绍遗传算法与其它方法计算成果的比较情况如图4.27
所示
118 土质边坡稳定分析?原理
方法
程序图 4,27 ACAD考题3(a)计算示意图在使用遗传算法进行分析时选用了以下参数进化次数为300群体大小为5交叉概率为0.75变异概率为0.02图4.28示其进化过程在第290次时计算收敛计算次数为290×5=1450其中5为自由度数目获得的最小安全系数为1.286 Goh通过表4.4 比较了多种优化方法的成果其中涉及本书作者(Chen & Shao,1988)的成果系作者为ACAD
所作计算成果之误参见第11章表11.15
图 4,28 进化次数与安全系数的关系表 4,4 对[例4,15]的计算结果比较方法楔体数目安全系数迭代次数
GWEDGEM 6 1.288 1820
GWEDGEM 7 1.288 1740
GWEDGEM 8 1.483 2490
GWEDGEM 9 2.099 3750
GAWEDGE (群体大小为5) 6 1.286 290×5=1450
SLOPE (Fredlund & Krahn,1977) 1.364-1.378
Spencer (1973) 1.24
Chen & Shao (1988) 1.242
Donald & Giam (1989) 1.27
第4章 确定最小安全系数的最优化方法 119
4,8 本章附录
4,8,1 本章数据文件表4.5
表 4,5 本 章 数 据 文 件有关章节系列号数据文件名内容
04-01-01 CPT4-3.DAT
[例4.1] 光滑滑裂面
4.2.1
04-01-02 CTP4-4.DAT
[例4.1] 折线滑裂面
04-02-01 EX41.DAT [例4.2] 两个自由度的算例用单形法
4.3
04-02-02 EX42.DAT [例4.3] 使用单形法计算图4.6例
4.5
04-03-01 CPT4-5.DAT [例4.10] 随机搜索
4.6
04-04-01 TIAN.DAT [例4.13] 天生桥滑坡例
4.7 04-05-01 CPT4-6.DAT [例4.14] 模拟退火参考文献
1 Boutrop,A,W,and Lovell,C,W,Search technique in slope stability analysis,Engineering Geology,1980,16:
51-61
2 Chen,Z,and Shao,C,Evaluation of minunmum factor of safety in slope stability analysis,Canadian
Geotechnical Journal,1988,25(4),735-748
3 陈祖煜,邵长明,最优化方法在确定边坡最小安全系数方面的应用,岩土工程学报,1988,10(4)
4 陈祖煜,邵长明,用最优化方法搜索最小安全系数的一个改进,岩石力学在工程中的应用,第二次全国石力学与工程学术会议论文集知识出版社,1989,71-75
5 陈祖煜,用DFP法计算边坡稳定最小安全系数时导数的计算方法.岩土力学数值方法的工程应用,第三届全国岩石力学数值计算与模型实验学术研讨会论文集,上海同济大学出版社,1990,797-803
6 Chen,Z,Experience with the search of minimum factors of safety of slopes,Proceedings,6th Australian New
Zealand Conference on Geomechanics,1992,426-431
7 Chen,Z Random trials used in determining global minimum factors of safety of slopes,Canadian
Geotechnical Journal,1992,29(2),225-233
8 Celestino,Y,B,and Duncan,J,M,Simplified search for noncircular slip surface,10
th
Int,Conf,On Soil
Mech,And Found,Engr.1981,Vol,391-394
9 Li,K,S,and White,W,Rapid evaluation of of the critical slip surface in slope stability problems.
International Journal for Numerical and Analytical Methods in Geomechanics,1987,11,449-473
10 Fox,R,Optimisation methods in engineering designs,Addison-Wesley Publication Company,1971
11 华东水利学院,土工原理与计算,北京水利出版社,1979
12 潘家铮,建筑物的抗滑稳定和滑坡分析,北京水利出版社,1980
13 孙君实,条分法的数值分析岩土工程学报,1984,6(2),1-12
14 Shoup,T,E,and Mistroup,F,Optimisation methods with alications for personal computers,Prentice-Hall Inc.
Eaglewood Cliffs,1987
15 Siegel,R,A.,Kovacs,W,D,and Lovell,C,W,Random surface generation in stability analysis,ASCE
Journal of Geotechnical Engineering,1981,107,996-1002
16 Stewart,G,W,III,A modification of Davidon’s minimization method to accept difference approximations of
derivatives,Journal of Association for Computing Machinery,1966,1,72-83
17 唐立山谢云尤矢勇罗祖华,非数值并行算法,第一册,模拟退火算法,科学出版社,1994
120 土质边坡稳定分析?原理
方法
程序
18 肖专文,张奇志,梁力,林韵梅,遗传进化算法在边坡稳定性分析中的应用,岩土工程学报 1998,1:44-
46
19 闫中华,均质土坝与非均质土坝稳定安全系数极值分布规律及电算程序简介,水利水电技术,1983年第
7期
20 张天宝,四边形块方法和沉抗土坝的稳定分析兼论复合土坡的最危险滑弧位置,水利学报,1981
21 邹广电,边坡稳定分析条分法的一个全局优化算法,岩土工程学报,2002,24(3),309-312
22 周文通,最优化方法在土坝稳定分析中的应用,土石坝工程,1984
23 Nguyen,V,U,Determination of critical slip surface,Journal of Geotechnical Engineering,ASCE,1985,111:
238-251
24 Goh,A.T.C,Genetic algorithm search for critical slip surface in multi-wedge stability analysis,Candian
Geotechnical Journal,36,(2)383-391,1999.
25 Greco,Y.R,Efficient Monte-Carlo technique for locating critical slip surface,ASCE J,Geotechnical
Engineerinf,122,(7)517-525,1997.
26 朱伯芳黎展盾张壁诚,结构优化设计原理与应用,北京水利电力出出版社,1984
4,1 概述
4,1,1 应用最优化方法确定最小安全系数边坡稳定分析应包含下面两个步骤
1) 对滑坡体内某一滑裂面按第2或第3章介绍的方法确定其抗滑稳定安全系数
2) 在所有可能的滑裂面中重复上述步骤找出相应最小安全系数的临界滑裂面在讨论了计算一个滑裂面的安全系数的方法后本章介绍边坡稳定分析极限平衡法的第二步就是要寻找最小安全系数的方法如果滑裂面曲线为y(x)那么这个问题具体化为寻找下列泛函的极值
)(yFF = (4.1)
岩土工程中的边坡的几何形状各异材料通常是非均质性纯解析的变分原理很难进行极值计算用最优化方法通过数值方法求解是一个比较现实可行的途径我国学者较早开展应用数值规划方法求解安全系数极值的问题闫中华,1983张天宝,
1981周文通,1984 20世纪80年代初期孙君实(1983)和Nguyen (1985)分别提出了使用复形法和单形法搜索任意形状和圆弧滑裂面的最小安全系数的方法80年代中期有更多的学者发表了有关研究工作(Celestino and Duncan,1981; Li,and White,1987) Chen &
Shao(1988)采用单形法和牛顿法进行任意形状滑裂面搜索的研究成果这篇论文所附的几道例题后来在国外多篇论文中被引用作为检验新的优化方法的考题借此书的机会作者将这些例题和原数据文件在本章介绍国内外有关的研究成果和STAB软件十几年在全国推广的实践证明应用计算机自动搜索临界滑裂面是可行的同时作者也发现相对于三维或二维斜分条的极限平衡法垂直条分法的极小值搜索问题比较简单采用任何一种优化计算方法配合本章介绍的随机搜索法即可快速地找到临界滑裂面因此目前在STAB程序中向用户提供的只是单形法这个功能以后作者还将把本章中介绍的各种方法补充进程序供用户选用
4,1,2 最优化方法最优化方法是近代数学规划中十分活跃的一个领域目前已有许多十分成熟的计算方法这些计算方法总的来看可以分为以下三大类
1,枚举法枚举法的基本思想是根据一定的模式比较不同自变量的目标函数经过筛选最
88 土质边坡稳定分析?原理
方法
程序终找到最小值这是最原始简单的方法如图4.1中示任一圆弧可用其圆心坐标(x
0
,y
0
)和半径r确定其相应的安全系数F可表达为
),,(
00 s
DyxfF = (4.2)
式中D
s
为滑弧深度即圆弧最低点的坐标可知
0
yrD
s
+= (4.3)
显然这是三个自由度的问题采用枚举法不断地改变x
0
,y
0
和D
s
的数值逐一比较相应的安全系数最终找到最小的安全系数在具体操作中先固定一个D
s
然后在圆心可能的位置中布置一个网格见第12章图12.1设网格的中心坐标为x
c
和y
c
在左右方向各布置了N
x
格在上下方向各布置了N
y
格则共计有(2N
x
+1)×(2N
y
+1)个网格点分别以该网格点为圆心以D
s
为滑弧深度计算相应安全系数找出最小的安全系数然后改变一个D
s
值重复相同的步骤在这一过程中有可能出现圆弧和边坡不相交的情况则应令程序自动抛弃该圆弧同样D
s
也是以一个中心值起算在其上下各布置N
d
层这样总计计算(2N
x
+1)×(2N
y
+1)×(2N
d
+1)个圆弧在枚举法的基础上用0.618法或其他方法提高搜索最小安全系数工作的效率这方面的工作可参阅有关文献图 4,1 圆弧滑裂面
2,数值分析方法随着计算机的发展数值分析方法逐步形成一门完整的学科统称为最优化方法(Method
of optimization)这一方法又可分为两大类第一类称模式搜索法(Pattern search method)其基本思想是根据一定的模式比较不同自变量的目标函数确定最优的搜索方向最终找到最小值另一类称牛顿法它要通过解析手段寻找使目标函数F对自变量z
i
的偏导数为零的极值点(?F/?z
i
=
0,i
=
1,2,...,n)同时从理论上讲还需要满足由二阶导数形成的Hessian矩阵正定这个达到极小值的充分条件此类方法中以导数为研究的主要对象因此也称为以导数为基础的方法(Gradient?based?method)一般认为当自由度较多时直接搜索法效率较低此时需要考虑牛顿法体系的分析方法第4章 确定最小安全系数的最优化方法 89
本节介绍直接搜索法中的单纯形法和Powell法以及牛顿法中负梯度法和DFP法
(Davidon?Fletcher?Powell method)由于这些方法的原理在众多的文献及教科书中都有所介绍这里不拟全面阐述其原理而着重讨论在边坡稳定分析领域需解决的特殊问题
3,非数值分析方法第三类方法是在近期计算机发展基础上形成的称为非数值方法这一方法被广泛应用于管理科学计算机科学分子物理学及超大规模集成电路设计中用于解决组合优化问题非数值分析利用计算机具有容量大计算速度快的优点通过大量随机采样来找到目标函数的最优值近期涌现了诸如模拟退火遗传算法神经网络和蚂蚁算法等均属此类在边坡稳定分析中上述算法均有人尝试本章4.7 节将对此类方法作一简要介绍
4,2 任意形状滑裂面的模拟和目标函数的确立
4,2,1 任意形状滑裂面的模拟最优化问题的提法是对于一个具有n个自变量的向量 z
=
(z
1
,z
2
,...,z
n
)确定其目标函数F的最小安全系数F
m
相应的自变量为z
m
为此需要对式(4.1)中的曲线y(x)用若干参数来模拟也就是说需要将任意形状滑裂面y(x)用z来近似表达将滑裂面曲线用m个点A
1
,A
2
,…,A
m
离散见图4.2也就是将此m个点用直线或光滑的曲线连起来以近似模拟此曲线此m个点的坐标用z
i
(i
=
1,2,,m)
表示
=
i
i
i
y
x
z (4.4)
图 4,2 任意形状滑裂面一旦这种连接的模式确定安全系数F即可表达成此m个点的坐标x
1
,y
1
,x
2
,y
2,
...,x
n
,y
n
的函数
()
mm
yxyxyxFF,,,,,
2211
L= (4.5)
90 土质边坡稳定分析?原理
方法
程序在进行最优化搜索过程中A
1
,A
2
,…,A
m
将移到临界滑裂面的位置
m21
B,...,B,B ′′见图
4.2此处m
=
6其中端点A
1
,A
m
原来在边坡线上有可能移到边坡线外或内如图4.2中的
m1
B,B ′′为此需要通过一定的处理方式分别找到它们和边坡线的交点B
1
和B
m
仍以
m21
B,...,B,B作为研究的对象对均匀的土质边坡通常希望滑裂面比较光滑此时可采用三次或更高次的样条函数连接这些点第11章将详细介绍使用样条函数构筑光滑曲线滑裂面的方法当然也可以采用直线和光滑曲线的组合构筑滑裂面例如图4.2,A
3
,A
4
,A
5
,
A
6
用曲线相联A
2
,A
3
用直线相联应用样条函数构筑光滑滑裂面对于减少自由度提高数值计算效率具有重要意义
[例4.1] 折线和光滑滑裂面比较如图4.3所示例如果用四个点模拟滑裂面那么折线滑裂面和曲线滑裂面相应的最小安全系数分别为1.489和1.364差别颇大换句话说如果使用折线模式要得到与曲线模式相同的精度就要更多的离散点这意味着更多的自由度而数值分析的机时和收敛难度是随自由度按指数增加的值得注意的是在笔者所知的所有这方面的文献中还没有其他研究者使用过曲线模拟任意形状滑裂面在稳定分析中对土质比较均匀的边坡常采用圆弧滑裂面此时自变量可以是圆弧的圆心的x,y坐标和半径r由于自变量仅三个最优化方法搜索最小安全系数收敛性能很好因此将不作为重点在本章中讨论图 4,3 说明用不同模式模拟滑裂面的差别滑裂面1?直线模式滑裂面2?曲线模式
4,1,2 目标函数的确定当滑裂面z的离散模型确定后安全系数便是m个控制点A
1
,A
2
,...,A
m
的函数在优化计算过程中这m个点中有n个点各沿某一设定方向 β
i
向临界滑裂面移动或者不规定方向任其自由移动其余m?n个点由于问题本身的要求可以固定滑裂面上任意一点的z
i
可以用相应于一个初始滑裂面的相对坐标来代表见图4.2
=
o
i
o
io
i
y
x
z (4.6)
+=?+=
i
i
i
o
it
o
ii
d
β
β
sin
cos
zzzz (4.7)
式中i
=
1,2,…,n,d
i
为第i个点沿β
i
移动的距离第4章 确定最小安全系数的最优化方法 91
对于固定的点其自由度为零如图4.2中的A
2
,A
3
当某点需沿一个软弱夹层移动时则其自由度为1如对该点的移动方向无特殊要求如图4.2中的其它点该点的自由度为
2问题的自由度为各点自由度的总和搜索最小安全系数的问题具体化为求下列函数的最小值问题
),...,,,,,,(
2121 nn
dddFF βββL= (4.8)
如前所述滑裂面的一部分可能是光滑曲线另一部分受软弱夹层的控制为直线STAB
程序是这样规定这一定义格式的设滑裂面被m个点A
1
,A
2
,……,A
m
分为m?1段从上交点向下交点编号为1,2,…,m?1在STAB程序中用LNO代表此m?1段中为直线段的线段总数在数组LOO(I),I
=
1,2,…,LNO中存入这些直线段的编号其它线段则默认为曲线当LNO为零时滑裂面为一光滑曲线当LNO为一大于m?1的数时则程序默认为全部为直线段这两种情况都无须填写LOO(I)
4,3 模式搜索法
4,3,1 单形法朱伯芳等,1984
1,建立初始单形对某一初始向量z
0
按下面模式构筑n个向量z
i
(i
=
1,2,..,n)组成n+1个顶点
T
m
n
T
m
T
m
pzqzqz
qzpzqz
qzqzpz
],,,[
],,,[
],,,[
00
2
0
1
00
2
0
1
2
00
2
0
1
1
+++=
+++=
+++=
LL
LL
LL
LL
z
z
z
(4.9)
其中
a
n
nn
p
2
1)1(?++
= (4.10)
a
n
n
q
2
1)1(?+
=
(4.11)
式中a为初始步长可根据实际情况进行调试
2,确定搜索方向计算单纯形n+1个顶点的目标函数值比较其大小从中找出目标函数最大点z
H
和次大点z
G
和最小点z
L
然后按下式确定下一步的搜索方向
∑=
+
v
H
v
i
n
i
v
n
n
zzz
1
2
1
(4.12)
式中上标v表示迭代次数
92 土质边坡稳定分析?原理
方法
程序显然从
ν
H
z到
ν
2+n
z目标函数是下降的因此以
ν
H
z和
ν
2+n
z两点的联线作为搜索方向
3,优化计算
(1) 反射沿搜索方向
ν
H
z和
ν
2+n
z点向前再走一步步长为)(
2
v
H
v
n
zz?
+
α到达
ν
3+n
z点
)(
223
v
H
v
n
v
n
v
n
zzzz?+=
+++
α (4.13)
ν
3+n
z称为反射点α>0称为反射系数计算反射点的目标函数并根据它的大小来决定下一步的走法
(2) 扩张如果反射点的函数值小于z
L
点的函数值即
)()(
v
L
FF zz
v
3n
<
+
(4.14)
则表明反射后情况有所改善沿搜索方向
v
n
v
H 2+
z,z还可以试探一下是否可以走得更远一些即是否可以扩张到
v
n 4+
z点
)(
2324
v
n
v
n
v
n
v
n
v
++++
+= zzzz (4.15)
式中v>1为扩张系数如果)()(
4
v
L
v
n
FF zz <
+
则以
v
n 4+
z点替换原来的最坏点
H
z构成新的单纯形转入第4
步进行收敛判断如果)()(
4
v
L
v
n
VV zz ≥
+
则以
v
n 3+
z点替换
v
H
z构成新的单纯形也转入第(4)步进行收敛判断
(3) 收缩如果反射点的函数值大于次坏点的函数值即
)()(
3
v
G
v
n
VV zz >
+
(4.16)
表明反射点走得太远需要缩回去按下式计算
)(
225
v
n
v
H
v
n
v
n +++
+= zzzz β (4.17)
式中0
<
β
<
1是收缩系数用
v
n 5+
z代替
v
H
z构成新的单纯形并转入第(4)步
(4) 缩小边长如果反射点函数值大于最坏点函数值即)()(
3
v
H
v
n
zVzV >
+
则缩小单纯形的边长以最好点
L
z为顶点其它各顶点向
L
z移近一半距离即按下式计算
ni
v
L
v
i
v
i
,...,2,1,0 ),(5.0 =?+= zzzz
v
L
(4.18)
得到新的单纯形转入第(4)步重复计算
4,收敛判断按照一定的方式通过反射扩充和收缩使单形不断更新逼近极值点收敛准则为
[]ε≤?∑
+
+
2
2
0
)()(
1
1
v
n
v
i
n
i
vv
n
zz (4.19)
第4章 确定最小安全系数的最优化方法 93
式中ε为要求的计算精度如式(4.19)满足则结束计算并以
v
L
z作为极小点否则置v
=
v+1转至第3步重复计算据已有经验可取α=1; 0.4
≤
β
≤
0.6; 2.0
≤
γ
≤
3.0经Nelder和Mead论证为使单纯形适应函数的性态及便于收敛α不宜比1大很多而α
<
1的计算次数比α=1要多故折衷取α
=
1另外认为β的数值变化对搜索效率的影响比γ要大他们推荐采用
6.04.0 ≤≤β (4.20)
0.38.2 ≤≤γ (4.21)
[例4.2] 说明单形法计算过程例为了形象直观地了解单形法在搜索最小安全系数时的工作状况我们考察图4.4所示的一个有两个自由度的例子滑裂面它由ABC组成计算时令C点固定不动A,B两点沿水平线移动则该滑裂面的安全系数F由A点的x坐标x
1
和B点的x坐标x
2
决定图4.5
示F相应x
1
,x
2
的等值线图根据枚举法可以发现在x
1
=
92.0,x
2
=
143.0时安全系数获得最小值1.257相应临界滑裂面如图4.4中标5的那个滑裂面图 4,4 具有两个自由度的算例
1,2,3,4?初始滑裂面5?临界滑裂面如果使用单形法则按式(4.9)初始生成的单形即三个滑裂面如图4.5在左下角三角形所示最大次大和最小分别为A B C三点D为按式(4.12)所得的新的顶点AD代表了两点联线的方向D代表了第一次反射和扩张后达到的点第一次迭代后形成了B
C D构成的新的单形开始新的一轮迭代依次循环最终达到安全系数的极值E搜索过程如图4.5中折线1所示最终收敛到F
m
=
1.257相应的z
m
=
(92.00,143.00)
T
94 土质边坡稳定分析?原理
方法
程序图 4,5 对图4.4示例用单形法计算过程图 4,6 具有软弱夹层的算例
[例4.3] 使用单形法计算最小安全系数图4.6具有分层剖面和一个倾斜的软弱夹层图中示4个不同的初始滑裂面点A和
B可以不受限制地自由移动故自由度为2同时此两点有意未设在边坡面上以说明程序在搜索过程中将进行判断和调整自动找到滑面与坡面的交点C和B点位于夹层上因此自由度为l该两点在搜索过程中将沿一定方向β移动β为该夹层的倾角图4.7中滑裂面3为采用单纯形法获得的临界滑裂面相应安全系数为1.025
第4章 确定最小安全系数的最优化方法 95
图 4,7 对图4.6示例用最优化方法计算过程临界滑裂面1?改良后的DFP法F
m
=1.025 2?负梯度法F
m
=1.025
3?单形法F
m
=1.025 4?DFP法F
m
=1.025
4,1,2 Powell法
1,基本原理有关此法的详细原理可参阅文献朱伯芳等1984其基本算法如下
(1) 给定初始点z
0
并令
Ti
)0,...,1,...,0,0(=S (4.22)
式(4.22)右端的1处于第i位i
=
1,2,…,n
(2) 相应迭代步i
=
1进行一次一维搜索第一次搜索方向p
i
=
S
i
使极小=+
)(
1 iii
V pz α (4.23)
并令
iiii
pzz α+=
1
(4.24)
然后令i
=
i+1重复上述过程直至i
=
n为止
(3) 取p
i
=
p
i+1
,i
=
1,2,…,n?1并令
0
zzp?=
nn
(4.25)
通过这一步增加
0
zz?
n
丢弃原来的坐标方向
1
p
(4) 沿z
n
z
0
作一次一维搜索求α
n
使V(z
n
+α
n
p
n
)
=极小并令
nnn
pzz α+=
0
(4.26)
然后转向第(2)步
(5) 重复进行上述运算n次即沿n个互相共轭的方向都进行了一维搜索然后作如下判断
εα <
nn
p (4.27)
96 土质边坡稳定分析?原理
方法
程序如上述判据满足则结束计算此时z
0
即为所求极小点否则转向第(1)步
2,应用Powell法时所做的改进在采用Powell法进行计算时第一步需要确定n个独立的搜索方向在大多数情况下都采用与n个自变量z
i
坐标轴方向一致的单位向量即由式(4.22)定义的S
i
这种处理方式并不适于边坡稳定分析参见图4.8(a)所示一个有两个自由度例子按图 4,8 在边坡稳定分析中使用Powell法作的改进
(a) 未作改进(b) 改进后
Powell法第一步固定B点不同让A
点沿水平方向移动到A′使安全系数达到极小值显然能够移动的范围是很有限的第二步固定A′不动B在其移动方向到达B′再次得一个极小值其移动范围更为有限这样的搜索效率很低当自由度较多时事实上无法实现有效的搜索设想如果在A点移动时让B点也向前移动一定值第一次一维搜索就可以使初值z
0
向极值z
m
迈进一大步如图4.8 (b)
因此在构筑S
i
时我们的思路是当某一控制点在移动方向搜索最小值时要求其邻近的点也按一定的比例移动使滑裂面整体而不是局部在优化计算中不断地调整以迅速地达到极值点为此按下式构筑n个向量( ),...,3,2,1 ni
i
=S
Τ
= ),...,,(
21
i
n
iii
eeeS (4.28)
按下式确定( ),...,3,2,1 nie
i
j
=
>
≤
=
时当时当
ji
xx
xx
ji
xx
xx
e
in
jn
ij
i
j
ij
i
j
)cos(
)cos(
1
1
ββ
ββ
(4.29)
当i
=
1时总是用(4.29)的下式i
=
n时总是用式(4.29)的上式
i
j
e的物理意义是当第i个变量沿其搜索方向移动一个单位值时邻近的点j都将在其移动方向获得一个增量
i
j
e这个增量值取决于第j点距端点与j同侧的那个的水平距离在i
=
0和i
=
n时等于零i
=
j时等于1其余点按线性分布原则确定内插参见图4.9
第4章 确定最小安全系数的最优化方法 97
图 4,9 Powell法搜索方向
i
j
e的改进
(a) 1<i<n情况(b) i=1情况(c) i=n情况
[例4.4] 说明对Powell法所做的改进例本例和[例4.2]具有相同的计算条件采用折线滑裂面参见图4.4如果初始滑裂面设为1图4.5中线2和3分别代表未经改进和改进后的计算过程可见改进后搜索效率大大提高采用这种改进的方法要求各控制点都按同一方向朝临界滑裂面逼近例如图4.2所示例可以将初始滑裂面布置在离估计的临界滑裂面稍远的地方但A,B,C点都采用水平向左的移动方向即式(4.29)中的β均为零通常可以很快地找到极值
4,4 牛顿法
4,4,1 负梯度法负梯度法的基本思想是对一个初始滑裂面寻找一个使安全系数减少速率最大的方向数学上就是
T
n
z
F
z
F
z
F
,,,
21
L这个向量在这个方向上进行一维搜索找到这一方向安全系数的低谷点完成了这第一次迭代后再在这个新的起点即上述低谷区重复这样的运算直到收敛至极值点
[例4.5] 对[例4.2]使用负梯度法计算最小安全系数在使用负梯度法时分别以(84.0,160.0),(70.0,145.0),(112.0,150.0)作为起点相应滑裂面如图4.4中1,2,3示使用负梯度法的搜索路径如图4.10中1,2,3三条折线所示可见每一次搜索均是沿着等值线的法线方向即下降速率最大的方向进行的对这样一个比较简单的算例负梯度法相应不同初始输入的滑裂面都可以成功地找到极值
98 土质边坡稳定分析?原理
方法
程序图 4,10 对[例4.2]使用负梯度法和DFP法计算最小安全系数
4,1,2 DFP法
1,基本原理
DFP法为Davidon?Flefchen?Powell法的简称其原理和详细分析方法参见文献朱伯芳等1984
数学上的极值点是目标函数对各自变量的一次导数为零的点也就是说在极值点
),,,(
21 n
z
F
z
F
z
F
= LG均为零同样要求由二阶导数构成的海色(Hessian)矩阵正定海色矩阵由下式定义
=
2
2
2
2
2
2
2
1
2
21
2
2
1
2
.,,
.,,
n
n
n
z
F
zz
F
z
F
zz
F
zz
F
z
F
称对MO
H (4.30)
第4章 确定最小安全系数的最优化方法 99
DFP法试图通过一系列迭代步骤使假定的海色矩阵的逆阵A
0
逐渐过渡到其真值每一次迭代的方向均按下式决定
vvv
GAS?= (4.31)
在迭代步 v +1通过下式确定A
ν+1
vvvv
DCAA?+=
+1
(4.32)
其中 v
=
0,1,2,…
v
T
v
T
vvv
Yz
zz
C
= (4.33)
v
vT
v
T
v
v
v
v
v
YAY
YAYA
D
)(
= (4.34)
vv
v
GGY?=
+1
(4.35)
vv
v
zzz?=?
+1
(4.36)
DFP法即是基于这一基本思路的一种效率较高的方法DFP法的基本思路是通过迭代使A
0
,A
1
,…,A
n
逼近极值点处Hessian阵的逆阵A通常取A
0
为单位矩阵
2,对DFP法的改进本书作者在20世纪80年代开始研究使用DFP法计算安全系数的极值问题发现在使用这一方法时还需要针对边坡问题的特点解决量纲问题和导数的准确计算问题下面分述此两个问题
(1) 关于量纲问题如前所述DFP法的基本思路是通过迭代使A
0
,A
1
,…,A
n
逼近极值点处Hessian阵的逆阵在边坡稳定分析中?
2
F/?z
i
z
j
的量纲通常在10
2
～10
5
的范围内H
1
的各元素数值的量纲大约在10
2
～10
5
范围内常规的DFP法采用A
0
为单位矩阵那么要从A
0
过渡到10
2
~10
5
的量级就变得十分困难为此改进的方法取
p
101
0
×=A (4.37)
式中p为对H
1
量纲的平均值的一个估值可按下式确定
1
1
122
10
|/|
10
+
=
<
<
∑
p
n
i
ip
n
zF
(4.38)
[例4.6] 对[例4.2]使用改进后的DFP法计算最小安全系数同样我们通过[例4.2]计算以图4.4中的第4个滑裂面为初始滑裂面如果不采用本节引入的改进则根据式(4.32)第一次迭代获得搜索方向为
T1
)7.1,3.19(=S以图4.10中
BC
1
代表沿着该方向显然无法到达极值点D如果看一下计算这个方向的细节
100 土质边坡稳定分析?原理
方法
程序
=
+
=
+=
24.3 217.7-
217.7- 2303.8
0.408 0.492-
0.492- 5917.0
23.7 218.2-
218.2- 4.2003
1 0
0 1
0001
DCAA
可以看到(A
0
D
0
)与C不是同一量级的量(A
0
D
0
)对A
1
的贡献非常小如果使用本节提出的改进那么根据式(4.38)可得到p
=
3,A
1
计算细节为
=
+
=
+=
615.7 273.8
273.8 2411.7
408.0 492.0-
492.0- 5917
23.7 218.2-
218.2- 4.2003
1000.0 0.0
0.0 1000.0
0001
DCAA
可以看到使用了改进的方法后(A
0
D
0
)和C
0
具有相同的量级据此算得的
,)45.3,07.3(
1 T
=S为图4.10中BC所示方向这就保证了第二次迭代快速地向极值点逼进并顺利地达到极值点
(2) 关于导数的计算问题牛顿法体系的最优化方法均是建立在寻找目标函数的导数的基础上的准确地确定?F/?z
i
的数值是个十分重要的问题我们遇到的困难是
(a)?F/?z
i
不能通过一个公式直接算得因此只能通过差分法计算?F/?z
i
的近似值?F/?z
i
(b) 由于F本身只能通过迭代求解获得因此F的有效位很少在差分计算时步长?z
i
也就不能取得很小这就使算得的导数值的误差较大另一方向当趋近极值点时导数值本身却趋近于零这两个因素迭加在一起就使计算所处得的导数在极值点附近相对误差急剧增加导致数值计算无法顺利地收敛到极值点大部分数值分析问题都不可能通过解析方法算得导数值在使用牛顿法进行非线性分析时解决在差分运算中导数计算精度问题是一个具有普遍意义的问题因此本节结合边坡稳定分析的特点讨论使用中心差分时步长的合理选定问题同时提出了一个计算导数的半解析法使用这个方法可以获得较精确的导数值
(c) 用差分法确定导数的步骤首先简单回顾Stewart (1967)讨论的有关差分计算中步长确定的原则对于某一函数)(zφφ =用向前差分计算dφ / dz在z
=
0处值的公式是
δ
φδφφφ )0()(
d
d?
=
≈
zz
(4.39)
中心差分计算公式是第4章 确定最小安全系数的最优化方法 101
δ
δφδφφ
2
)()(
d
d
≈
z
(4.40)
式中δ为选用的步长使用式(4.39)计算dφ/dz其误差来自两个方面
(a) 截断误差?φ的准确值为
...
2
1
=
)0()(
2
++
=?
ρδγδ
φδφφ
(4.41)
式中γ,ρ分别为?φ/?z,?
2
φ/?z
2
在z
=
0处的值因此用式(4.39)计算γ实际上忽略了ρδ
2
/2
由此引起的相对截断误差η
t
近似为
γ
ρδ
η
2
2
1
=
t
(4.42)
(b) 舍入误差任何一个数值在计算机中都是一个有限字长的数字近似代表的如果δ
选取得太小计算机将无法识别φ(δ)和φ(0)之间的区别由于φ的有限精度造成的按式(4.39)
计算的γ的误差称舍入误差舍入误差η
c
按下式估算
η
φ
φ
η
≤
)0(
2
c
(4.43)
式中η为φ(0)的相对误差如果φ(z)的数值是通过某一公式对自变量z进行若干算术运算得到的那么η值首先取决于z本身的相对误差在计算机中z是被一个在区间[ )1(),1(
mm
zz εε +? ]的某个数值近似代表的其中ε
m
称为机器误差它的大小取决于计算机的种类和精度级别(单双精度)
可以按图4.11所示框图编一个简单的程序来测试ε
m
的数值参见Dennis and Shanabl,1983
年对于16位的计算机可以知道ε
m
=
5.96×10
8
图 4,11 确定计算机精度的计算框图
η值还取决于在进行φ(z)的算术运算时积累的误差它使φ(z)的相对误差η比z的相对误差ε
m
大几倍或几十倍如果φ(z)的数值不能用显式计算需要通过迭代求得则η取决于迭代精度例如稳定分析的安全系数F它是通过迭代求解式(2.31)和式(2.32)获得的在单精度运算时迭代收敛标准是10
5
由于安全系数的数值在1左右故可大致认定F的相对误差
η<10
4
Steward提出的确定步长δ的步骤如下
102 土质边坡稳定分析?原理
方法
程序
(a) 要求式(4.42)的和式(4.43)右端相等到即要求舍入误差和截断误差平衡据此计算
δ值
(b) 用式(4.39)计算γ的值并按式(4.42)估算η
t
如果η
t
超过了允许值ε(例如ε = 5%)
则要用中心差分即使用中心差分式(4.40)计算γ
(c) 如果使用中心差分则δ通过解下式求得
0)0(2
2
1
2
=?+φηδγεεδp (4.44)
图4.12示dφ/dz数值其中tanα
0
为准确值tanα
1
为中心差分计算值tanα
2
为向前差分计算值tanα
3
为向后差分计算值式(4.44)是令式(4.43)右端等于ε获得的其中?φ按式(4.41)确定在稳定分析中由于F的精度低故总是用中心差分因此我们不进行
Steward提出的第一步计算在用式(4.44)
计算δ时取ε
=
5%,η
=
10
4
ρ,γ的值近似采用在DFP法计算时上一次迭代已经算得的值假定φ在(?δ,δ)区间呈抛物线分布这样γ值仍按式(4.39)计算ρ的计算公式为图 4,12 向前差分可能获得相反符号的导数示意图
2
2
)0(2)()(
δ
φδφδφ
ρ
+
= (4.45)
采用中心差分的另一个理由是由于F的精度低δ的数值相对较大当快要接近极值时向前差分可能得到与真实值符号截然相反的导数值见图4.12中心差分可在很大程度上避开这一危险
3,半解析法计算导数的步骤通用条分法的特点是在计算F值时各项导数不通过差分而是通过解析公式获得的利用这一特点我们可以将计算?F/?z
i
所包含的差分运算减少到最低量设想由于z
i
的一个增长量?z
i
使式(2.31)和(2.32)中的F和λ得各有一增量?F和?λ则
()0,,=?+?+?+=
iinn
zzFFGG λλ (4.46)
0),,( =?+?+?+=
iinn
zzFFMM λλ (4.47)
全微分式(4.46)和式(4.47)可得第4章 确定最小安全系数的最优化方法 103
i
n
i
n
i
n
z
G
z
G
z
F
F
G
=
+
λ
λ
(4.48)
i
n
ii
n
zz
M
z
F
F
M
Μ?
=
+
λ
λ
(4.49)
由此可得
F
ΜGM
F
G
z
GΜ
z
ΜG
z
F
nnnn
i
nn
i
nn
i
=
λλ
λλ
(4.50)
式(4.50)给出了一个不直接计算F值确定?F/?z
i
的公式式(4.50)右边的?M
n
/?λ,?M
n
/?F,
G
n
/?λ,?G
n
/?F都是在计算一个滑裂面安全系数时已经算得的而且是通过精确的计算公式确定的即第2章中式(2.37)至式(2.40)?G
n
/?z
i
,?M
n
/?z
i
仍然需要通过差分求得但是由于计算G
n
和M
n
的公式(2.31)和式(2.32)不需迭代故其导数的精度远比对F差分精度高用中心差分计算?G
n
/?z
i
和?M
n
/?z
i
时δ的选用原则是要求按式(4.43)确定的舍入误差相对值η
c
小于允许值ε G
n
和M
n
是通过式积分算得的在两个相邻的控制点1,+ii之间被分成了
1+i
N个土条因此要求在进行差分运算时每个土条高度h的相对误差都小于规定值当第i个控制点在其移动方向z有一增量δ 时其(i,i+1),(i?1,i)两段的每个土条的高度
h都将获得一个增量h?其中h获得增量最小的那个土条相应值为δ
/N
m
这里
),max(
1+
=
iim
NNN,
1
,
+ii
NN分别为(i?1,i)和(i,i+1)段的土条数分别以η
c
=
ε,?φ
=
δ
/
N
m
,φ
=
z
m
代入式(4.43)可得一确定δ的近似公式
η
ε
δ
mm
Nz
2≥ (4.51)
式中η为z
m
的相对误差可取η
=
ε
m
选用5%~10% z
m
为x
i
,y
i
中的大值下面列举一个例题以验证本节建议的计算导数的方法的可行性
[例4.7] 测试计算导数精度例题图4.13示一垂直均匀边坡假设滑裂面为通过坡趾A与水平线倾角为α的直线那么沿该滑裂面滑动的安全系数F仅为B点的x坐标值的函数并且安全系数F和dF/dx均可通过理论推导获得:
αγ
αφ
2sin
4
cottan
h
c
F += (4.52)
x
h
c
dx
dF
/)2cot4cot(tan
γ
ααφ+= (4.53)
104 土质边坡稳定分析?原理
方法
程序图 4,13 测试计算导数精度例题当F获得极值时
s/21tan +=α (4.54)
φγ tan
4
h
c
s = (4.55)
对图4.13示边坡取h
=
25m,φ
=
35
°
,c
=
49kPa,γ
=
1.8g/cm
3
分别采用差分和半解析法计算dF/dx并与理论值对照列于表4.1
此例当x
=
12.283时获得安全系数极小值F
m
=
0.9056表4.1所列计算结果表明对于本例半解析解和差分法均可以获得和理论值完全一致的最小安全系数和临界滑裂面表 4,1 理论值和近似法结果对照
x 25 8 5 12.283
理论值式4.52) 1.1446 0.09896 1.2956 0.9056
F
程序计算1.1447 0.9896 1.2956 0.9055
理论值式4.53 0.02801 -0.0499 -0.1854 ≈0
δ 0.019 0.011 0.011 0.011
半解析法
dF/dx 0.02801 -0.0499 -0.1853 0.713×10
-4
δ 0.035 0.078 0.014 0.378
x
F
d
d
数值差分法
dF/dx 0.02806 -0.0499 -0.1854 0.367×10
-4
[例4.8] 对[例4.3]使用DFP法计算成果用经过改进的DFP法计算[例4.3]成果示于图4.14临界滑裂面1至4相应图4.6
中的四个不同的初始滑裂面可见基本上收敛在相同的位置F
m
在1.009~1.025之间
4,5 确定整体极值的随机搜索方法
4,5,1 概述数学中的最优化理论研究的目标函数在其自变量空间是光滑连续的同时要求在极值点附近目标函数是自变量的二次函数而岩土工程中的边坡一般为非均质的因此在实际应用时经常遇到多极值问题在使用数值分析方法进行边坡稳定的最大化分析时笔者发现几乎所有的方法都存在陷入局部极值陷井无法获得整体极值的问题第4章 确定最小安全系数的最优化方法 105
图 4,14 对[例4.3]使用负梯度法和DFP法计算成果
[例4.9] 说明确定整体极值困难一例如图4.15所示已知边坡的整体极值为1.364对于选定的那个初始滑裂面1,F
o
=
1.593
单形法和DFP法计算均终止在线2相应安全系数F
=
1.452其相应临界滑裂面和已知的实际位置相差甚远经验表明在下面的两种情况下数值计算的收敛问题变为极为严峻
(1) 当问题包括较多的自由度时可以设想此时目标函数的自变量空间存在许多局部极值点使搜索整体极值的任务变得十分困难
(2) 当边坡断面较复杂时此时也会出现许多安全系数的局部极值图 4,15 说明确定整体极值困难一例
1?初始滑裂面; 2?应用单形法获得的临界滑裂面如图4.15所示计算终止的那个滑裂面的一部分实际上和两种不同土层的交界线重合很可能是一个局部极值点这道题如果改成均质土坡同样的初始滑裂面位置则计算所得的临界滑裂面可很好地收敛到整体极值点的位置由此可见为了使最优化方法真正成为一个在边坡稳定分析中实用性强的工具还需作进一步的努力有效地解决确定整体极值问题在这方面一个十分简单又十分有效的思路是设法确定一个靠近整体极值的初值显然初值离整体极值愈近丢失整体极值的可能性越小不少教科书中都指出确定一个好的初值的重要性同时指出这个任务可以由随机搜索圆满完成(Fox,1971; Shoup,et,al,1987; Box,1965)
随机搜索的基本思想是应用随机数在研究域内构筑一系列自变量比较其相应的目标函数寻找最小的目标函数这样在进行常规的最优化计算以前对自变量空间进行一次均匀的高密度的扫瞄把整体极值的大体位置确定下来再将这个最小目标函数作为初值进行最优化计算
106 土质边坡稳定分析?原理
方法
程序根据Monte Carlo法的原理当随机搜索的次数趋于无穷大时目标函数的最小值就是整体极值因此随机搜索本身也可以独立地成为一种计算最小安全系数的方法已经有一些文献报导这方面的研究成果(Boutrup,1980; Siegel,1981)在笔者发表了有关使用随机搜索和常规最优化方法相结合的方法求解安全系数极值的论文后Greco (1997)针对笔者的论点作过讨论他认为在一定模式指导下纯随机搜索方法一样可以精确地确定临界滑裂面但是这种纯粹的随机方法通常需要几千几万次的随机搜索才能达到足够的精度而本节所研究的把随机方法和确定性方法结合起来的途径在数值计算效率和收敛性能两个方面均具有明显的优势
4,1,2 最小安全系数的随机搜索对于某一边坡根据问题特点确定一个搜索区域如图4.16所示其轴线用z
o
表示其宽度为D
i
半带宽为d
i
=
D
i
/2这个搜索区域左右边界分别用两个滑裂面z
1
,z
r
来代表即
=
i
i
i
o
i
l
i
d
β
β
sin
cos
zz (4.56)
+=
i
i
i
o
i
r
i
d
β
β
sin
cos
zz (4.57)
式中D
=
(D
1
,D
2
,,D
n
)
T
称为搜索区宽度图 4,16 生成随机滑裂面示意图在搜索区内任意一个滑裂面可用下式表示
+=
i
i
ii
o
ii
dr
β
β
sin
cos
)5.0(zz (4.58)
式中
t
n
rrr ),...,,(
21
=r为该滑裂面和各控制点相对于轴线距离的系数
n
rrr,...,,
21
为伪随机数其值均在(0,1)之间随机搜索的步骤如下第4章 确定最小安全系数的最优化方法 107
1) 计算相应于滑裂面z
1
的安全系数F
0
2) 使用计算机伪随机发生器产生n个随机数
n
rrr,...,,
21
应用式(4.58)确定一个滑裂面z计算其相应的安全系数F
1
3) 比较F
1
和F
0
如果F
1
小于F
0
则F
0
和z
0
用F
1
和z
1
更新否则直接转入步骤(4)
4) 重复步骤2和3直到比较的次数足够大获得的最小安全系数足够小作为最优化法的初值足够好为止由于计算机产生的伪随机数具有很好的均匀性可以认为搜索区域内的滑裂面空间的每个部分都机会均等地被扫描了一遍搜索次数越多扫描密度越高成果越佳应用上述步骤在求解某一边坡的具体问题时需要根据经验确定一个搜索区域和一个搜索次数现在需要解决的问题是如何确定一个合适的随机搜索次数N
设想所寻找的临界滑裂面存在着一个条带称为保证区间任何落在该区间内的滑裂面如果被用作是最优化方法的初值都可以保证计算收敛那么现在的问题就变成按上述的随机方法生成滑裂面要生成多少个即可保证至少有一个落在这个保证区间中我们进一步假定该保证区间的宽度d
i
和搜索区间宽度D
i
存在一个固定和关系
n
n
D
d
D
d
D
d
m,..
2
2
1
1
=== (4.59)
式中m称为置信系数(0
<
m
<
1)
那么当随机数发生器为滑裂面第一个控制点按式(4.58)生成一个相对坐标时该点落在保证区间中的概率为m根据概率论理论n个控制点全部落在保证区间中的概率为m
n
由二项定理可知在随机生成的N个滑裂面中有r个落在保证区间中有N?r个未落在保证区间中的概率是
rNnrnr
N
r
N
mmCp
= )1()()(
(4.60)
那么在N次中至少有一次落在保证区间中的概率称成功率是
Nn
Nnn
N
m
mmCP
)(11
)1()(1
00
=
=
(4.61)
举例在图4.15所示的例中假如用3个控制点并且规定了搜索区域如图4.16中条带部分所示如果假定m
=
0.4那么如果进行了50次随机搜索其成功率是
%96)5.01(1
253
==P
从另一角度说如果取m
=
0.3为了使成功率达到96%我们需要作122次随机搜索
%96)3.01(1
1223
==P
应用上述步骤在求解某一边坡的具体问题是需要根据经验确定一个搜索区域和一个相应此搜索区的置信系数m搜索区域划得越大说明对临界滑裂的位置推测的把握越小m值也越小这样根据式(4.61)所算出的需随机搜索次数也就越多一般来说根据经验是不难划定搜索区域并确定m值的也许反复试算调整几次才能找到较好的结果
108 土质边坡稳定分析?原理
方法
程序使用计算机这也不是一件很难的事因此在实际应用时也可以根据经验来确定搜索次数STAB程序默认的功能则自动设定一个搜索次数
4,1,3 提高随机搜索效率的途径在将这一套随机搜索的方法运用到实际解题中去以前为了提高随机搜索效率还需要解决两个问题
1,滑裂面的几何合理性问题在随机生成的滑裂面中应该删去那些几何形状明显不合理的滑裂面这些形状不合理的滑裂面的存在不仅耗费许多无谓的计算时间而且在计算这种滑裂面的安全系数时还会带来数值分析收敛的麻烦根据边坡稳定分析实际特点对滑裂面和几何形状作出如下限制
(1) 相邻两个控制点的x坐标值应满足:
ii
xx >
+1
(4.62)
(2) 相邻两个控制点所连的直线与x轴的夹角α
i
应控制在一个合理的范围内如
oo
8045 ≤≤?
i
α (4.63)
2,应用简化方法提高随机搜索的效率在自由度较多时为使搜索达到一定保证率需要的随机搜索次数很多因此有必要找到一个途径将随机搜索使用的时间尽量减少一个很有效的方法就是在随机搜索时采用简化方法计算安全系数因为随机搜索的目的是确定一个初始滑裂面过高的安全系数精度是没有必要的在第3.6节中详细介绍的简化法就是用于随机搜索的简化方法简化法1用于光滑形状的滑裂面简化法2用于折线形状滑裂面在一般情况下使用简化法可以使随机搜索的时间减少2/3至1/2
[例4.10] 采用随机搜索重新分析[例4.9]
回到本节一开始提到的那个搜索整体极值失败的[例4.9]如图4.15所示仍使用四个点并用光滑曲线相连来代表滑裂面搜索区间如图4.17所示采用简化法1进行45次随机搜索图 4,17 对[例4.9]采用随机搜索求解示意第4章 确定最小安全系数的最优化方法 109
图4.18示其中的一些典型的随机生成的滑裂面由于使用简化法很快即可在此45
个滑裂面中找到一个相应最小安全系数的滑裂面如图4.17线1所示以此滑裂面作为初始值任何一种最优化方法均很轻松地收敛到相应整体极值的临界滑裂面2上相应安全系数为1.364 Morgenstern?Price法图 4,18 对[例4.9]采用随机生成的滑裂面
4,6 应用实题本节介绍采用上述最优化方法进行边坡稳定分析的工程实例
110 土质边坡稳定分析?原理
方法
程序
[例4.11] 紫坪铺库区左岸堆积体在第2章[例2.1]中讨论了紫坪铺库区左岸堆积体沿夹泥层ABCEFGHIJ的稳定安全系数在进行详细核算时需分析各种可能的滑裂面的安全系数因此要对这一折线型的滑裂面分段核算如图4.19所示初始输入一个滑裂面LGHM相应F
0
=
1.476令G沿GH
方向移动优化计算获得的临界滑裂面为L′G′H′M′相应F
m
=
1.121这一计算确定了以GH
中的一部分为底滑面临界破坏模式下一步还要核算其它可能的破坏模式第13章13.7
节详细介绍了这一核算过程STAB程序为折线型的底滑面提供了分段进行优化计算的工具但不可能一次自动地就把整体极值找出来图 4,19 紫坪铺库区左岸堆积体例
[例4.12] 小浪底水库大坝沿软弱夹层滑动分析高167m的小浪底水库大坝地基中存在缓倾角软弱夹层需要复核大坝下游坝坡在8
度地震情况下沿此软弱夹层滑动的最小安全系数用五个点来模拟滑裂面图4.20中AB
BC均为曲线CD段为软弱夹层故为直线DE段较短采用直线首先划定一个搜索区域如阴影部分所示采用简化方法2作223次随机搜索作为研究作者也曾用严格的通用条分法计算这223个随机滑裂面发现使用简化法2和使用严格法发现的最小安全系数的滑裂面都是第193个但使用时间减少一半根据随机搜索获得的初值最后获得的上下游边坡临界面滑裂面如图4.20中ABCDE所示分析此临界滑裂面图中的点位置是软弱夹层上部的终止点在这样的位置滑裂面中软弱夹层的长度最大因而安全系数最小故计算成果和常识判断是一致的图4.20
中E点的位置肯定了这样一个见识滑裂面通常从坝趾处逸出此处的抗滑力最小表4.2对比了以下三种情况的数值计算效率和精度采用纯粹的随机搜索不采用随机搜索单纯应用最优化方法本节介绍的将两者结合起来的方法第4章 确定最小安全系数的最优化方法 111
可以看到单纯依靠随机搜索当达到2000次后安全系数改善得很慢机时达11170
秒所得到的最小安全系数为1.253与精确值1.199相差颇远说明Boutrup (1980)和Siegel
(1981)建议的纯随机搜索的方法不甚有效单纯使用最优化方法虽然得到较小的最小安全系数1.234但无论从所用的机时和还是从得到的安全系数的精度来衡量都不如经过随机搜索确定初值的方法有效需要说明的是上述讨论系1988年的成果耗费机时系相应
IBM?PC?XT微机图 4,20 随机搜索在小浪底水库边坡稳定分析中的应用表 4,2 三种情况下数值计算效率与精度对比状况单纯随机搜索次数单纯最优化方法两者结合的方法
223 605 1336 2356 0 223
机时秒699 3431 6437 11170 2940 2661
最小安全系数F
m
1.293 1.270 1.261 1.253 1.234 1.199
注 两者结合的方法先进行223次随机搜索然后再以获得的最小安全系数为初值使用单形法确定F
m
耗费机
时系相应IBM?PC?XT微机
[例4.13] 天生桥二级闸首滑坡反演分析该滑坡在水电工程中为一个值得吸取教训的实例发生于1985年12月24日48人丧生由于边坡中存在一层饱和软粘土见图4.21中第5层当边坡开挖将此层粘土揭露出后即发生了滑面滑面下部沿堆积土和基岩的接触面各土层的岩土力学参数如表4.3
所示初始滑裂面如图4.22中A
1
B
1
C
1
D
1
E
1
F
1
G
1
H
1
所示最终获得的临界滑裂面A
2
B
2
…G
2
H
2
临界滑裂面的位置在上部与实测的偏离较大下部则与实测的吻合也是沿基岩接触面下滑
112 土质边坡稳定分析?原理
方法
程序图 4,21 天生桥二级闸首滑坡工程地质图表 4,3 各层土的物理力学特性指标图4.21
抗剪强度指标土层特性密度ρ (g/cm
3
)
φ (°) c (kPa)
1新填土1.85 21.8 19.6
2旧填土1.85 21.8 19.6
3第四系人工堆积1.85 21.8 0.0
4中细砂1.85 20.8 29.4
5砂质淤泥1.81 10.2 34.3
6卵砾石及砂1.90 24.2 0.0
7灰岩2.40 45.0 39.2
第4章 确定最小安全系数的最优化方法 113
图 4,22 天生桥二级闸首滑坡计算
4,7 非数值方法
4,7,1 模拟退火法模拟退火法是模拟金属退火过程的一种最优化计算方法唐立山等1994在金属退火时温度徐徐降低在每一个温度阶段系统的能量都要达到平衡状态使其达到最小值假设在某一随机过程中获得了一个相应某一随机变量的安全系数为F(z
i+1
)如果它与已往已经得到的最优解F(z
i
)的差值?F
=
F(z
i+1
)-F(z
i
)小于零那么自然要用F(z
i+1
)来代替
F(z
i
)但是如果F大于等于零则传统的随机搜索方法就要抛弃这一选择在模拟退火法中根据Monte Carlo法的Metropolis准则还要作一次抽签试验也就是要让计算机产生一个在(0,1)之间的随机数r
i
将其与exp[F/T]比较其中T就是这一阶段的温度在这里T是一个具有安全系数F特征的量如果r
i
<
]/exp[ TF那么仍然要用z
i+1
来代替z
i
这一处理有利于跳过局部极限形成的陷井但是这一做法终究抛弃了一个比现有值更优的自变量因此是带有一定的风险的我们用承担一定的风险作为代价使计算跳出局部极限形成的陷井可以看出在这一决策过程中以下做法是明智的
(1) 在退火过程中温度是逐渐缓慢减少的初期T较大用z
i+1
代替z
i
的可能性也较大因而允许自变量在较大的范围内变动以利于跳过陷阱确定一个较好的方向越到计算接近终止这一替代的可能性越小
(2) 只有在?F的绝对值足够少时即F(z
i+1
)与F(z
i
)十分接近时F(z
i+1
)这个表面上欠优的安全系数才有可能入选
(3) 随机数r
i
足够的小时F(z
i+1
)才有可能入选某种意义上说这意味我们把一部分成功的希望寄托在运气上
114 土质边坡稳定分析?原理
方法
程序这里关键的问题是T的数值在开始的时候取得很大其变小的过程是十分缓慢因此即使这一选择错了以后仍然有机会纠正过来在实际应用这一过程时还采用了一种记忆功能如果在某一温度计算步结束时实践证明把搜索方向从F(z
i
)改为F(z
i+1
)
并没有带来更好的结果我们还是要将原来被抛弃的F(z
i
)取回来进入下一温度计算步根据上述原则模拟退火法的计算步骤如下
(1) 设置初始温度T
0
这一温度应足够高以防止计算落入局部极值的陷井同时过高的初始温度可能降低搜索的效率建议先进行一定数量比如100次的随机搜索找到安全系数的最大值和最小值F
max
和F
min
将初始温度设为
minmax
FFT
o
= (4.64)
将优化计算过程分成m段第k温度段的温度按下式确定
k
Tk
TT )(
0
α?= (4.65)
式中T
k
为第k步迭代时的温度α
T
为温度冷却的阻尼系数α
r
应在[0,1]区间内通过试算确定
(2) 对于一个具有n个自由度的自变量z产生n个随机数v
=
(v
1
,v
2
,…,v
n
)每个ν
i
都在(?1,1)区间并将其换成方向导数u
=
(u
1
,u
2
,…,u
n
)
2/1
1
2
)(
∑
=
=
n
i
i
i
i
v
v
u (4.66)
为了减少自变量空间的粗糙度我们还引入一个系数
)....,,2,1(
)}(,),(,),(max{
)(
11
ni
zzzzzz
zz
nnii
ii
i
=
=
+?+?+
+
κ (4.67)
式中z
和z
+
分别代替第i个自变量在探索区域的最大值和最小值这样相应某一迭代步的自变量搜索为
1
0
)(
=?
k
rk
rr α (4.68)
式中?r
k
为第k步的带宽同样要随着迭代步的增加逐步减少这样自变量的增量为
nidurkz
iikii
,,2,1 L==? (4.69)
第4章 确定最小安全系数的最优化方法 115
[例4.14] 采用模拟退火法重新分析[例4.9]
对于[例4.9]第4.5.3节曾以随机搜索找到整体极值现使用模拟退火法进行分析采用同样的初始滑裂面同样用4个点联成一条光滑的曲线位于坡趾的点在优化过程中不动其它3个点沿水平方向移动因此共有3个自由度计算中采用α
T
=
0.9,α
r
=
0.98首先通过100次随机搜索生成初始温度T
0
随后进行m
k
=
300次退火模拟相应?r
o
=
0.27,T
k
=
1.006和1.001
从图4.23中可见在安全系数从初值1.427开始经过21676次迭代趋于极小值1.364
在向最终值过渡的过程中有时F值反而增加但最终顺利地收敛到整体极值图 4,23 采用模拟退火法分析[例4.9]安全系数收敛到整体极值过程
4,1,2 遗传算法肖专文等(1998)和Goh (1999)分别用遗传算法计算圆弧和任意形状临界滑裂面这里只对基本遗传算法作一概括的介绍有关遗传算法更为详细的内容请参见有关专著遗传算法基于自然界优胜劣汰的生物进化机制是一种自适应搜索技术在特定的竟争环境中依据这种规则适应性强的个体生存下来的可能性大Holland (1975)在其著作中介绍了GA的这一基本原理近来遗传算法作为一种很有应用潜力的优化方法倍受关注原因是它能有效地搜索一个复杂问题的整个定义域全局搜索遗传算法是一个迭代操作过程该迭代操作针对多个待求问题的候选解进行的这些候选解目标函数组成一个固定规模的群体候选解是模拟生物体的染色体将待求问题编码而形成的最初候选解的群体是随机生成的每一个染色体代表给定优化问题的一个可能的解每一个候选解由一串基因组成通常表示为一由0和1组成的二进制串组成染色体的每一个基因二进制子串代表一个待优化的参数整个二进制串称为染色体或个体由代表各个变量的基因子串集合而成使用目标函数可计算一个染色体对应的目标函数值安全系数进而可以
116 土质边坡稳定分析?原理
方法
程序确定与每一个染色体相对应的适应值安全系数的函数染色体通过迭代而进化在每一个迭代步骤中由多个染色体组成的群体称为一代下一代群体中的每一个染色体是由上一代群体中的两个染色体相互结合交叉操作形成也可以通过直接改变上一代群体中的某个染色体变异操作而形成从父代和子代中选择某些适应值大的染色体而淘汰其它适应值小的染色体选择操作可以形成新一代的群体而且通过这种优胜劣汰机制每一代群体中染色体个体的个数保持不变染色体所对应的适应值决定了其生存与产生后代的能力适应值最大最小安全系数最小的染色体最有可能被选择并用于产生下一代染色体这一迭代过程延缓下去直到寻找到最优解或者达到某一个迭代终止的标准图4.24所示的流程图给出了遗传算法的主要操作步骤为了实现遗传算法的一个优化循环需要选择交叉和变异三个基本操作下面对这些基本操作作一简单的介绍
1,遗传算法的基本步骤
(1) 选择操作选择操作可以确定选用当代群体中哪些个体生成下一代群体中新的个体选择操作的基本思想是将适应值超出平均水平的染色体以概率的方式挑选出来比如在轮盘选择策略(Holland,1975)中一个父代个体被选择的概率与其对应的适应值成正比
(2) 交叉操作交叉操作用于确定父代的染色体怎样结合形成子代的染色体交叉操作针对于两个经选择操作选定的父代染色体这两个父代染色体根据给定的交叉概率相互交换基因值二进制数单点交叉是最基本的交叉操作在这一操作中随机选取二进制串上的交叉点然后两个染色体在该点相互交换基因值在双点交叉中随机选取两个交叉点然后两个染色体相互交换该两点间的染色体部分图4.25 (a)和(b)分别给出了单点交叉操作和双点交叉操作的例子图 4,24 遗传算法的流程图图 4,25 交叉操作
(a) 单点交叉(b) 双点交叉第4章 确定最小安全系数的最优化方法 117
(3) 变异操作变异操作是按给定的概率随机改变某个染色体中一个二进制位上的数值例如在图4.26中由0变为1
变异操作的主要目的是提供新的基因以防止迭代落入一个局部峰值区域
(4) 参数初始化群体规模交叉概率变异概率的选取依不同的问题而异通常需要做一些试验确定图 4,26 变异操作
2,与传统方法的比较遗传算法的实现是对一个群体的多方位操作该群体是由许多潜在解组成的而不是仅对单一的解进行操作这一并行性质决定了遗传算法全局优化的特点遗传算法使用适应值而不使用适应函数的梯度或更高次的微分及其它信息所以适用于解决各种不同的优化问题具有鲁棒性强的特点遗传算法属于非确定性算法但是它不同于纯随机的算法遗传算法将随机搜索与方向性搜索结合在一起因此可以平衡优化搜索的两个相互矛盾的策略搜寻最优解与全局搜索爬山法是搜寻最优解但忽视全局搜索的典型例子而随机搜索是全局搜索的典型方法但它忽视了可能存在最优解的局部空间正如其它优化算法一样遗传算法也有其局限性主要缺点是要化费大量的计算时间但随着计算机计算速度的迅速提高这一缺点将被克服
3,遗传算法的实现
Goh (1999)的算例使用楔体分析的程序GWEDGEM该程序使用多楔体方法确定边坡的安全系数GWEDGEM的当时版本将楔体的数量限制在6块以内该程序引入遗传算法搜索临界滑裂面由于遗传算法通常用于搜索最大值故将目标函数转化为适应函数
2
min
max
)(
1
F
f = (4.70)
使用二进制编码表示参数变量二进制串的长度取决于所需的精度例如对于一个四位的二进制编码其精度大约为搜索区间的1/2
4
(1/16)当然根据所需的精度不同的变量可以使用不同长度的二进制串表示在GAWEDGE中每一个变量采用9位长度的二进制串表示这样对于一个含6个楔体的问题共17个变量染色体的总长度为
153位
[例4.15] ACAD考题3(a)
Goh在其论文中介绍的3个考题恰为将在第11章第11.4.2节中ACAD程序中的Exl(e)
Ex3(a)和Ex4现择Ex3(a)为例介绍遗传算法与其它方法计算成果的比较情况如图4.27
所示
118 土质边坡稳定分析?原理
方法
程序图 4,27 ACAD考题3(a)计算示意图在使用遗传算法进行分析时选用了以下参数进化次数为300群体大小为5交叉概率为0.75变异概率为0.02图4.28示其进化过程在第290次时计算收敛计算次数为290×5=1450其中5为自由度数目获得的最小安全系数为1.286 Goh通过表4.4 比较了多种优化方法的成果其中涉及本书作者(Chen & Shao,1988)的成果系作者为ACAD
所作计算成果之误参见第11章表11.15
图 4,28 进化次数与安全系数的关系表 4,4 对[例4,15]的计算结果比较方法楔体数目安全系数迭代次数
GWEDGEM 6 1.288 1820
GWEDGEM 7 1.288 1740
GWEDGEM 8 1.483 2490
GWEDGEM 9 2.099 3750
GAWEDGE (群体大小为5) 6 1.286 290×5=1450
SLOPE (Fredlund & Krahn,1977) 1.364-1.378
Spencer (1973) 1.24
Chen & Shao (1988) 1.242
Donald & Giam (1989) 1.27
第4章 确定最小安全系数的最优化方法 119
4,8 本章附录
4,8,1 本章数据文件表4.5
表 4,5 本 章 数 据 文 件有关章节系列号数据文件名内容
04-01-01 CPT4-3.DAT
[例4.1] 光滑滑裂面
4.2.1
04-01-02 CTP4-4.DAT
[例4.1] 折线滑裂面
04-02-01 EX41.DAT [例4.2] 两个自由度的算例用单形法
4.3
04-02-02 EX42.DAT [例4.3] 使用单形法计算图4.6例
4.5
04-03-01 CPT4-5.DAT [例4.10] 随机搜索
4.6
04-04-01 TIAN.DAT [例4.13] 天生桥滑坡例
4.7 04-05-01 CPT4-6.DAT [例4.14] 模拟退火参考文献
1 Boutrop,A,W,and Lovell,C,W,Search technique in slope stability analysis,Engineering Geology,1980,16:
51-61
2 Chen,Z,and Shao,C,Evaluation of minunmum factor of safety in slope stability analysis,Canadian
Geotechnical Journal,1988,25(4),735-748
3 陈祖煜,邵长明,最优化方法在确定边坡最小安全系数方面的应用,岩土工程学报,1988,10(4)
4 陈祖煜,邵长明,用最优化方法搜索最小安全系数的一个改进,岩石力学在工程中的应用,第二次全国石力学与工程学术会议论文集知识出版社,1989,71-75
5 陈祖煜,用DFP法计算边坡稳定最小安全系数时导数的计算方法.岩土力学数值方法的工程应用,第三届全国岩石力学数值计算与模型实验学术研讨会论文集,上海同济大学出版社,1990,797-803
6 Chen,Z,Experience with the search of minimum factors of safety of slopes,Proceedings,6th Australian New
Zealand Conference on Geomechanics,1992,426-431
7 Chen,Z Random trials used in determining global minimum factors of safety of slopes,Canadian
Geotechnical Journal,1992,29(2),225-233
8 Celestino,Y,B,and Duncan,J,M,Simplified search for noncircular slip surface,10
th
Int,Conf,On Soil
Mech,And Found,Engr.1981,Vol,391-394
9 Li,K,S,and White,W,Rapid evaluation of of the critical slip surface in slope stability problems.
International Journal for Numerical and Analytical Methods in Geomechanics,1987,11,449-473
10 Fox,R,Optimisation methods in engineering designs,Addison-Wesley Publication Company,1971
11 华东水利学院,土工原理与计算,北京水利出版社,1979
12 潘家铮,建筑物的抗滑稳定和滑坡分析,北京水利出版社,1980
13 孙君实,条分法的数值分析岩土工程学报,1984,6(2),1-12
14 Shoup,T,E,and Mistroup,F,Optimisation methods with alications for personal computers,Prentice-Hall Inc.
Eaglewood Cliffs,1987
15 Siegel,R,A.,Kovacs,W,D,and Lovell,C,W,Random surface generation in stability analysis,ASCE
Journal of Geotechnical Engineering,1981,107,996-1002
16 Stewart,G,W,III,A modification of Davidon’s minimization method to accept difference approximations of
derivatives,Journal of Association for Computing Machinery,1966,1,72-83
17 唐立山谢云尤矢勇罗祖华,非数值并行算法,第一册,模拟退火算法,科学出版社,1994
120 土质边坡稳定分析?原理
方法
程序
18 肖专文,张奇志,梁力,林韵梅,遗传进化算法在边坡稳定性分析中的应用,岩土工程学报 1998,1:44-
46
19 闫中华,均质土坝与非均质土坝稳定安全系数极值分布规律及电算程序简介,水利水电技术,1983年第
7期
20 张天宝,四边形块方法和沉抗土坝的稳定分析兼论复合土坡的最危险滑弧位置,水利学报,1981
21 邹广电,边坡稳定分析条分法的一个全局优化算法,岩土工程学报,2002,24(3),309-312
22 周文通,最优化方法在土坝稳定分析中的应用,土石坝工程,1984
23 Nguyen,V,U,Determination of critical slip surface,Journal of Geotechnical Engineering,ASCE,1985,111:
238-251
24 Goh,A.T.C,Genetic algorithm search for critical slip surface in multi-wedge stability analysis,Candian
Geotechnical Journal,36,(2)383-391,1999.
25 Greco,Y.R,Efficient Monte-Carlo technique for locating critical slip surface,ASCE J,Geotechnical
Engineerinf,122,(7)517-525,1997.
26 朱伯芳黎展盾张壁诚,结构优化设计原理与应用,北京水利电力出出版社,1984
