0.   引言

为保障通风安全，深部矿井通常在地面配置2台主要通风机。其中一台用于井下通风，另一台作为备用通风机。主要通风机切换过程是选择备用通风机作为运行通风机，以提高矿井通风系统的可靠性^[1]。然而切换过程中存在供给风量大范围波动甚至中断的问题，导致瓦斯浓度极易超出安全允许界值^[2]。因此有必要控制主要通风机切换过程，以保障切换过程的安全性与平稳性。

目前常用的主要通风机切换控制方法为基于模型的方法。文献[3]建立了切换过程的非线性约束规划模型，采用逐步调节风门叶片角度的策略对供给风量进行控制。文献[4]在文献[3]的基础上，采用粒子群优化算法对风门角度进行调整，以实现供给风量优化控制。上述方法基于静态模型，当干扰出现时，只有在切换过程达到新的稳态时才能进行优化，因此造成优化滞后^[5]。文献[6]通过静态工况下的数据建立了切换过程的RBF神经网络模型，并采用顺序控制策略进行供给风量控制。文献[7]建立了切换过程的动态数学模型，并基于李雅普诺夫稳定性理论，提出了理想工况下的分散控制器设计方法。但上述方法难以很好地对约束进行处理。

另外，一些学者将智能算法应用于主要通风机切换过程控制。针对切换过程难以用精确模型描述的问题，文献[8-9]采用模糊PID技术实现供给风量控制。文献[10-11]将基于神经网络和多模型切换的智能解耦控制器与PID算法相结合，提出了智能解耦控制方法，用于通风机切换过程供给风量优化。然而智能算法依赖模拟领域专家的知识和经验，具有一定的主观性和随意性。

本文以中国平煤能源化工集团有限责任公司（以下称平煤集团）二矿主要通风机切换过程为研究背景，探求动态特性具有强非线性、系统状态受约束的切换过程供给风量优化控制问题。首先，基于流体动力学方程和图论概念，得到切换过程的动态模型；然后，通过泰勒展开式对其进行线性化处理，降低计算复杂度，从而将基于模型预测控制(Model Predictive Control，MPC)的供给风量优化问题描述为一个二次规划问题；最后，采用原对偶神经网络对所建优化问题进行在线求解，实现对切换过程中供给风量的实时优化控制。

1.   主要通风机切换过程模型

1.1   主要通风机切换过程动态模型

主要通风机切换过程如图1所示。该过程涉及3个部分：地下矿井、2条主通风管道、2台主要通风机。从图1可看出，2条主通风管道的风量相互耦合、风门的协同动作、地下矿井的生产活动等多种因素导致主要通风机切换过程特性随运行工况动态变化。因此，主要通风机切换过程本质上是一个复杂的非线性动态系统。

图  1  矿井主要通风机切换过程结构

Figure  1.  Structure of switchover process of mine main ventilators

下载: 全尺寸图片 幻灯片

根据文献[12]，将主要通风机切换过程描述为如图2所示的网络拓扑结构。节点①为1号出风口；节点②为2号出风口；节点③为2条主通风管道连接处；节点④为地下矿井； $ {{\rm{R}}_1} $ ， $ {{\rm{R}}_2} $ 分别为1号主要通风机水平风门、垂直风门的风阻； $ {{\rm{R}}_3} $ ， $ {{\rm{R}}_4} $ 分别为2号主要通风机垂直风门、水平风门的风阻； $ {{\rm{R}}_5} $ 为地下矿井风阻； $ {{\rm{M}}_1} $ ， $ {{\rm{M}}_2} $ 分别为1号、2号主要通风机；1—5为风量分支号，分别对应1号水平风门、1号垂直风门、2号垂直风门、2号水平风门、总供风风量分支。

将通过主通风管道的风量视作流体，其运动方程可描述为

式中： $S$ 为管道截面积；L为管道长度； $\rho $ 为管道内气体密度； $v$ 为管道内气体流速；t为时间；H为管道2个端点间的压降； $\tau $ 为流体表面的切应力； $D$ 为管道直径。

图  2  矿井主要通风机切换过程网络拓扑结构

Figure  2.  Topology structure of switchover process of mine main ventilators

下载: 全尺寸图片 幻灯片

由Darcy−Weisbach方程可得

式中 ${f}$ 为摩擦因数。

风量可表示为^[13]

由式（3）求导可得

将式（2）和式（4）代入式（1）可得

对式（5）进行数学变换，可得

