计算机建模和解算

以上所给出的通风网络方程组及其牛顿法求解过程涉及大量的矩阵运算，尤其是大规模的线性代数运算。这些运算如果只用编程语言的基础功能（如用 Python 的列表表示矩阵）来建模和求解，是极其困难的，而需要借助某些科学计算软件包（如 Python + NumPy + SciPy 或 MATLAB）。以下文档叙述基于 Python 的 NumPy 和 SciPy 包进行矿井通风网络的计算机建模与解算。

配置开发环境#

NumPy 是使用 Python 进行科学计算的基础包。它为 Python 提供了一个强大多维数组对象，各种派生对象（如掩码数组和矩阵），以及一系列对数组进行快速操作的例程，包括数学、逻辑、形状操作、排序、选择、I/O、离散傅立叶变换、基本线性代数、基本统计操作、随机模拟，等等。当将 Python 用于科学和工程领域时，几乎都要用到此包。

SciPy 是构建在 NumPy 之上的数学算法和常用函数的集合。其不同的子包对应不同的应用，如插值、积分、优化、图像处理、统计、特殊函数等。这里我们主要使用 SciPy 中最优化子包 optimize 中用于非线性方程组求解的 fsolve 函数。

像 NumPy 和 SciPy 这类科学计算软件包一般包含 FORTRAN、C、C++ 等语言编写的二进制编译文件依赖，因此直接用 Python 自带的包安装器 pip 安装并不容易。NumPy 的官网提供了该包的多种安装方法，对于 Windows 操作系统，其中最简单的方法就是直接安装 Anaconda 科学计算平台，这是一个包管理器、一个环境管理器、一个 Python 数据科学发行版以及众多开源包的集合，其中当然也包含 NumPy 和 SciPy。本次课程设计将使用这个开发平台，请从官网下载其安装包（比较大，大约 800MB，而 Python 安装包才 25MB 左右）并按照提示安装。在安装完成后，Windows 操作系统的开采菜单将会出现“Anaconda3 (64-bit)”组，其他包含多个启动项，我们这里只需要选择“Anaconda Prompt”或“Anaconda Powershell Prompt”启动项，就可以像正常安装 Python 那样通过命令行使用 Python 解释器。另外，还需要设置 VS Code 或 PyCharm 开发环境，使其使用的 Python 解释器正好就是刚刚安装 Anaconda 的位置（例如：anaconda3\python.exe）。

相关技术介绍#

前面已经说过，在通风网络方程组（式(9)）和其相应的雅克比矩阵（式(10)）中，存在大量的向量和矩阵。这些向量和矩阵并不是用 Python 基础的一维列表、二维列表来表示，而是用 NumPy 中定义的一维数组、二维数组（矩阵）来表示。NumPy 还提供了操作这些数组的基础函数，详情请见 NumPy 的文档。

在求解非线性方程组时，则需要使用 scipy.optimize.fsolve 函数，该函数实际上是对著名 MINPACK 库的 hybrd 和 hybrj 算法的封装。该函数返回一个由 func(x) = 0 定义的（非线性）方程组在一个初始估计值 x0 附近的根。封装后 fsolve 函数的签名如下：

scipy.optimize.fsolve(func, x0, args=(), fprime=None, full_output=0, col_deriv=0, xtol=1.49012e-08, maxfev=0, band=None, epsfcn=None, factor=100, diag=None)

可以看出，fsolve 的签名非常复杂。鉴于该函数的重要性，现将其英文文档翻译为中文：

参数：

  • func：可调用的 f(x, *args)

    一个函数，该函数接受至少一个（可能为向量）参数，并返回同样长度的值。

  • x0：ndarray 类型

    func(x) = 0 根的初值（初始估计值）。

  • args：元组类型，可选的

    func 的额外实参。

  • fprime：可调用的 f(x, *args)，可选的

    一个函数，用于计算 func 的雅可比矩阵，其导数跨越各行。默认没有提供此函数时，将自动估算雅可比矩阵。

  • full_output：布尔类型，可选的

    当此参数为 True 时，将反正可选的输出。

  • col_deriv：布尔类型，可选的

    指定雅可比函数是否沿列计算导数，因为没有转置操作，沿列计算会更快。

  • xtol：浮点类型，可选的

    如果两次连续迭代之间的相对误差达到最大值 xtol，计算将终止。

  • maxfev：整数类型，可选的

    对函数的最大调用次数。如果为零，则最大次数为 100*(N+1) ，其中 N 是 x0 中元素的数量。

  • band：元组类型，可选的

    如果设置为包含雅克比矩阵带内的子对角线和超对角线数量的两个序列，则雅克比矩阵被认为是带状的（仅适用于 fprime=None）。

  • epsfcn：浮点类型，可选的

    雅可比矩阵前向差分近似的合适步长（fprime=None）。如果 epsfcn 小于机器精度，则假设函数中的相对误差与机器精度相当。

  • factor：浮点类型，可选的

    确定初始步长界限（factor * || diag * x||）的参数 。应该在区间 (0.1，100) 内。

  • diag：序列类型，可选的

    N 个正条目，用作变量的比例因子。

