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.