PDE Toolbox和最小二乘法在导热问题中的应用

朱秋旻 陈晓平

(江苏大学电气信息工程学院，江苏 镇江 212013)

0 引言

热传导方程式描述了一个区域内的温度如何随时间变化。由于许多模型没有解析解，因此必须以数值方法计算模型给出的定价［1］。这里直接用Matlab软件求解，并给出图形。

法国科学家勒让德于1806年独立发现了最小二乘法(least square，LS)［2］。从此，这个方法被广泛应用于各行各业的工程计算中。在Matlab中可以编写最小二乘法的应用程序进行系统参数辨识［3］。

偏微分方程工具箱提供了研究和求解空间二维偏微分方程问题的一个强大而又灵活实用的环境［4］，用它可方便地求解热方程，得出金属体的内部特性。

1 热传导方程的解法

1.1 热传导方程介绍

热传导方程由热力学第一定律和傅里叶定律推导得到。为简单起见，现只讨论单纯热传导现象，并忽略体积膨胀，同时假设金属板静止，可得到如下方程:

式中:T为某时某点温度;t为时间;λ为导热率，通常为常数;ρ为密度;cv为等体比热容;▽为拉普拉斯算子。此公式表示某点温度随时间的变化与这一点温度的散度成正比。

由于抛物型偏微分方程标准形式为:

令 d=1、c=k、a=0、f=0，式(2)就成了所要求的热方程式(1)了。

方程边界条件一般有Dirichlet条件和Neumann条件两种，即:

式中:h为空间采样步长;u为方程的未知函数;n为边界切面的法向量;h、r、c、q、g 都为已知函数或常数。

1.2 解金属板导热问题

考虑一个带有矩形孔的金属板上的热传导问题。板左边的温度保持在100℃，板右边热量从板向环境空气定常流动，其他边及内孔边界保持绝缘。

当初始板温t=0℃时，则可概括为如下定解问题:

上述问题中金属板域的外边界顶点坐标分别为(－0.5，－0.8)、(0.5，－ 0.8)、(0.5，0.8)、(－ 0.5，0.8);内边界 顶点坐 标分别为(－ 0.05，－ 0.4)、(0.05，－0.4)、(0.05，0.4)、(－ 0.05，0.4)。接着使用Matlab图形用户界面(GUI)求解这一问题。在PDE Toolbox窗口的工具栏中选择Generic Scalar模式，然后经过区域设置、边界条件设置、方程类型设置、网格剖分、初值和误差设置，可以得到此问题的数值解和解的图形。

部分Matlab程序如下。

1.3 热传导方程的差分格式

热方程为偏微分方程，求解它可用求解微分方程的数值方法。

考虑上述二维常系数热传导方程为:

要求解这个方程，常用差分格式的递推公式，其分为古典显格式和古典隐格式两种。

1.3.1 古典显格式

取网格的空间步长为h=1/N，时间步长为τ。在(xj，tk)处时间用向前差商，空间用中心差商近似，则可得到古典显格式，即:

式中:局部截断误差为O(τ+h2)，r=τ/h2(网格比)，其中，O表示高阶无穷小。当且仅当r＜1/4时，格式是稳定的。

1.3.2 古典隐格式

在(xj，tk+1)处时间用向后差商，空间用中心差商近似，则可得古典隐格式，即:

式中:局部截断误差为 O(τ+h2)，r=τ/h2，隐格式对任何网格比都是稳定的。

用差分方法同样可以在Matlab上求解热传导方程，在键入程序后，就可得到板上各点的温度值，对于更复杂的问题，此方法有明显的优势。

2 最小二乘原理与算法

2.1 最小二乘基本原理

最小二乘法(LS)是用于参数估计的数学方法，它使数学模型在误差平方和最小的意义上拟合试验数据，也是一种涉及较少数学基础而又被大量应用的一种基本方法。

最小二乘法提供一个估算方法，用来得到一个在最小方差意义上与试验数据最好拟合的数学模型。由最小二乘法获得的估计在一定的条件下有最佳的统计特性，即估计的结果是无偏的、一致的和有效的。

