前言

??濕空氣物性是我們制冷低溫專業(yè)非常常用的工具，目前大家要計(jì)算相應(yīng)工況的參數(shù)，通常是查表或是套用他人擬合的公式。目前大多數(shù)的擬合公式是基于最小二乘法來擬合實(shí)驗(yàn)數(shù)據(jù)得到的。

??不過剛好這學(xué)期在《最優(yōu)化方法》課上學(xué)習(xí)了共軛梯度下降法，剛好可以結(jié)合之前在機(jī)器學(xué)習(xí)課上學(xué)到的回歸預(yù)測方法來自己擬合個(gè)公式試試看。同時(shí)自己寫一個(gè)共軛梯度下降的優(yōu)化器，替代先前普通梯度下降。

??整體程序是基于MATLAB實(shí)現(xiàn)的。不過最終的擬合結(jié)果是挺讓人驚喜的，意外的精度非常高，達(dá)到了0.1%的級別。

濕空氣物性曲線

訓(xùn)練集選取與前處理

??經(jīng)典的有監(jiān)督學(xué)習(xí)流程如下

??因此第一部是訓(xùn)練集的選取，具體參數(shù)來源我們選擇的是制冷界認(rèn)可度最高的Ashrae Handbook，這是Ashrae團(tuán)隊(duì)在1988年做的一系列關(guān)于空氣物性的實(shí)驗(yàn)，并將實(shí)驗(yàn)數(shù)據(jù)記錄在其中的手冊。制冷界在對物性擬合時(shí)大多是參考其中的數(shù)據(jù)。

Ashrae handbook

??書中的原始數(shù)據(jù)中溫度t采用的是攝氏溫度，壓力單位是pa。對于一個(gè)合理的熱力學(xué)公式模型，溫度用開氏溫度（K）表示是最為合適的，因此我們公式中的溫度T都是基于原攝氏溫度t上，增加了273.15，以轉(zhuǎn)化為開式溫度表示，即
??????T=t+273.15 （單位：K）
??其中t是攝氏溫度℃，T是與之對應(yīng)的開氏溫度°K
??對于壓力，用帕斯卡(pa)來表示的話，容易造成低溫下，壓力僅有個(gè)位數(shù)的數(shù)值；而高溫時(shí)，壓力達(dá)到成千上萬。這種情況對于數(shù)值分析來說是應(yīng)該極力避免的，因?yàn)閿?shù)據(jù)的跨度過大的話，一點(diǎn)微小的計(jì)算機(jī)誤差都會(huì)對結(jié)果造成巨大的影響。因此我們需要對壓力做特征標(biāo)定（Feature scaling）,具體做法是對原壓力取對數(shù)，因而有
??????Ps=lg(p)
??其中Ps是新標(biāo)度下的壓力，p是手冊上記載的原始壓力

最后將數(shù)據(jù)集按4:1的比例分割為訓(xùn)練集與交叉驗(yàn)證集。至此，數(shù)據(jù)來源與數(shù)據(jù)預(yù)處理部分完成。

模型選取

??選擇函數(shù)模型是監(jiān)督學(xué)習(xí)中相當(dāng)重要的一步，若沒能選擇合適的模型，不管有多少訓(xùn)練樣本，不管迭代多少步都不會(huì)得到理想的預(yù)測結(jié)果。模型過于簡單會(huì)欠擬合；而若用過于復(fù)雜的模型去訓(xùn)練，則容易發(fā)生過擬合。

因此為了較好地?cái)M合數(shù)據(jù)，需要選擇一個(gè)適合的模型，我們根據(jù)濕空氣的物性曲線，同時(shí)結(jié)合前人在濕空氣物性公式中，所發(fā)現(xiàn)的擬合度比較好的模型，綜合考慮選擇的模型如下：
（1）四次多項(xiàng)式模型，經(jīng)典的高復(fù)雜度模型[3]

經(jīng)典四次多項(xiàng)式模型

(2)參考Riedel-Plank-Miller對臨界壓力擬合方程[4]：

