《土质边坡稳定分析》教学资料（原理、方法、程序）第四章 确定最小安全系数的最优化方法

第4章确定最小安全系数的最优化方法 4.1概述 4.1.1应用最优化方法确定最小安全系数 边坡稳定分析应包含下面两个步骤 )对滑坡体内某一滑裂面按第2或第3章介绍的方法,确定其抗滑稳定安全系数: 2)在所有可能的滑裂面中,重复上述步骤,找出相应最小安全系数的临界滑裂面 在讨论了计算一个滑裂面的安全系数的方法后,本章介绍边坡稳定分析极限平衡法的 第二步,就是要寻找最小安全系数的方法。 如果滑裂面曲线为y(x)。那么,这个问题具体化为寻找下列泛函的极值 F=F) 岩土工程中的边坡的几何形状各异,材料通常是非均质性,纯解析的变分原理很难进 行极值计算。用最优化方法通过数值方法求解,是一个比较现实可行的途径。 我国学者较早开展应用数值规划方法求解安全系数极值的问题(闫中华,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.枚举法 枚举法的基本思想是,根据一定的模式,比较不同自变量的目标函数,经过筛选,最

第4章 确定最小安全系数的最优化方法 4. 1 概述 4. 1. 1 应用最优化方法确定最小安全系数 边坡稳定分析应包含下面两个步骤 1) 对滑坡体内某一滑裂面按第 2 或第 3 章介绍的方法 确定其抗滑稳定安全系数 2) 在所有可能的滑裂面中 重复上述步骤 找出相应最小安全系数的临界滑裂面 在讨论了计算一个滑裂面的安全系数的方法后 本章介绍边坡稳定分析极限平衡法的 第二步 就是要寻找最小安全系数的方法 如果滑裂面曲线为 y(x) 那么 这个问题具体化为寻找下列泛函的极值 F = F( y) (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土质边坡稳定分析一原理·方法程序 终找到最小值。这是最原始、简单的方法 如图41中示,任一圆弧可用其圆心坐标x,y)和半径r确定。其相应的安全系数F可 表达为 F=f(x0,y0,D3) 式中:D为滑弧深度,即圆弧最低点的坐标,可知 D Vo 显然这是三个自由度的问题。采用枚举法,不断地改变x。y和D的数值,逐一比较 相应的安全系数,最终找到最小的安全系数。在具体操作中,先固定一个D,然后在圆 可能的位置中布置一个网格,见第12章图121。设网格的中心坐标为x和y2。在左、右方 向,各布置了N格,在上、下方向各布置了N格,则共计有(2N2+1)×(2N+1)个网格点,分 别以该网格点为圆心,以D,为滑弧深度计算相应安全系数,找出最小的安全系数。然后改 变一个D值,重复相同的步骤。在这一过程中,有可能出现圆弧和边坡不相交的情况,则 应令程序自动抛弃该圆弧。同样,D也是以一个中心值起算,在其上、下各布置N层,这 样总计计算(2N2+)×(2N+1)x(2N+1)个圆弧。 在枚举法的基础上,用0618法或其他方法,提高搜索最小安全系数工作的效率。这 方面的工作可参阅有关文献。 图4.1圆弧滑裂面 2.数值分析方法 随着计算机的发展,数值分析方法逐步形成一门完整的学科,统称为最优化方法( Method of optimization)。这一方法又可分为两大类 第一类称模式搜索法( Pattern search method)。其基本思想是,根据一定的模式,比较不 同自变量的目标函数,确定最优的搜索方向,最终找到最小值。 另一类称牛顿法。它要通过解析手段寻找使目标函数F对自变量κ的偏导数为零的极 值点(aFOx2=0,i=1,2…,m)。同时,从理论上讲,还需要满足由二阶导数形成的 Hessian矩 阵正定这个达到极小值的充分条件。此类方法中以导数为研究的主要对象,因此,也称为 以导数为基础的方法( Gradient- based- method)。一般认为,当自由度较多时,直接搜索法效 率较低。此时需要考虑牛顿法体系的分析方法

