Bouc-Wen 模型参数识别的非线性自适应遗传算法和试验验证

章红梅，胡 帆，段元锋

(浙江大学建筑工程学院，浙江，杭州 310058)

结构动力行为中滞回模型是表征结构动力响应行为的重要规则[1]。Bouc-Wen 模型是一种通用的非线性光滑滞回模型，可以表达结构构件滞回特征中的刚度退化、强度退化等多种力学特征，在表征和模拟结构复杂滞回行为上具有强大的适应性[2]。然而，Bouc-Wen 模型本身参数众多，且其中的多数参数不具有明确的物理意义，这使得获得恰当的参数组较为困难，Bouc-Wen 模型的应用因此也受到极大的阻碍。针对该问题，Ismail等[3]直接用解析的方法对模型进行参数求解，但是其中的大量假设，使得所得到的Bouc-Wen 模型失去了广泛性。相比解析方法，利用数值方法从试验中识别Bouc-Wen 模型参数具有更广泛的适用性。

目前，从已知试验数据中识别Bouc-Wen 模型参数的算法研究大致可以分为两类：一类为传统的数值算法，如Kalman 滤波器法[4-6]、最小二乘法[7]、梯度下降法[8]等，对于模型的参数识别效果不算优异；另一类是各种智能算法，如遗传算法[9]、粒子群算法[10]、微分进化算法[11]、GSO 算法[12]等，该类算法虽然可以达到较好的效果，但是实现过程较复杂，还有巨大的改进空间。

其中应用较多的遗传算法是1975 年由Holland[13]教授首次提出的，目前已广泛应用于遗传合成、VLSI 技术、土木工程和机器学习。刘永强等[14]采用遗传算法对MR 阻尼器进行了参数识别，通过缩小参数取值范围的方法取得了较好的拟合效果；李得民等[15]通过比较由遗传算法识别MRE隔振器的Bouc-Wen 模型和遗传算法优化对MRE力学行为仿真建模的BP 神经网络两者的结果，验证了两种模型的正确性；Negash 等[16]提出了一种针对环节优化的新型遗传算法。然而，现阶段应用遗传算法对Bouc-Wen 模型进行参数识别仍然有许多不足，主要在于由局部收敛导致的拟合精度不足和由耗时较长导致的实用性下降。

本文提出了一种非线性自适应遗传算法(NAGA)，通过考虑染色体库策略、最优保存策略、改良的适应度函数和选择算子以及改进的自适应交叉算子和变异算子，对变换轴压比的钢筋混凝土剪力墙低周反复加载试验进行Bouc-Wen 参数识别并和标准遗传算法进行对比，结果表明该方法显著提高了识别的有效性。

1 遗传算法识别模型

Bouc[17]提出了一种多变量、多用途的光滑滞回模型，Wen[18]在1976 年对其进行了改进，称为Bouc-Wen 模型，并首先用于模拟单自由度系统的非线性滞回特性。如图1 所示的单自由度系统的运动方程可以用式(1)表示。

图1 单自由度Bouc-Wen 模型Fig. 1 Single-degree-of-freedom Bouc-Wen model

式中：u(t)为质量块m相对于地面的位移；c为线性粘滞阻尼系数；F(t) 为 外部激励荷载；z(t)为滞回位移； α为每个滞回环峰值处切线刚度与初始刚度之比，其范围为0～1；k为初始刚度。

Bouc-Wen 模型采用式(2)来描述滞回位移。式中，n是控制滞回曲线屈服的尖锐程度参数；A、β、γ三个滞回参数决定了滞回环的基本形状。

2 标准遗传算法策略设置

标准遗传算法(SGA)是由Goldberg[22]在研究控制优化应用中提出的，这为研究人员应用遗传算法解决各类实际问题奠定了基础。本文先按照SGA 流程依次进行策略设置，然后，针对其中出现的搜索方向不稳定、收敛效率低的问题进行改进，提出一种识别效果显著改善的识别方法(NAGA)。主要算法策略介绍如下。