将式（3）代入式（6）可得

因惯性系数 $ K =\dfrac{S}{L} $ ，风阻 $ R = \dfrac{{{f}L}}{{{\text{2}}\rho D{S^2}}} $ ，则可得管道内风量模型：

根据式（8），可将主要通风机切换过程的分支动态模型描述为

式中 $ {Q_j}，{K_j}，{R_j}，{H_j} $ 分别为分支 $j$ （j=1，2，…，5）对应的风量、惯性系数、风阻和压降。

根据基尔霍夫风压定律和风量定律，可得主要通风机切换过程的动态模型：

式中： $ {\boldsymbol{Q}}_{\rm{w}} $ 为切换过程中4个风门的风量， ${\boldsymbol{Q}}_{\rm{w}} = {[{Q_1}\;\;{Q_2}\;\;{Q_3}\;\;{Q_4}]^{\text{T}}}$ ； ${{\boldsymbol{Z}}_1} = - {\left( {{\boldsymbol{K}}_{\rm{c}}^{ - 1} + {\boldsymbol{V}}_{\rm{c}}^{\rm{T}}{\boldsymbol{K}}_{\rm{a}}^{ - 1}{{\boldsymbol{V}}_{\rm{c}}}} \right)^{ - 1}}$ ， ${{\boldsymbol{K}}_{\rm{c}}} = {\rm{diag}} \left( {K_1}， {K_2}，{K_3}，{K_4} \right)$ ， $ {{\boldsymbol{V}}_{\rm{c}}} $ 为风量平衡方程的常数矩阵， $ {{\boldsymbol{K}}_{\rm{a}}} = [{K_5}] $ ； $ {{\boldsymbol{Q}}_{\rm{a}}} = {\rm{diag}} \left( {{Q_1},\;{Q_2},\;{Q_3},\;{Q_4}} \right) $ ； $ {\boldsymbol{u}} $ 为控制量， ${\boldsymbol{u}} = [{R_1}\;\;\; {R_{2}}\;\;\; {R_3}\;\;\; {R_4}]^{\text{T}}$ ； ${{\boldsymbol{Z}}_2} = {{\boldsymbol{K}}_{\rm{c}}}{\boldsymbol{V}}_{\rm{c}}^{\rm{T}}{\left( {{{\boldsymbol{V}}_{\rm{c}}}{{\boldsymbol{K}}_{\rm{c}}}{\boldsymbol{V}}_{\rm{c}}^{\rm{T}} + {{\boldsymbol{K}}_{\rm{a}}}} \right)^{ - 1}}{{\boldsymbol{K}}_{\rm{a}}}$ ； ${{\boldsymbol{Z}}_3} = \left( {\boldsymbol{K}}_{\rm{c}}^{ - 1} + {\boldsymbol{V}}_{\rm{c}}^{\rm{T}}{\boldsymbol{K}}_{\rm{a}}^{ - 1} {{\boldsymbol{V}}_{\rm{c}}} \right)^{ - 1}$ ； ${\boldsymbol{d}} = {[{P_1}\;\;{P_1}\;\;{P_2}\;\;{P_2}]^{\rm{T}}}$ ， $ {P_1} $ ， $ {P_2} $ 分别为1号和2号主要通风机风压。

P₁，P₂可通过风量来描述^[14]：

式中 $ {a_0},{a_1},{a_2} $ 为待定常数，由主要通风机特性曲线决定。

1.2   线性误差模型

将式（10）写为以下形式：

主要通风机切换过程需要4个风门的风量按照参考值变化，因此可将其视为跟踪问题。将参考轨迹描述为

式中： $ {{\boldsymbol{Q}}_{\rm{r}}} $ 为风门参考风量， $ {{\boldsymbol{Q}}_{\rm{r}}} = {[{Q_{1{\rm{r}}}}\;\;{Q_{2{\rm{r}}}}\;\;{Q_{3{\rm{r}}}}\;\;{Q_{4{\rm{r}}}}]^{\text{T}}} $ ， $ {Q_{1{\rm{r}}}} $ ， $ {Q_{2{\rm{r}}}} $ ， $ {Q_{3{\rm{r}}}} $ ， $ {Q_{4{\rm{r}}}} $ 分别为1号水平风门、1号垂直风门、2号垂直风门、2号水平风门的参考风量； $ {{\boldsymbol{u}}_{\rm{r}}} $ 为参考控制量。

