CN102855368A - 一种炼油装置高温部位腐蚀剩余寿命预测方法 - Google Patents

Abstract

本发明涉及一种炼油装置高温部位腐蚀剩余寿命预测方法；通过模拟试验获得炼油装置高温部位金属材料在不同的油品酸值、硫含量、温度和流速条件下的腐蚀速率；运用最小二乘支持向量机方法，建立油品酸值、硫含量、温度、流速与腐蚀速率之间对应关系的数学模型；测量炼油装置高温部位工艺介质的酸值、硫含量、温度和流速，根据建立的数学模型，进行腐蚀速率的预测；测量炼油装置高温部位的平均剩余壁厚，根据设备设计条件获得最小要求壁厚，利用得到的平均剩余壁厚、最小要求壁厚和预测出的腐蚀速率，计算出腐蚀剩余寿命；本方法综合考虑多种因素的协同腐蚀作用，能够快速、准确的预测出加工高硫高酸原油装置高温部位的腐蚀剩余寿命。

Description

一种炼油装置高温部位腐蚀剩余寿命预测方法

技术领域

本发明涉及一种加工高硫高酸原油炼油装置高温部位(酸值范围0～4.0wt％，硫含量0～6.0mgKOH/g，温度范围180～420℃)腐蚀剩余寿命预测方法。

背景技术

炼油工业是以石油为原料生产汽油、煤油、柴油等石油产品的能源行业，是国民经济最重要的支柱产业之一。随着世界各国对石油的需求量剧增，石油的深度开采，使得劣质、重质原油的产量逐年增加，这些劣质、重质原油呈现出高酸、高硫的变化趋势，而在炼制过程中，这些原油中的酸和硫会腐蚀炼油装置，影响装置的安全运行，影响产品的收率和质量，增加企业检维修成本，造成巨大的经济损失，甚至还有可能造成极其严重的环境污染、灾难性的人身及设备安全事故。特别是高温部位的腐蚀已经成为威胁炼油装置安全平稳运行的主要隐患之一。

现有研究发现，高硫高酸原油中的环烷酸和硫是造成炼油装置高温设备腐蚀失效的最主要原因，而影响高温环烷酸和硫腐蚀的因素又很多，包括材料种类、油品酸值、硫含量、温度条件、流速条件等，这些因素之间复杂的协同作用，导致设备发生腐蚀，目前没有定量的物理模型可以对其进行准确的描述。因此如何有效的预测受环烷酸和硫腐蚀影响的炼油设备的剩余寿命，保障炼油装置的安全运行是摆在我们面前亟待解决的一个问题。

针对炼油装置高温部位腐蚀剩余寿命的预测问题，目前采用的方法主要有定期定点测厚；现场腐蚀挂片；基于灰色预测理论或者神经网络方法。对于定期定点测厚来说，检测位置和检测周期的确定比较困难，往往依据工程经验制定检测规范，缺乏理论依据。现场腐蚀挂片周期过长，安装和拆卸极不方便。基于灰色预测理论方法所需样本数量少，计算简单，但由于在实际情况中腐蚀速率具有较大波动特性，其时间序列不够光滑，采用灰色预测模型获得的预测精度偏低。而基于神经网络的方法容易陷入局部最优解，训练效率低，网络结构与模型参数难以确保最优，且由于以经验风险最小化为训练准则，存在着过学习和低泛化的问题，难以保证预测精度，限制了其在实际过程中的应用。

发明内容

本发明的目的在于提供一种用于加工高硫高酸原油炼油装置高温部位腐蚀剩余寿命预测方法，解决现有技术中难以精确预测腐蚀剩余寿命的问题。

本发明所述的加工高硫高酸原油炼油装置高温部位腐蚀剩余寿命预测方法，步骤如下：

1)通过模拟试验获得炼油装置高温部位金属材料在不同的油品酸值、硫含量、温度和流速条件下的腐蚀速率；

2)依据实验数据，运用最小二乘支持向量机方法(LS-SVM)，建立油品酸值、硫含量、温度、流速等影响因素与腐蚀速率之间对应关系的数学模型；