返回值：

  • x：ndarray 类型

    解（或不成功调用时的最后一次迭代的结果）。

  • infodict：字典类型

    可选输出的字典，带有键：

    • nfev：函数调用的次数
    • njev：雅克比函数调用的次数
    • fvec：在输出时计算所得的函数值
    • fjac：由最终近似雅可比矩阵 QR 分解产生的正交矩阵 q，按列存储
    • r：由相同矩阵 QR 分解产生的上三角矩阵
    • qtf：向量 (transpose(q) * fvec)

  • ier：整数类型

    一个整数标志。如果找到解，则设置为 1，否则请参考 mesg 了解更多信息。

  • mesg：字符串类型

    如果没有找到解决方案，mesg 会详细说明失败的原因。

核心代码分析#

为了降低编程的难度，这里给出了程序的核心代码，并对代码给出了大量的注释，请从此链接下载，并进行解压（如果解压到 D 盘根目录，则不用对代码进行任何修改）。本小节进一步对源代码进行解释。

首先，mine_vertilation_network.py 代码文件的开头是如下 3 行：

import os
import numpy as np
from scipy.optimize import fsolve

导入 os 模块是为了使用 os.chdir 函数；接下来一行不仅导入了 numpy 包，并将其重命名为 np；scipy 包的内容很多，我们只需要使用其中的 optimize 子包的 fsolve 函数。

os.chdir('D:/mine_vertilation_network')
incidence_matrix_path = 'incidence_matrix.csv'  # 关联矩阵 M
edge_resistances_path = 'edge_resistances.csv'  # 分支风阻 R
incidence_matrix = np.loadtxt(incidence_matrix_path, dtype=np.int64, delimiter=',')
edge_resistances = np.loadtxt(edge_resistances_path, dtype=np.float64, delimiter=',')

我们并不想把建模数据写死在源代码中，而是把重要的建模数据，如关联矩阵、分支风阻分别放到 incidence_matrix.csv 和 edge_resistances.csv 文件中，这使我们可以使用 Excel 或 VS Code 编辑其中的数据。然后在代码中使用 np.loadtxt 分别将这两个文件中的数据载入到 incidence_matrix 和 edge_resistances 变量中，这两个变量分别是矩阵和向量。

(node_num, edge_num) = incidence_matrix.shape
x_num = node_num + edge_num

这两行分别计算了节点数 node_num（关联矩阵的行数，前面的 m 值）和分支数 edge_num（关联矩阵的列数，前面的 n 值），这时通过调用 incidence_matrix 的 shape 属性得到的，该属性返回一个整型元组，其各个元素分别为多维数组各个维度的元素个数。node_num 和 edge_num 之和正好就是方程未知数的个数 x_num（m + n）。我们提前把这些数值算出来，以方便后续调用。

surface_pressure = 0.0
ref_node_index = 0

在前面得出的数学模型中，我们假定第 r 个节点的风压 Pt r 为参考节点的风压，并一般以矿井地表空气对应的节点作为 r 节点。在实际建模中，一般习惯将大气节点作为第一个节点，其索引 ref_node_index 为 0，其对应的全压 surface_pressure 为 0.0Pa。当然，这两个值都是可以改变的，只是我们还没有来得及不将他们编死在代码中。

fans = [{'ni': 3, 'a': 0.0, 'b': 0.0327, 'c': -18.464, 'd': 1146.3},
{'ni': 4, 'a': 0.0, 'b': -0.1971, 'c': 18.75, 'd': -18.322}]

这两行关于通风机的参数也是暂时编死在代码中的。一般可使用一元二次、一元三次方程组描述通风机个体特性曲线，形如：

$$ h\_f=aq^3+bq^2+cq+d $$

矿井中只有少数巷道存在通风机，因此描述通风机时，只需要给出其对应的节点索引 ni，以及各个通风机个体特性曲线的参数 a、b、c 和 d 即可。本模型中有两个通风机，其所在分支的索引分别为 3 和 4。这两个通风机的三次项参数 a 都为 0.0，说明实际上用一元二次方程组就可以描述他们的个体特性曲线。

flowrates_init = 30.0 * np.ones(edge_num, dtype=np.float64)  # 风量初值
pressures_init = -50.0 * np.ones(node_num, dtype=np.float64)  # 全压初值
pressures_init[ref_node_index] = surface_pressure  # 设置地面节点的全压
x0 = np.concatenate((flowrates_init, pressures_init))

这 4 行设置了未知数的初始估计值。这种批量设置初值的方法极为粗糙，实际设置初值时有更多复杂的技巧，不过暂时能用就行。在实际计算中，我们设置的初值可能并不能得出满意的计算结果，如对于抽出式通风，计算所得的节点相对全压存在大量大于 0 的值显然是不合理的，这里只能手工尝试设置不同的初值，直到满意为止。风量初值和风压初值都属于未知数 x0 的初值，我们通过 np.concatenate 将他们拼接在一起。