88 土质边坡稳定分析 原理 ⋅ 方法 ⋅ 程序 终找到最小值 这是最原始 简单的方法 如图 4.1 中示 任一圆弧可用其圆心坐标(x0, y0)和半径 r 确定 其相应的安全系数 F 可 表达为 ( , , ) 0 0 Ds F = f x y (4.2) 式中 Ds为滑弧深度 即圆弧最低点的坐标 可知 0 D r y s = + (4.3) 显然这是三个自由度的问题 采用枚举法 不断地改变 x0, y0 和 Ds 的数值 逐一比较 相应的安全系数 最终找到最小的安全系数 在具体操作中 先固定一个 Ds 然后在圆心 可能的位置中布置一个网格 见第 12 章图 12.1 设网格的中心坐标为 xc和 yc 在左 右方 向 各布置了 Nx 格 在上 下方向各布置了 Ny 格 则共计有(2Nx+1)×(2Ny+1)个网格点 分 别以该网格点为圆心 以 Ds 为滑弧深度计算相应安全系数 找出最小的安全系数 然后改 变一个 Ds 值 重复相同的步骤 在这一过程中 有可能出现圆弧和边坡不相交的情况 则 应令程序自动抛弃该圆弧 同样 Ds 也是以一个中心值起算 在其上 下各布置 Nd 层 这 样总计计算(2Nx+1)×(2Ny+1)×(2Nd+1)个圆弧 在枚举法的基础上 用 0.618 法或其他方法 提高搜索最小安全系数工作的效率 这 方面的工作可参阅有关文献 图 4. 1 圆弧滑裂面 2. 数值分析方法 随着计算机的发展 数值分析方法逐步形成一门完整的学科 统称为最优化方法(Method of optimization) 这一方法又可分为两大类 第一类称模式搜索法(Pattern search method) 其基本思想是 根据一定的模式 比较不 同自变量的目标函数 确定最优的搜索方向 最终找到最小值 另一类称牛顿法 它要通过解析手段寻找使目标函数 F 对自变量 zi 的偏导数为零的极 值点(∂F/∂zi = 0, i = 1,2,...,n) 同时 从理论上讲 还需要满足由二阶导数形成的 Hessian 矩 阵正定这个达到极小值的充分条件 此类方法中以导数为研究的主要对象 因此 也称为 以导数为基础的方法(Gradient−based−method) 一般认为 当自由度较多时 直接搜索法效 率较低 此时需要考虑牛顿法体系的分析方法