??這兩個(gè)公式都有研究者曾做過相應(yīng)的物性擬合工作，在2013版的Ashrae handbook中以模型一，通過數(shù)值分析中的插值算法實(shí)現(xiàn)的擬合。本文也將采用模型一對濕空氣物性進(jìn)行擬合，然后與Ashrae handbook中的擬合結(jié)果作對比，觀察我們的方式擬合的結(jié)果是否具有更高的精度。

??將模型一改寫為通用機(jī)器學(xué)習(xí)的符號形式，即把ln Ps改為換為hypothesis(假設(shè))來代表我們的計(jì)算值，將t改寫為x，模型參數(shù)b改寫為θ，從而得到如下形式

損失函數(shù)

??對于代價(jià)函數(shù)，沒有使用目前在機(jī)器學(xué)習(xí)中較為主流的交叉熵形式，采用的是均方誤差，主要出于兩方面的考量：
1、我們的需要擬合的曲線本身就是很好的二次凸函數(shù)，交叉熵往往是為了彌補(bǔ)原函數(shù)非凸而采取的方式，因此也沒有必要使用。
2、方便我們在后面使用共軛梯度下降法，共軛梯度下降法不適用于交叉熵函數(shù)。

其中，m：訓(xùn)練樣本的個(gè)數(shù)

??hθ(x)：通過模型預(yù)測的結(jié)果值

??Y：原訓(xùn)練樣本中的標(biāo)簽值

均方誤差的表現(xiàn)形式較為直觀，即模型預(yù)測值與真實(shí)值的差的平方的累加。

優(yōu)化算法——FR共軛梯度下降法

??由于機(jī)器學(xué)習(xí)的目標(biāo)是調(diào)整模型中的參數(shù)，使得代價(jià)函數(shù)最小。因此可以將該機(jī)器學(xué)習(xí)問題轉(zhuǎn)化為無約束最優(yōu)化問題的表述方式，F(xiàn)R共軛梯度下降法能很好地處理這類問題。

??FR共軛梯度法，其基本思想是把共軛性于最速下降法相結(jié)合，利用已知點(diǎn)處的梯度構(gòu)造一組共軛方向，并沿這組方向進(jìn)行搜索，求出函數(shù)的極小點(diǎn)，對于二次凸函數(shù)，沿著共軛方向經(jīng)過有限次迭代必定到達(dá)極小點(diǎn)。

梯度計(jì)算與海森矩陣

對于我們的代價(jià)函數(shù)

其中h表示預(yù)測值

繼續(xù)求二階導(dǎo)，得到海森矩陣（5X5）

程序?qū)崿F(xiàn)

主要使用的是Matlab 2018a進(jìn)行程序的編寫，程序部分主要分為代價(jià)函數(shù)與共軛梯度法這兩部分的編寫，并著重介紹編程過程中采用的數(shù)據(jù)結(jié)構(gòu)與算法思想，并展示部分關(guān)鍵代碼。

1、代價(jià)函數(shù)

代價(jià)函數(shù)部分同時(shí)包括了代價(jià)函數(shù)值、梯度與海森陣的計(jì)算，使用了matlab的函數(shù)功能，輸入為theta，同時(shí)輸出代價(jià)函數(shù)值、梯度值與海森陣。

function [jval,gradient,hess] = costfunction(theta)

函數(shù)需要先提取訓(xùn)練數(shù)據(jù)，因此第一步是讀取儲存在excel文件中的數(shù)據(jù)集，由于數(shù)據(jù)集容量較大，直接用循環(huán)結(jié)構(gòu)會(huì)造成運(yùn)行效率低下。因此本程序中將主要使用矩陣的形式進(jìn)行計(jì)算，矩陣運(yùn)算會(huì)采用并行計(jì)算的模式，相較于循環(huán)的單線結(jié)構(gòu)，效率會(huì)大大地提升。

t=xlsread('t2.xlsx');
p=xlsread('ps2');
t2=t.*t;
t3=t2.*t;
t4=t3.*t;
t_1=t.^-1;%數(shù)據(jù)預(yù)處理
t_log=log(t);
m=91;
r=ones(91,1);

