轧制力能参数计算精度直接影响带钢板形和板厚的控制质量，很多专家学者对此进行了广泛的研究。针对神经网络无法提升小批量、多品种轧制模式下的轧制力设定精度的问题，文献[1]采用数据统计分析的方法对预设定模型的修正系数重新计算。文献[2]采用有限元法与人工神经网络相结合的方式进行轧制力预报，精度有所提高且满足在线使用需求。文献[3]提出了一种基于贝叶斯方法的神经网络预测冷轧轧制力，在收敛速度和预测精度上均优于传统的BP神经网络。文献[4-5]采用神经网络对变形抗力和摩擦系数等非线性因素进行系数校正，再结合轧制力数学解析模型，建立了轧制力在线计算模型。

主电机功率是冷连轧生产中需要确定的重要力能参数之一，其主要作用包括：为主电机容量选取提供合理依据；在工艺参数设定时，必须保证主电机的输出功率不超过电机本身容许的最大功率^[6]；在轧制负荷分配时，通常使各机架按照功率成比例或等功率裕度的方式进行分配，以充分发挥各机架主电机的能力^[7-10]。由此可见，提升电机功率模型的计算精度在提高产品控制精度、优化轧制规程和发挥设备能力并保证生产安全等方面具有重要意义。

主电机功率的确定必须先计算传动轧辊的轧制力矩、附加摩擦力矩、空转力矩和动力矩，其中轧制力矩是最大和有效部分，剩余部分为有害力矩^[11]。轧制力矩理论模型主要有2种：传统模型和上界法。由于传统模型形式简单、物理意义清晰，因此工程上经常使用^[12]。上界法是以上界定理为依据，通过假定满足运动许可条件的速度场，计算上界功率并对其求得最小值，最后解得力能和变形参数的方法^[13]。此方法明确了各参数间的物理关系，且力能参数结果往往比实际值略高，对设备的选择和保护十分有利。近些年不少专家学者从速度场设定角度出发进行研究，如抛物线速度场^[14]、余弦速度场^[15]、正切速度场^[16]以及流函数法^[17]。采用上述2种方法计算主电机功率都不可避免的需要确定摩擦系数、力臂系数以及损失力矩，但这些模型参数常常难以确定^[18-19]。例如：力臂系数影响因素众多，目前还没有公认的解析公式^[20]；损失力矩会受到轧制速度和轧制力等参数变化的影响；摩擦系数会随着轧制长度的增大而减小。因此，这些参数对轧制力矩和电机功率的计算精度有很大影响。

鉴于上述问题，本文提出一种机理模型与数据驱动相结合的冷轧电机功率预测模型。其中，机理模型部分采用上界法得到解析解，然后将解析解与相关变量输入长短时记忆(long short-term memory, LSTM)网络中训练，以提取变量数据之间的深度特征和时序相关性。实验验证表明：所提出的时序耦合预测模型与上界法或不考虑数据时序特性的人工神经网络相比预测精度均有所提高，相对误差分别从不超过±8.5%和±6.3%提升到±3.9%；而且因有传统机理模型参与，预测结果更符合物理规律。

1 轧制力矩解析模型

冷轧过程带钢宽厚比远大于10，因此可以忽略宽展并认为带钢在变形区处于平面应变状态。再结合沿轧制线的对称性，便可得到其1/2变形区示意见图 1。

图 1中，以入口处和轧制线的交点为原点建立坐标系，H和h为带钢入口和出口厚度，v_H和v_h为带钢入口和出口速度，h_x为变形区任意位置带钢厚度，v_r为轧制速度，α为接触角，θ为咬入角，R为工作辊半径，l为变形区长度。根据图中几何关系，接触弧方程及其一阶导数可以表示为:

$ h_x=h+2 R-2 \sqrt{R^2-(l-x)^2}=h+2 R-2 R \cos \alpha $

(1)

$ h_x^{\prime}=-\frac{2(l-x)}{\sqrt{R^2-(l-x)^2}}=-2 \tan \alpha $

(2)

1.1 二维速度场与应变速率场

考虑到冷连轧过程中宽展量可以忽略，因此认为轧制过程是平面应变问题，则ε_z=0；采用平断面假设，即任意x截面上的水平速度相同，并假设截面上的垂直速度呈线性分布，得到如下二维速度场：