第4章确定最小安全系数的最优化方法89 本节介绍直接搜索法中的单纯形法和 Powell法以及牛顿法中负梯度法和DFP法 ( Davidon- Fletcher-Powell method)由于这些方法的原理在众多的文献及教科书中都有所介 绍,这里不拟全面阐述其原理,而着重讨论在边坡稳定分析领域需解决的特殊问题。 3.非数值分析方法 第三类方法是在近期计算机发展基础上形成的,称为非数值方法。这一方法被广泛应 用于管理科学、计算机科学、分子物理学及超大规模集成电路设计中,用于解决组合优化 问题。非数值分析利用计算机具有容量大、计算速度快的优点,通过大量随机采样来找到 目标函数的最优值。近期涌现了诸如模拟退火,遗传算法,神经网络和蚂蚁算法等。均属 此类。在边坡稳定分析中,上述算法均有人尝试。本章47节将对此类方法作一简要介绍 4.2任意形状滑裂面的模拟和目标函数的确立 4.2.1任意形状滑裂面的模拟 最优化问题的提法是:对于一个具有n个自变量的向量z=(1,z2,…,m),确定其目标 函数F的最小安全系数Fm,相应的自变量为zn 为此,需要对式(41)中的曲线y(x)用若干参数来模拟。也就是说,需要将任意形状滑 裂面yx)用z来近似表达。将滑裂面曲线用m个点A1,A2,…,Am离散,见图42。也就是将 此m个点用直线或光滑的曲线连起来,以近似模拟此曲线。此m个点的坐标用x(i=1,2,…m) 表示: 次弱夹层 图4.2任意形状滑裂面 旦这种连接的模式确定,安全系数F即可表达成此m个点的坐标x,y2x2y2…xmyn 的函数 F=F(x,yi, x2, y2,. xm,y

第 4 章 确定最小安全系数的最优化方法 89 本节介绍直接搜索法中的单纯形法和 Powell 法以及牛顿法中负梯度法和 DFP 法 (Davidon−Fletcher−Powell method) 由于这些方法的原理在众多的文献及教科书中都有所介 绍 这里不拟全面阐述其原理 而着重讨论在边坡稳定分析领域需解决的特殊问题 3. 非数值分析方法 第三类方法是在近期计算机发展基础上形成的 称为非数值方法 这一方法被广泛应 用于管理科学 计算机科学 分子物理学及超大规模集成电路设计中 用于解决组合优化 问题 非数值分析利用计算机具有容量大 计算速度快的优点 通过大量随机采样来找到 目标函数的最优值 近期涌现了诸如模拟退火 遗传算法 神经网络和蚂蚁算法等 均属 此类 在边坡稳定分析中 上述算法均有人尝试 本章 4.7 节将对此类方法作一简要介绍 4. 2 任意形状滑裂面的模拟和目标函数的确立 4. 2. 1 任意形状滑裂面的模拟 最优化问题的提法是 对于一个具有 n 个自变量的向量 z = (z1, z2, ..., zn) 确定其目标 函数 F 的最小安全系数 Fm 相应的自变量为 zm 为此 需要对式(4.1)中的曲线 y(x)用若干参数来模拟 也就是说 需要将任意形状滑 裂面 y(x)用 z 来近似表达 将滑裂面曲线用 m 个点 A1, A2, …, Am离散 见图 4.2 也就是将 此 m 个点用直线或光滑的曲线连起来 以近似模拟此曲线 此 m 个点的坐标用 zi (i=1,2, ,m) 表示       = i i i y x z (4.4) 图 4. 2 任意形状滑裂面 一旦这种连接的模式确定 安全系数 F 即可表达成此 m 个点的坐标 x1, y1, x2, y2, .... xn, yn 的函数 F = F( ) x1, y1, x2 , y2 ,Lxm , ym (4.5)

90土质边坡德定分析一原理方法程序 在进行最优化搜索过程中,A1,A2…,Am将移到临界滑裂面的位置B,B2,Bn,见图 42,此处m=6,其中端点A,Am原来在边坡线上,有可能移到边坡线外或内,如图42中 的B1,Bn。为此,需要通过一定的处理方式,分别找到它们和边坡线的交点B1和Bn,仍 以B1,B2,…,Ba作为硏究的对象。对均匀的土质边坡,通常希望滑裂面比较光滑,此时,可 采用三次或更高次的样条函数连接这些点。第11章将详细介绍使用样条函数构筑光滑曲线 滑裂面的方法。当然,也可以采用直线和光滑曲线的组合构筑滑裂面。例如图42,A3,A,A3 A用曲线相联;A2,A3用直线相联。 应用样条函数构筑光滑滑裂面,对于减少自由度,提高数值计算效率具有重要意义 例41]折线和光滑滑裂面比较 如图43所示例,如果用四个点模拟滑裂面,那么,折线滑裂面和曲线滑裂面相应的 最小安全系数分别为1489和1.364,差别颇大。换句话说,如果使用折线模式,要得到与 曲线模式相同的精度,就要更多的离散点,这意味着更多的自由度。而数值分析的机时和 收敛难度是随自由度按指数增加的。值得注意的是,在笔者所知的所有这方面的文献中, 还没有其他研究者使用过曲线模拟任意形状滑裂面 在稳定分析中,对土质比较均匀的边坡,常采用圆弧滑裂面。此时,自变量可以是圆 弧的圆心的x,y坐标和半径r,由于自变量仅三个,最优化方法搜索最小安全系数收敛性能 很好,因此将不作为重点在本章中讨论。 图4.3说明用不同模式模拟滑裂面的差别 滑裂面1-直线模式;滑裂面2-曲线模式 4.1.2目标函数的确定 当滑裂面z的离散模型确定后,安全系数便是m个控制点A1,A2,,A的函数。在优 化计算过程中,这m个点中有n个点各沿某一设定方向β向临界滑裂面移动,或者不规定 方向任其自由移动。其余m-n个点由于问题本身的要求可以固定。滑裂面上任意一点的z 可以用相应于一个初始滑裂面的相对坐标来代表,见图42 =+ 式中:=1,2,n,d为第个点沿β移动的距离