将主要通风机切换过程动态模型的一般形式在点 $ \left( {{{\boldsymbol{Q}}_{\rm{r}}}，{{\boldsymbol{u}}_{\rm{r}}}} \right) $ 处用泰勒展开式展开，忽略其高阶项，得

则误差模型可表示为

式中： $ {\boldsymbol{e}} $ 为风量误差； ${\boldsymbol{A}}_{\rm{e}}=\dfrac{{\partial g\left( {\boldsymbol{Q}_{\rm{w}},\boldsymbol{u}} \right)}}{{\partial \boldsymbol{Q}_{\rm{w}}}} \left| {_{\begin{array}{*{20}{l}} {\boldsymbol{Q}_{\rm{w}} = {\boldsymbol{Q}_{\rm{r}}}} \\[-2pt] {{\boldsymbol{u}} = {{\boldsymbol{u}}_{\rm{r}}}} \end{array}}} \right.$ ； ${{\boldsymbol{B}}}_{{\rm{e}}}= \dfrac{{\partial g\left( {{\boldsymbol{Q}}_{\rm{w}},{\boldsymbol{u}}} \right)}}{{\partial {\boldsymbol{u}}}}\left| {_{\begin{array}{*{20}{l}} {{\boldsymbol{Q}}_{\rm{w}} = {{\boldsymbol{Q}}_{\rm{r}}}} \\[-2pt] {{\boldsymbol{u}} = {{\boldsymbol{u}}_{\rm{r}}}} \end{array}}} \right.$ ；u_e为控制量误差。

为了基于上述模型来设计反馈控制器，需对式（16）进行离散化处理，得到离散误差模型：

式中： $ k $ 为采样时刻；e（k）为k时刻风量误差； $ {\boldsymbol{A}}(k) = {\boldsymbol{I}} + T {{\boldsymbol{A}}_{\rm{e}}} $ ， $ {\boldsymbol{I}} $ 为单位矩阵，T为采样时间； ${\boldsymbol{B}}(k) = T {{\boldsymbol{B}}_{\rm{e}}} $ ；u_e（k）为k时刻控制量误差。

2.   主要通风机切换过程MPC系统设计

2.1   MPC系统结构

主要通风机切换过程是一种离散的多输入多输出系统，采用MPC算法可很好地满足各种约束，且采用滚动优化策略，使风量较好地跟踪其参考值，保证主要通风机切换过程中通风稳定。

MPC算法由预测模型、反馈校正、滚动优化3个部分组成^[15-16]。本文设计的主要通风机切换过程MPC系统结构如图3所示。 $ {{\boldsymbol{Q}}_{\rm{r}}}(k) $ 为k时刻参考风量； $ {\boldsymbol{Q}}_{\rm{w}}(k) $ 为k时刻实际风量； $ {{\boldsymbol{u}}_{\rm{r}}}(k) $ 为k时刻参考控制量； $ {\boldsymbol{u}}(k) $ 为k时刻最优控制量。则可将系统控制器设计为 $ {\boldsymbol{u}}(k) = {{\boldsymbol{u}}_{\rm{e}}}(k) + {{\boldsymbol{u}}_{\rm{r}}}(k) $ 。

图  3  主要通风机切换过程MPC系统结构

Figure  3.  Model predictive control(MPC) system structure of switchover process of main ventilators

下载: 全尺寸图片 幻灯片

预测模型的构建是实现MPC系统的基础，可通过离散数学模型及当前时刻状态量和控制量对未来时刻的状态量进行预测。反馈校正环节用于获得系统实际风量，可通过采集实时数据避免部分干扰与误差，从而对预测模型进行修正。滚动优化描述的是通过移动的优化区间实施反复在线优化的过程。在每个采样时刻，系统会在优化区间内进行优化。当系统运行到下一时刻时，优化区间会随之前移1个时刻，并利用相同的优化方法求解得到系统当前最优控制量。

2.2   预测模型

在当前采样时刻 $ k $ 预测的k+N时刻风量误差 $ {\boldsymbol{e}}(k + N\left| k \right.) $ 可由式（17）迭代得到。定义以下变量：

式中 ${\boldsymbol{u}}_{\rm{e}}(k+1\left|k\right.) $ 为当前采样时刻k预测的k+1时刻控制量误差。

则风量误差预测模型可描述为

2.3   目标函数

根据MPC系统迭代优化特点，需在采样时刻 $ k $ 测量或估计系统当前的状态变量，并通过优化定义的目标函数来获取当前时刻最优控制量。目标函数的作用是确保风量快速准确地跟踪到拟定的参考轨迹，因此考虑将系统状态偏差及控制量写入目标函数。目标函数可表示为

