基于有限差分法的多介质高温作业服设计的优化数学模型

An Optimization Mathematical Model for Design of Multi-Medium High Temperature Operation Suit Based on Finite Difference Method

摘要: 本文主要研究一种典型高温作业专用服的一维温度分布和厚度优化。首先，运用有限差分法对连续的传热模型进行离散化，利用有约束非线性规划对未定参数进行估计，仿真75摄氏度高温环境下服装内部温度分布。其次，通过减少约束条件简化单目标规划问题，利用Fibonacci法、黄金分割法、二分法和二次插值法进行求解，最优值为17.6毫米。最后，本文提出一种有效缩短区间长度求解可行边界点的算法，采用简单完全分层序列法对多目标规划模型进行求解，求得II、IV层最优厚度分别为19.2毫米与6.4毫米。

Abstract: This paper mainly studies the one-dimensional temperature distribution and thickness optimiza-tion of a typical high temperature operation special suit. Firstly, the finite difference method is used to discretize the continuous heat transfer model, and the undetermined parameters are estimated by using the constrained nonlinear programming, and the internal temperature distribution of the garment in 75 degrees Celsius high temperature environment is simulated. Secondly, by reducing the constraint conditions to simplify the single-objective planning problem, the optimal value is 17.6 mm by using the Fibonacci method, the Golden section method, the dichotomy method and the two-interpolation method. Finally, an algorithm for effectively shortening the interval length to solve the feasible boundary points is proposed, and a simple complete layered sequence method is used to solve the multi-objective programming model, and the optimal thickness of the II and IV layers is 19.2 mm and 6.4 mm.

1. 引言

高温作业专用服装通常由三层织物材料构成，记为I、II、III层，其中I层与外界环境接触，III层与皮肤之间还存在空隙，将此空隙记为IV层。性能良好的服装既可以保证工作人员的生命安全，又可以最大限度地提高工作效率。高温作业专用服装设计问题，概括的说是“在织物材料性能参数已知的前提下，优化各织物层的厚度”，使人在高温环境下工作一段时间后皮肤外侧的温度在合理范围内，从而达到降低生产成本、穿着轻便的意图。为了解决此问题，本文基于2018年全国大学生数学建模竞赛A题所提供的数据进行研究。

数据背景：将体内温度控制在37℃的假人放置在75℃高温环境的实验室中，II层厚度为6 mm、IV层厚度为5 mm、工作时间为90分钟，测量假人皮肤外侧的温度。

数据1：专用服装材料的参数值(见表1)；数据2：假人皮肤外侧的测量温度。

Table 1. Parameter values for special garment materials

表1. 专用服装材料的参数值

本文通过建立相关数学模型来解决以下问题：

问题1：还原数据背景下的高温专业服温度分布。

问题2：优选II层厚度，保证当环境温度为65℃、IV层的厚度为5.5 mm、工作60分钟时，假人皮肤外侧温度不超过47℃，且超过44℃的时间不超过5分钟。

问题3：优选II层和IV层厚度，确保当环境温度为80℃时、工作30分钟时，假人皮肤外侧温度不超过47℃，且超过44℃的时间不超过5分钟。

2. 问题分析

关于问题1，本文首先建立基本假设，合理简化实际问题中的三维模型。根据分析专业服装和空气层的主要热传递方式，结合Fourier实验定律和能量守恒定律，推导I至IV层织物材料各个层内的热传导方程，由数据2可以判断出此方程为非稳态的。依据各个层间的交界处“温度连续变化”和“热流密度连续变化”，给出各个织物材料层间的交界条件；依据牛顿冷却定律，得到传热系统边界条件；通过赋予初值条件，建立“作业服传热模型”(下述简称“模型一”)。运用前人对有限差分法的研究对“模型一”进行离散化，通过建立有约束非线性规划问题估计对流换热系数，最后选取恰当的空间步长、时间步长，利用追赶法求解三对角形线性方程组，仿真75℃高温环境下作业服的温度分布。

关于问题2，以第II织物材料层厚度为目标函数，建立单目标规划模型(下述简称“模型二”)以求得满足约束条件的最小值。沿用问题1的温度离散模型，首先对65℃高温环境、37℃恒温、工作60分钟后的人体皮肤表面温度的上界进行估计，从而简化“模型二”的约束条件，降低计算复杂度。其次利用4种一维搜索方法，即“Fibonacci法”、“黄金分割法”、“二分法”和“二次插值法”，对第II织物材料层厚度求解最优值。最后综合考虑“最优解的精度”、“程序运行的时长”两方面因素，选择最优计算方法，并给出最优解。