3)测量炼油装置高温部位工艺介质的酸值、硫含量、温度和流速，根据建立的数学模型，进行腐蚀速率的预测；

4)测量炼油装置高温部位的平均剩余壁厚，根据设备设计条件获得最小要求壁厚，利用得到的平均剩余壁厚、最小要求壁厚和预测出的腐蚀速率，计算出腐蚀剩余寿命。

与现有技术相比，本发明具有如下优点和效果：综合考虑多种因素的协同腐蚀作用，运用基于结构最小化原则的最小二乘支持向量机方法，建立了腐蚀影响因素与腐蚀速率之间相互关系的数学模型，并编制成计算机程序，能够快速、准确的预测出加工高硫高酸原油装置高温部位的腐蚀剩余寿命，为炼油装置的安全运行提供技术支持。

附图说明

图1是本发明方法的流程图。

具体实施方式

以下结合附图和实施例对本发明的技术方案做进一步的描述。

1)通过模拟试验获得炼油装置高温部位金属材料在不同的油品酸值、硫含量、温度和流速条件下的腐蚀速率。

试验装置采用高温反应釜，采用白油为基油，以环烷酸模拟原油酸值，以十二硫醇模拟原油中硫含量。酸值范围0～4.0wt％，硫含量0～6.0mgKOH/g，温度范围180～420℃，流速范围0～30m/s。试验的金属材料包括碳钢、低合金钢、不锈钢、镍基合金、铜、铝、钛或其他金属材料。

试验的金属材料依据SY/T 5273-91《油田采出水用缓蚀剂性能平价方法》中的规定的试片规格进行加工，试片使用前进行丙酮和乙醇除油、干燥处理后称重。试验完毕后，从高温反应釜中取出试片，用丙酮洗去试片上的油渍和腐蚀产物，再用乙醇洗去试片上的丙酮，稍微用冷风干燥，置于干燥器24小时后称重，并按下式计算平均腐蚀速率V_corr：

$V_{\mathrm{corr}}(\mathrm{mm}/a)=\frac{8.76\×{10}^4}{7.8\mathrm{St}}\stackrel{\‾}{m}---\left(1\right)$

其中，

为平均失重，单位为g；S为试片面积，单位为cm²；t为试验周期，单位为h。

2)依据实验数据，运用最小二乘支持向量机方法，建立油品酸值、硫含量、温度、流速等影响因素与腐蚀速率之间对应关系的数学模型。

(1)数据预处理和样本集划分

将每一组试验数据作为一个样本，把油品酸值、硫含量、温度、流速作为样本的输入数据，对应的腐蚀速率作为输出数据。

对于给定的全部样本集，采用下面公式对样本的输入和输出数据进行归一化处理，处理之后的数值按比例缩放到[0，1]范围内。

$x_{i,j}^{\′}=\frac{x_{i,j}-\underset{i}{\min }\left(x_{i,j}\right)}{\underset{i}{\max }\left(x_{i,j}\right)-\underset{i}{\min }\left(x_{i,j}\right)}---\left(2\right)$

$y_i^{\′}=\frac{y_i-\underset{i}{\min }\left(y_i\right)}{\underset{i}{\max }\left(y_i\right)-\underset{i}{\min }\left(y_i\right)}---\left(3\right)$

式中：i＝1，Λ，m，m为全部样本的个数；j＝1，Λ，n，n为输入数据中变量的个数；x_i，j和x′_i，j分别为归一化前后的第i个输入数据中第j个变量的值；

和

分别为归一化前的输入数据中第j个变量的最大值和最小值；y_i和y′_i分别为归一化前后的第i个输出数据的值；和

分别为归一化前的输出数据的最大值和最小值。

将预处理之后的样本集划分为训练样本和测试样本。其中，训练样本用于模型训练，测试样本用来评价训练好之后的模型预测性能。

(2)建立数学模型

对于给定训练样本集：

                 T＝{(x_i，y_i)}，i＝1，L，l    (4)

式中：x_i为输入数据，y_i为输出数据，

