岩石力学强度指标的测定是岩石力学与工程研究的基础,岩石力学试验机是测量岩石力学性质最常见的测试仪器。从系统原理来看,岩石压力、拉力试验机是典型的液压系统,试验机已经从原来的人工手动调节阀门流量控制油缸动作、输出压力,普遍发展为电液伺服数字控制系统[1]。国外数字控制岩石试验机处于领先地位,以美国MTS公司、英国GCTS公司、法国TOP公司为代表占据了大量市场,近年来,国内试验机技术进展迅速,东北大学深部金属矿山安全开采教育部重点实验室联合长春朝阳、科意、中机试验等厂商研发了具有自主知识产权的国际领先系列真三轴设备[2,3]。目前国内对控制精度要求高的试验机大多采用德国DOLI公司的EDC控制器或美国DELTA公司的RMC控制器,控制器厂商设计的PID算法内嵌在控制器中。国外控制器性能优异,但价格昂贵,且控制技术封闭。国内试验机控制器需要软件和硬件的双重突破。

试验机控制器的PID参数通常由技术人员在出厂时设置好就不再更改,然而岩石力学实验是一个复杂的过程,不同种类岩石强度差异甚远,同一岩石试件在加载过程的不同阶段也呈现出不同的性质[4],另外还会出现设备维修老化等造成的实验环境变化。岩石变形控制对控制精度要求极高,需要变换PID参数以满足多变的应用环境,然而试验机PID参数调节过程复杂,即使有经验的技术人员也需要反复调节很长时间。当实验环境发生变化,旧的参数不满足控制需求时,只能耗费大量时间和财力通过仪器生产厂家再次调节。使用一种具有自适应功能的控制算法及时调节控制策略,满足不同条件下的控制需求是必要的。GPC广义预测控制是在自适应控制方法的研究中发展起来的一类预测控制算法。与传统自校正控制算法相比,GPC采用实时辨识的建模策略,对控制对象的数学模型精度要求更低[5],GPC已经在多种复杂工业过程获得了应用[6,7]。GPC实时辨识岩石试验机加载过程模型参数,根据模型预测值和未来期望值求算控制量,可以有效规避繁杂的PID参数调节过程并智能适应实验条件变化。

岩石压力试验机是一个液压控制系统,包含液压油源、控制器和作动器。液压油源输出动力油源,伺服阀通过流量的方向和大小调节作动器的位移速度,作动器通过压头加载在岩样上,使岩样产生变形。岩石试验机不同于金属万能试验机[8],其吨位大、试样刚度变化大、控制精度要求更高,根据不同科研需要,加载过程复杂多变,尤其是测定岩石峰后承载能力时,是对非线性时变系统的高精度控制。

图1岩石试验机控制系统框图

图1为试验机伺服控制系统的框图。控制器接收期望变形值ysp和变形传感器反馈的实际值y的误差e,经控制算法计算后输出控制电压u给伺服阀,进行方向和流量调节。

图2为岩样加载初始状态图,变形传感器为LVDT差动位移传感器,作动器在试验机最底部,加载时带动底部垫块、岩样和球头向上移动,球头与上压头接触,最终使岩样产生变形。岩样加载过程存在多个非线性环节,一方面是由于试验机本身造成的,如由电液伺服阀和系统摩擦产生的非线性环节[9,10],作动器与压头刚接触时的受力状态变化等;另一方面是由于岩样本身性质造成的,同种岩石从施加应力到破坏完毕大多会呈现出弹性阶段、塑性阶段、脆性阶段以及峰后阶段等多种性质变化[11],存在严重的时变特性,不仅如此,不同岩石加载过程中所呈现的力学性质也存在巨大差异,且实验过程常要求岩石恒速率变形,对控制精度要求极高,这都为试验机控制算法带来了巨大的挑战。

图2岩样加载初始状态图

1、岩石试验机变形控制模型辨识

系统存在诸多不确定非线性环节,为简化问题,将系统辨识为空载和加载时的2个线性模型,辨识方法采用最小二乘法。

考虑式(1)单输入输出系统:

A(z−1)y(k)=B(z−1)u(t−nb)+e(k){A(z−1)=1+a1z−1+a2z−2+⋯+anaz−naB(z−1)=1+b1z−1+b2z−2+⋯+bnbz−nb   (1)A(z-1)y(k)=B(z-1)u(t-nb)+e(k){A(z-1)=1+a1z-1+a2z-2+⋯+anaz-naB(z-1)=1+b1z-1+b2z-2+⋯+bnbz-nb   (1)

式中,na——输出阶次

nb——输入阶次

y(k)——系统输出

u(k)——系统输入

e(k)——白噪声输入信号

z-1——后移算子

记:

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪h(k)=[−y(k−1),−y(k−2),⋯,−y(k−na),    u(k−1),u(k−2),⋯,u(k−nb)]Tθ=[a1,a2,⋯,an,b1,b2,⋯,bn]Tψ(k)=[v(t−1),v(t−2),v(t−na),0,⋯,0]T{h(k)=[-y(k-1),-y(k-2),⋯,-y(k-na),    u(k-1),u(k-2),⋯,u(k-nb)]Τθ=[a1,a2,⋯,an,b1,b2,⋯,bn]Τψ(k)=[v(t-1),v(t-2),v(t-na),0,⋯,0]Τ

可得式(2):

y(k)=hT(k)θ+e(k)e(k)=ψT(k)θ+v(k)         (2)y(k)=hΤ(k)θ+e(k)e(k)=ψΤ(k)θ+v(k)         (2)

最小二乘的辨识的准则函数为:

J(θ)=∑k=1∞[e(k)]2=∑k=1∞[y(k)−hT(k)θ]2 (3)J(θ)=∑k=1∞[e(k)]2=∑k=1∞[y(k)-hΤ(k)θ]2 (3)

对准则函数J(θ)求导数,并令其一阶导数为0,能够得到参数θ的估计值θˆθ^,通常θˆθ^称作θ的最小二乘估计。为减小单次计算量,加强新采集数据在系统辨识中所占权重,防止数据饱和,实时控制系统编程时常采用式(4),带遗忘因子的递推最小二乘形式[12,13]:

⎧⎩⎨⎪⎪⎪⎪θˆ(k)=θˆ(k−1)+K(k)[y(k)−hT(k)θˆ(k−1)]K(k)=P(k−1)h(k)[hT(k)P(k−1)h(k)+μ]−1P(k)=1μ[I−K(k)hT(k)]P(k−1)         (4){θ^(k)=θ^(k-1)+Κ(k)[y(k)-hΤ(k)θ^(k-1)]Κ(k)=Ρ(k-1)h(k)[hΤ(k)Ρ(k-1)h(k)+μ]-1Ρ(k)=1μ[Ι-Κ(k)hΤ(k)]Ρ(k-1)         (4)

式中,μ为遗忘因子,0<μ<1,一般取0.95～0.98之间;μ为1时是标准递推最小二乘形式;P(k)为正定的协方差矩阵。

分别在空载状态和加载高刚度垫块的情况下,使伺服阀控制电压均匀变化,记录变形传感器的读数。假定不同系统阶次,选取部分输入输出数据进行参数辨识和误差检验,另选取部分数据对辨识结果进行验证,使用均方误差值来表征模型辨识准确度。

空载试验下,辨识均方误差随阶次变化如图3所示;垫块试验下,辨识均方误差随阶次变化如图4所示。

图3系统阶次—均方误差图(空载)

从空载实验和垫块加载实验的辨识结果来看,系统取4阶和5阶时均方误差较小,总是小于取更低阶次时的误差,加载试验系统取5阶时误差比4阶时更大。空载实验的辨识误差和检验误差都小于垫块加载实验的辨识误差,系统取4阶时,垫块加载实验检验误差约是空载实验的15倍,说明垫块加载的过程相对复杂,开环加载实验使用的是质地均匀的实心铁块,非均质岩石加载过程模型会更加复杂,非线性和时变特征更强。最终,取均方差较小的阶次,认为岩石验机变形控制系统为4阶系统。

图4系统阶次—均方误差图(加载)

2、GPC仿真实验

2.1 GPC预测控制

在GPC预测控制策略中,即使被控对象具有非线性或时变特性,依旧认为预测模型可以采用局部线性化和实时辨识的方式表征被控对象。被控对象当前时刻模型可通过式(4)遗忘因子最小二乘辨识。系统当前时刻的辨识结果形式符合式(1)差分方程形式。

图5GPC预测控制结构框图

图5为GPC预测控制结构框图,其中,ysp为期望输出,yr为参考轨迹,u为输入,y为实际输出值,θˆθ^为反馈校正后模型系数,ym为模型多步预测输出。为了利用辨识模型式(1)导出j步后输出y(k+j|k)y(k+j|k)的预测值,引入式(5)关于z-1的多项式方程:

1=Ej(z−1)A(z−1)Δ+z−jFj(z−1)         (5)1=Ej(z-1)A(z-1)Δ+z-jFj(z-1)         (5)

用EjΔzj乘以式(1)两端,并忽略未来时刻噪声,得到k时刻为起点的未来j步的最合适预测输出值表达式:

yˆ(k+j|k)=Gj(z−1)Δu(k+j−1|k)+Fj(z−1)y(k)         (6)y^(k+j|k)=Gj(z-1)Δu(k+j-1|k)+Fj(z-1)y(k)         (6)

式中,Gj(z-1)=Ej(z-1)B(z-1)。结合式(1)可得:

Gj(z−1)=B(z−1)A(z−1)Δ[1−z−jFj(z−1)]=gj,0+gj,1z−1+⋯         (7)Gj(z-1)=B(z-1)A(z-1)Δ[1-z-jFj(z-1)]=gj,0+gj,1z-1+⋯         (7)

因此,多项式Gj(z-1)中前j项的系数与对象阶跃响应前j项的采样值一一对应,记做g1,…,gj,对于j=1,2,…,N,式(8)可以写成矩阵方程:

yˆ=GΔu+f         (8)y^=GΔu+f         (8)

式中,yˆ=(yˆ(k+1),yˆ(k+2),⋯,yˆ(k+N)),Ny^=(y^(k+1),y^(k+2),⋯,y^(k+Ν)),Ν为建模时域;Δu=(Δu(k),Δu(k),Δu(k+M-1)),M为控制时域;f=(f1(k),f2(k),…,fN(k))T,fi(k)=zi-1[Gi(z-1)-z-(i-1)gi,i-1,…,gi,0]Δu(k)+Fi(z-1)y(k),i=1,2,…,N)

G=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪g1g2⋮gp⋮gN0g1⋮gp−1⋮gN−1⋯⋯⋱⋯⋱⋯00⋮g1⋮gN−M+1⎫⎭⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪G={g10⋯0g2g1⋯0⋮⋮⋱⋮gpgp-1⋯g1⋮⋮⋱⋮gΝgΝ-1⋯gΝ-Μ+1}

其中多项式E,F,f可通过递推求解丢番图得到,参考文献[14,15]。式(9)为GPC在k时刻优化性能指标。

J(k)=E⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪∑j=N1N2[y(k+j|k)−yr(k+j)]2+∑j=1Mλ(j)[Δu(k+j−1|k)]2⎫⎭⎬⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪  (9)J(k)=E{∑j=Ν1Ν2[y(k+j|k)-yr(k+j)]2+∑j=1Μλ(j)[Δu(k+j-1|k)]2}  (9)

式中,E(·)——{·}的数学期望

yr——系统的期望输出

N1——优化时域起始时刻

N2——优化时域终止时刻

M——控制时域,在M步后控制增量0

λ(j)——控制指标的权系数

式(9)中系统输出值用预测值式(8)代替,则有当λI+GTG非奇异时,得到使性能指标式(9)最优的解为:

Δu(k)=(λI+GTG)−1GT(yr(k)−f(k))  (10) Δu(k)=(λΙ+GΤG)-1GΤ(yr(k)-f(k))  (10)

此时最优控制量则可由式(11)给出:

u(k)=u(k−1)+dT(yr(k)−f(k))         (11)u(k)=u(k-1)+dΤ(yr(k)-f(k))         (11)

式中,dT为矩阵(λI+GTG)-1GT的第一行。

2.2 仿真实验

考虑到GPC控制算法中存在大量矩阵运算,并且模型辨识环节实时进行,可以对控制对象局部线性化,因此在控制算法中使用2阶模型在线辨识、预测输出并计算控制量来降低计算量。

岩石变形实验多为不同速率下匀速加载过程,因此仿真实验期望输出采用斜坡函数。GPC仿真效果图如图6所示,系统初始变形为1mm,在第2～10秒,10～14秒,14～22秒的各时段之间采用不同速率匀速加载,分别为0.005,0.003,0.001mm/s,期望变形值在22s后保持2.5mm不变。

图6GPC仿真过程曲线

根据系统在第1章模型辨识结果,对2个4阶被控对象模型仿真实验。仿真开始阶段,被控对象使用的是空载辨识模型,系统输出跟随性很好,变形速率与期望速率吻合。在第6秒,被控对象由空载模型突变为垫块加载模型,此时控制算法模型辨识参数处于调整阶段,故在6～7秒时,虽然期望速率未发生变化,但原来已经稳定的跟踪误差发生了变化,这是控制器适应被控对象模型变化的过程。

对于仿真过程,空载模型和垫块加载模型使用了完全相同的控制参数,包括预测时域、控制时域、遗忘因子和性能指标均衡参数,体现了控制算法的自适应能力。GPC参数调节过程十分简单,根据控制效果简单调整后即可获得较为理想的控制效果。从整体控制效果来看,基于GPC的变形控制系统具有很好的动态性能和稳态性能,说明2阶模型不仅当前时刻可以很好地表征被控对象,在预测时域的局部范围内对被控对象的预测输出偏差也是满足控制需求的。

3、实验结果

通过前两章的系统辨识与仿真,将GPC应用到试验机变形控制的可行性已经得到验证,并且控制结构和相关参数已经基本确定。使用CSharp语言编写控制算法,并开展岩石力学单轴实验,验证实际控制效果。软件编程完毕后,采用高刚度垫块的加载实验进行设备调试,得到如图7所示仿真结果。初始期望加载速率为0.01mm/s,此时试验机上压头还未与球头接触,直到第13.52秒时,垫块开始受力,实验力显著增加,但是此过程参数变化平滑,且变化幅值较小,在控制量的及时调整下,系统输出依旧保持了预定加载速率。随后分别进行了0.004,0.001mm/s和-0.001,-0.005mm/s的加载、保载和卸载试验。图7b显示跟踪误差在速率变换时较大,但也在保持在微米级,且控制器调节迅速,整体上取得了满意的控制效果。

图7GPC岩石试验机实验过程曲线(铁块)

从图8辨识模型参数变化过程来看,辨识模型参数实时变化,并与加载速率呈现出一定的相关性,变形速率为0时,模型参数幅值变化较小,趋于稳定。

图8GPC岩石试验机模型参数变化曲线(铁块)

设备调试完毕后,使用设备开展岩石力学实验,测试红砂岩和花岗岩的单轴抗压强度。图9为花岗岩强度测试过程曲线。花岗岩恒速率变形加载实验使用0.001mm/s的变形速率进行加载,并在实验过程中保载了117s。在GPC控制策略作用下,岩石变形量与期望值基本重合,控制效果较好。岩石在1113秒时进入峰后破坏阶段,强度变小,内部开裂严重[16],导致在峰后阶段出现了2处较大的跟踪误差,GPC及时调整了控制量,能够在岩石峰后较弱的承载能力下继续按期望加载而不使岩石发生立即完全破坏。

图9GPC岩石试验机实验过程曲线(花岗岩)

GPC辨识的模型参数变化如图10所示,花岗岩加载过程中模型参数实时调整,在恒加载速率下,参数变化幅值较小,加载速率为0时,参数趋于稳定。最终测得0.001mm/s变形速率下,花岗岩试件实验力峰值为114.76kN,抗压强度为73.76MPa,同时得到了岩样峰后压力和变形值,取得了十分满意的实验效果。

图10GPC岩石试验机模型参数变化曲线(花岗岩)

图11为砂岩测试的实验过程曲线,岩石变形按照设定的0.004mm/s速率稳定加载。由于砂岩样品内部存在裂隙,在58.52秒时引起了压力值波动,出现了允许范围内的较大跟踪误差,控制器快速调整了控制策略,继续按期望加载,完成了砂岩强度测试。图12为砂岩实验过程模型参数变化曲线,模型参数整体稳定,在58.52秒时有小幅度突变。砂岩样品试验力峰值为127.25kN,强度为64.84MPa。

图11GPC岩石试验机实验过程曲线(砂岩)

图12GPC岩石试验机模型参数变化曲线(砂岩)

图13为岩石破坏后形态,可以看出样品中存在较大节理,出现提前脱落现象,极大降低了岩石强度[17]。2种岩石破碎效果差异明显,花岗岩内部存在更多节理裂隙,破碎后存在大量块状形态,砂岩均质性更好,多以粉末状破碎。虽然岩石性质差异巨大,但没有影响控制效果,验证了GPC对于不同种类岩石较好的适应性。

图13岩石破坏后形态

4、结论

通过对试验机电液伺服控制系统的模型辨识和GPC广义预测控制策略的研究,确认了系统阶次,明确了算法可行性。岩石轴向变形控制系统是一个4阶系统,可以采用局部线性化的方式在GPC中使用2阶模型实时辨识,减小控制算法计算量。

用CSharp语言编写的GPC控制策略可以应用于岩石试验机电液伺服系统,完成不同速率下的岩石力学恒速率加载实验。GPC广义预测控制,参数调节过程简单,能够有效规避PID参数调节的复杂过程,可针对不同岩石种类和岩石加载过程的非线性、时变性质实时调整控制策略,自适应和智能化程度高,具有较好的应用推广价值。

参考文献：

[1]牛学超,张庆喜,岳中文.岩石三轴试验机的现状及发展趋势[J].岩土力学,2013,34(2):

[2]张希巍,冯夏庭,孔瑞,等.硬岩应力-应变曲线真三轴仪研制关键技术研究[J].岩石力学与工程学报,2017,36(11):2629-2640.

[4]付小敏,邓荣贵.室内岩石力学试验[M].成都:西南交通大学出版社,2012.

[5]关守平,尤富强,徐林,等.计算机控制理论与设计[M].北京:机械工业出版社,2012.

[6]何伟,鲁明,李国强,等.SNCR脱硝系统的广义预测控制[J].石油化工自动化,2016,52(1):38-41.

[7]马草原,朱信尚,韩永刚,等.基于PSO的自适应广义预测微燃机控制[J].控制工程,2019,26(2):179-184.

[8]苏世杰,游有鹏,齐继阳,等.电液伺服试验机力控系统负载刚度自适应控制[J].控制理论与应用,2018,35(4):11-19.

[9]毛尾,纪朱珂,韦海利,等.电液比例伺服系统模糊PID复合控制应用研究[J].液压与气动,2019,(1):95-99.

[10]郭庆.电液伺服非线性控制技术研究进展综述[J].液压与气动,2018,(3):1-9.

[11]郭建强,黄质宏.循环荷载作用下岩石疲劳本构模型初探[J].岩土工程学报,2015,(9):156-162.

[12]荀倩,王培良,李祖欣,等.基于递推最小二乘法的永磁伺服系统参数辨识[J].电工技术学报,2016,31(17):161-169.

[13]鲍海静,张韬,张静.基于遗忘因子递推最小二乘法的伺服系统转动惯量辨识方法[J].上海电气技术,2017,(3):9-13,18.

[15]钱积新,赵军,徐祖华.预测控制[M].北京:化学工业出版社,2007.

[16]谢璨,李树忱,平洋,等.峰后裂隙岩石非线性损伤特性与数值模拟研究[J].岩土力学,2017,(7):2128-2136.

[17]胡伟,邬爱清,陈胜宏.含隐裂隙柱状节理玄武岩单轴力学特性研究[J].岩石力学与工程学报,2017,36(8):1880-1888.

高继开,田军,姚志宾,何本国.基于广义预测控制的智能岩石试验机研究[J].液压与气动,2020(11):74-80.

基金:国家自然科学基金(51809038).