之后直接按照前文推導(dǎo)公式，輸入代價(jià)函數(shù)值、梯度值與海森陣的計(jì)算公式即可，計(jì)算過程中調(diào)用了很多matlab自帶的矩陣計(jì)算函數(shù)。
代價(jià)函數(shù)代碼：

jval=(1/2*m)*sum(([r,t,t2,t3,t4]*theta-p).^2);

梯度代碼（部分）：

gradient=zeros(5,1);
gradient(1)=(1/m)*sum([r,t,t2,t3,t4]*theta-p);
gradient(2)=(1/m)*sum(([r,t,t2,t3,t4]*theta-p).*t);
gradient(3)=(1/m)*sum(([r,t,t2,t3,t4]*theta-p).*t2);
gradient(4)=(1/m)*sum(([r,t,t2,t3,t4]*theta-p).*t3);
gradient(5)=(1/m)*sum(([r,t,t2,t3,t4]*theta-p).*t4);

海森陣代碼（部分）：

hess=zeros(5,5);
hess(1,1)=(1/m)*m;
hess(1,2)=(1/m)*sum(t);
hess(1,3)=(1/m)*sum(t2);
hess(1,4)=(1/m)*sum(t3);
hess(1,5)=(1/m)*sum(t4);

至此代價(jià)函數(shù)的程序編寫完畢。

2、共軛梯度迭代法

為了提升程序效率以及方便查看每個(gè)參數(shù)整體迭代結(jié)果，同樣采用了矩陣與向量的數(shù)據(jù)結(jié)構(gòu)

beta=zeros(11,1);%為循環(huán)內(nèi)的變量事先開辟內(nèi)存，提高運(yùn)行速度
theta=zeros(5,50);
cost=zeros(1,50);

之后按照共軛梯度法的步驟進(jìn)行，F(xiàn)R法第一部需要先按照最速下降法進(jìn)行

theta(:,1)=zeros(5,1);%為θ選擇初始值，我們把初始θ全部設(shè)為0
[cost(1),g(:,1),hess]=costfunction(initialtheta); %計(jì)算代價(jià)函數(shù)、梯度、海森陣
d(:,1)=-g(:,1);
lambda(1)=-(g(:,1)'*d(:,1))/(d(:,1)'*hess*d(:,1));
theta(:,2)=initialtheta+lambda(1)*d(:,1);   

在第一部迭代后，將采用共軛梯度法進(jìn)行之后的迭代步驟。因?yàn)樾枰ＷC最終訓(xùn)練結(jié)果在極小值點(diǎn)，本程序一共進(jìn)行了50次迭代操作

for i=2:50
    [cost(i),g(:,i),hess]=costfunction(theta(:,i));
    beta(i-1)=(norm(g(:,i)))^2/(norm(g(:,i-1)))^2;
    d(:,i)=-g(:,i)+beta(i-1)*d(:,i-1); 
    lambda(i)=-(g(:,i)'*d(:,i))/(d(:,i)'*hess*d(:,i));
    theta(:,i+1)=theta(:,i)+lambda(i)*d(:,i);
end

至此共軛梯度迭代法的編程全部完成。
除此之外還有測試部分以及誤差計(jì)算部分的程序代碼。將各個(gè)部分的程序都完成后，以函數(shù)調(diào)用的方式在主程序中將各部分串聯(lián)起來，從而完成我們的訓(xùn)練過程。

擬合結(jié)果

按上文描述完成編程工作后，開始運(yùn)行程序，經(jīng)過大約一小時(shí)的訓(xùn)練與迭代。以下是最終的訓(xùn)練結(jié)果。

1、代價(jià)函數(shù)

衡量一次訓(xùn)練結(jié)果是否成功，一個(gè)重要的指標(biāo)是代價(jià)函數(shù)是否減小了；一般而言，代價(jià)函數(shù)下降的越顯著，說明模型的預(yù)測值與真實(shí)值越接近，從而訓(xùn)練的結(jié)果就越理想。

可以看到，代價(jià)函數(shù)值從初始的10^4數(shù)量級，經(jīng)過數(shù)次迭代最終下降到203.8887，同時(shí)函數(shù)值沒有出現(xiàn)回升的情況，因此訓(xùn)練確實(shí)使得模型向正確的反向前進(jìn)了。