$ \left\{\begin{array}{l} v_x=\frac{U}{h_x b} \\ v_y=\frac{U}{b} \frac{h_x^{\prime}}{h_x^2} y=v_x \frac{h_x^{\prime}}{h_x} y \\ v_z=0 \end{array}\right. $

(3)

式中：$U=v_{\mathrm{H}} H b=v_{\mathrm{h}} h b=v_{\mathrm{n}} h_{\mathrm{n}} b=v_{\mathrm{r}} \cos \alpha_{\mathrm{n}} b(h+2 R-\left.2 R \cos \alpha_{\mathrm{n}}\right)$，h_n为中性面厚度，α_n为中性角。

由Cauchy方程可知，对应的应变速率场表示为

$ \left\{\begin{array}{l} \dot{\varepsilon}_x=\frac{\partial v_x}{\partial x}=-v_x \frac{h_x^{\prime}}{h_x} \\ \dot{\varepsilon}_y=\frac{\partial v_y}{\partial y}=v_x \frac{h_x^{\prime}}{h_x} \\ \dot{\varepsilon}_z=0 \end{array}\right. $

(4)

由式(3)、(4)可得：$\dot{\varepsilon}_x+\dot{\varepsilon}_y+\dot{\varepsilon}_z=0$；当x=0时，v_x=v_H；当y=0时，v_y=0；在轧辊与轧件的接触面处v_y=－v_xtan α。因此速度场满足体积不变条件和速度边界条件。

1.2 冷轧总功率计算 1.2.1 塑性变形功率

由于Mises屈服准则是非线性的，会提高积分计算难度，因此本文采用内接十二边形(inscribed dodecagon, ID)线性屈服准则计算塑性变形功率^[21]。由式(4)可知，$\dot{\varepsilon}_{\max }=\dot{\varepsilon}_x, \dot{\varepsilon}_{\min }=\dot{\varepsilon}_y$，则内部塑性变形功率可表示为

$ \begin{aligned} \dot{W}_i= & \int_V D\left(\dot{\varepsilon}_{i j}\right) \mathrm{d} V= \\ & \frac{2 b}{\sqrt{3}} \sigma_s \int_0^l \int_0^{h_x / 2}\left(\dot{\varepsilon}_{\max }-\dot{\varepsilon}_{\min }\right) \mathrm{d} x \mathrm{~d} y=\frac{2 U}{\sqrt{3}} \sigma_\mathrm{s} \ln \left(\frac{H}{h}\right) \end{aligned} $

(5)

1.2.2 摩擦功率

摩擦功率可由轧辊和带钢表面沿x和y方向的相对滑动速度与摩擦应力乘积并积分得到，摩擦应力采用摩擦因子模型计算。需要注意的是，前后滑区相对滑动速度方向发生改变，因此积分过程中要以中性面为界进行分段积分。令A₁=1+h/2R，则摩擦功率可表示为

$ \begin{aligned} \dot{W}_{\mathrm{f}}= & 2 m k b\left(\int_0^l\left|v_{\mathrm{r}} \cos \alpha-v_x\right| \mathrm{d} x+\int_0^l\left|v_{\mathrm{r}} \sin \alpha-v_x \tan \alpha\right| \tan \alpha \mathrm{d}_x\right)= \\ & 2mkb\left\{ {R{v_{\rm{r}}}\left( {\theta - 2{\alpha _{\rm{n}}}} \right) - \frac{U}{{2b}}\frac{{\sec \bar \alpha }}{{\sqrt {{A_1} - 1} }}} \right.\\ & \left. {\left[ {\arctan \left( {\sqrt {\frac{{{A_1} + 1}}{{{A_1} - 1}}} \tan \frac{\theta }{2}} \right) - 2\arctan \left( {\sqrt {\frac{{{A_1} + 1}}{{{A_1} - 1}}} \tan \frac{{{\alpha _{\rm{n}}}}}{2}} \right)} \right]} \right\} \end{aligned} $

(6)

式中：k为带钢的剪切强度，m为摩擦因子。

1.2.3 剪切功率

由式(3)可知，在变形区入口处$v_y=-v_{\mathrm{H}} \frac{2 \tan \theta}{H} y$，垂直速度不为0，入口处存在速度不连续现象；在变形区出口处v_y=0，出口处垂直速度为0，不存在速度不连续面。因此，仅在入口处有由于速度不连续导致的剪切功率，其表达式为

$ \dot{W}_{\mathrm{s}}=2 b \int_0^{H / 2} k\left|\Delta v_y\right| \mathrm{d} y=\frac{U k}{2} \tan \theta $

(7)

1.2.4 张力功率

在冷连轧生产中，为了防止带钢跑偏、降低轧件变形抗力和提高生产效率，常常在带钢两端施加张力。张力改变会引起变形区各参数的变化，因此需要额外考虑张力引发的功率变化，其表达式如下：

$ \dot{W}_{\mathrm{t}}=2 b\left(\int_0^{H / 2} \sigma_{\mathrm{b}} v_{\mathrm{H}} \mathrm{d} y-\int_0^{h / 2} \sigma_{\mathrm{f}} v_{\mathrm{h}} \mathrm{d} y\right)=U\left(\sigma_{\mathrm{b}}-\sigma_{\mathrm{f}}\right) $

(8)

式中：σ_f为前张力，σ_b为后张力。

1.2.5 总功率泛函最小化

将式(5)~(8)代入$\varPhi=\dot{W}_{\mathrm{i}}+\dot{W}_{\mathrm{f}}+\dot{W}_{\mathrm{s}}+\dot{W}_{\mathrm{t}}$即可得到总功率泛函。令总功率泛函对中性角α_n求偏导并等于0。

$ \frac{\partial \varPhi}{\partial \alpha_\mathrm{n}}=\frac{\partial \dot{W}_{\mathrm{i}}}{\partial \alpha_\mathrm{n}}+\frac{\partial \dot{W}_\mathrm{f}}{\partial \alpha_\mathrm{n}}+\frac{\partial \dot{W}_\mathrm{s}}{\partial \alpha_\mathrm{n}}+\frac{\partial \dot{W}_1}{\partial \alpha_\mathrm{n}}=0 $

(9)

式中各功率偏导表达式如下：

$ \frac{\partial \dot{W}_{\mathrm{i}}}{\partial \alpha_{\mathrm{n}}}=\frac{2}{\sqrt{3}} \sigma_{\mathrm{s}} \ln \left(\frac{H}{h}\right) \frac{\partial U}{\partial \alpha_{\mathrm{n}}} $

(10)

$ \begin{array}{l} \frac{{\partial {{\dot W}_{\rm{f}}}}}{{\partial {\alpha _{\rm{n}}}}} = 2mkb\left\{ { - 2R{v_{\rm{r}}} + \frac{{U\quad \sec {\alpha _{\rm{n}}}}}{{2b{A_1} - \cos {\alpha _{\rm{n}}}}} - \frac{{\sec \bar \alpha }}{{2b\sqrt {{A_1} - 1} }}} \right.\\ \left. {\left[ {\arctan \left( {\sqrt {\frac{{{A_1} + 1}}{{{A_1} - 1}}} \tan \frac{\theta }{2}} \right) - 2\arctan \left( {\sqrt {\frac{{{A_1} + 1}}{{{A_1} - 1}}} \tan \frac{{{\alpha _{\rm{n}}}}}{2}} \right)} \right]\frac{{\partial U}}{{\partial {\alpha _{\rm{n}}}}}} \right\} \end{array} $

(11)

$ \frac{\partial \dot{W}_{\mathrm{s}}}{\partial \alpha_{\mathrm{n}}}=\frac{k}{2} \tan \theta \frac{\partial U}{\partial \alpha_{\mathrm{n}}} $

(12)

$ \frac{\partial \dot{W}_{\mathrm{t}}}{\partial \alpha_\mathrm{n}}=\left(\sigma_{\mathrm{b}}-\sigma_{\mathrm{f}}\right) \frac{\partial U}{\partial \alpha_\mathrm{n}} $

(13)

$ \frac{\partial U}{\partial \alpha_{\mathrm{n}}}=v_{\mathrm{r}} b R \sin 2 \alpha_{\mathrm{n}}-v_{\mathrm{r}} b(R+h) \sin \alpha_{\mathrm{n}} $

(14)

将解得的α_n代入总功率泛函便可得总功率最小值Φ_min。于是，轧制力矩可以表示为

$ M_{\mathrm{z}}=R \varPhi_{\min } / 2 v_{\mathrm{r}} $

(15)

上界法既无法考虑因轧制长度造成的摩擦系数时变性，又无法计算损失功率，导致其计算结果虽然符合工程需求，但精度却不够高。因此将塑性变形功率、张力功率和剪切功率解析解作为数据驱动模型的输入量，以提高模型精度。

2 LSTM网络模型

人工神经网络(artificial neural network, ANN)是通过大量简单的神经元之间的相互连接来构造复杂的网络结构，但同层神经元之间没有信息传递导致无法处理时序数据。传统的循环神经网络(recurrent neural network, RNN)允许信息在同层之间传递，但在训练数据样本量过大时又会出现梯度爆炸或梯度消失的现象^[22]。针对这一问题，学者们采用一种特殊的RNN模型成功的解决了长期依赖问题^[23]，即LSTM网络，其原理见图 2。

LSTM网络由多个单元构成，每个单元包含3个门结构并以特殊的方式组合，遗忘门控制哪些历史信息将被遗忘，记忆门控制哪些新输入的信息将被保留，输出门则决定哪些信息会被输出。图 2中各输出向量的计算如下所示：

$ \boldsymbol{f}_t=\sigma\left(\left[\boldsymbol{X}_t, \boldsymbol{h}_{t-1}\right] \boldsymbol{w}_{\mathrm{f}}+\boldsymbol{b}_{\mathrm{f}}\right) $

(16)

$ \boldsymbol{i}_t=\sigma\left(\left[\boldsymbol{X}_t, \boldsymbol{h}_{t-1}\right] \boldsymbol{w}_{\mathrm{i}}+\boldsymbol{b}_{\mathrm{i}}\right) $

(17)

$ \widetilde{\boldsymbol{C}}_t=\tanh \left(\left[\boldsymbol{X}_t, \boldsymbol{h}_{t-1}\right] \boldsymbol{w}_{\mathrm{C}}+\boldsymbol{b}_{\mathrm{C}}\right) $

(18)

$ \boldsymbol{o}_t=\sigma\left(\left[\boldsymbol{X}_t, \boldsymbol{h}_{t-1}\right] \boldsymbol{w}_{\mathrm{o}}+\boldsymbol{b}_{\mathrm{o}}\right) $

(19)

$ \boldsymbol{h}_t=\boldsymbol{o}_t \tanh \boldsymbol{C}_t $

(20)

$ \boldsymbol{C}_t=\boldsymbol{f}_t \boldsymbol{C}_{t-1}+\boldsymbol{i}_t \widetilde{\boldsymbol{C}}_t $

(21)

式中：$\boldsymbol{f}_t 、\boldsymbol{i}_t 、\widetilde{\boldsymbol{C}}_t$和o_t分别为遗忘门、记忆门和输出门的输出向量，C_t为细胞状态，σ和tanh为激活函数，w为权重矩阵，b为偏置向量。

3 上界法与LSTM时序耦合预测模型

本文所有数据均来自某1 450 mm六辊五机架冷连轧产线。考虑到各工艺参数的采集来自不同装置，它们的采样频率并不统一，因此以最大采样间隔0.24 s为基准进行等间隔降采样。从工作辊刚上机开始，提取一个换辊周期内的稳态轧制数据共计32 062条、17维特征参数，并可根据特征属性分为4类，其中轧件信息6条，工艺参数6条，轧辊数据2条以及功率解析解3条，见表 1。

时序耦合预测模型结构见图 3。首先从PDA及二级数据库中导出所需的14维特征参数，然后将轧件信息v₁~v₆及工艺参数v₇~v₁₂代入上界法求解塑性变形功率v₁₅、剪切功率v₁₆和张力功率v₁₇的解析解，再将所有特征输入LSTM网络中学习训练摩擦功率和损失功率等具有时序特点部分，最后将训练结果用于实际电机功率预测。

经过结构优化并反复测试，最终确定本文的LSTM网络由1个输入层、2个隐含层和1个输出层构成，每个隐含层包括25个LSTM单元节点，层与层之间采用全连接。图 4为不同学习率下的训练集均方根误差(root mean square error, RMSE)，可以发现均方根误差随学习率的增大先下降后上升，由此确定最优学习率为0.001 5。此外，时间步长也是影响LSTM网络预测效果的重要超参数，不同步长下的训练集均方根误差见图 5，最终选取时间步长为20。

在训练之前，首先要对数据进行缺失值处理及数据归一化。实际生产中常常会因传感器故障等问题导致数据缺失，本文中由于缺失值数据不超过20条，故采用上条数据进行填充。对于不同的工艺参数，其数值上会有成百上千倍的差异。若不进行归一化处理，模型的训练精度就会大打折扣，甚至无法得出正确的结果。因此采用最常见的归一化方法，如下式所示：

$ x_i^*=\frac{x_i-x_{\min }}{x_{\max }-x_{\min }} $

(22)

式中：x_i^*为归一化后数据，x_i为原始数据，x_min为特征中的最小值，x_max为特征中的最大值。

最后将所有的数据分为训练集和测试集，数据集占比为8∶ 2。训练集用于模型权重训练，测试集用于模型的泛化能力评估。模型训练过程中的损失函数见图 6，可以看到，当迭代1 000轮次左右后基本达到稳定状态，最优平均绝对误差(mean absolute error, MAE)可达到0.041 7。

4 结果与讨论

为了验证预测精度，本文对上界法、ANN和时序耦合预测模型的预测结果进行对比。其中，上界法采用本文建立的速度场得到总功率解析解；ANN和时序耦合预测模型则将上界法中部分功率的解析解(塑性变形功率、张力功率和剪切功率)和其他14维参数作为输入量，按照顺序输入模型学习训练。

不同模型的测试集预测效果对比见图 7。从整体上，上界法、ANN和时序耦合预测模型均能实现电机功率预测，但ANN和耦合模型的效果明显优于上界法；三者的最大相对误差分别为±8.5%、±6.3%和±3.9%，可见耦合模型的预测效果最好。值得注意的是，当相邻卷带钢轧制规程变化较大时(如采样序列800~1000和2400~2500)，耦合模型的预测效果明显优于其他2种，体现了考虑数据时序特性的优势。从局部上，耦合模型的误差绝对值更小且分布更均匀。

图 8为3种模型的残差分布直方图，可以看出：上界法残差集中在±50 kW附近且不符合正态分布。其主要原因是没有捕获到变量之间的强耦合关系，如：摩擦系数和损失功率与其他参数间的关系。而ANN和耦合模型则均基本遵循正态分布，表明这2种模型更加合理，其中耦合模型的残差在0附近的频数更高、方差更小，说明耦合模型的拟合精度要高于ANN。

为进一步对比模型预测效果，选取广泛应用于拟合模型性能评估的平均绝对误差、均方根误差和相关系数对文中模型进行评价，结果见表 2。可以发现：上界法、ANN和耦合模型的平均绝对误差和均方根误差依次递减，但相关系数值逐步增大，这是由于基于数据驱动的ANN和耦合模型可以更好地提取数据之间的特征，相较于传统机理模型具有较大的提升；LSTM在ANN的基础上进一步考虑了轧制长度与摩擦系数的时序变化特点，因此耦合模型具有更好的综合预测性能，平均绝对误差、均方根误差和相关系数分别可达到28.0、35.9和0.936。

应用时序耦合预测模型分析了电机功率与轧制长度、轧制速度的关系，见图 9。电机功率与轧制长度呈反比，其原因是：工作辊刚上机时粗糙度较大，此时摩擦系数也会变大，但粗糙度会随轧制长度增加而衰减，因此电机功率减小。电机功率与轧制速度呈正比，其原因是：轧制速度越高，单位时间通过辊缝的带钢越多，电机功率也随之增加。

5 结论

1) 提出了一种基于上界法和LSTM网络耦合的电机功率预测模型，其中塑性变形功率、剪切功率和张力功率采用上界法计算解析解，摩擦功率和损失功率采用LSTM网络预测学习数据的时序特性。

2) 上界法、ANN和时序耦合预测模型的预测结果对比表明：考虑数据时序特性的耦合模型明显优于其他2种方法，残差分布更合理，计算误差不超过±3.9%，满足在线控制要求。

3) 时序耦合预测模型中考虑轧制长度的时序特性对预测精度有很大提升，主要是由于轧辊粗糙度会随轧制长度增加而衰减，并进一步影响轧制力矩和电机功率。