因为LS的原则是希望某个量S(0)和a的值能使观测值和由模型的计算值之间的误差为最小，所以各次观测误差可表示为:

式中:i=1，2，…，l。

整个观测过程的误差是由各次观测误差所组成的，采用每个误差的平方和作为总误差:

所选的误差平方和函数J就是估计参数时所采用的性能指标。显然，J越小越好，即所选的S(0)和a的值能使每个误差的平方和J的值最小。由于平方运算也称二乘运算，因此称上述方法为最小二乘估计法。要使J达到极小值，只需分别对S(0)和a求偏导，令它们等于0。这样就可以得到关于S(0)和a的两个估计表达式。

2.2 最小二乘的基本关系式

对于单输入单输出离散时间动态系统，设u(k)为输入、y(k)为输出、z(k)为量测输出、v(k)为白噪声，其系统数学关系可用如下随机差分方程描述，即:

利用数据序列{u(k)}、{z(k)}，极小化下列准则函数。

在文献［3］中，根据所求解的问题的J不同，在不同场合下，函数J往往有不同的名称，它是一个标量。对J求导，令其等于0可得:

式(14)就是所要求的参数估计式，从统计学的角度考虑。观测总次数l必须大大超过设定的未知参数的数目，此时由观测所提供的方程式的数目才能超过确定出方程组的唯一解所需的数目。

3 基于最小二乘算法的仿真

3.1 用积矩矩阵估计阶次

考虑干扰是白噪声的情况，令:

当n从1开始逐一增加时，若上式值有明显的增加，则在显著增加处的n值就是模型的阶次。由仿真得出本文所讨论的问题的阶次大约在2的附近，所以本问题的差分方程模型的阶次选为2。

3.2 模型的参数辨识

程序框图如图1所示。

图1 最小二乘算法程序框图Fig.1 Block diagram of the LS

根据试验得到的数据，将本问题的模型阶次选为2，并用线性差分方程来描述，可得:

式中:v(k)为服从正态分布的白噪声N(0，1);z(k)为金属板中心温度与初始温度(高于室温)之差;u(k)为1，表示金属板两端加热源，为－1，表示不加热源。输入信号采用4阶M序列，幅度为1。则按式(17)构造zl和HL，数据长度 L=14，加权阵取 I，并利用式(14)计算参数估计值。

根据金属板的ARMX模块估计结构，由程序运行可得如图2所示的温度输出曲线和输入输出误差曲线图。

图2 ARMX模块估计输出曲线和误差曲线Fig.2 Output curve and error curve of ARMX estimated module

由图2可以看出，温床系统的估计输出曲线与量测温度值相符，且误差曲线也符合要求。

图3给出了一定点温度在阶跃功率输入下的输出变化，随时间输出温度趋于稳定，反映了金属板的热传导性质。从图4所示的频率响应曲线可以看出，在低频信号下输出有较好的稳定性，当外界有高频噪声时，采用ARMX模块估计就会出现较大的偏差。

4 结束语

从仿真结果看，本文所采用的方法能较好地描述金属板的内、外部特性，得到了所需系统模型，较好地解决了上面提出的问题，可以满足实际的要求。在更复杂的环境中，还应对噪声和信号进行分析，以便更精确地获得系统模型的其他数据。

［1］王竹溪.热力学［M］.2版.北京:北京大学出版社，2005.

［2］John F.Partial differential equations［M］.Springer，1982.

［3］陆君安，尚涛，谢进，等.偏微分方程的 MATLAB解法［M］.武汉:武汉大学出版社，2001.

［4］侯媛彬，汪海，王立琦.系统辨识及其MATLAB仿真［M］.北京:科学出版社，2004.

［5］Iserles A.微分方程数值分析基础教程［M］.刘晓艳，刘学深，译.北京:清华大学出版社，2005.

［6］Hou Yuanbin.A decoupling control method with improving genetic algorithm［C］∥Proceedings of 2002 International Conference on Machine Learning and Cybernetics，China，2002:2112 －2115.

［7］周彤.含有不确定因素的模型检验及其非伪概率估计［J］.控制理论与应用，1996(2):145－152.