2、梯度

對于FR共軛梯度法，判定是否最終達(dá)到極小點(diǎn)的依據(jù)是，梯度是否為0。

本文中由于模型有五個(gè)參數(shù)，因此對應(yīng)梯度為5X1的向量，對每次迭代的梯度取無窮范數(shù)，繪制梯度趨勢圖如下

可以看到，梯度值最開始非常大，經(jīng)過幾次迭代后迅速收斂，這就是共軛梯度法二次收斂性的好處，在第29步后，梯度變?yōu)?.03x10-7，可以認(rèn)為此時(shí)梯度為0。并且之后保持這個(gè)數(shù)據(jù)直到第50步為止都不變化，因此可以認(rèn)為最終達(dá)到了極小值點(diǎn)。

3、擬合公式與誤差分析

對于我們選取的模型

擬合得到參數(shù)，為

從0℃到90℃，每隔5℃選取一濕空氣樣本，這是之前設(shè)定的交叉驗(yàn)證集，沒有被拿來訓(xùn)練，來驗(yàn)證我們擬合得到的最終曲線精度情況，對比情況如下表

通過誤差情況可以看出，曲線在低溫段誤差較大，在t=5℃時(shí)，誤差最大達(dá)到了70%，但是當(dāng)溫度提高后，誤差顯著縮小。在t=20℃到40℃段，誤差開始明顯的減小，小于10%，并一直在下降；在t=40℃到90℃段，誤差最大不超過1%，對于制冷研究計(jì)算而言，1%的誤差是足夠精準(zhǔn)的。

與標(biāo)準(zhǔn)Ashrea擬合公式對比

??Ashrae相當(dāng)于制冷領(lǐng)域的圣經(jīng)一般的書本了，幾乎大家查公式，或要查某些實(shí)驗(yàn)參數(shù)都會(huì)以它為基準(zhǔn)。
??不過Ashrae上擬合的公式都是基于最小二乘法的，它的公式也是目前在大多數(shù)空調(diào)系統(tǒng)數(shù)字化設(shè)計(jì)中所采用的公式。我們用同樣的公式模型看看機(jī)器學(xué)習(xí)的方法擬合到的公式有什么差別

??對比發(fā)現(xiàn)Ashrae handbook公式對于各個(gè)溫度段的誤差都在1%以上，可以看到本文擬合得到的公式，在t=40℃到90℃段誤差均小于1%，因此精度是明顯優(yōu)于Ashrae handbook公式。、

總結(jié)

??對比傳統(tǒng)的利用數(shù)值分析中的插值方法對參數(shù)進(jìn)行曲線擬合的方法，本文通過利用機(jī)器學(xué)習(xí)中的監(jiān)督學(xué)習(xí)法與共軛梯度下降法得到的濕空氣物性公式是具有優(yōu)越性與實(shí)際意義的。
??對于低溫度下的誤差過大現(xiàn)象產(chǎn)生的原因，可能是因?yàn)榈蜏夭糠謽颖緦?yīng)的壓力值也比較低，同時(shí)在總樣本數(shù)中占據(jù)比例小，因此在代價(jià)函數(shù)中的權(quán)重不高，從而使得機(jī)器在學(xué)習(xí)過程中更傾向于去“接近”高溫度情況下的數(shù)據(jù)，而忽略的低溫度的數(shù)據(jù)。

程序源代碼請移步我的github:
https://github.com/huangchuhccc/Conjugate-Gradient-Method-Moist-air-fitting

參考文獻(xiàn)

[1]ASHRAE(2013). ASHRAE Handbook of fundamentals Atlanta,GA. American Sociert of Heating Refrigeration and Air---conditioning Engineers, inc.
[2]韓學(xué)孟.一種實(shí)用的水飽和蒸氣壓擬合方程[J]. 山西農(nóng)業(yè)大學(xué)學(xué)報(bào)，1996,16(3):278 -280.
[3]王婷，谷波，邱峰.濕空氣飽和水蒸氣曲線計(jì)算模型的建立與分析.制冷技術(shù)，2010,38（5）：47-49