2.1 编码策略

本文采用二进制编码，使用基于字符集{0,1}构成的基因串表示Bouc-Wen 模型中的实际参数，下文中由MATLAB 程序所计算出的值均为二进制数转化而来的实数值，不代表实际控制精度，同时令 β/|γ|≥1[19-20]。Bouc-Wen 模型中，[α0,k,γ,β,δν,δη,A0,δA]这8 个参数都是有一定范围约束的，并且部分参数不具备明确的物理意义，具体取值也难以标准化，Ma 等[23]用理论推导的方式证明了这8 个参数是有冗余的，因此部分参数的范围应该是动态的。如表1 所示，本文半定量半定性地选取了较大的解空间，因此，可以在广泛的8 维参数空间中求解最优解。本文的染色体是这8 个参数的组合，根据每个参数的上下限和精度的要求，本文中染色体长度设置为121，各参数具体分布见图2。

表1 Bouc-Wen 模型中各参数的上下限和精度Table 1 Limit and precision of Bouc-Wen model parameters

图2 染色体上各参数对应基因的分布Fig. 2 Distribution of genes corresponding to parameters on chromosome

2.2 初始化种群

由于对于Bouc-Wen 模型的最优解数量及在解空间的分布不明确，为了有效地在实际参数空间中均匀选取样本，本文利用MATLAB 自带的随机函数rand，在染色体基因位上随机生成0 或1，形成一个初始种群。种群规模N的大小对遗传算法的搜索质量有着非常大的影响，一般建议值N[24]为20～100。考虑到整个参数空间比较大，达到2121的规模，为了避免算法发生“早熟”现象，本文的种群规模设定为500。

2.3 适应度函数

对于规模为N 的种群P={I1,I2,···,IN}，染色体Ii∈P的目标函数值E(Ii)可以通过式(9)计算，这个函数通过计算试验数据和Bouc-Wen 模型拟合数据的均方根误差(RMSE)值来评估该条染色体的优劣程度。适应度函数可以通过式(10)进行计算。

2.4 选择算子

2.5 交叉算子和变异算子

本文采用均匀变异策略。作为能够增加种群多样性的有效手段，变异概率直接决定算法收敛速度。变异概率如果设置的太小，算法可能陷入早熟，如果过高的话，算法收敛速度过慢，难以生成一个稳定的优解。一般取pm=0.001 ～0.05[25]，本文采用Srinivas[25]推荐的0.005。

3 非线性自适应遗传算法的提出

在上述SGA 流程的基础上，本文提出染色体库策略、最优保存策略，设计了一种双选择算子，引入自适应交叉变异策略和有界多点交叉方法来对其进行改进。方法介绍如下。

3.1 染色体库策略

本文种群规模N=500，遗传代数为G=200。在本文第4 节中4 片试验剪力墙在3 次SGA 过程中的平均染色体重复率见图3，由此可以看出重复出现的个体占30%左右，甚至还有重复出现14 次的个体，如果能够记录重复个体的适应度，整个优化过程将会节省大量算力。

图3 3 次SGA 过程中各片剪力墙的染色体平均重复率Fig. 3 Chromosome average repetition rate of each shear wall in three SGA processes

本文采用牛顿-拉夫逊迭代求解Bouc-Wen 模型，要提高算法运算效率，就需要尽量减少耗时占总时长99%以上的数值求解函数的调用次数，由此提出染色体库策略。染色体库是将运算过程中出现过的染色体储存到一个数组中，同时也将该染色体对应的RMSE(E)值储存到同一个数组中，这样形成了一个对于单个问题的染色体库。在每次求解Bouc-Wen 方程之前，先查询染色体库，如果有相同染色体的结果，直接调用；如果没有，则进行计算，同时将这些还未收录的染色体连带其RMSE(E)值收录进染色体库。若染色体数量达到染色体库可收录上限，则通过比较替换染色体库中RMSE(E)值最大的，也就是适应度最差的染色体。这是由于随着遗传进化的推进，种群逐步收敛到适应度较高的个体上，那么这些适应度高的染色体出现地概率更大，如此染色体库能够得到更有效地调用。染色体库大小可设计为种群规模的5 倍～10 倍，本文中由于N和G值较大，因此取10。