def mvn_func(x):
    """该函数根据迭代过程中临时得出的 x 向量自动计算通风网络方程中的向量函数值 F(x)，该方程用于 fslove() 函数中的 func 参数"""

    # 未知数 x 的前 edge_num 项为分支风量 flowrates，后面各项为节点全压 pressures。
    # pressures 中第 ref_node_index 项为地面节点的全压，该项风压总是不变的，因此需要在每个循环中重新设定
    flowrates = x[:edge_num]  # 分支风量
    pressures = x[edge_num:]  # 节点全压
    pressures[ref_node_index] = surface_pressure  # 地面节点全压

    # 计算 M×Q
    f_nodes = incidence_matrix @ flowrates
    # 计算 diag(R)×diag(|Q|)×Q - M'×P
    f_edges = (np.diag(edge_resistances) @ np.diag(np.abs(flowrates)) @ flowrates) - (incidence_matrix.T @ pressures)
    # 计算通风机风压，并附加的特定节点的 f_edges 上
    for fan in fans:
        ni = fan.get('ni')
        ri = abs(flowrates[ni])
        f_edges[ni] -= fan.get('a') * pow(ri, 3) + fan.get('b') * pow(ri, 2) + fan.get('c') * ri + fan.get('d')

    # 拼接通风网络方程组中左侧函数值的两项，得到函数值向量
    f = np.concatenate((f_nodes, f_edges))
    return f

mvn_func(x) 函数与式(9)所示的通风网络方程组左侧向量函数相对应。该函数的参数 x 即为待求解的变量，它由 fsolve 在迭代过程中自动求得，我们需要根据 x 计算式(9)左侧部分的值。由于该式分为多个部分，为了计算方便，我们也把 x 切分为风量向量 flowrates 和风压向量 pressures，分别计算各个节点对应的风量守恒方程左侧向量 f_nodes，以及各个分支对应的能量守恒方程左侧向量 f_edges。后者的计算要更复杂一点，因为我们要对通风机 fans 进行变量，将其产生的风压附加的分支上。最后，将这两个向量的值拼接在一起并返回。

np.diag() 函数的作用是将一个向量转换为对角矩阵，np.abs() 函数的作用则是计算数组中各个元素的绝对值，并重新构造一个数组。尤其需要注意的是，在 NumPy 中，a * b 并不是执行矩阵乘法，而是将两个等尺寸的矩阵的对应元素相乘，或者将一个矩阵的每个元素乘以一个标量（即所谓的广播），a @ b 执行的才是矩阵乘法。

def mvn_jac_func(x):
    """该函数根据迭代过程中临时得出的 x 向量自动计算通风网络方程中的向量函数值 F(x) 的雅克比矩阵，该方程用于 fslove() 函数中的 fprime 参数"""

    flowrates = x[:edge_num]
    # 拼接 M 和 0 (m×m），得到雅克比矩阵的第一行，axis=1 参数表示对列进行拼接
    f_jac_nodes = np.concatenate((incidence_matrix, np.zeros((node_num, node_num))), axis=1)

    # 计算 2×diag(R)×diag(|Q|) 项
    f_jac_edges = 2.0 * (np.diag(edge_resistances) @ np.diag(np.abs(flowrates)))
    # 计算 diag(Hf) 项
    for fan in fans:
        ni = fan.get('ni')
        ri = abs(flowrates[ni])
        f_jac_edges[ni] -= 3 * fan.get('a') * pow(ri, 2) + 2 * fan.get('b') * ri + fan.get('c')
    # 将雅克比矩阵中第二行的两项进行拼接，得到雅克比矩阵的第二行，axis=1 参数表示对列进行拼接
    f_jac_edges = np.concatenate((f_jac_edges, -1.0*incidence_matrix.T), axis=1)

    # 拼接雅克比矩阵中的两行，得到整个雅克比矩阵
    f_jac = np.concatenate((f_jac_nodes, f_jac_edges))
    return f_jac

mvn_jac_func(x) 函数与式(10)所示的雅克比矩阵相对应。该函数的实现方式和前面的 mvn_func(x) 整体类似。

(root, infodict, ier, mesg) = fsolve(mvn_func, x0, fprime=mvn_jac_func, full_output=True)
# (root, infodict, ier, mesg) = fsolve(mvn_func, x0, full_output=True)

前面已经铺垫好了，在此处进行方程组的求解。这两行值需要使用其中一个就行。前者显式给出了雅克比矩阵，后者由系统自动估算雅克比矩阵值。由于前者能精确计算雅克比矩阵，因此其速度要更快一点。如果使用后者，则完全不必定义前面的 mvn_jac_func(x) 函数。通过为 fsolve() 函数提供 full_output=True，使返回值变为 4 个（否则只有一个 root），这样可以基于这些返回值作出更多的判断。

if ier:
    print('求解成功！\n各分支风量：')
    for i in range(edge_num):
        print(f'{i}:\t{root[i]:.3f}')

    print('各节点风压：')
    for i in range(node_num):
        print(f'{i}:\t{root[edge_num + i]:.3f}')
else:
    print('求解失败！\n错误信息为：\n', mesg)

最后，当我们判断 fsolve 执行成功时（ier 的值为 1），打印求解结果，否则打印错误信息。

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