l为训练样本的个数。

构造下面的回归函数：

                     y(x)＝w^Tφ(x)+b          (5)

式中：w为权值向量；b为偏置值；φ(·)为非线性映射函数，将输入空间映射到一个高维特征空间。

根据结构风险最小化原则，可以得到如下式的优化问题：

$\min J(w,\mathrm{\ξ})=\frac{1}{2}{w}^{T}w+\frac{1}{2}C\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\ξ}}_i^2---\left(6\right)$

等式约束条件为：

                  y_i＝w^Tφ(x_i)+b+ξ_i，i＝1，L，l    (7)

式中：C和ξ_i分别为正则化参数和误差变量，C是用于确定在训练误差和模型复杂度之间的权衡。

为了求解上述优化问题，构造Lagrange函数：

$L(w,b,\mathrm{\ξ},\mathrm{\α})=J(w,\mathrm{\ξ})-\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\α}}_i[w^T\mathrm{\φ}\left(x_i\right)+b+{\mathrm{\ξ}}_i-y_i]---\left(8\right)$

式中：α_i为引入Lagrange乘子。

根据Karush-Kuhn-Tucker(KKT)条件，可以得到如下等式：

$\left\{\begin{array}{cc}\frac{\∂L}{\∂w}=0\⇒w=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\α}}_i\mathrm{\φ}\left(x_i\right)& \\ \frac{\∂L}{\∂b}=0\⇒\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\α}}_i=0& \\ \frac{\∂L}{\∂{\mathrm{\ξ}}_i}=0\⇒{\mathrm{\α}}_i=C{\mathrm{\ξ}}_i,& i=1,L,l\\ \frac{\∂L}{\∂{\mathrm{\α}}_i}=0\⇒w^T\mathrm{\φ}\left(x_i\right)+b+{\mathrm{\ξ}}_i-y_i=0,& i=1,L,l\end{array}\right.---\left(9\right)$

通过消去w和ξ，将式(9)中的等式合并，可以得到下面的矩阵方程：

$\left[\begin{array}{cccc}0& 1& L& 1\\ 1& \mathrm{\φ}{\left(x_1\right)}^T\mathrm{\φ}\left(x_1\right)+\frac{1}{C}& L& \mathrm{\φ}{\left(x_1\right)}^T\mathrm{\φ}\left(x_l\right)\\ M& M& M& M\\ 1& \mathrm{\φ}{\left(x_l\right)}^T\mathrm{\φ}\left(x_1\right)& L& \mathrm{\φ}{\left(x_l\right)}^T\mathrm{\φ}\left(x_l\right)+\frac{1}{C}\end{array}\right]\left[\begin{array}{c}b\\ {\mathrm{\α}}_1\\ M\\ {\mathrm{\α}}_l\end{array}\right]=\left[\begin{array}{c}0\\ y_1\\ M\\ y_l\end{array}\right]---\left(10\right)$

根据Mercer条件，定义核函数：

                 K(x_i，x_j)＝φ(x_i)^Tφ(x_j)    (11)

式(10)可转化为：

$\left[\begin{array}{cc}0& 1_l^T\\ 1_l& \mathrm{\Ω}+C^{-1}I\end{array}\right]\left[\begin{array}{c}b\\ \mathrm{\α}\end{array}\right]=\left[\begin{array}{c}0\\ y\end{array}\right]---\left(12\right)$

式中：1_l和I分别代表l×1的单位向量和l×l的单位矩阵，1_l＝[1，L，1]^T；α＝[α₁，L，α_l]^T；y＝[y₁，L，y_l]^T；Ω为核矩阵，其第i行j列的元素为Ω_i，j＝K(x_i，x_j)，i，j＝1，L，l。

通过式(12)求出α和b，最后由此可以得到LS-SVM模型的回归函数：

$y\left(x\right)=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\α}}_{i}K(x,x_i)+b---\left(13\right)$

式中：核函数K(x_i，x_j)选取径向基函数K(x_i，x_j)＝exp[||x_i-x_j||²/(-2σ²)]。

(3)模型参数优化