3.2 最优保存策略

Rudolph[26]通过齐次有限马尔科夫链分析证明，由于SGA 的不可化归性，无论如何初始化，目标函数、交叉变异算子如何选取，其都不会收敛到全局最优。因此，本文采用Negash 等[16]使用的最优保存策略，即在每一代中选取表现最好的染色体直接复制到下一代种群中，以此来保证每代中的最优个体不会在遗传过程中遭到破坏。同时，Rudolph[26]证明了采用此策略后算法可以收敛到全局最优解。

3.3 适应度函数改进

图4 SGA 进化过程中最差染色体RMSE(E)的常用对数Fig. 4 Common logarithm of RMSE(E) of worst chromosome during SGA evolution process

为了研究a值对变换后适应度值离散度的影响并进一步确定a的值，本文选择剪力墙SW1-1在SGA 过程中第200 代染色体的适应度值作为分析对象。因为当a的值不同时，变换后的适应度的量级通常是不同的，所以不可能直接使用标准差来进行判断。本文提出了离散度的概念来进行判断，见式(14)。

SGA 进化过程后期的种群，染色体适应度值较为集中，因此 Υ较小。由表2 可以看出，a值从0.5 到1.0 时，种群的 Υ增加，而当a从1.0 到2.5 时，种群 Υ逐渐减小。当取1.0 时，分析对象的离散度是最大的，因此，本文中a取1.0。

表2 不同a 值对应的离散度Table 2 Measure of dispersion at different values of a

3.4 选择算子改进

SGA 进化后期，由于染色体的适应度都非常接近，每个染色体被选择的概率也相近，这会使得轮盘赌选择过程变成随机选择过程，如此算法整体收敛性会大打折扣。因此本文提出一种双选择算子对此进行改进。在遗传进化前中期，使用轮盘赌策略在种群中保持较小的选择压力，使得各种基因模式都能得到保留；而到进化后期，为加快算法收敛速度，提高选择压力，采用基于排序的确定性选择算子。双选择算子的流程图见图5。引入遗传代数g作为控制变量，当g≤0.6G时，采用轮盘赌算子；当g＞0.6G时，采用基于排序的确定性选择算子，即根据适应度对每一代的染色体进行排序，将排在前25%的染色体复制两次到下一代，排在中间50%的个体直接复制到下一代，排在末尾25%的个体直接淘汰。这样后期通过确定性的筛选，种群会快速收敛。

图5 双选择算子Fig. 5 Double-choice operator

3.5 交叉算子和变异算子改进

遗传算法的全局搜索能力由交叉算子控制，而局部搜索能力由变异算子控制。对于SGA 来说，其交叉概率和变异概率都是恒定的，交叉概率和变异概率的合理设定一般需要通过大量的试验，而且很难找到一个放在遗传进化各个阶段都合适的常数。同时固定的交叉概率和变异概率使得算法在遗传进化后期，有一定可能性破坏表现较好的染色体，引起种群的“退化”现象。交叉概率和变异概率如果能够依照每个染色体自身适应度值进行调整，使得较高适应度个体的基因组合能够有较大概率保留，而较低适应度个体能够通过大交叉概率和大变异概率突破现有的种群模式，这种方法使得整个算法能够跳出局部收敛，将充分改善算法的运行效率和鲁棒性。