关于问题3，通过构建刻画生产高温作业服所需成本的目标函数 $v\left(d_1,d_2,d_3\right)$ 、刻画高温作业服重量大小的目标函数 $w\left(d_1,d_2,d_3,d_4\right)$ 和薄厚程度的目标函数 $h\left(d_1,d_2,d_3,d_4\right)$ ，其中 $d_1,d_2,d_3,d_4$ 为织物材料层的厚度(包括空气层)，建立多目标规划模型(下述简称“模型三”)以求得目标函数的最小值。沿用问题1的温度离散模型，首先对“模型三”的可行域进行定性分析；其次根据可行域边界特点，提出一种有效缩短区间长度求解可行边界点的算法；最后运用简单完全分层序列法对模型三进行求解，并给出最优解。

3. 模型的建立

在建立数学模型之前，本文首先做出如下假设：

假设1：忽略衣服褶皱，将织物层视为多层平行材料；

假设2：假设各层材料质地均匀，且保持各向同性，本文选定热量从外界环境垂直于皮肤表面进行传递；

Figure 1. Model Abstraction diagram

图1. 模型抽象示意图

假设3：系统传热过程忽略水汽、汗液的影响 [1] ；

假设4：假设每层材料介质的热传导率在各个方向相同。

其次，对一类三织物材料层的作业服抽象为一维模型，见图1。将一名身着高温作业服的工作人员抽象化为一个多均匀质地层的圆柱体，在空间上，高温作业服和人体是三维的。但对于本传热问题，基于模型假设2且无其他不均匀热源及传热过程，因而可将三维模型简化为一维传热模型，仅研究热量由作业服外层到皮肤表面的传热过程，并建立坐标系。

3.1. “模型一”的建立

对于热传导、热辐射、热对流三种基本传热方式，由表1可获知空气层厚度最大为6.4 mm (<8 mm)，间隙太小无法形成对流运动，这时空气层以热传导为主 [2] 。由于作业服第I层织物材料阻挡了大部分外界高温环境的辐射传热，因此在I至III层织物材料中的传热过程可以忽略热辐射 [1] 。作业服的温度不会超过外界温度，此时辐射传递的热量较对流和传导的热量小，在考虑第IV层的传热方式中可以忽略热辐射 [3] 。人体皮肤下毛细血管中存在大量流动的血液，所以在空气层与人体皮肤的交界 ${\tau }_4$ 以热对流为主 [1] [2] 。根据上述对每一层织物材料(包括空气层和交界面)传热方式的分析(见表2)，本文所需要建立的作业服传热模型为一维分段常系数热传导方程，在热扩散系数为间断点处添加交界条件使温度分布函数 $u\left(x,t\right)$ 连续。

Table 2. Main heat transfer mode of fabric material and air layer

表2. 织物材料与空气层的主要传热方式

Figure 2. Diagram of the Thin bar

图2. 细杆示意图

设温度分布函数 $u\left(x,t\right)$ 关于时间t有一阶连续的导数，关于空间x有二阶连续的导数，在推导一维非稳态热传导方程之前，本文选择了一个单一介质、质地均匀的细杆(见图2)作为研究对象，考虑其热量传播过程。根据热量计算公式

$\Delta Q=cm\Delta u$ ,

其中c为细杆的比热容，m为体积微元的质量， $\Delta u$ 为温度变化，在温度场 $u\left(x,t\right)$ 中，时间间隔 $\Delta t$ 内，体积微元 $\Delta V$ (长度微元为 $\Delta x$ )吸收的热量为

$\Delta Q=c\rho S\Delta x\left(u\left(x,t+\Delta t\right)-u\left(x,t\right)\right)$ ,

当 $\Delta t\to 0_{+}$ 时，

$\Delta Q=c\rho S\Delta x\frac{\partial u}{\partial t}\left(x,t\right)\Delta t$ .(1)

根据Fourier实验定律可知，当物质内外存在温度差时，热量从高温侧流向低温侧，时间间隔 内，流过一个面积为 的截面的热量与截面外法线方向的温度变化成正比 [4] ，即

$\Delta Q=-\lambda \frac{\partial u}{\partial e_n}\left(x,t\right)S\Delta t$ ,

其中 $e_n$ 为截面的单位外法向量， $\lambda$ 为细杆在x处的热传导系数，负号表示热量从温度高的一侧流向温度低的一侧。参照假设2，本文研究的一维热传导问题中 $\frac{\partial u}{\partial e_n}=\frac{\partial u}{\partial x}$ ，由此可知热流密度(单位时间流过单位面积的热量)为

. (2)

结合上述分析，可以得出以下两条结论：

其一，在时间间隔 $\Delta t$ 内，细杆x处流入的热量为

$Q_{in}=-\lambda \frac{\partial u}{\partial x}\left(x,t\right)S\Delta t$ ; (3)

其二，在时间间隔 $\Delta t$ 内，细杆 $x+\Delta x$ 处流出的热量为

