本发明属于工业过程预测与控制领域，涉及一种基于分布估计局部优化的最小二乘支持向量机软测量建模方法。
背景技术：
软测量建模技术即为在实际的工业生产过程中，一些过程变量和质量变量难以使用传感器进行直接测量或者测量成本过高，因此人们往往会利用较为容易测量的一些过程变量，通过建立数学模型的方法来估计那些难以测量的过程或质量变量，进而能够很好地控制产品质量，提高生产效率。
在常用的软测量数学模型中，最小二乘支持向量机模型(以下简称为lssvm)以其使用训练样本小，求解简单，运算速度快等优点被广泛地用于各类工业过程变量的软测量。但该模型本身仍具有较大的缺陷，其模型中的核函数参数和惩罚因子往往需要人为进行设定，而这两个参数会直接影响模型的训练效果，若选择不当则会造成较大的预测偏差。
为了解决这个问题，目前较为常用的方法为基于参数优化的lssvm，即基于训练样本，利用优化的方式对lssvm中的核函数参数及惩罚因子进行寻优的方法，常用的寻优方法主要有遗传算法、粒子群算法、差分进化算法等。该方法能够在lssvm进行预测之前，以最小化训练集的偏差为目标，搜索出一组较佳的模型参数，能够有效地提高lssvm模型的预测精度，该种方法由于在优化过程中使用全部训练集的预测偏差作为目标函数因此可称其为全局优化。
但该方法仍然具有一定的缺陷，且主要分为两个部分，一方面如果只以最小化训练样本集的偏差为优化目标会使算法的泛化能力变差，这种情况主要发生在训练集和测试集的数据分布差异较大时，即利用训练集优化出的模型参数会对训练集有较好的拟合效果，而测试集可能差异较大；另一方面，该种方法对寻优算法的要求也较高，即优化算法应该具有更好的收敛性和寻优效果，而目前很多常用的智能优化算法如(遗传算法)仍不具备这些特性。
技术实现要素：
针对现有基于参数优化lssvm的不足，本发明提出了一种基于分布估计局部优化的最小二乘支持向量机软测量建模方法，该方法通过局部模型建立优化目标，并采用了具有较快收敛速度的分布估计算法进行寻优，使算法即使在训练和测试样本集差异较大的情况下仍能以较快的速度对测试样本进行精准预测。
本发明的具体技术方案如下：
一种基于分布估计局部优化的最小二乘支持向量机软测量建模方法，其特征在于，该方法包括以下步骤：
步骤一：收集历史工业过程数据建立训练集，训练集中共有n个样本，每个样本具有m个易于测量的过程变量，以及一个需要被估计的质量变量，该训练样本集表示为s_train＝{(xi,yi)|i＝1,2…,n}，其中xi∈r1×m，yi∈r1，并用s_train(x)表示样本中的过程变量，s_train(y)表示样本中的质量变量，其余样本集合表示方法相同；模型的测试集与训练集具有相同的形式，表示为s_test＝{(xi,yi)|i＝1,2…,nt}，nt为测试集中的样本个数；
步骤二：将步骤一中的所得到的训练样本集和测试样本集中的xi进行标准化，即将过程变量的均值化为0，方差化为1，所得到标准化后的训练样本集和测试样本集为s_std_train和s_std_test；
步骤三：针对标准化后测试样本集中的每个测试样本在训练集中找到与之欧式距离最近的样本构建用于最小二乘支持向量机模型参数局部优化的目标集s_obj，并在标准化后的训练样本集中s_std_train中删除所选出的用于作为目标集的训练样本，得到用于局部优化的训练集
步骤四：以目标集s_obj中每个样本的预测值与实际值偏差绝对值的和作为用于最小二乘支持向量机模型参数局部优化的适应值函数f，同时以径向基核函数参数σ2及惩罚因子γ两个待优化参数的取值范围作为局部优化的约束条件；
步骤五：初始化单变量边缘分布算法的种群大小、优化变量个数、最大迭代次数及截断选择率等相关参数，并根据待优化参数的取值范围随机采样建立初始种群pop_init；
步骤六：将初始种群中的每一个体，即为pop_init的每一行，带入步骤四中的适应值函数，即利用局部优化训练集对模型训练，并对局部优化目标集进行预测，并与真实值代入适应值函数中计算适应值；
步骤七：将步骤五得到的初始种群，以及步骤六每个个体对应的适应值带入单变量边缘分布算法中进行迭代，直至满足最大迭代次数为止，返回一组优化后的最小二乘支持向量机模型参数(γbest,σ2best)；
步骤八：将步骤七所得到的优化后的模型参数带回最小二乘支持向量机模型中并利用全体训练集重新对模型进行训练，所得模型如下式所示
y(x)＝lssvms_std_train(x,γbest,σ2best)
步骤九：将测试集带入步骤八训练的模型中，具体如下式所示，即可实现对质量变量的软测量
y＝lssvms_std_train(s_std_test,γbest,σ2best)。
进一步地，所述步骤三中选择用于局部优化的训练集和目标集的具体步骤如下：
(a)针对标准化后的测试集s_std_test中每一个样本在标准化后的训练集s_std_train中找到与之欧式距离最近的一个样本，并记录所找到的样本在训练集中的位置；
(b)将所找到的nt个样本组建成用于局部优化的目标集s_obj，将记录目标集在原训练集位置的索引组建成目标集的索引集s_obj_index；
(c)根据目标集的索引集s_obj_index，在s_std_train中删除所选出用于构建目标集的样本，得到用于局部优化的训练集为
进一步地，所述步骤四具体为：
(a)使用
(b)计算局部优化目标集各个样本的预测值与实际值偏差绝对值的和作为适应值函数f，其表达式如下式所示：
(c)以模型参数γ和σ2的最大与最小取值范围来建立约束条件，如下式所示；
进一步地，所述步骤五具体为：
(a)设定单变量边缘分布算法参数，包括优化的变量个数pop_id，种群大小pop_size，最大迭代次数max_iteration以及截断选择率selectrate；
(b)在步骤四中所设定的模型参数约束条件范围内进行随机采样，构建大小为pop_size×pop_id的初始种群pop_init。
进一步地，所述的步骤三中的样本间欧式距离的计算公式为
dij＝||s_std_test(xi)-s_std_train(xj)||2
式中i＝1…nt，j＝1…n，dij为测试集中第i个样本与训练集中第j个样本之间的欧式距离。
本发明的有益效果如下：
本发明利用测试集建立用于lssvm模型参数优化的局部目标集，使得优化的过程更贴合与测试集，有效地防止了传统使用全体训练集作为优化目标而产生的过拟合现象，对分布变化程度较高的数据具有很强的适应能力，并且使用分布估计算法进行参数优化能够有效提高收敛精度，因此基于本发明所构建的软测量模型可以较为精准地实现关键质量变量的预测和控制。
附图说明
图1为基于分布估计局部优化的最小二乘支持向量机软测量模型的预测曲线图；
图2为基于分布估计全局优化的最小二乘支持向量机软测量模型的预测曲线图；
图3为分别基于分布估计和标准粒子群算法的局部优化的最小二乘支持向量机软测量模型的收敛曲线。
具体实施方式
下面根据附图和优选实施例详细描述本发明，本发明的目的和效果将变得更加明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
本发明的基于分布估计局部优化的最小二乘支持向量机软测量建模方法中，所使用的最小二乘支持向量机算法建模过程如下，最小二乘支持向量机回归模型如下式所示
其中k(x,xi)为核函数，αi和b为模型的参数，l为训练样本个数，xi为训练样本，y(x)为模型输出即为预测值。
模型中的参数αi和b是通过下式所示的线性方程组进行求解的
其中y＝[y1,y2,...,yl]t，z＝[1,...,1]t，α＝[α1,...,αl]t，d＝diag(γ-1,...,γ-1)，
所述最小二乘支持向量机模型使用的核函数为径向基核函数，如下式所示
因此在lssvm模型的训练过程中，参数αi和b可以通过训练集求得，而惩罚因子γ和径向基核函数参数σ2却需要人为设定，故也为变量，因此lssvm模型由训练样本集s，惩罚因子γ和径向基核函数参数σ2确定，因此lssvm模型可写为
y(x)＝lssvms(x,γ,σ2)
本发明所使用的分布估计算法(eda)为单变量边缘分布算法(下述称umda)，其优化的过程如下7个步骤所示。
(a)在n个变量的阈值范围内随机产生n×m维初始种群pop(0)，其中m为种群大小；
(b)将初始种群带入适应值函数中计算适应值并按照适应值由低到高排序；
(c)利用ms＝floor(selectrate×m)计算截断选择的精英种群个体数ms，并取原种群适应值较低的ms个体组成精英种群popmaster(0)，selectrate为截断选择率；
(d)计算精英种群的均值和标准差μi和σi(i＝1…n)；
(e)由μi和σi(i＝1…n)构建n个一维累积高斯概率密度函数，并随机产生(m-ms)维概率矩阵，带入n个一维累积高斯分布函数中进行采样，获得n×(m-ms)维下一代部分种群poppart(t+1)；
(f)再将本代精英种群popmaster(t)和采样出的下一代种群poppart(t+1)进行组合得到下一代种群pop(t+1)；
(g)检查是否到达迭代终止条件，若到达则输出当代最优个体和最优适应值，若未达到则将pop(t+1)返回步骤2进行迭代运算。
所述一维累积高斯概率密度函数为
其中i＝1…n，xi为变量i的随机变量，xi为变量i的取值，f(xi＝xi)为随机变量取值为xi的概率，μi和σi分别表示变量i在种群中的均值和标准差。
因此，本发明的基于分布估计局部优化的最小二乘支持向量机软测量建模方法具体步骤如下：
步骤一：收集历史工业过程数据建立训练集，训练集中共有n个样本，每个样本具有m个易于测量的过程变量，以及一个需要被估计的质量变量，该训练样本集表示为s_train＝{(xi,yi)|i＝1,2…,n}，其中xi∈r1×m，yi∈r1，并用s_train(x)表示样本中的过程变量，s_train(y)表示样本中的质量变量，其余样本集合表示方法相同；模型的测试集与训练集具有相同的形式，表示为s_test＝{(xi,yi)|i＝1,2…,nt}，nt为测试集中的样本个数；
步骤二：将步骤一中的所得到的训练样本集和测试样本集中的xi进行标准化，即将过程变量的均值化为0，方差化为1，所得到标准化后的训练样本集和测试样本集为s_std_train和s_std_test；
步骤三：构建局部优化的训练集和目标集
(a)针对标准化后的测试集s_std_test中每一个样本在标准化后的训练集s_std_train中找到与之欧式距离最近的一个样本，并记录所找到的样本在训练集中的位置，所述样本间欧式距离的计算公式为
dij＝||s_std_test(xi)-s_std_train(xj)||2
式中i＝1…nt，j＝1…n，dij为测试集中第i个样本与训练集中第j个样本之间的欧式距离；
(b)将所找到的nt个样本组建成用于局部优化的目标集s_obj，将记录目标集在原训练集位置的索引组建成目标集的索引集s_obj_index；显然s_obj是s_std_train的子集；
(c)根据目标集的索引集s_obj_index，在s_std_train中删除所选出用于构建目标集的样本，得到用于局部优化的训练集为
步骤四：建立优化的适应值函数及约束条件
(a)使用
(b)计算局部优化目标集各个样本的预测值与实际值偏差绝对值的和作为适应值函数f，其表达式如下式所示：
从适应值函数中可知，f显含模型参数γ和σ2；
(c)以模型参数γ和σ2的最大与最小取值范围来建立约束条件，如下式所示；
步骤五：初始化单变量边缘分布算法的参数及建立初始种群
(a)设定单变量边缘分布算法参数，包括优化的变量个数pop_id，种群大小pop_size，最大迭代次数max_iteration以及截断选择率selectrate；
(b)在步骤四中所设定的模型参数约束条件范围内进行随机采样，构建大小为pop_size×pop_id的初始种群pop_init，由于待优化的变量为最小二乘支持向量机模型的两个参数，故pop_id＝2；
步骤六：将初始种群中的每一个体，即为pop_init的每一行，带入步骤四中的适应值函数f，即利用局部优化训练集对模型训练，并对局部优化目标集进行预测，并与真实值代入适应值函数中计算适应值；
步骤七：将步骤五得到的初始种群，以及步骤六每个个体对应的适应值带入单变量边缘分布算法中进行迭代，直至满足最大迭代次数为止，返回一组优化后的最小二乘支持向量机模型参数(γbest,σ2best)；
步骤八：将步骤七所得到的优化后的模型参数带回最小二乘支持向量机模型中，并使用全体训练集对最小二乘支持向量机进行训练，所得的模型如下式所示
y(x)＝lssvms_std_train(x,γbest,σ2best)
步骤九：将测试集带入步骤八训练的模型中，具体如下式所示，即可实现对质量变量的软测量
y＝lssvms_std_train(s_std_test,γbest,σ2best)。
以下结合一个具体的工业过程实例来验证本发明的有效性。co2吸收塔是工业合成氨的一个子单元，来自前一单元的工艺气在经过初步降温后送入工艺冷凝分离罐进行再次降温，并送入吸收塔；经过吸收塔后的工艺气，再送入除雾分离罐中，残余的co2含量使用仪器进行记录；吸收塔吸收co2后，吸收液由贫液，半贫液变成富液，富液从吸收塔底部经富液闪蒸槽送入再生塔，进行溶液的再生操作；之后再生的溶液又被抽回吸收塔内。
co2吸收所发生的化学反应方程式为
co2+k2co3+h2o←→2khco3+q
从上述工艺可知，为了能够使工艺气中的残余co2含量尽可能地减少，因此需要对其含量进行实时检测，而残余的co2含量即可使用本发明进行软测量。
因此选择co2吸收塔中易于测得的11个过程变量来预测工艺气中残余的co2含量，变量的编号及说明如表1所示。
表1吸收塔工艺中变量说明
从上述数据中进行采样，取1500个样本数据构建训练集，500个样本数据构建测试集，且测试样本的采样点距离训练样本采样点较远，即训练样本与测试样本在数据分布上会存在一定差异；
设置umda算法参数为种群大小pop_size＝50，最大迭代次数max_iteration＝10以及截断选择率selectrate＝0.3；
设置待优化的lssvm模型参数范围，惩罚因子γmin＝0.001，γmax＝1000，γmin≤γ≤γmax，径向基核函数参数σ2min＝0.001，σ2max＝1000，σ2min≤σ2≤σ2max
利用本发明对500个测试样本的残余co2含量进行预测，同时也使用基于umda算法全局优化(优化的训练集和目标集均为原始的训练集)的lssvm算法与之进行对比，二者优化出的模型参数及rmse，max_error值如表2所示，二者测试样本的预测曲线分别如图1、图2所示。
其中误差标准均方根误差(rmse)，计算公式如下：
其中nt为测试样本个数，yj为残余co2含量的真实值，ypj为模型预测值。
其中预测的最大误差(max_error)，计算公式如下
表2基于umda局部优化和全局优化的lssvm预测效果对比
从表2及图1-2即可看出，基于umda局部优化的lssvm模型比全局优化的lssvm模型具有更高的预测精度，且基于全局优化的lssvm在训练集上的表现几乎没有偏差，这更验证了该种方法易发生对训练集的过拟合，而基于局部优化的lssvm则可以解决这个问题。
接下来验证umda算法具有更好地收敛效果，利用标准粒子群算法(pso)局部优化lssvm模型参数，与之进行对比，设定umda和pso的种群大小均为50，最大迭代次数均为10，lssvm参数取值范围不变，并设定pso算法学习因子分别为1和2，所得两种方法优化出的参数及偏差如表3所示，两种优化方法的收敛曲线如图3所示。
表3基于umda局部优化和基于pso局部优化的lssvm预测效果对比
从表3中基于umda和pso局部优化的lssvm模型的预测误差可知，umda算法能够更好地控制lssvm的预测偏差，从图3的收敛曲线也可以看出umda的最终适应值要小于pso；从多次实验结果来看，umda的最终适应值总是小于pso，验证了基于umda局部优化lssvm模型参数会有更高的精准度和更好的收敛性。
本领域普通技术人员可以理解，以上所述仅为发明的优选实例而已，并不用于限制发明，尽管参照前述实例对发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实例记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在发明的精神和原则之内，所做的修改、等同替换等均应包含在发明的保护范围之内。