Srinvivas 等[25]首先提出自适应遗传算法，其交叉概率和变异概率只考虑了与适应度的关系，且是线性关系。本文在此基础上引入遗传代数和种群规模作为控制变量来调整交叉概率和变异概率。在遗传初期，交叉变异概率维持在基本水平，到进化后期，抑制表现较好的个体交叉变异，加剧表现较差的个体交叉变异，如此在加快算法在当前最优解域收敛的同时，尽可能探索更多的基因模式。交叉概率和变异概率的具体形式见式(15)～式(16)[27]。

考虑到一点交叉有一定概率破坏较好的基因模式，不利于长距离模式的保留和重组。因此本文借鉴建筑模块假说思想[13]，提出一种有界多点交叉方法，定向对每个参数所在编码区域进行交叉操作，如此使得算法搜索具有一定方向性，详细的操作过程见图6。首先对于每个参数生成一个随机数，若该随机数小于交叉概率则该参数对应的基因片段要进行交叉，再根据各随机数确定各交叉点。

图6 有界多点交叉Fig. 6 Bounded multipoint crossover

4 RC 剪力墙低周反复加载试验验证

4.1 试验简介

作为一种重要的竖向承载和抗侧力构件，钢筋混凝土剪力墙被广泛地应用在工业与民用建筑中。但是由于剪力墙在服役过程中受到弯、剪、扭的共同作用，和混凝土材料的复杂非线性行为等，使得剪力墙的滞回特性难以准确预测。本文中用于验证的4 片变换轴压比剪力墙[28]参数见表3。

表3 剪力墙构件参数Table 3 Shear wall component parameters

RC 剪力墙低周反复加载试验见图7。试验采用拟静力方式进行加载。正式试验前，先在墙体顶部施加大小为预定荷载40%的竖向压力，重复加载2 次～3 次，以消除试件内部组织不均匀性，然后将竖向压力加至预定荷载并在试验中保持不变。

图7 RC 剪力墙低周反复加载试验Fig. 7 Low-cycle cyclic loading test of RC shear wall

侧向加载分两个阶段。第一阶段采用力控制，先单调逐级加载至开裂，循环一次，以后各级荷载循环一次，直至屈服；第二阶段采用位移控制，分别按屈服时顶点位移的倍数逐级加载，每级循环三次，直至构件的承载力下降到峰值承载力的80%左右为止。在整个加载的过程中，虽有加载方式的区别，但是无论是力还是位移控制加载，本文都是取用最后测量得到的不同荷载步上的一组力-位移值，大量力-位移值就构成了本文的识别数据组，在这样的数据组空间中，本文采用了位移作为控制变量的目标函数式(9)，值得一提的是，该控制变量与试验加载的控制方式的概念并不相同。SW1-1 的具体侧向加载历程如图8所示。

图8 SW1-1 的侧向加载历程Fig. 8 Specific lateral loading history of SW1-1

4.2 识别结果

本文半定量半定性地给8 个参数选择了较大的解域，SGA 和NAGA 识别得到的多组8 参数未聚集到到所设定范围的边界，说明所选解域已包括全局最优解。NAGA 克服了SGA 易局部收敛的问题，可以不断向全局最优解接近，从而提高识别效果。

4.2.1 参数识别结果

由于遗传算法有一定的随机性，本文对每片墙进行10 次识别，引入变异系数cv来衡量其结果的波动情况，具体表达式见式(17)。表4 为四片剪力墙平均RMSE(E)值μ及其变异系数，总体来说，NAGA 参数识别效果要明显优于SGA，并且变异系数更小，表明算法收敛性更好。RMSE(E)值最小的SGA、NAGA 参数识别结果分别见表5、表6，并将基于该参数组合进行后续分析。

表4 4 片剪力墙参数识别波动情况Table 4 Parameter identification fluctuation of four shear walls

表5 SGA 参数识别结果Table 5 Parameter identification results based on SGA

表6 NAGA 参数识别结果Table 6 Parameter identification results based on NAGA

4.2.2 滞回曲线

