朱秋旻 陳曉平

(江蘇大學(xué)電氣信息工程學(xué)院，江蘇 鎮(zhèn)江 212013)

0 引言

熱傳導(dǎo)方程式描述了一個(gè)區(qū)域內(nèi)的溫度如何隨時(shí)間變化。由于許多模型沒有解析解，因此必須以數(shù)值方法計(jì)算模型給出的定價(jià)［1］。這里直接用Matlab軟件求解，并給出圖形。

法國(guó)科學(xué)家勒讓德于1806年獨(dú)立發(fā)現(xiàn)了最小二乘法(least square，LS)［2］。從此，這個(gè)方法被廣泛應(yīng)用于各行各業(yè)的工程計(jì)算中。在Matlab中可以編寫最小二乘法的應(yīng)用程序進(jìn)行系統(tǒng)參數(shù)辨識(shí)［3］。

偏微分方程工具箱提供了研究和求解空間二維偏微分方程問題的一個(gè)強(qiáng)大而又靈活實(shí)用的環(huán)境［4］，用它可方便地求解熱方程，得出金屬體的內(nèi)部特性。

1 熱傳導(dǎo)方程的解法

1.1 熱傳導(dǎo)方程介紹

熱傳導(dǎo)方程由熱力學(xué)第一定律和傅里葉定律推導(dǎo)得到。為簡(jiǎn)單起見，現(xiàn)只討論單純熱傳導(dǎo)現(xiàn)象，并忽略體積膨脹，同時(shí)假設(shè)金屬板靜止，可得到如下方程:

式中:T為某時(shí)某點(diǎn)溫度;t為時(shí)間;λ為導(dǎo)熱率，通常為常數(shù);ρ為密度;cv為等體比熱容;▽為拉普拉斯算子。此公式表示某點(diǎn)溫度隨時(shí)間的變化與這一點(diǎn)溫度的散度成正比。

由于拋物型偏微分方程標(biāo)準(zhǔn)形式為:

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

方程邊界條件一般有Dirichlet條件和Neumann條件兩種，即:

式中:h為空間采樣步長(zhǎng);u為方程的未知函數(shù);n為邊界切面的法向量;h、r、c、q、g 都為已知函數(shù)或常數(shù)。

1.2 解金屬板導(dǎo)熱問題

考慮一個(gè)帶有矩形孔的金屬板上的熱傳導(dǎo)問題。板左邊的溫度保持在100℃，板右邊熱量從板向環(huán)境空氣定常流動(dòng)，其他邊及內(nèi)孔邊界保持絕緣。

當(dāng)初始板溫t=0℃時(shí)，則可概括為如下定解問題:

上述問題中金屬板域的外邊界頂點(diǎn)坐標(biāo)分別為(－0.5，－0.8)、(0.5，－ 0.8)、(0.5，0.8)、(－ 0.5，0.8);內(nèi)邊界 頂點(diǎn)坐 標(biāo)分別為(－ 0.05，－ 0.4)、(0.05，－0.4)、(0.05，0.4)、(－ 0.05，0.4)。接著使用Matlab圖形用戶界面(GUI)求解這一問題。在PDE Toolbox窗口的工具欄中選擇Generic Scalar模式，然后經(jīng)過區(qū)域設(shè)置、邊界條件設(shè)置、方程類型設(shè)置、網(wǎng)格剖分、初值和誤差設(shè)置，可以得到此問題的數(shù)值解和解的圖形。

部分Matlab程序如下。

1.3 熱傳導(dǎo)方程的差分格式

熱方程為偏微分方程，求解它可用求解微分方程的數(shù)值方法。

考慮上述二維常系數(shù)熱傳導(dǎo)方程為:

要求解這個(gè)方程，常用差分格式的遞推公式，其分為古典顯格式和古典隱格式兩種。

1.3.1 古典顯格式

取網(wǎng)格的空間步長(zhǎng)為h=1/N，時(shí)間步長(zhǎng)為τ。在(xj，tk)處時(shí)間用向前差商，空間用中心差商近似，則可得到古典顯格式，即:

式中:局部截?cái)嗾`差為O(τ+h2)，r=τ/h2(網(wǎng)格比)，其中，O表示高階無窮小。當(dāng)且僅當(dāng)r＜1/4時(shí)，格式是穩(wěn)定的。

1.3.2 古典隱格式

在(xj，tk+1)處時(shí)間用向后差商，空間用中心差商近似，則可得古典隱格式，即:

式中:局部截?cái)嗾`差為 O(τ+h2)，r=τ/h2，隱格式對(duì)任何網(wǎng)格比都是穩(wěn)定的。

用差分方法同樣可以在Matlab上求解熱傳導(dǎo)方程，在鍵入程序后，就可得到板上各點(diǎn)的溫度值，對(duì)于更復(fù)雜的問題，此方法有明顯的優(yōu)勢(shì)。

2 最小二乘原理與算法

2.1 最小二乘基本原理

最小二乘法(LS)是用于參數(shù)估計(jì)的數(shù)學(xué)方法，它使數(shù)學(xué)模型在誤差平方和最小的意義上擬合試驗(yàn)數(shù)據(jù)，也是一種涉及較少數(shù)學(xué)基礎(chǔ)而又被大量應(yīng)用的一種基本方法。

最小二乘法提供一個(gè)估算方法，用來得到一個(gè)在最小方差意義上與試驗(yàn)數(shù)據(jù)最好擬合的數(shù)學(xué)模型。由最小二乘法獲得的估計(jì)在一定的條件下有最佳的統(tǒng)計(jì)特性，即估計(jì)的結(jié)果是無偏的、一致的和有效的。