$Q_{out}=-\lambda \frac{\partial u}{\partial x}\left(x+\Delta x,t\right)S\Delta t$ . (4)

由(3)式与(4)式可得

$\Delta Q=Q_{in}-Q_{out}=-\lambda \left(\frac{\partial u}{\partial x}\left(x,t\right)-\frac{\partial u}{\partial x}\left(x+\Delta x,t\right)\right)S\Delta t$ . (5)

根据能量守恒定律，联立(1)式和(5)式可得

$c\rho \Delta x\frac{\partial u}{\partial t}\left(x,t\right)=\lambda \left(\frac{\partial u}{\partial x}\left(x+\Delta x,t\right)-\frac{\partial u}{\partial x}\left(x,t\right)\right)$ ,

移项后有

$\frac{\partial u}{\partial t}\left(x,t\right)=\frac{\lambda }{c\rho }\frac{\left(\frac{\partial u}{\partial x}\left(x+\Delta x,t\right)-\frac{\partial u}{\partial x}\left(x,t\right)\right)}{\Delta x}\triangleq a\frac{\left(\frac{\partial u}{\partial x}\left(x+\Delta x,t\right)-\frac{\partial u}{\partial x}\left(x,t\right)\right)}{\Delta x}$ ,

其中a称为热扩散率，当 $\Delta x\to 0_{+}$ 时

$\frac{\partial u}{\partial t}\left(x,t\right)=a\frac{{\partial }^{2}u}{\partial x^2}\left(x,t\right)$ .

本文的研究对象为3层织物材料的隔热服，其一维热传导方程具体为

$\frac{\partial u}{\partial t}=a_i\frac{{\partial }^{2}u}{\partial x^2},\left(x,t\right)\in \left(X_i,X_{i+1}\right)\times (0,T],a_i=\frac{{\lambda }_i}{c_i\Delta {\rho }_i}\left(i=1,2,3,4\right)$ ,

其中 ${\lambda }_i,c_i,{\rho }_i$ 分别为第i层织物材料的热传导率、比热和密度，T为总时长。根据表2给出的各层传热分析，对方程做出如下限定条件：

(A1)初始条件：假设进入高温环境时，人体与作业服已达到稳定状态，作业服温度分布的初始值为假人温度37℃，即 $u\left(x,0\right)=37\left(x\in \left[X_1,X_5\right]\right)$ 。

(A2)边界条件：牛顿冷却定律表明，当物体表面与周围存在温度差时，单位时间从单位面积散失的热量与温度差成正比 [5] 。在传热系统左边界 $x=X_1$ 和右边界 $x=X_5$ 分别应用之，结合式(2)得到两个第三类边界条件，即

$-{\lambda }_1\frac{\partial u}{\partial x}\left(X_1,t\right)=h_1\left(u_{env}-u\left(X_1,t\right)\right)$ ,

${\lambda }_4\frac{\partial u}{\partial x}\left(X_5,t\right)=h_2\left(u\left(X_5,t\right)-37\right)$ ,

其中 $u_{env}$ 为外界高温环境温度， $h_1$ 与 $h_2$ 为对流换热系数。

(A3)交界条件：假设材料间接触良好，忽略接触热阻，在交界面上温度与热流密度连续，结合式(2)，可得