式中： $ {{\boldsymbol{Z}}_{\rm{Q}}} $ ， $ {{\boldsymbol{Z}}_{\rm{R}}} $ 为权重矩阵； $ {{{N}}_{\rm{p}}} $ 为预测时域； $ {N_{\rm{c}}} $ 为控制时域， $ {N_{\rm{p}}} \geqslant {N_{\rm{c}}} \geqslant 1 $ ； $ {{\boldsymbol{u}}_{\rm{e}}}_{\min } (k)$ ， $ {{\boldsymbol{u}}_{\rm{e}}}_{\max } (k)$ 分别为控制量误差最小值和最大值； $ {{\boldsymbol{Q}}_{{\rm{w}} \min }}(k) $ ， $ {{\boldsymbol{Q}}_{{\rm{w}} \max }}(k) $ 分别为风量最小值和最大值。

目标函数中的第1项表示系统状态与参考量之间的差值，可通过该项判断系统跟踪参考轨迹的能力；第2项为系统对控制量的约束。

结合预测模型，可将目标函数重写为

式中： ${\boldsymbol{U}}_{\rm{min}} (k)$ ， ${\boldsymbol{U}}_{\rm{max}} (k) $ 分别为 ${\boldsymbol{U}}(k) $ 最小值和最大值； ${\boldsymbol{E}}_{\rm{min}} (k)$ ， ${\boldsymbol{E}}_{\rm{max}} (k) $ 分别为 ${\boldsymbol{E}}(k) $ 最小值和最大值。

将上述目标函数转换为以下二次规划问题：

式中： ${\boldsymbol{G}} = {\overline {\boldsymbol{B}}^{\rm{T}}}{{\boldsymbol{Z}}_{\rm{Q}}}\overline {\boldsymbol{B}} + {{\boldsymbol{Z}}_{\rm{R}}}$ ； ${\boldsymbol{F}} = 2{\boldsymbol{e}}^{\rm{T}}{(k)}{\overline {\boldsymbol{A}}^{\rm{T}}}{{\boldsymbol{Z}}_{\rm{Q}}}\overline {\boldsymbol{B}}$ ； ${\boldsymbol{X}} = [ {\boldsymbol{I}} \;\; { - {\boldsymbol{I}}} \;\; {\overline {\boldsymbol{B}}} \;\; { - \overline {\boldsymbol{B}}}]^{\rm{T}}$ ， ${\boldsymbol{b}}(k) = [{{\boldsymbol{U}}_{\max }}(k)\;\;\;\;{{\boldsymbol{U}}_{\min }}(k)\;\;\;\;{{\boldsymbol{E}}_{\max }}(k) - {\overline {\boldsymbol{A}}}{\boldsymbol{e}}(k)\;\;\;\;{{\boldsymbol{E}}_{\min }}(k) + {\overline {\boldsymbol{A}}}{\boldsymbol{e}}(k)]^{\rm{T}}$ 。

通过原对偶神经网络对上述优化问题进行在线求解，得到控制时域内的最优控制序列。将该序列中的第1个向量作用于系统，得到控制器 $ {\boldsymbol{u}}(k) = {{\boldsymbol{u}}_{\rm{e}}}({{k}}) + {{\boldsymbol{u}}_{\rm{r}}}(k) $ 。

3.   试验及结果分析

3.1   主要通风机切换过程参数配置

平煤集团二矿采用2台GAF28−16−1FB型轴流式通风机构建主要通风机切换系统。通过试验测量和工程经验，设置主要通风机切换过程参数： $ {R_5} = 0.035\;{\rm{kg}}/{{\rm{m}}^7} $ ，K₁=K₂=K₃=K₄=0.3，K₅=0.017， $ {{\boldsymbol{V}}_{\rm{c}}} = [0 \;\; - 1 - 1\;\; 0] $ ， $ {a_0} = 6\;847 $ ， $ {a_1} = 19.6 $ ， $ {a_2} = - 0.070\;94 $ ，1号水平风门、1号垂直风门、2号垂直风门、2号水平风门的风量初始值分别为1，299，1，150 m³/s。

为了更好地观察风量跟踪情况，需事先拟定4个风门风量的参考风量。考虑实际矿井的供给风量需求，给定参考风量 $Q_{5{\rm{r}}}=300 \;{{\rm{m}}^3}/{\rm{s}}$ ，4个风门参考风量为