由图9 可知，所提出的非线性自适应遗传算法识别效果更好，模拟滞回曲线的饱和度和试验滞回曲线更接近。从4 片剪力墙的滞回曲线可以看到，加载过程中没有出现明显的捏拢现象，说明第二节中对于Bouc-Wen 模型的简化是合理的。

图9 基于SGA 和NAGA 的模拟滞回曲线比较Fig. 9 Comparison of simulated hysteretic curves based on SGA and NAGA

决定系数是用来评价模型的拟合程度的，一般来说决定系数越大，说明模拟的效果越好，其值不大于1，具体可由式(18)计算。由表7 可知改进之后算法识别的效果更好。

表7 SGA 和NAGA 的决定系数Table 7 Coefficient of determination of SGA and NAGA

4.2.3 力与位移响应过程

从图10 可以看到，对于SGA 由于识别出的初始刚度过大使得初期计算出的侧向力较大，而刚度退化和强度退化过大使得在加载中后期侧向力计算值小于试验值；而对于NAGA 来说，虽然由于识别出的初始刚度偏小使得初期模拟曲线的侧向力小于试验力，但当数据量足够多时识别和拟合效果优异。

图10 基于SGA 和NAGA 的力与位移响应过程比较Fig. 10 Comparison of force and displacement response process based on SGA and NAGA

4.2.4 算法效率分析

由图11 可知，引入最优保存策略后，种群的最优个体稳定地朝着使RMSE(E)更小的方向进化，最优个体不会出现波动甚至“退化”现象。NAGA 在100 代左右就收敛到当前最优解，收敛速度较为快速，也表明G的取值是合理的。

图11 SGA 与NAGA 进化过程中最优染色体的RMSE(E)变化对比Fig. 11 Comparison between change of RMSE(E) of best chromosome during evolution of NAGA and SGA

各片剪力墙最优的RMSE(E)及其识别的时间见表8。由此可见，在识别效果方面，NAGA 最优染色体的RMSE(E)比SGA 的值小了40%；在识别速度方面，由于染色体库策略的引入使得数值求解函数的调用次数减少，同时算法收敛性的增强也使得数值求解过程中可以更快地满足各荷载步滞回位移的精度要求，如此NAGA 的速度相较SGA 有了一定的优势。

表8 SGA 和NAGA 实用性比较Table 8 Practicality comparison between SGA with NAGA

5 结论

本文在标准遗传算法的基础上，通过引入染色体库策略、最优保存策略，改良适应度函数和选择算子，并在现有自适应交叉算子和变异算子的基础上进一步改良，提出了一种非线性自适应遗传算法。通过4 组变换轴压比的RC 剪力墙低周反复加载试验数据验证了该方法的有效性，与已有的SGA 方法相比，本方法在Bouc-Wen 模型参数识别上具有明显的优越性。本文主要结论如下：

(1) 本文所提出的非线性自适应遗传算法(NAGA)对于8 参数Bouc-Wen 模型参数的识别效果良好，通过所识别的参数能够较好地模拟剪力墙这种结构构件在低周反复加载中的滞回行为；

(2) 与试验数据相比，在本文识别表征RC 剪力墙滞回行为的Bouc-Wen 模型参数的拟合效果显著优于现有的标准遗传算法(SGA)；

(3) 研究中染色体库策略的引入和算法收敛性的增强使得参数识别速度有了一定提升，算法实用性得到有效增强。

由于遗传算法的随机性，要想通过改变其超参数，用分析滞回曲线、RMSE(E)或时间等定量指标来分析是困难的，目前只能做到定性分析来得到一组合适的遗传算法参数。本文所述方法的提出可用于混合模拟试验的模型参数识别中，通过试验子结构的滞回性能参数实时识别，并及时反馈给计算子结构进行整体性能模拟分析，也同时适用于大型复杂结构在服役期间根据所识别得到的关键结构构件的滞回特性预测整体结构的抗灾能力。