$\{\begin{array}{l}u\left(X_i-0,t\right)=u\left(X_i+0,t\right)\\ {\lambda }_i\frac{\partial u}{\partial x}\left(X_{i+1}-0,t\right)={\lambda }_{i+1}\frac{\partial u}{\partial x}\left(X_{i+1}+0,t\right)\end{array},i=1,2,3$ .

综合上述的讨论，本文建立“一维非稳态作业服传热模型”，即

$\{\begin{array}{l}\frac{\partial u}{\partial t}=a_i\frac{{\partial }^{2}u}{\partial x^2},\left(x,t\right)\in \left(X_i,X_{i+1}\right)\times \left(0,T\right),a_i=\frac{{\lambda }_i}{c_i\Delta {\rho }_i}\left(i=1,2,3,4\right)\\ \left(初值条件\right)u\left(x,0\right)=37\left(x\in \left[X_1,X_5\right]\right)\\ \left(边界条件\right)\{\begin{array}{l}-{\lambda }_1\frac{\partial u}{\partial x}\left(X_1,t\right)=h_1\left(u_{env}-u\left(X_1,t\right)\right)\\ {\lambda }_4\frac{\partial u}{\partial x}\left(X_5,t\right)=h_2\left(u\left(X_5,t\right)-37\right)\end{array}\\ \left(交界条件\right)\{\begin{array}{l}u\left(X_i-0,t\right)=u\left(X_i+0,t\right)\\ {\lambda }_i\frac{\partial u}{\partial x}\left(X_{i+1}-0,t\right)={\lambda }_{i+1}\frac{\partial u}{\partial x}\left(X_{i+1}+0,t\right)\end{array}\left(i=1,2,3\right)\end{array}$ .

3.2. “模型二”、“模型三”的建立

在n层织物材料厚度未知的情况下，不妨设每一层织物材料的厚度为 $d_i\left(i=1,2,\cdots ,n\right)$ ，记空气层厚度为 $d_{n+1}$ ，根据表1所给数据可知，作业服的部分织物材料层厚度是固定的，为此构造示性函数

${\delta }_j=\{\begin{array}{l}0,j\notin A\\ 1,j\in A\end{array}$ ,

其中A为指标集，对于 $\forall j\in A$ ，织物材料层厚度 $d_j$ 是可变的。在构建目标函数时，可以只考虑集合 中的层指标，构建刻画生产高温作业服所需成本的目标函数 $v\left(d_1,d_2,\cdots ,d_{n+1}\right)$ 为

$v=p_1{\delta }_1{d}_1+p_2{\delta }_2{d}_2+\cdots +p_n{\delta }_n{d}_n$ ,

其中 $p_j$ 为第j层织物材料的单价(元/mm)且 $p_j>0$ ；构建刻画高温作业服重量大小的目标函数 $w\left(d_1,d_2,\cdots ,d_{n+1}\right)$ 为

$w={\rho }_1{\delta }_1{d}_1+{\rho }_2{\delta }_2{d}_2+\cdots +{\rho }_n{\delta }_n{d}_n+{\rho }_{n+1}{\delta }_{n+1}{d}_{n+1}$ ,

其中 ${\rho }_i$ 为第i织物材料层的密度；构建刻画高温作业服薄厚程度的目标函数 $h\left(d_1,d_2,\cdots ,d_{n+1}\right)$ 为

$h={\delta }_1{d}_1+{\delta }_2{d}_2+\cdots +{\delta }_n{d}_n+{\delta }_{n+1}{d}_{n+1}$ .

本文研究的问题2与问题3中，指标集A是不同的，表3中对具体问题进行具体分析。

Table 3. Objective function of establishing multi-objective programming model

表3. 建立多目标规划模型的目标函数

根据表3，建立“模型二”如下

$\begin{array}{l}\min p_2{d}_2\\ \min {\rho }_2{d}_2\\ \min d_2\\ \text{s}.\text{t}.\{\begin{array}{l}u\left(X_5,3600\right)\le 47\\ u\left(X_5,3300\right)\le 44\\ 0.6\le d_2\le 25\end{array}\end{array}$ $\iff$ $\begin{array}{l}\min d_2\\ \text{s}.\text{t}.\{\begin{array}{l}u\left(X_5,3600\right)\le 47\\ u\left(X_5,3300\right)\le 44\\ 0.6\le d_2\le 25\end{array}\end{array}$

同理建立“模型三”如下

$\begin{array}{l}\min p_2{d}_2\\ \min {\rho }_2{d}_2+{\rho }_4{d}_4\\ \min d_2+d_4\\ \text{s}.\text{t}.\{\begin{array}{l}u\left(X_5,3600\right)\le 47\\ u\left(X_5,3300\right)\le 44\\ 0.6\le d_2\le 25\\ 0.6\le d_4\le 6.4\end{array}\end{array}$ $\iff$ $\begin{array}{l}\min d_2\\ \min {\rho }_2{d}_2+{\rho }_4{d}_4\\ \min d_2+d_4\\ \text{s}.\text{t}.\{\begin{array}{l}u\left(X_5,3600\right)\le 47\\ u\left(X_5,3300\right)\le 44\\ 0.6\le d_2\le 25\\ 0.6\le d_4\le 6.4\end{array}\end{array}$

4. 模型的求解

4.1. 问题1求解

常见的热传导方程数值解法有有限差分法、蒙特卡罗法、有限元分析法和分子动力学模拟等 [4] ，本文采用有限差分算法进行求解。首先对非稳态传热模型求解区域进行网络剖分，可以剖分成正交网络、光滑网络、Kershaw网络、随机网络等 [6] 。本文将求解区域剖分成正交网络，在t方向上取剖分 $0=t_1<t_2<\cdots <t_n=T$ ，在x方向上取剖分 $X_1=x_1<x_2<\cdots <x_m=X_5$ ，其中存在正整数 $z,s,l$ 满足 $x_z=X_2,x_s=X_3,x_l=X_4$ 。其次对“模型一”进行离散化，在文献 [7] [8] 中，作者系统地介绍了一维常系

数热传导方程 $\frac{\partial u\left(x,t\right)}{\partial t}=a\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}\left(x\in R,t>0\right)$ 的差分格式，在表4中列举了几种经典的差分格式解的稳

定性条件和截断误差。不同的差分格式具有不同的优缺点，如“跳点格式”可以节省存储和计算量，“Saul’yev不对称格式”采用取平均值的方法构造差分格式可以抵消截断误差 [7] 。

Table 4. Stability and truncation error of several differential format solutions, of which r = Δ t Δ x 2

表4. 几种差分格式解的稳定性及截断误差，其中 $r=\frac{\Delta t}{\Delta x^2}$

运用差分格式近似求解偏微分方程时，综合考虑解的稳定性和截断误差精度是必要的，结合表4，本文选用Crank-Nicolson格式，即权重为0.5的加权隐式差分格式

$\frac{u_j^k-u_j^{k-1}}{\Delta t}=\frac{a}{2}\left(\frac{u_{j+1}^k-2u_j^k+u_{j-1}^k}{\Delta x^2}+\frac{u_{j+1}^{k-1}-2u_j^{k-1}+u_{j-1}^{k-1}}{\Delta x^2}\right)$

对热传导方程离散化，其中 $u_j^k$ 表示 $u\left(x_j,t_k\right)$ 。定义时间一阶向后差商算子 ${\partial }_{t,\Delta t}f\left(t\right)=\frac{f\left(t+\Delta t\right)-f\left(t\right)}{\Delta t}$ ，空间一阶向后差商算子 ${\partial }_{x,\Delta x}f\left(x\right)=\frac{f\left(x+\Delta x\right)-f\left(x\right)}{\Delta x}$ ，对“模型一”的边界条件与交界条件离散化，得到离散温度分布模型如下

$\{\begin{array}{l}-a_{i}ru\left(x_{j-1},t_k\right)+\left(2a_{i}r+2\right)u\left(x_j,t_k\right)-a_{i}ru\left(x_{j+1},t_k\right)=a_{i}ru\left(x_{j-1},t_{k-1}\right)+\left(2-2a_{i}r\right)u\left(x_j,t_{k-1}\right)+a_{i}ru\left(x_{j+1},t_{k-1}\right)\\ \left(初值条件\right)u\left(x_j,0\right)=37\\ \left(边界条件\right)\{\begin{array}{l}\left(\frac{{\lambda }_1}{\Delta x}+h_1\right)u\left(x_1,t_k\right)-\frac{{\lambda }_1}{\Delta x}u\left(x_2,t_k\right)=u_{env}{h}_1\\ \left(\frac{{\lambda }_4}{\Delta x}+h_2\right)u\left(x_m,t_k\right)-\frac{{\lambda }_4}{\Delta x}u\left(x_{m-1},t_k\right)=37h_2\end{array}\\ \left(交界条件\right)-\frac{{\lambda }_i}{\Delta x}u\left(x_{B_i-1},t_k\right)+\frac{{\lambda }_i+{\lambda }_{i+1}}{\Delta x}u\left(x_{B_i},t_k\right)-\frac{{\lambda }_{i+1}}{\Delta x}u\left(x_{B_i+1},t_k\right)=0\left(i=1,2,3\right)\end{array}$

其中 $x_{B_i}$ 为织物材料i与 $i+1\left(i=1,2,3\right)$ 层交界处， $x_m$ 为右边界点。综上，求解隔热服的温度分布可以归结为求解系数矩阵A为

$\left(\begin{array}{ccc}\begin{array}{ccc}\frac{{\lambda }_1}{\Delta x}+h_1& -\frac{{\lambda }_1}{\Delta x}& \\ -a_{1}r& 2a_{1}r+2& -a_{1}r\\ & -a_{1}r& 2a_{1}r+2\end{array}& \begin{array}{ccc}& & \\ & & \\ & & \end{array}& \begin{array}{ccc}& & \\ & & \\ & & \end{array}\\ \begin{array}{ccc}& & \ddots \\ & & \\ & & \end{array}& \begin{array}{ccc}\ddots & & \\ -\frac{{\lambda }_1}{\Delta x}& \frac{{\lambda }_1+{\lambda }_2}{\Delta x}& -\frac{{\lambda }_2}{\Delta x}\\ & \ddots & \ddots \end{array}& \begin{array}{ccc}& & \\ & & \\ & & \end{array}\\ \begin{array}{ccc}& & \\ & & \\ & & \end{array}& \begin{array}{ccc}& & \\ & & \\ & & \end{array}& \begin{array}{ccc}2a_{4}r+2& -a_{4}r& \\ -a_{4}r& 2a_{4}r+2& -a_{4}r\\ & -\frac{{\lambda }_4}{\Delta x}& \frac{{\lambda }_4}{\Delta x}+h_2\end{array}\end{array}\right)$

的三对角形线性方程组 $A\gamma \left(u\left(t_{k+1}\right)\right)=\gamma \left(u\left(t_k\right)\right)$ ，其中 $\gamma \left(u\left(t_k\right)\right)$ 表示一个关于 $u\left(t_k\right)$ 的函数，求解过程可以表示为如下示意过程

$\left(\begin{array}{c}\begin{array}{c}37\\ 37\end{array}\\ \begin{array}{c}⋮\\ 37\end{array}\end{array}\right)\stackrel{解方程}{\to }\left(\begin{array}{c}\begin{array}{c}u\left(x_1,t_1\right)\\ u\left(x_2,t_1\right)\end{array}\\ \begin{array}{c}⋮\\ u\left(x_m,t_1\right)\end{array}\end{array}\right)\stackrel{解方程}{\to }\cdots \stackrel{解方程}{\to }\left(\begin{array}{c}\begin{array}{c}u\left(x_1,t_n\right)\\ u\left(x_2,t_n\right)\end{array}\\ \begin{array}{c}⋮\\ u\left(x_m,t_n\right)\end{array}\end{array}\right)$

由于方程组系数矩阵中含有未知参数 $h_1$ 与 $h_2$ ，在求解温度分布之前需要对其进行估计。

设 $e_r\left(h_1,h_2\right)$ 为在75℃高温环境下工作90分钟后人体皮肤外侧温度相对误差向量，本文首先对未知参数 $h_1$ 与 $h_2$ 有约束非线性规划问题

，

利用Matlab优化工具箱fmincon函数进行求解 [9] ，得到表5所示结果。

其次，选择对取值区间 $\left(h_1,h_2\right)\in \left[100,140\right]\times \left[8,8.6\right]$ 尝试遍历求解，本文分两步进行探索。

第一步，对参数 $h_1$ 以1为步长在 $\left[100,140\right]$ 取值、对参数 $$ 以0.01为步长在 $\left[8.1,8.6\right]$ 取值进行精确搜索，绘制 $h_1-h_2-e_r$ 图(见图3)可知 $h_1=114$ ， $h_2=8.35$ 为最优解，且在求解区域内最优值存在于子区域 $\left(h_1,h_2\right)\in \left(113,115\right)\times \left(8.34,8.36\right)$ 。

第二步，对参数 $h_1$ 以0.01为步长、对参数 $h_2$ 以0.001为步长进行搜索，并绘制 $h_1-h_2-e_r$ 图(见图3)可知 $h_1=113.69$ ， $h_2=8.350$ 为最优解，此时 ${\Vert e_r\left(h_1,h_2\right)\Vert }_2=0.006978$ 。

通过遍历方法得到高精度最优值所需计算量是庞大的，根据表5所示信息可得知，建立有约束非线性规划求解未知参数 $h_1$ 与 $h_2$ 既高效又精确，最优解 $h_1=113.62$ ， $h_2=8.35$ (保留小数点后两位)。

对于三对角形可逆系数矩阵的线性方程组 $Ax=b$ 求解可直接得到 $x=A^{-1}b$ ，但实际问题中的数据规模较大且系数矩阵A中除了主对角线和次对角线上元素非零，其余元素均为零，为了节省程序运行耗时和内存，本文采用追赶法对三对角形线性方程组 $Ax=b$ 求解 [10] ，步骤如下。

Table 5. Optimal solution h 1 and h 2 obtained by constrained nonlinear programming

表5. 有约束非线性规划求解所得最优解 $h_1$ 与 $h_2$

(a) (b)(c) (d)

Figure 3. Contour map obtained from the traversal search and $h_1-h_2-{\Vert e_r\left(h_1,h_2\right)\Vert }_2$ figure

图3. 遍历搜索所得等值线图和 $h_1-h_2-{\Vert e_r\left(h_1,h_2\right)\Vert }_2$ 图

第一步，将系数矩阵A进行LU分解使得 $A=L\cdot U$ ，其中L为上三角矩阵，U为下三角矩阵，从而方程组 $Ax=b$ 转化为 $LUx=b$ ；

第二步，求解方程 $Ly=b$ ；

第三步，求解方程 $Ux=y$ 。

求解过程中对 $\Delta x$ 与 $\Delta t$ 选取不同值计算所得人体皮肤外侧温度相对误差存在差异，当 $\Delta x$ 和 $\Delta t$ 较大时，截断误差较大；当 $\Delta x$ 和 $\Delta t$ 较小时，计算次数增多，累积的舍入误差较大。表6列举了几种步长选择的具体情况，发现选取 $\Delta x=0.1\text{mm}$ 与 $\Delta t=1\text{s}$ 时所得相对误差和程序运行耗时较为理想。图4(a)描绘了 $h_1=113.62$ ， $h_2=8.35$ 、 $\Delta x=0.1\text{mm}$ 、 $\Delta t=1\text{s}$ 时求解所得相对误差与人体皮肤外侧温度计算结果，图4(b)展示了75℃高温环境下、工作时长为90分钟时作业服与空气层的温度分布。

Table 6. Where Δ x units are in millimeters and Δ t units are in seconds

表6. 其中 $\Delta x$ 的单位为毫米， $\Delta t$ 的单位为秒

Figure 4. (a) Calculation results and errors of the lateral temperature of human skin. (b) Calculation results of temperature distribution

图4. (a)人体皮肤外侧温度计算结果与误差；(b)温度分布计算结果

4.2. 问题2求解

设第II织物材料层厚度为 $d_2$ ，第3300、3360秒人体皮肤外侧温度为 $u_{3300}\left(d_2\right)$ 和 $u_{3600}\left(d_2\right)$ ，根据表1可得 $d_2\in \left[0.6,25\right]$ ，易知 $u_{3300}\left(d_2\right)$ 与 $u_{3600}\left(d_2\right)$ 关于 $d_2$ 单调递减，即

$u_{3600}\left(d_2\right)\le u_{3600}\left(0.6\right)=45.101124$

图6(a)展示了以 $\Delta d_2=0.01\text{mm}$ 为步长遍历区间 $\left[0.6,25\right]$ 的计算结果，描绘了 $d_2$ 与 $u_{3300}\left(d_2\right)$ 和 $u_{3600}\left(d_2\right)$ 的关系，验证了上述理论分析。由此可知 $d_2$ 在区间 $\left[0.6,25\right]$ 上有 $u_{3600}\left(d_2\right)\le 47$ 恒成立，将“模型二”转化为

$\begin{array}{l}\min d_2\\ \text{s}\text{.t}\text{.}\{\begin{array}{l}u\left(X_5,3300\right)\le 44\\ 0.6\le d_2\le 25\end{array}\end{array}$

依据问题1的离散温度分布模型，本节以 $e=0.01$ 为搜索区间长度限制，采用Fibonacci法、黄金分割法、二分法和二次插值法4种方法逐一对“模型二”进行求解 [11] ，表7给出一些主要数值计算结果，图5展示了各种方法迭代区间变化情况。

Figure 5. Diagram of iterative interval of 4 search methods

图5. 4种搜索方法迭代区间示意图

Table 7. Calculation results of Model Two

表7. 模型二计算结果

从表7中的计算结果可知，单层织物材料厚度优化问题中，既要保证解的可靠性，又要尽可能节省运行时间，在对 $d_2$ 精确到小数点后1位时，优先选用二次插值法；在对 $d_2$ 精确到小数点后2位时，优先选用二分法。在与表1数据精度一致，即保留小数点后1位时，最优解 $d_2=17.6\text{mm}$ ，当外界温度为65℃时，高温作业时长在60分钟时人体皮肤外侧温度变化如图6(b)所示，超过44℃的时长为273秒，且最高温度为44.071775℃。

Figure 6. (a) The thickness of the II fabric material layer is $d_2$ at 0.01 for 3,300 seconds and 3,600 seconds. (b) The temperature change of the lateral temperature of the human skin within 60 minutes of the operation of the II fabric material layer thickness $d_2=17.55\text{mm}$

图6. (a)第II织物材料层厚度 $d_2$ 以0.01 mm为步长进行逐一求解3300秒和3600秒温度。(b)第II织物材料层厚度 $d_2=17.55\text{mm}$ 时人体皮肤外侧温度在作业60分钟内的温度变化

4.3. 问题3求解

为了求解多目标规划问题

$\begin{array}{l}\min d_2\\ \min {\rho }_2{d}_2+{\rho }_4{d}_4\\ \min d_2+d_4\\ \text{s}.\text{t}.\{\begin{array}{l}u\left(X_5,3600\right)\le 47\\ u\left(X_5,3300\right)\le 44\\ 0.6\le d_2\le 25\\ 0.6\le d_4\le 6.4\end{array}\end{array}$ ，

基于问题2的求解过程，本节首先从“减少约束条件”入手。记 $u_{1800}\left(d_2,d_4\right)$ 与 $u_{1500}\left(d_2,d_4\right)$ 分别表示第1800秒和第1500秒时的人体皮肤外侧温度，因为 $\frac{\partial u_{1800}}{\partial d_2}\le 0$ 、 $\frac{\partial u_{1800}}{\partial d_4}\le 0$ ，故有

$u_{1800}\left(d_2,d_4\right)\le u_{1800}\left(0.6,0.6\right)=58.553266$ (6)

成立，同理可得

$u_{1500}\left(d_2,d_4\right)\le u_{1500}\left(0.6,0.6\right)=58.553266$ (7)

成立。综合式(6)与(7)可知在 $\left(d_2,d_4\right)\in \left[0.6,25\right]\times \left[0.6,6.4\right]$ 中， $u_{1800}\left(d_2,d_4\right)\le 47$ 或 $u_{1500}\left(d_2,d_4\right)\le 44$ 均不恒成立。

在实际问题中，对成本、重量与厚度目标函数的优先考虑顺序是不同的。对于完全分层多目标规划问题 ${\min }_{x\in D}{\left[P_i\left(F_i\left(x\right)\right)\right]}_{i=1}^s$ ，其中 $s\ge 2$ ， $F_i\left(x\right)\left(i=1,2,\cdots ,s\right)$ 是目标函数， $P_i\left(i=1,2,\cdots ,s\right)$ 是优先层次的符

号，表示其后面括号内的目标函数 $F_i\left(x\right)$ 属于第i优先层次，且满足 $P_i\gg P_{i+1}\left(i=1,2,\cdots ,s-1\right)$ ，符号“ $\gg$ ”表示第i优先层“优先于”第 $i+1$ 优先层。在文献 [11] 中，作者提出了一种完全分层序列法，并证明了简单完全分层序列法可以得到一个有效最优解，宽容完全分层序列法只能得到一个弱有效最优解。本文依据简单完全分层序列法对模型三进行求解，其算法流程如下。

第一步，确定初始可行域，将原问题 ${\min }_{x\in D}{\left[P_i\left(F_i\left(x\right)\right)\right]}_{i=1}^s$ 的可行域作为第一优先层的可行域，即 $D^1=D$ ，令 $k=1$ ；

第二步，在第k优先层的可行域 $D^k$ 上求解第k优先层目标函数 $F_k\left(x\right)$ 的优化问题 ${\min }_{x\in D^k}{F}_k\left(x\right)$ ，得到最优解 $x^k$ 和最优值 $F_k\left(x^k\right)$ ；

第三步，若 $k=m$ ，输出 $x^k$ ；若 $k<m$ ，转第四步；

第四步，建立下一层的可行域 $D^{k+1}=\left\{x\in D^k|F_k\left(x\right)\le F_k\left(x^k\right)\right\}$ ， $k=k+1$ ，转第二步。

在问题三中， $k=3$ ，本文分6种分类情况进行求解模型三。确定初始可行域时，本文首先定性分析。当 $d_4=6.4\text{mm}$ 时，一定存在一个最小的 $d_2$ 满足约束条件，若不然，模型三的解集为空集，同理可得 $d_2=25\text{mm}$ 也符合上述分析。对于 $\forall d_4\in \left[0.6,6.4\right]$ ，若存在满足模型三的 $d_2$ ，则可以确定唯一一个最小的 $d_{2\min }\in \left[0.6,25\right]$ ；若不存在，则 $d_4$ 为初始可行域 $D=\left\{\left(d_4,d_2\right)\right\}$ 中 $d_4$ 分量的一个下界，从而可以确定 $d_4$ 分量存在一个下确界 $\overline{d_4}$ ，同理可知 $d_2$ 分量存在下确界 $\overline{d_2}$ 。不妨设在可行域上存在边界曲线 $S\left(x\right)$ 使得 $d_{2\min }=S\left(d_4\right)$ ，易知 $S\left(x\right)$ 是单调递减的，进一步可以得到可行域D的具体范围为 $D=\left\{\left(d_4,d_2\right)|\left[\overline{d_4},6.4\right]\times \left[S\left(x\right),25\right]\right\}$ 。通过以0.1 mm为步长对区间 $\left[0.6,6.4\right]\times \left[0.6,25\right]$ 遍历，计算初始可行域D的范围，图7(a)将计算结果可视化，符合上述分析。

对于计算 $S\left(x\right)$ 上的可行点，本文提出了一种有效缩短求解区间的可行边界点计算方法，具体流程见图8。对 $d_4$ 在区间 $\left[\overline{d_4},6.4\right]$ 上以0.1 mm为步长遍历，求解可行点耗时节约80%。结合简单完全分层序列法，在表8中给出了不同分类情况时模型三的最优解。为了与表1数据精度一致，所得结果保留小数点后1位，即最优解为 $d_2=19.2\text{mm}$ ， $d_4=6.4\text{mm}$ 。

Figure 7. (a) Initial feasible domain; (b) Number of iterations of new algorithms in the application of different one-dimensional search methods

图7. (a)初始可行域(b)新算法在不同一维搜索方法的应用下迭代次数

Table 8. Model Three Solution results

表8. 模型三求解结果

Figure 8. Algorithm flowchart

图8. 算法流程图

致谢

感谢此次研究的指导老师吴宜均老师，对我的研究提出中肯的建议。对所有给予转载和引用权的资料、图片、文献、研究思想和设想的所有者表示感谢！