最小二乘支持向量机模型需要优化的参数有2个，分别为正则化参数C和径向基函数的核宽度σ。采用遗传算法来优化模型参数。

①参数实数编码和种群初始化

采用实数编码，正则化参数C和核宽度σ的优化区间分别在[0，10000]和[0，10]范围内，根据这两个参数取值范围，随机产生代表LS-SVM模型参数值的初始化种群。

②适应度函数设计和计算

选择负均方根误差(-RMSE)作为适应度函数，其定义如下：

${\mathrm{Fitness}}_{\mathrm{LS}-\mathrm{SVM}}=-\mathrm{RMSE}=-\sqrt{\frac{1}{l}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{(y_{i,\mathrm{fit}}-y_{i,\mathrm{act}})}^2}---\left(14\right)$

式中：y_i，fit和y_i，act分别为第i个样本的拟合值和真实值。

③遗传操作

遗传操作需要进行选择操作、交叉操作和变异操作来产生下一代新的种群。选择操作根据个体适应度的值从种群中选择优良的个体，并淘汰掉劣质的个体，使得适应度越高的个体在子代中被选择的机会更大。交叉操作从种群内按照预先设定的交叉概率选出要进行交叉的个体，搭配成对，然后对每对个体进行交叉操作，产生子代新个体。变异操作以变异概率对交叉后产生新个体进行变异操作。

④终止条件判断

设定最大进化代数为终止条件。重复步骤②至③，直到满足终止条件。

(4)评价模型预测性能

为了验证模型的有效性，将训练好之后的LS-SVM模型使用测试样本来检验预测效果。

采用平均绝对百分误差(MAPE)和均方根误差(RMSE)作为模型预测性能的评价指标：

$\mathrm{MAPE}=100\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\left|\frac{y_{i,\mathrm{pre}}-y_{i,\mathrm{act}}}{y_{i,\mathrm{act}}}\right|---\left(15\right)$

$\mathrm{RMSE}=\sqrt{\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{(y_{i,\mathrm{pre}}-y_{i,\mathrm{act}})}^2}---\left(16\right)$

式中：M为测试样本的个数；y_i，pre和y_i，act分别为第i个样本的预测值和真实值。

3)测量炼油装置高温部位工艺介质的酸值、硫含量、温度和流速，根据建立的数学模型，进行腐蚀速率的预测。

按照建立的最小二乘支持向量机模型编写计算机程序，然后将工业现场测定的酸值、硫含量、温度和流速参数输入计算机程序，计算得到腐蚀速率V。

4)测量炼油装置高温部位的平均剩余壁厚，根据设备设计条件获得最小要求壁厚，利用得到的平均剩余壁厚、最小要求壁厚和预测的腐蚀速率，计算出腐蚀剩余寿命。

腐蚀剩余寿命的计算方法为，用当前测量的平均剩余壁厚δ_c(单位：mm)减去最小要求壁厚δ_min(单位：mm)的差值除以腐蚀速率V(单位：mm/year)，即为剩余寿命R(单位：year)，即

$R=\frac{{\mathrm{\δ}}_c-{\mathrm{\δ}}_{\min }}{V}---\left(17\right)$

实施例

下面以5％Cr钢的腐蚀寿命预测为实施例对本发明进行详细说明，但本发明并不被这些实施例所限定。

以5％Cr钢为试验材料，通过模拟试验获得炼油装置高温部位金属材料在不同的油品酸值、硫含量、温度和流速条件下的腐蚀速率，具体的试验数据见表1。

表1

序号	酸值(KOHmg/g)	硫含量(wt％)	温度(℃)	流速(m/s)	腐蚀速率(mm/year)
						1	0.5	0.1	220	0	0.0254
2	0.5	0.8	270	5	0.1016
						3	0.5	1.5	330	10	0.381
4	0.5	2.5	385	20	0.9144
						5	0.5	3.5	412	30	1.0668
6	1.0	0.1	270	10	0.127
						7	1.0	0.8	330	20	0.4318
8	1.0	1.5	385	30	0.9398
						9	1.0	2.5	412	0	1.143
