范 鹏
(辽宁省水利水电勘测设计研究院有限责任公司,辽宁 沈阳 110006)
重力坝结构设计较为复杂,需要满足经济和安全的要求。目前常用的重力坝设计方法为数值模拟[1]和理论计算[2]的方法,两者计算量均较大,设计效率较低。因此,有人提出使用代理模型的方法,代替数值模拟模型进行重力坝结构优化设计。目前,较为常用的代理模型方法为克里金插值法;该方法在结构设计方面已经有了较多的使用,如液压阀结构设计[3]、土地景观利用稳定性评价[4]、钢桥结构优化[5]等方面有较多的使用并取得了较好的效果,使用代理模型进行结构优化设计是可靠的。粒子群算法是一种优化算法,可以较快的速度获取全局最优值,在无人机航程规划[6]、电磁态势预测[7]、水库优化调度[8]、植物冠层图像分割[9]等方面有较多的使用,优化方案准确可靠。
结合前人研究,提出使用PCE近似模型代替数值模拟计算,结合改进粒子群算法对重力坝结构优化设计进行研究。
在数学上已经证明任何函数尽可以使用多项式进行展开,且展开精度是可靠的。使用随机展开是一个随机计算的过程,可以将系统响应近似为随机变量的展开,从而对系统的响应规律进行概率统计。该方法在不确定分析中具有较多的使用。
建立概率空间Ω内的一个不确定模型:
y=Y(x,t)
(1)
式中,x—不确定因素参数,取值范围为Ω;定义一个N组独立的不确定因素输入变量:x=(x1,x2,…,xN)T,第i个变量边缘部位的概率密度为:fi(xi);则联合概率密度如下:f(x)=f1(x1)f2(x2)…fN(xN)。因此,原模型式(1)可以使用多项式进行展开:
(2)
式(2)可简化为:
(3)
式中,ξ—满足随机分布的变量形式;a—确定系数;Φ—P阶多项式。
在进行多项式展开时,基函数是所有变量产生的正交多项式的乘积,满足正交关系:
(4)
式中,E—数学期望;δ—Kronecher算式。
(5)
前P阶多项式阶段形式为:
(6)
保留P阶多项式计算方法如下:
(7)
根据式(4)计算每一项的展开式
ai(t)=E[Y(x,t)Φi(ξ)]
(8)
可根据多项式的混沌系数求取各阶矩。
E[y′(x,t)]=a0
(9)
(10)
使用高斯投影法进行混沌系数分配
(11)
式中,ξi1—各个变量的高斯积分点;ωiN—各变量的权重。
以实际工程中遇到的重力坝设计为例,选取非溢流断面作为优化设计断面。重力坝设计高为25m,顶宽为4.5m。设计断面如图1所示,浇筑材料见表1。
图1 重力坝坝体断面图
浇筑位置材料号弹性模量E/GPa泊松比容重/(KN/m3)抗压强度/MPa抗拉强度/MPa坝体125.50.172013.41.54坝基4200.2525//
根据表1的浇筑材料参数属性,建立关于重力坝安全系数、重力坝断面面积和坝体断面最大宽度的近似模型。近似模型建立方法如下:使用数值模拟计算的方法,确定输入变量和输出结果之间的关系,选取一种方法进行数学关系拟合,对近似模型的计算精度进行精度检验,满足设计要求,则可对新输入变量所的输出结果进行预测分析。
建立近似模型之前,首先需要选取足够代表性的、均匀分布的样本点数据。目前,选取样本点常用的方法为拉丁超立方法,使用该方法选取的样本点具有均匀分布、具有代表性和计算快速收敛的特点。样本点选取个数的要求如下:当变量个数为n个时,样本点选取个数应该大于(n+1)(n+2)/2个。选取上游边坡系数(C1)、下游边坡系数(C2)、上游折点度与坝高比值(C3)、下游折点高度与坝高比值(C4)一共4个变量作为重力坝结构设计变量。使用拉丁超立方的方法选取样本点55个,用来建立近似模型,使用15组数据进行精度检验。
3.2 重力坝安全系数、断面面积、断面最大宽度近似模型建立
使用ANSYS作为数值模拟计算软件,该软件是目前较为常用的一种有限元数值模拟计算软件,在各行业内有较为广泛的使用[10- 11]。使用MATLAB将ANSYS计算结果进行调用,获取坝体安全系数和输入变量之间的关系,根据上述方法建立近似模型,并进行精度检验。
通常使用下式来判断近似模型精度。
δ1=Max{|fi-f|}(i=1,…,N)
(12)
(13)
式中,δ1、δ2—分别代表最大误差和均方根误差;N—检测点个数;i—第i个检测点。
式(12)、式(13)的值越小表示代理模型精度越高。
精度评价结果见表2。
表2 KRIGING代理计算模型精度评价
根据表2的计算精度检验结果可以得出最大误差仅为4.13%,均方根误差为0.23%,近似型精度满足工程精度要求。可以使用所建立的近似模型代替ANSYS数值模拟计算结果。
基础的粒子群算法具有局部最优值的局限性,因此引入交叉因子、学习因子和动态惯性因子克服粒子群算法的局限性。
(14)
式中,randj(0,1)—在[0,1]范围内随机的一个数;jrand—[1,D]范围内随机的一个整数;p—两个粒子之间产生交叉的概率。
计算完成后对粒子的最优位置进行更新,更新结果如下:
(15)
(2)动态惯性因子
对于复杂的问题,惯性权重并非简单的按照线性衰减。因此,在进行粒子群算法获取最优值时可以使用动态惯性因子对粒子群算法进行改进。计算动态惯性权重,需要明确粒子进化度和聚合度两个含义:
(16)
式中,m—粒子进化度;n—粒子聚合度;Pg,i、Pg,i-1—当前时刻全局最优值和上一时刻的历史全局最优值;Ps—粒子群的规模,即总数。
惯性权重是一个与粒子进化度和粒子聚合度密切相关的参数,计算方法如下:
ω=f(m,l)=ω0-0.6m+0.15n
(17)
式中,ω0—初始惯性权重。
设计变量采用向量的形式表示,Z=[C1,C2,C3,C4],各参数意义如下:上游边坡系数(C1)、下游边坡系数(C2)、上游折点高度与坝高比值(C3)、下游折点高度与坝高比值(C4)。
Z1=[0,0.60,0.32,0.86];Z1=[0.21,0.80,0.66,0.93];K≥3;D≥25m
(18)
式中,Z1、Z2—分别表示各设计变量的最小值和最大值;K—坝体安全系数;D—坝体横断面最大宽度。
设计目标函数时需要引入惩罚函数,惩罚函数的引入可以加快运算速率提高计算精度。建立新的目标函数如下:
F=f(A)·e[f(a)+f(b)]
(19)
式中,f(a)、f(b)—惩罚函数;F—建立的新的目标函数。
使用MATLAB编写粒子群算法基于近似模型(PCE)的模型,使用电脑进行多次计算,获取重力坝结构最优设计方案,初始设计变量值初始和优化后的变量取值结果见3。
表3 结构设计优化前后各变量参数
通过重力坝结构优化设计,在满足安全系数的要求下,断面面积减小率为(337.55-313.6)/474=7.1%。根据近似模型优化设计结果建立数值模拟模型,计算得到安全系数K=6.51,满足重力坝安全设计要求。
(1)使用PCE方法建立的数值模拟近似模型具有较高的计算精度,可以满足工程设计要求,同时可以有效减少数值模拟在重力坝设计中的计算量。
(2)使用粒子进化度和粒子聚合度对惯性权重进行优化设计可以提高粒子群基本算法的稳定性;通过重力坝结构优化,断面面积减小7.1%,具有较好的经济效益,同时安全系数为6.51满足工程安全要求。