90 土质边坡稳定分析 原理 ⋅ 方法 ⋅ 程序 在进行最优化搜索过程中 A1, A2, …, Am将移到临界滑裂面的位置B1,B2,...,Bm ′ ′ 见图 4.2 此处 m = 6 其中端点 A1, Am原来在边坡线上 有可能移到边坡线外或内 如图 4.2 中 的B1, Bm ′ ′ 为此 需要通过一定的处理方式 分别找到它们和边坡线的交点 B1 和 Bm 仍 以B1,B2,...,Bm 作为研究的对象 对均匀的土质边坡 通常希望滑裂面比较光滑 此时 可 采用三次或更高次的样条函数连接这些点 第 11 章将详细介绍使用样条函数构筑光滑曲线 滑裂面的方法 当然 也可以采用直线和光滑曲线的组合构筑滑裂面 例如图 4.2, A3, A4, A5, A6用曲线相联 A2, A3用直线相联 应用样条函数构筑光滑滑裂面 对于减少自由度 提高数值计算效率具有重要意义 [例 4.1] 折线和光滑滑裂面比较 如图 4.3 所示例 如果用四个点模拟滑裂面 那么 折线滑裂面和曲线滑裂面相应的 最小安全系数分别为 1.489 和 1.364 差别颇大 换句话说 如果使用折线模式 要得到与 曲线模式相同的精度 就要更多的离散点 这意味着更多的自由度 而数值分析的机时和 收敛难度是随自由度按指数增加的 值得注意的是 在笔者所知的所有这方面的文献中 还没有其他研究者使用过曲线模拟任意形状滑裂面 在稳定分析中 对土质比较均匀的边坡 常采用圆弧滑裂面 此时 自变量可以是圆 弧的圆心的 x, y 坐标和半径 r 由于自变量仅三个 最优化方法搜索最小安全系数收敛性能 很好 因此将不作为重点在本章中讨论 图 4. 3 说明用不同模式模拟滑裂面的差别 滑裂面 1−直线模式 滑裂面 2−曲线模式 4. 1. 2 目标函数的确定 当滑裂面 z 的离散模型确定后 安全系数便是 m 个控制点 A1, A2, ..., Am的函数 在优 化计算过程中 这 m 个点中有 n 个点各沿某一设定方向 βi向临界滑裂面移动 或者不规定 方向任其自由移动 其余 m−n 个点由于问题本身的要求可以固定 滑裂面上任意一点的 zi 可以用相应于一个初始滑裂面的相对坐标来代表 见图 4.2         = o i o o i i y x z (4.6)       = + ∆ = + i i i o t i o i i d β β sin cos z z z z (4.7) 式中 i=1,2,…,n, di为第 i 个点沿βi移动的距离

第4章确定最小安全系数的最优化方法 对于固定的点,其自由度为零,如图42中的A2,A3,当某点需沿一个软弱夹层移动时, 则其自由度为1。如对该点的移动方向无特殊要求,如图42中的其它点,该点的自由度为 2。问题的自由度为各点自由度的总和,搜索最小安全系数的问题具体化为求下列函数的最 小值问题 F=F(d1,d2,…,dn,阝1,阝2,,阝n) 如前所述,滑裂面的一部分可能是光滑曲线,另一部分受软弱夹层的控制为直线。STAB 程序是这样规定这一定义格式的。设滑裂面被m个点A,A2 Am分为m-1段,从上 交点向下交点编号为1,2,…,m-1。在STAB程序中,用LNO代表此m-1段中为直线段的 线段总数,在数组LOO①,I=12,,LNO中存入这些直线段的编号。其它线段则默认为曲 线。当LNO为零时,滑裂面为一光滑曲线;当LNO为一大于m-1的数时,则程序默认为 全部为直线段。这两种情况都无须填写LOO(D 4.3模式搜索法 4.3.1单形法(朱伯芳等1984) 1.建立初始单形 对某一初始向量z,按下面模式构筑n个向量z(i=1,2,,n,组成m+1个顶点 z=[=+p,=2+q,…2+q z2=[=8+q-2+p…,z0 z=[1+g,22+q,……,m+p 其中 (4.10) (n+1) (4.11) 式中:a为初始步长,可根据实际情况进行调试 2.确定搜索方向 计算单纯形n+1个顶点的目标函数值,比较其大小,从中找出目标函数最大点z和次 大点和最小点z,然后按下式确定下一步的搜索方向 式中上标ν表示迭代次数

第 4 章 确定最小安全系数的最优化方法 91 对于固定的点 其自由度为零 如图 4.2 中的 A2, A3 当某点需沿一个软弱夹层移动时 则其自由度为 1 如对该点的移动方向无特殊要求 如图 4.2 中的其它点 该点的自由度为 2 问题的自由度为各点自由度的总和 搜索最小安全系数的问题具体化为求下列函数的最 小值问题 ( , , , , , ,..., ) F = F d1 d2 L dn β1 β2 βn (4.8) 如前所述 滑裂面的一部分可能是光滑曲线 另一部分受软弱夹层的控制为直线 STAB 程序是这样规定这一定义格式的 设滑裂面被 m 个点 A1, A2, ……, Am分为 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 z q z q z p z q z p z q z p z q z q [ , , , ] [ , , , ] [ , , , ] 0 0 2 0 1 0 0 2 0 1 2 0 0 2 0 1 1 = + + + = + + + = + + + LL LL LL LL z z z (4.9) 其中 a n n n p 2 ( +1) + −1 = (4.10) a n n q 2 ( +1) −1 = (4.11) 式中 a 为初始步长 可根据实际情况进行调试 2. 确定搜索方向 计算单纯形 n+1 个顶点的目标函数值 比较其大小 从中找出目标函数最大点 zH 和次 大点 zG 和最小点 zL 然后按下式确定下一步的搜索方向       = ∑ − − + v H v i n i v n n z z z 1 2 1 (4.12) 式中上标 v 表示迭代次数