10	1.0	3.5	220	5	0.1574
						11	2.0	0.1	330	30	0.7366
12	2.0	0.8	385	0	0.8892
						13	2.0	1.5	412	5	1.1684
14	2.0	2.5	220	10	0.1778
						15	2.0	3.5	270	20	0.4064
16	3.0	0.1	385	5	1.2954
						17	3.0	0.8	412	10	1.0668
18	3.0	1.5	220	20	0.254
						19	3.0	2.5	270	30	0.508
20	3.0	3.5	330	0	1.016
						21	5.0	0.1	412	20	2.160
22	5.0	0.8	220	30	0.2286

23	5.0	1.5	270	0	0.762
						24	5.0	2.5	330	5	1.321
25	5.0	3.5	385	10	1.778
						26	0.35	0.1	365	0	0.2032
27	3.1	0.55	250	5	0.1524
						28	4.8	0.96	410	10	1.27
29	1.7	1.1	265	20	0.381
						30	1.9	2.8	370	25	0.889

以上共30个数据样本，前25个数据为训练样本，后5个为测试样本。每组样本包括输入数据(酸值x₁、硫含量x₂、温度x₃、流速x₄)和输出数据(腐蚀速率y)。利用最小二乘支持向量机方法建立酸值、硫含量、温度、流速等影响因素与腐蚀速率之间相互关系的数学模型。

模型的最优参数为C＝9978.4，σ＝0.0817。拟合的MAPE为0.0124，RMSE为0.0014，MAPE和RMSE这两个评价指标的值都较小，表明该模型具有很好的拟合效果。

根据以上模型，编制计算机程序，输入测试样本的所有测试样本的输入数据(酸值、硫含量、温度、流速)，将模型的预测值和实验比较，发现二者吻合得较好，其绝对百分误差的值都比较小，满足工程需要(见表2)，这表明，本发明方法的预测效果良好，精度较高。

表2

现测得某一加工高硫高酸原油装置高温部位的油品酸值为2.0KOHmg/g，硫含量为1.0wt％，温度为300℃，流速为5m/s，该部位的材质为5％Cr钢。将这些参数输入编制的计算程序中，计算得到腐蚀速率为0.4418mm/year。

该部位的最小要求壁厚为5mm，利用超声波测厚仪测得当前的平均剩余壁厚为10.2mm，利用式(17)计算得到该部位的腐蚀剩余寿命为11.77年。

Claims (2)

1.一种炼油装置高温部位腐蚀剩余寿命预测方法，温度范围180～420℃，其特征在于：包括以下步骤：

1)通过模拟试验获得炼油装置高温部位金属材料在不同的油品酸值、硫含量、温度和流速条件下的腐蚀速率；按下式计算平均腐蚀速率V_corr：

$V_{\mathrm{corr}}(\mathrm{mm}/a)=\frac{8.76\×{10}^4}{7.8\mathrm{St}}\stackrel{\‾}{m}---\left(1\right)$

其中，

为平均失重，单位为g；S为试片面积，单位为cm²；t为试验周期，单位为h；

2)根据实验数据，运用最小二乘支持向量机方法，建立油品酸值、硫含量、温度、流速与腐蚀速率之间对应关系的数学模型；

3)测量炼油装置高温部位工艺介质的酸值、硫含量、温度和流速，根据建立的数学模型，进行腐蚀速率的预测；

4)测量炼油装置高温部位的平均剩余壁厚，根据设备设计条件获得最小要求壁厚，利用得到的平均剩余壁厚、最小要求壁厚和预测出的腐蚀速率，计算出腐蚀剩余寿命；

腐蚀剩余寿命的计算方法为，用当前测量的平均剩余壁厚δ_c，单位：mm减去最小要求壁厚δ_min，单位：mm的差值除以腐蚀速率V，单位：mm/year，即为剩余寿命R，单位：year)，即

$R=\frac{{\mathrm{\δ}}_c-{\mathrm{\δ}}_{\min }}{V}---\left(17\right)$