因?yàn)長(zhǎng)S的原則是希望某個(gè)量S(0)和a的值能使觀測(cè)值和由模型的計(jì)算值之間的誤差為最小，所以各次觀測(cè)誤差可表示為:

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

整個(gè)觀測(cè)過程的誤差是由各次觀測(cè)誤差所組成的，采用每個(gè)誤差的平方和作為總誤差:

所選的誤差平方和函數(shù)J就是估計(jì)參數(shù)時(shí)所采用的性能指標(biāo)。顯然，J越小越好，即所選的S(0)和a的值能使每個(gè)誤差的平方和J的值最小。由于平方運(yùn)算也稱二乘運(yùn)算，因此稱上述方法為最小二乘估計(jì)法。要使J達(dá)到極小值，只需分別對(duì)S(0)和a求偏導(dǎo)，令它們等于0。這樣就可以得到關(guān)于S(0)和a的兩個(gè)估計(jì)表達(dá)式。

2.2 最小二乘的基本關(guān)系式

對(duì)于單輸入單輸出離散時(shí)間動(dòng)態(tài)系統(tǒng)，設(shè)u(k)為輸入、y(k)為輸出、z(k)為量測(cè)輸出、v(k)為白噪聲，其系統(tǒng)數(shù)學(xué)關(guān)系可用如下隨機(jī)差分方程描述，即:

利用數(shù)據(jù)序列{u(k)}、{z(k)}，極小化下列準(zhǔn)則函數(shù)。

在文獻(xiàn)［3］中，根據(jù)所求解的問題的J不同，在不同場(chǎng)合下，函數(shù)J往往有不同的名稱，它是一個(gè)標(biāo)量。對(duì)J求導(dǎo)，令其等于0可得:

式(14)就是所要求的參數(shù)估計(jì)式，從統(tǒng)計(jì)學(xué)的角度考慮。觀測(cè)總次數(shù)l必須大大超過設(shè)定的未知參數(shù)的數(shù)目，此時(shí)由觀測(cè)所提供的方程式的數(shù)目才能超過確定出方程組的唯一解所需的數(shù)目。

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

3.1 用積矩矩陣估計(jì)階次

考慮干擾是白噪聲的情況，令:

當(dāng)n從1開始逐一增加時(shí)，若上式值有明顯的增加，則在顯著增加處的n值就是模型的階次。由仿真得出本文所討論的問題的階次大約在2的附近，所以本問題的差分方程模型的階次選為2。

3.2 模型的參數(shù)辨識(shí)

程序框圖如圖1所示。

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

根據(jù)試驗(yàn)得到的數(shù)據(jù)，將本問題的模型階次選為2，并用線性差分方程來描述，可得:

式中:v(k)為服從正態(tài)分布的白噪聲N(0，1);z(k)為金屬板中心溫度與初始溫度(高于室溫)之差;u(k)為1，表示金屬板兩端加熱源，為－1，表示不加熱源。輸入信號(hào)采用4階M序列，幅度為1。則按式(17)構(gòu)造zl和HL，數(shù)據(jù)長(zhǎng)度 L=14，加權(quán)陣取 I，并利用式(14)計(jì)算參數(shù)估計(jì)值。

根據(jù)金屬板的ARMX模塊估計(jì)結(jié)構(gòu)，由程序運(yùn)行可得如圖2所示的溫度輸出曲線和輸入輸出誤差曲線圖。

圖2 ARMX模塊估計(jì)輸出曲線和誤差曲線Fig.2 Output curve and error curve of ARMX estimated module

由圖2可以看出，溫床系統(tǒng)的估計(jì)輸出曲線與量測(cè)溫度值相符，且誤差曲線也符合要求。

圖3給出了一定點(diǎn)溫度在階躍功率輸入下的輸出變化，隨時(shí)間輸出溫度趨于穩(wěn)定，反映了金屬板的熱傳導(dǎo)性質(zhì)。從圖4所示的頻率響應(yīng)曲線可以看出，在低頻信號(hào)下輸出有較好的穩(wěn)定性，當(dāng)外界有高頻噪聲時(shí)，采用ARMX模塊估計(jì)就會(huì)出現(xiàn)較大的偏差。

4 結(jié)束語

從仿真結(jié)果看，本文所采用的方法能較好地描述金屬板的內(nèi)、外部特性，得到了所需系統(tǒng)模型，較好地解決了上面提出的問題，可以滿足實(shí)際的要求。在更復(fù)雜的環(huán)境中，還應(yīng)對(duì)噪聲和信號(hào)進(jìn)行分析，以便更精確地獲得系統(tǒng)模型的其他數(shù)據(jù)。

［1］王竹溪.熱力學(xué)［M］.2版.北京:北京大學(xué)出版社，2005.

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

［3］陸君安，尚濤，謝進(jìn)，等.偏微分方程的 MATLAB解法［M］.武漢:武漢大學(xué)出版社，2001.

［4］侯媛彬，汪海，王立琦.系統(tǒng)辨識(shí)及其MATLAB仿真［M］.北京:科學(xué)出版社，2004.

［5］Iserles A.微分方程數(shù)值分析基礎(chǔ)教程［M］.劉曉艷，劉學(xué)深，譯.北京:清華大學(xué)出版社，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］周彤.含有不確定因素的模型檢驗(yàn)及其非偽概率估計(jì)［J］.控制理論與應(yīng)用，1996(2):145－152.