前言
湿空气物性是我们制冷低温专业非常常用的工具，目前大家要计算相应工况的参数，通常是查表或是套用他人拟合的公式。目前大多数的拟合公式是基于最小二乘法来拟合实验数据得到的。
不过刚好这学期在《最优化方法》课上学习了共轭梯度下降法，刚好可以结合之前在机器学习课上学到的回归预测方法来自己拟合个公式试试看。同时自己写一个共轭梯度下降的优化器，替代先前普通梯度下降。
整体程序是基于MATLAB实现的。不过最终的拟合结果是挺让人惊喜的，意外的精度非常高，达到了0.1%的级别。
湿空气物性曲线
训练集选取与前处理
经典的有监督学习流程如下
因此第一部是训练集的选取，具体参数来源我们选择的是制冷界认可度最高的Ashrae Handbook，这是Ashrae团队在1988年做的一系列关于空气物性的实验，并将实验数据记录在其中的手册。制冷界在对物性拟合时大多是参考其中的数据。
Ashrae handbook
书中的原始数据中温度t采用的是摄氏温度，压力单位是pa。对于一个合理的热力学公式模型，温度用开氏温度（K）表示是最为合适的，因此我们公式中的温度T都是基于原摄氏温度t上，增加了273.15，以转化为开式温度表示，即
T=t+273.15 （单位：K）
其中t是摄氏温度℃，T是与之对应的开氏温度°K
对于压力，用帕斯卡(pa)来表示的话，容易造成低温下，压力仅有个位数的数值；而高温时，压力达到成千上万。这种情况对于数值分析来说是应该极力避免的，因为数据的跨度过大的话，一点微小的计算机误差都会对结果造成巨大的影响。因此我们需要对压力做特征标定（Feature scaling）,具体做法是对原压力取对数，因而有
Ps=lg(p)
其中Ps是新标度下的压力，p是手册上记载的原始压力
最后将数据集按4:1的比例分割为训练集与交叉验证集。至此，数据来源与数据预处理部分完成。
模型选取
选择函数模型是监督学习中相当重要的一步，若没能选择合适的模型，不管有多少训练样本，不管迭代多少步都不会得到理想的预测结果。模型过于简单会欠拟合；而若用过于复杂的模型去训练，则容易发生过拟合。
因此为了较好地拟合数据，需要选择一个适合的模型，我们根据湿空气的物性曲线，同时结合前人在湿空气物性公式中，所发现的拟合度比较好的模型，综合考虑选择的模型如下：
（1）四次多项式模型，经典的高复杂度模型[3]
经典四次多项式模型
(2)参考Riedel-Plank-Miller对临界压力拟合方程[4]：
对上述方程改写为由温度预测饱和压力的形式，以下简称模型二
这两个公式都有研究者曾做过相应的物性拟合工作，在2013版的Ashrae handbook中以模型一，通过数值分析中的插值算法实现的拟合。本文也将采用模型一对湿空气物性进行拟合，然后与Ashrae handbook中的拟合结果作对比，观察我们的方式拟合的结果是否具有更高的精度。
将模型一改写为通用机器学习的符号形式，即把ln Ps改为换为hypothesis(假设)来代表我们的计算值，将t改写为x，模型参数b改写为θ，从而得到如下形式
损失函数
对于代价函数，没有使用目前在机器学习中较为主流的交叉熵形式，采用的是均方误差，主要出于两方面的考量：
1、我们的需要拟合的曲线本身就是很好的二次凸函数，交叉熵往往是为了弥补原函数非凸而采取的方式，因此也没有必要使用。
2、方便我们在后面使用共轭梯度下降法，共轭梯度下降法不适用于交叉熵函数。
其中，m：训练样本的个数
hθ(x)：通过模型预测的结果值
Y：原训练样本中的标签值
均方误差的表现形式较为直观，即模型预测值与真实值的差的平方的累加。
优化算法——FR共轭梯度下降法
由于机器学习的目标是调整模型中的参数，使得代价函数最小。因此可以将该机器学习问题转化为无约束最优化问题的表述方式，FR共轭梯度下降法能很好地处理这类问题。
FR共轭梯度法，其基本思想是把共轭性于最速下降法相结合，利用已知点处的梯度构造一组共轭方向，并沿这组方向进行搜索，求出函数的极小点，对于二次凸函数，沿着共轭方向经过有限次迭代必定到达极小点。
具体算法原理就不赘述了，最优化课本上写了十几页呢··介绍下大概的步骤
与adams优化法，FR法在解决代价函数为2次凸函数时具有速度与精度的双重准确性，是更为合适的。
梯度计算与海森矩阵
对于我们的代价函数
其中h表示预测值
对于J(θ)依次θ0到θ4求导，得到梯度为
继续求二阶导，得到海森矩阵（5X5）
在确定梯度与海森阵的形式后，只需要带入先前预处理好的数据集即可。
至此我们已经完成了对代价函数的梯度与海森阵的计算，理论推导部分全部完成，可以开始编写程序开始拟合工作了
程序实现
主要使用的是Matlab 2018a进行程序的编写，程序部分主要分为代价函数与共轭梯度法这两部分的编写，并着重介绍编程过程中采用的数据结构与算法思想，并展示部分关键代码。
1、代价函数
代价函数部分同时包括了代价函数值、梯度与海森阵的计算，使用了matlab的函数功能，输入为theta，同时输出代价函数值、梯度值与海森阵。
function [jval,gradient,hess] = costfunction(theta)
函数需要先提取训练数据，因此第一步是读取储存在excel文件中的数据集，由于数据集容量较大，直接用循环结构会造成运行效率低下。因此本程序中将主要使用矩阵的形式进行计算，矩阵运算会采用并行计算的模式，相较于循环的单线结构，效率会大大地提升。
t=xlsread('t2.xlsx');
p=xlsread('ps2');
t2=t.*t;
t3=t2.*t;
t4=t3.*t;
t_1=t.^-1;%数据预处理
t_log=log(t);
m=91;
r=ones(91,1);
之后直接按照前文推导公式，输入代价函数值、梯度值与海森阵的计算公式即可，计算过程中调用了很多matlab自带的矩阵计算函数。
代价函数代码：
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);
至此代价函数的程序编写完毕。
2、共轭梯度迭代法
为了提升程序效率以及方便查看每个参数整体迭代结果，同样采用了矩阵与向量的数据结构
beta=zeros(11,1);%为循环内的变量事先开辟内存，提高运行速度
theta=zeros(5,50);
cost=zeros(1,50);
之后按照共轭梯度法的步骤进行，FR法第一部需要先按照最速下降法进行
theta(:,1)=zeros(5,1);%为θ选择初始值，我们把初始θ全部设为0
[cost(1),g(:,1),hess]=costfunction(initialtheta); %计算代价函数、梯度、海森阵
d(:,1)=-g(:,1);
lambda(1)=-(g(:,1)'*d(:,1))/(d(:,1)'*hess*d(:,1));
theta(:,2)=initialtheta+lambda(1)*d(:,1);
在第一部迭代后，将采用共轭梯度法进行之后的迭代步骤。因为需要保证最终训练结果在极小值点，本程序一共进行了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
至此共轭梯度迭代法的编程全部完成。
除此之外还有测试部分以及误差计算部分的程序代码。将各个部分的程序都完成后，以函数调用的方式在主程序中将各部分串联起来，从而完成我们的训练过程。
拟合结果
按上文描述完成编程工作后，开始运行程序，经过大约一小时的训练与迭代。以下是最终的训练结果。
1、代价函数
衡量一次训练结果是否成功，一个重要的指标是代价函数是否减小了；一般而言，代价函数下降的越显著，说明模型的预测值与真实值越接近，从而训练的结果就越理想。
可以看到，代价函数值从初始的10^4数量级，经过数次迭代最终下降到203.8887，同时函数值没有出现回升的情况，因此训练确实使得模型向正确的反向前进了。
2、梯度
对于FR共轭梯度法，判定是否最终达到极小点的依据是，梯度是否为0。
本文中由于模型有五个参数，因此对应梯度为5X1的向量，对每次迭代的梯度取无穷范数，绘制梯度趋势图如下
可以看到，梯度值最开始非常大，经过几次迭代后迅速收敛，这就是共轭梯度法二次收敛性的好处，在第29步后，梯度变为5.03x10-7，可以认为此时梯度为0。并且之后保持这个数据直到第50步为止都不变化，因此可以认为最终达到了极小值点。
3、拟合公式与误差分析
对于我们选取的模型
拟合得到参数，为
从0℃到90℃，每隔5℃选取一湿空气样本，这是之前设定的交叉验证集，没有被拿来训练，来验证我们拟合得到的最终曲线精度情况，对比情况如下表
通过误差情况可以看出，曲线在低温段误差较大，在t=5℃时，误差最大达到了70%，但是当温度提高后，误差显著缩小。在t=20℃到40℃段，误差开始明显的减小，小于10%，并一直在下降；在t=40℃到90℃段，误差最大不超过1%，对于制冷研究计算而言，1%的误差是足够精准的。
与标准Ashrea拟合公式对比
Ashrae相当于制冷领域的圣经一般的书本了，几乎大家查公式，或要查某些实验参数都会以它为基准。
不过Ashrae上拟合的公式都是基于最小二乘法的，它的公式也是目前在大多数空调系统数字化设计中所采用的公式。我们用同样的公式模型看看机器学习的方法拟合到的公式有什么差别
对比发现Ashrae handbook公式对于各个温度段的误差都在1%以上，可以看到本文拟合得到的公式，在t=40℃到90℃段误差均小于1%，因此精度是明显优于Ashrae handbook公式。、
总结
对比传统的利用数值分析中的插值方法对参数进行曲线拟合的方法，本文通过利用机器学习中的监督学习法与共轭梯度下降法得到的湿空气物性公式是具有优越性与实际意义的。
对于低温度下的误差过大现象产生的原因，可能是因为低温部分样本对应的压力值也比较低，同时在总样本数中占据比例小，因此在代价函数中的权重不高，从而使得机器在学习过程中更倾向于去“接近”高温度情况下的数据，而忽略的低温度的数据。
程序源代码请移步我的github:
https://github.com/huangchuhccc/Conjugate-Gradient-Method-Moist-air-fitting
参考文献
[1]ASHRAE(2013). ASHRAE Handbook of fundamentals Atlanta,GA. American Sociert of Heating Refrigeration and Air---conditioning Engineers, inc.
[2]韩学孟.一种实用的水饱和蒸气压拟合方程[J]. 山西农业大学学报，1996,16(3):278 -280.
[3]王婷，谷波，邱峰.湿空气饱和水蒸气曲线计算模型的建立与分析.制冷技术，2010,38（5）：47-49
网友评论