2.根据权利要求1所述的炼油装置高温部位腐蚀剩余寿命预测方法，其特征在于：步骤2)由如下方法实现：

(1)数据预处理和样本集划分

将每一组试验数据作为一个样本，把油品酸值、硫含量、温度、流速作为样本的输入数据，对应的腐蚀速率作为输出数据；

对于给定的全部样本集，采用下面公式对样本的输入和输出数据进行归一化处理，处理之后的数值按比例缩放到[0，1]范围内；

$x_{i,j}^{\′}=\frac{x_{i,j}-\underset{i}{\min }\left(x_{i,j}\right)}{\underset{i}{\max }\left(x_{i,j}\right)-\underset{i}{\min }\left(x_{i,j}\right)}---\left(2\right)$

$y_i^{\′}=\frac{y_i-\underset{i}{\min }\left(y_i\right)}{\underset{i}{\max }\left(y_i\right)-\underset{i}{\min }\left(y_i\right)}---\left(3\right)$

式中：i＝1，Λ，m，m为全部样本的个数；j＝1，Λ，n，n为输入数据中变量的个数；x_i，j和x′_i，j分别为归一化前后的第i个输入数据中第j个变量的值；

和

分别为归一化前的输入数据中第j个变量的最大值和最小值；y_i和y′_i分别为归一化前后的第i个输出数据的值；

和

分别为归一化前的输出数据的最大值和最小值；

将预处理之后的样本集划分为训练样本和测试样本；其中，训练样本用于模型训练，测试样本用来评价训练好之后的模型预测性能。；

(2)建立数学模型

对于给定训练样本集：

                  T＝{(x_i，y_i)}，i＝1，L，l    (4)

式中：x_i为输入数据，

y_i为输出数据，

l为训练样本的个数；

构造下面的回归函数：

                   y(x)＝w^Tφ(x)+b    (5)

式中：w为权值向量；b为偏置值；φ(·)为非线性映射函数，将输入空间映射到一个高维特征空间；

根据结构风险最小化原则，可以得到如下式的优化问题：

$\min J(w,\mathrm{\ξ})=\frac{1}{2}{w}^{T}w+\frac{1}{2}C\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\ξ}}_i^2---\left(6\right)$

等式约束条件为：

                y_i＝w^Tφ(x_i)+b+ξ_i，i＝1，L，l    (7)

式中：C和ξ_i分别为正则化参数和误差变量，C是用于确定在训练误差和模型复杂度之间的权衡；

为了求解上述优化问题，构造Lagrange函数：

$L(w,b,\mathrm{\ξ},\mathrm{\α})=J(w,\mathrm{\ξ})-\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\α}}_i[w^T\mathrm{\φ}\left(x_i\right)+b+{\mathrm{\ξ}}_i-y_i]---\left(8\right)$

式中：L为Lagrange函数；α_i为引入Lagrange乘子；

根据Karush-Kuhn-Tucker(KKT)条件，可以得到如下等式：

$\left\{\begin{array}{cc}\frac{\∂L}{\∂w}=0\⇒w=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\α}}_i\mathrm{\φ}\left(x_i\right)& \\ \frac{\∂L}{\∂b}=0\⇒\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\α}}_i=0& \\ \frac{\∂L}{\∂{\mathrm{\ξ}}_i}=0\⇒{\mathrm{\α}}_i=C{\mathrm{\ξ}}_i,& i=1,L,l\\ \frac{\∂L}{\∂{\mathrm{\α}}_i}=0\⇒w^T\mathrm{\φ}\left(x_i\right)+b+{\mathrm{\ξ}}_i-y_i=0,& i=1,L,l\end{array}\right.---\left(9\right)$

通过消去w和ξ，将式(9)中的等式合并，可以得到下面的矩阵方程：