3.2   MPC参数设置

设置预测时域N_p和控制时域N_c均为3，采样时间T=0.1 s，运行时间为50 s，权重矩阵 $ {{\boldsymbol{Z}}_{\rm{Q}}} = 5{\boldsymbol{I}} $ ， $ {{\boldsymbol{Z}}_{\rm{R}}} = {\boldsymbol{I}} $ ，控制量最大值为120 kg/m⁷，最小值为0，风量最大值为350 m³/s，最小值为0。

根据切换过程动态特性，将参考控制量设计为 ${{\boldsymbol{u}}_{\rm{r}}} $ = $ [30.85\exp ( - 0.01t)\;\;\;0.02\exp (0.02t)\;\;\;16.54\exp ( - 0.02t) 0.3\exp (0.01t)]^{\rm{T}}$ 。

3.3   试验结果分析

为验证设计的MPC系统在主要通风机切换过程中风量跟踪控制方面的性能，以Matlab为工具，编写MPC算法进行仿真，得到4个风门风量的跟踪控制效果和供给风量波动信息，如图4—图6所示。

图  4  4个风门的控制量

Figure  4.  Control values of four air doors

下载: 全尺寸图片 幻灯片

图  5  4个风门风量跟踪控制效果

Figure  5.  Tracking control effect of air volume through four drampers

下载: 全尺寸图片 幻灯片

图  6  供给风量

Figure  6.  Supply air volume

下载: 全尺寸图片 幻灯片

从图4可看出，4个风门的控制量满足所设置的约束且较小，对风门起到很好的保护作用。另外，1号风门和3号风门风阻从大到小，实现了风门由关到开的切换，2号和4号风门风阻由小到大，实现了风门由开到关的切换，表明系统实现了主要通风机切换。

根据矿井通风要求，供给风量在通风机切换过程中需保持较小的波动。从图5可看出，4个风门在50 s内均完成了风量调整，1号水平风门风量 $ {Q_1} $ 和2号垂直风门风量 $ {Q_3} $ 逐渐增大，1号垂直风门风量 $ {Q_2} $ 和2号水平风门风量 $ {Q_4} $ 逐渐减小， $ {Q_2} $ 减小与 $ {Q_3} $ 增大保证了供给风量平稳。主要通风机切换过程中，4个风门风量与给定参考风量的最大偏差分别为11.37，16.49，14.71，7.29 m³/s，差值较小，说明4个风门风量可很好地按照参考风量变化。此外，每个采样时刻的运算时长为0.027 s，小于采样时间0.1 s，满足主要通风机切换过程的实时性要求。

为了避免出现瓦斯超限现象，供给风量必须保持稳定，即 $ {Q_5} $ 的波动在主要通风机切换过程中必须维持在指定范围内。矿井实际需求风量为 $ 300\;{{\rm{m}}^3}/{\rm{s}} $ ，根据《煤矿安全规程》规定，主要通风机切换过程中供给风量波动不得超过10%，则供给风量上限为 $ 330\;{{\rm{m}}^3}/{\rm{s}} $ ，下限为 $ 270\;{{\rm{m}}^3}/{\rm{s}} $ 。从图6可看出，供给风量最大值为 $ 302.704\;{{\rm{m}}^3}/{\rm{s}} $ ，最小值为 $ 298.552\;{{\rm{m}}^3}/{\rm{s}} $ ，风量波动最大值仅为0.9%，满足矿井实际通风需求。

将MPC系统控制结果与常用的PID控制结果（图7）进行对比，可知PID控制下主要通风机切换过程的供给风量波动较大，最大值为 $ 322.924\;{{\rm{m}}^3}/{\rm{s}} $ ，风量波动达7.64%，虽符合《煤矿安全规程》要求，但控制效果明显劣于MPC。

图  7  PID控制结果

Figure  7.  PID control results

下载: 全尺寸图片 幻灯片

4.   结论

（1） 通过建立主要通风机切换过程动态模型，设计了一种主要通风机切换过程MPC系统，可解决复杂非线性和约束下供给风量的优化控制问题。

（2） 采用泰勒展开法对主要通风机切换过程动态模型进行线性化处理，并将原对偶神经网络用于二次规划问题的在线优化求解，大大降低了计算复杂度。

（3） 试验结果表明，主要通风机切换过程MPC系统能够保证切换过程中供给风量的平稳性，且满足实时控制要求。