用牛顿法解算矿井通风网络方程组

式(9)这样的通风网络方程组的规模通常很大，因此一般只能进行数值解算。当方程组是零维时，一些通用的多变量非线性方程组的数值求解方法都同时适用于该方程的求解，最常用的是牛顿法以及各种拟牛顿法。相对于以往常用的通风网络求解的斯科特-恒斯雷法、京大二试法等，牛顿法具有数学逻辑较好、收敛速度快、对各种非线性方程（组）通用的优点。牛顿法的缺点是所得的解严重依赖所提供的初值，并且不能找出方程组所有的解。下面将推导用牛顿法求解式(9)的过程。

牛顿法介绍#

牛顿法（Newton’s method）又称为牛顿-拉弗森方法（Newton-Raphson method），它是一种在实数域和复数域上近似求解方程的方法。方法使用函数 f(x) 的泰勒级数的前面几项来寻找方程 f(x) = 0 的根。该方法同样可以用于方程组求解。

图3 用牛顿法迭代求解方程的过程示意图

对于如下具有 k 个变量和 l 个方程的连续可微方程组：

用牛顿法求解方程组的迭代步骤为：

式中：

• s——迭代过程中临时产生的未知数，且 s = x``i + 1 - x``i；
• J``F——向量值函数 F(x) 的雅可比矩阵（可以将其理解为由一组偏导数构成的矩阵），其计算方法为：

由以上介绍可以看出，在用牛顿法求解方程组时，首先解得方程组对应向量值函数 F(x) 的雅可比矩阵 J``F，然后给定未知数的初值 x0（将其代入 xi），每次迭代时通过式 $\boldsymbol{J}_F\left( \boldsymbol{x}_i \right) \boldsymbol{s}=-\boldsymbol{F}\left( \boldsymbol{x}_i \right)$ 解得 s，进一步再通过式 $\boldsymbol{x}_{i+1}=\boldsymbol{x}_i+\boldsymbol{s}$ 求得下一个迭代值 x``i + 1，如此反复迭代直到误差足够小。除了求解雅克比矩阵 J``F，其他步骤均可以数值计算软件代为完成。

通风网络方程组对应的雅可比矩阵#

对于式(9)，其分支风量 Q(e) 和节点风压 P(v) 按如下方式组成为未知数向量：

为了推导得出式(9)对应的雅克比矩阵，将该式中第 1 式左侧函数记为 S，将第 2 式左侧函数记为 E。即有：

S 的第 i 行所对应的雅可比矩阵行向量为：

式中：

• si——节点 i 所对应的风量守恒方程；
• qj——分支 j 的风量。

由于 $\frac{\partial s_i}{\partial q_j}=m_{ij}$，S 对应的雅克比矩阵为：

E 的第 i 行所对应的雅可比矩阵行向量为：

式中：

• Ei——分支 i 所对应的能量平衡方程；
• Pt k——节点 k 的全压。

E 中各项对 Q(e) 的偏导数具有如下规律：

1. E 中 diag(R(e))diag(|Q(e)|)Q(e) 项只是分支风量的函数，故该项中第 i 个元素对 Qj 的偏导数可分为两种情况：当 i ≠ j 时，其偏导数为 0；当 i = j 时，偏导数为 2Rj|Qj|。
2. E 中 HN(e) 项基本不受风量和全压的影响，故其偏导数总为 0。
3. E 中 Hf(e) 项中第 i 个元素对 Qj 的偏导数可分为两种情况：当 i ≠ j 时，其偏导数为 0；当 i = j 时，其偏导数为 $\frac{\partial H_{fi}}{\partial Q_i}$（具体形式根据风机个体特性曲线回归方程的形式求导确定）。

因此，E 中各行对 Q(e) 的偏导数为：

E 中仅有 MTPt(v) 项对 Pt(v) 的偏导数不为零，从而 E 中各行对 Pt(v) 的偏导数为：

E 对应的雅克比矩阵为：

注意，上式的矩阵由两项构成（不要将其看作是一项），其中 MT 是末项，它和前面一项没有减法关系。

拼接 JS 和 JE，得到式(9)左侧的向量值函数对应的雅可比矩阵为：

(10)