$\left[\begin{array}{cccc}0& 1& L& 1\\ 1& \mathrm{\φ}{\left(x_1\right)}^T\mathrm{\φ}\left(x_1\right)+\frac{1}{C}& L& \mathrm{\φ}{\left(x_1\right)}^T\mathrm{\φ}\left(x_l\right)\\ M& M& M& M\\ 1& \mathrm{\φ}{\left(x_l\right)}^T\mathrm{\φ}\left(x_1\right)& L& \mathrm{\φ}{\left(x_l\right)}^T\mathrm{\φ}\left(x_l\right)+\frac{1}{C}\end{array}\right]\left[\begin{array}{c}b\\ {\mathrm{\α}}_1\\ M\\ {\mathrm{\α}}_l\end{array}\right]=\left[\begin{array}{c}0\\ y_1\\ M\\ y_l\end{array}\right]---\left(10\right)$

根据Mercer条件，定义核函数：

                K(x_i，x_j)＝φ(x_i)^Tφ(x_j)    (11)

式(10)可转化为：

$\left[\begin{array}{cc}0& 1_l^T\\ 1_l& \mathrm{\Ω}+C^{-1}I\end{array}\right]\left[\begin{array}{c}b\\ \mathrm{\α}\end{array}\right]=\left[\begin{array}{c}0\\ y\end{array}\right]---\left(12\right)$

式中：1_j和I分别代表l×1的单位向量和l×l的单位矩阵，1_l＝[1，L，1]^T；α＝[α₁，L，α_l]^T；y＝[y₁，L，y_l]^T；Ω为核矩阵，其第i行j列的元素为Ω_i，j＝K(x_i，x_j)，i，j＝1，L，l；

通过式(12)求出α和b，最后由此可以得到LS-SVM模型的回归函数：

$y\left(x\right)=\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{\mathrm{\α}}_{i}K(x,x_i)+b---\left(13\right)$

式中：核函数K(x_i，x_j)选取径向基函数K(x_i，x_j)＝exp[||x_i-x_j||²/(-2σ²)]。

(3)模型参数优化

最小二乘支持向量机模型需要优化的参数有2个，分别为正则化参数C和径向基函数的核宽度σ；采用遗传算法来优化模型参数；

①参数实数编码和种群初始化

采用实数编码，正则化参数C和核宽度σ的优化区间分别在[0，10000]和[0，10]范围内，根据这两个参数取值范围，随机产生代表LS-SVM模型参数值的初始化种群；

②适应度函数设计和计算

选择负均方根误差-RMSE作为适应度函数，其定义如下：

${\mathrm{Fitness}}_{\mathrm{LS}-\mathrm{SVM}}=-\mathrm{RMSE}=-\sqrt{\frac{1}{l}\underset{i=1}{\overset{l}{\mathrm{\Σ}}}{(y_{i,\mathrm{fit}}-y_{i,\mathrm{act}})}^2}---\left(14\right)$

式中：y_i，fit和y_i，act分别为第i个样本的拟合值和真实值；

③遗传操作

遗传操作需要进行选择操作、交叉操作和变异操作来产生下一代新的种群；选择操作根据个体适应度的值从种群中选择优良的个体，并淘汰掉劣质的个体，使得适应度越高的个体在子代中被选择的机会更大；交叉操作从种群内按照预先设定的交叉概率选出要进行交叉的个体，搭配成对，然后对每对个体进行交叉操作，产生子代新个体；变异操作以变异概率对交叉后产生新个体进行变异操作；

④终止条件判断

设定最大进化代数为终止条件；重复步骤②至③，直到满足终止条件；

(4)评价模型预测性能

为了验证模型的有效性，将训练好之后的LS-SVM模型使用测试样本来检验预测效果；

采用平均绝对百分误差(MAPE)和均方根误差(RMSE)作为模型预测性能的评价指标：

$\mathrm{MAPE}=100\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\left|\frac{y_{i,\mathrm{pre}}-y_{i,\mathrm{act}}}{y_{i,\mathrm{act}}}\right|---\left(15\right)$

$\mathrm{RMSE}=\sqrt{\frac{1}{M}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{(y_{i,\mathrm{pre}}-y_{i,\mathrm{act}})}^2}---\left(16\right)$

式中：M为测试样本的个数；y_i，pre和y_i，act分别为第i个样本的预测值和真实值。

Images