一种VERNA-GA的汽油调合优化调度方法

摘要

本发明公开了一种VERNA-GA的汽油调合优化调度方法，该方法首先在高预测精度的调合模型和生产调度约束条件下，根据调合组分油、产品油和市场供求的信息，构成以利润和满足约束为目标的汽油调合调度模型；然后采用VERNA-GA方法，其计算结构包括主群体选择、交叉、变异操作和病毒群体感染、复制、删减操作，依方法中提供的算法参数取值范围进行参数设置并搜索优化配方；按所得配方进行汽油调合优化调度生产。本发明提供了一天和多天的汽油调合优化调度方案，无需配方经验，在生产约束范围内随机搜索配方，可快速找出质量卡边、利润高的配方，易于操作及优化调度规模的扩展。

说明书

技术领域

本发明涉及化工过程优化领域，尤其涉及一种病毒进化RNA遗传算法（Viral Evolution based RNA Genetic Algorithm(VERNA-GA)）的汽油调合优化调度方法。 

背景技术

随着社会的发展，石油炼制工业呈现出规模大型化、技术现代化和品种多样化的特点，其生产能力、产品质量和品种持续稳定地增长。出于技术经济的综合考虑，加上炼油装置工艺的局限性，各炼油装置生产的一次产品油性能一般都不能直接满足各种油品质量的要求，如汽油、柴油、润滑油类产品质量的要求。一次产品油常常称为半成品油或基础油等。为了降低成本、节约能源、提高效率、优化工艺，常常需要在一次产品油中加入添加剂，或通过双组分、多组分半产品油按不同比例的调合，充分利用不同组分油的物化性质，发挥各自的优良性能，相互取长补短，以达到用户要求的产品质量。汽油、柴油的质量升级和润滑油的高质量要求，使炼油厂为满足新的质量要求而付出高昂的代价。为此，应该通过油品调合手段，在满足汽油、柴油和润滑油指标的条件下，最大限度地将生产过程中产生的各种组分汽油、柴油及其他基础油，按一定的配方进行调合而生产出成本最低、质量合格的高品质汽油、柴油。 

汽油调合过程是炼油厂装置生产的汽油组分混合的过程，同时加入少量添加剂和抗腐蚀剂成为满足一定质量标准的混合物。需要满足包括调合组分量和产品存储设备的限制等大量操作约束，当满足所有产品质量和调合操作约束时，汽油调合配方优化和控制的目标是最大化调合利润。在实际生产过程中，由于产品的指标多为各组分油性质的非线性函数，同时又要满足各组分油的物料平衡，并且在很多情况下必须要考虑到价格、产品性能指标等参数的不确定性的问题，因此汽油调合过程是一个包含不确定性的复杂的具有多个线性和非线性不等式约束的非线性规划问题。 

作为炼油厂重要的产品，调合效益在生产企业的经济效益中占有举足轻重的地位，是影响炼油厂经济效益的重要因素之一。一个优化的汽油调合方案，不仅能够减少原油加工成本，降低库存水平，更重要的是能提高资源的利用率，增加产品收率，达到满足需求和最大化利润的目的，因此油品调合问题的解决对提高炼油厂效益具有极其深远的意义。 

发明内容

本发明的目的在于克服现有技术的不足，提供一种VERNA-GA汽油调合优化调度方法。 

本发明的目的是通过以下技术方案来实现的：一种VERNA-GA汽油调合优化调度方法，它包括如下步骤： 

（1）根据调合组分油、产品油和市场供求的信息，选定油品质量稳定的调合组分油，调合组分油包括催化重整油、直馏石脑油、催化裂化汽油、丁烷、烷基化油、加氢裂化汽油、高效添加剂； 

（2）根据选定组分油将调合组分油的质量指标、成本和供应量，与产品汽油的成本、价格和质量指标构成目标函数参数与组分油供应量和产品汽油产量不等式约束，调合模型预测值和实测值之间的裕度设为0-0.3； 

调合模型中辛烷值计算如下式所示： 

$\left(\begin{array}{c}{P}_{\mathrm{RON}}=r^{T}x+{\alpha }_1({\mathrm{cr}}^T\mathrm{diag}\left(s\right)x-\frac{\left(r^{T}x\right)\left(s^{T}x\right)}{e^{T}x})\\ +{\alpha }_2(o_s^{T}x-\frac{{\left(o^{T}x\right)}^2}{e^{T}x})+{\alpha }_3(a_s^{T}x-\frac{{\left(a^{T}x\right)}^2}{e^{T}x})\end{array}\right)$

$P_{\mathrm{MON}}=m^{T}x+{\alpha }_4(m^T\mathrm{diag}\left(s\right)x-\frac{\left(m^{T}x\right)\left(s^{T}x\right)}{e^{T}x})+$

${\alpha }_5(o_s^{T}x-\frac{{\left(o^{T}x\right)}^2}{e^{T}x})+\frac{{\alpha }_6}{10000\left(e^{T}x\right)}{(a_s^{T}x-\frac{{\left(a^{T}x\right)}^2}{e^{T}x})}^2$

式中，x为组分流率，r为每种组分的研究法辛烷值（RON），m为每种组分的马达法辛烷值（MON），s=r-m，o是每种组分的蜡含量，o_s是蜡含量的平方，e为单位列向量，a为参调组分的芳烃含量，a_s是芳烃含量的平方，α_i为模型参数（其中，α₁=0.03224，α₂=0.00101，α₃=0，α₄=0.0445，α₅=0.00081，α₆=-0.00645），P_RON为产品的RON值，P_MON为产品的MON值； 

雷德蒸气压（RVP）计算如下所示： 

$P_{\mathrm{RVP}}={\left(\underset{i=1}{\overset{n}{\Sigma }}{u}_i{\left({\mathrm{RVP}}_i\right)}^{1.25}\right)}^{0.8}$

式中，u_i为i组分在产品中的比例含量，n为参调组分数量,P_RVP为产品的RVP； 

模型中产品质量指标和产品组分两方面的约束条件如下所述： 

产品辛烷值指标约束A：P_RON,n≥RON_n，min； 

产品辛烷值指标约束B：P_MON,n≥MON_n，min； 

产品蒸汽压指标约束：P_RVP，n≤RVP_n，max； 

产品和调合组分之间的物料平衡约束:

调合产品的物料平衡约束：${\mathrm{Vp}}_{n,j}={\mathrm{Vp}}_{n,j-1}+\underset{m=1}{\overset{{\mathrm{Nc}}_n}{\Sigma }}{\mathrm{Vc}}_{n,m,j};$

组分储量约束：${\mathrm{Vc}}_{m,L}\le \underset{n=1}{\overset{N_p}{\Sigma }}{\mathrm{Vc}}_{n,m}\le {\mathrm{Vc}}_{m,U};$

产品需求量约束：Vp_n=Mp_n； 

上述各式中，P_RON,n为产品n的RON值，RON_n，min为产品n的RON最小值，P_MON,n为产品n的MON值，MON_n,min产品n的MON最小值，P_RVP,n为产品n的RVP值，RVP_n，max为产品n的RVP最大值，Vp_n，j为产品n第j天的储量，Vc_n，m，j为第j天参与调合产品n的组分m的体积，Vc_m，L为组分m的罐存储量下限，Vc_m，U为组分m的罐存储量上限，Mp_n为产品n的市场需求量； 

在满足组分油物料平衡和产品指标的条件下，使产品的利润目标函数由下式表示： 

$\max f\left(x\right)=\underset{t=1}{\overset{T}{\Sigma }}\underset{n=1}{\overset{N_p}{\Sigma }}({\mathrm{Cp}}_{n,t}·{\mathrm{Vp}}_{n,t}-\underset{m=1}{\overset{{\mathrm{Nc}}_n}{\Sigma }}{\mathrm{Cc}}_{m,t}·{\mathrm{Vc}}_{n,m,t});$

式中，等式右边第一项为调合产品的价值，第二项为参调组分的成本；T为调合调度的时间尺度，Nc_n是产品n的组分种类，Np为产品种类，Vc_n，m,t为产品n中组分m的含量，Vp_n，t是产品n的调合量，Cc_m，t为组分m的成本，Cp_n，t为产品n的价格； 

利润目标函数加上惩罚项后总目标函数为： 

$f^{\prime }(x,\mu )=\max f\left(x\right)+\mu \underset{i=1}{\overset{m}{\Sigma }}\min {(0,g_i\left(x\right))}^2;$

式中，μ是惩罚系数，g_i(x)是约束函数； 

（3）根据产品需求和实际数据设置参数：运行的最大代数G_max，参数编码长度l，个体编码长度L为每个参数编码长度l乘以参数的个数n，种群数N，病毒数M，自交叉概率p_adr，置换交叉概率p_C,位置自适应变异概率p_m，感染操作概率p_infect、复制操作概率p_copy、删减操作概率p_cut，衰减系数decay，及以终止准则； 

（4）运行VERNA-GA，对化工动力学模型中的总生产利润进行求解，通过最大化目标函数，得到满足物料平衡和产品指标要求的一定配方比例的组分油使用量。 

进一步地，所述步骤（3）中，所述终止准则为：运行代数达最大代数G_max或者得到的目标函数值达到实际要求。 

进一步地，所述的步骤（4）包括以下子步骤： 

（4.1）随机生成包含N个长度为L的RNA序列的初始种群，每一个RNA序列代表组分油使用量参数的可能解；通过无感染能力的病毒个体去复制主个体部分基因操作来产生M个病毒个体；组分油使用量参数均由字符集{0,1,2,3}编码为一个长度为l的RNA子序列，设未知参数个数为n，则一个RNA序列的编码产度为L=ln； 

（4.2）将每一个RNA序列解码为组分油使用量参数并计算其对应的目标函数值，再将目标函数值转化为对应的适应度值，适应度值最大的个体定义为这一代的最优个体； 

（4.3）根据适应度值优劣程度进行排序，将种群分为两类，一半为优质群体SuG，另一半为劣质群体InG，排名第一的个体为最优个体，定义Ncnew为经交叉操作产生新个体数目，令Ncnew=0；如果最优解连续g₁代无变化，则切换到SQP算法进行寻优。 

（4.4）执行自交叉操作、置换交叉操作、位置自适应变异操作，以及病毒感染操作； 

（4.5）重复步骤（4.2）至步骤（4.4）进行迭代，直到满足终止准则； 

（4.6）将最终种群的最优个体解码为组分油使用量。 

进一步地，所述步骤（4.4）中的RNA交叉操作为： 

（A）随机从劣质群体InG中选择一个RNA序列作为父代个体，以概率p_adr执行RNA自交叉操作，产生一个新个体,Ncnew=Ncnew+1； 

（B）随即从优质群体SuG中选择两条RNA序列作为父代个体，以概率p_c对它们置换交叉操作，产生两个子代个体，选择其中适应度高的一个个体Ncnew=Ncnew+1； 

（C）重复步骤A和B，直到产生的新个体数达到N/2，即Ncnew>0.5Size。 

进一步地，所述自交叉操作具体为：随机从劣质群体InG中选择一个RNA序列作为父代个体，然后在这个个体上随机选择一个子序列，子序列长度在[1,L]之间，对这个子序列进行重新排序，最后将子序列插回原来位置，得到一个子代个体； 

所述置换交叉操作具体为：随即从优质群体SuG中选择两个个RNA序列作为父代个体，再在各个个体中选取一段病毒或酶转座子序列进行交叉操作；操作中选取的一段序列，其位置和含有碱基数目都是随机的，但是两个父代个体 间序列的碱基数必须相同。 

进一步地，所述的步骤（4.4）中的自适应变异操作具体为：根据每个碱基在个体结构中的位置信息i计算p_m，即随机产生一个位置i，将种群个体序列分成两部分，高于i位部分和低于i位部分序列里的碱基分别以概率p_mh和p_ml进行变异操作，变异概率p_mh和p_ml依进化代数g改变，总变异概率p_m依碱基位置i变化： 

$\left(\begin{array}{c}{p}_{\mathrm{mh}}=a_1+\frac{b_1}{1+\exp (-a(g-g_0))}\\ p_{\mathrm{ml}}=a_2+\frac{b_2}{1+\exp (-a(g-g_0))}\\ p_m=\frac{i}{L}[\alpha ,\beta ]\left(\begin{array}{c}{p}_{\mathrm{mh}}\\ p_{\mathrm{ml}}\end{array}\right)\end{array}\right);\left(\begin{array}{c}\alpha =1\\ \beta =0,\end{array}\right)1=<i=<\frac{L}{2};\left(\begin{array}{c}\alpha =0\\ \beta =1,\end{array}\right)\frac{L}{2}<i=<L$

式中，i为在个体中随机产生的一个碱基位置，L为个体长度，g为当前进化代数，g₀为hotspot和coldspot的转折点，a为变化速率，a₁和a₂分最别为p_mh和p_ml终变异概率，b₁和b₂分别为p_mh和p_ml的变化范围。 

进一步地，所述的步骤（4.4）中的病毒感染操作具体为：对于病毒群体viruspop(t)中的M个病毒个体依次以概率p_infect对主群体hostpop(t)中的个体进行感染操作，包括如下步骤： 

（a）假设当前进化代数为t，病毒个体为i，令i＝1； 

（b）对于viruspop(t)中个体i以概率p_infect对hostpop(t)中每个个体执行感染操作，hostpop(t)中被病毒i感染后的群体为U，如果感染后群体中个体数为零，重复此步骤； 

（c）对U中每个个体k计算主个体k被病毒i感染前的适应度为计算如果则用个体k替换hostpop(t)中相应的个体； 

（d）计算如果则从U中随机选择一个主个体对病毒个体i执行概率为p_copy的复制操作；否则，对病毒个体i执行概率为p_cut的删减操作，将病毒个体i记入viruspop(t+1)； 

（e）根据病毒生命衰减系数decay，计算病毒生命力若life_i(t+1)<0，则从hostpop(t)中随机选择一主个体执行概率为p_copy的复制操作，产生新的病毒个体替换viruspop(t+1)中的病毒个体i且令λ_i(t+1)=0，life_i(t+1)=0； 

(f）令i＝i+1，如果i≤M，转至步骤2，否则结束循环，转至RNA遗传操作。 

本发明的有益效果是，本发明VERNA-GA汽油调合优化调度方法既实现了遗传操作在父子代群体间纵向继承信息进行全局搜索的功能，也实现了病毒感染操作在同一代群体中横向传播进化进行局部搜索的功能，因此能较快获得问题的满意解；同时在交叉和变异操作上的改进，使得收敛效率和精度上都达到较好的结果；提供的汽油调合优化调度方案实现了随机搜索，无需经验值来设置配方初始值寻找最优解，计算速度快，成品油质量卡边；方法中使用的约束处理方式，使方法使用直观、易于操作；用搜索得到的配方结果指导采购、生产和存储，可实现最小库存。 

附图说明

图1是复制操作图； 

图2是感染操作图； 

图3是剪切操作图； 

图4是汽油调合过程流程图； 

具体实施方式

下面根据附图和实施例详细描述本发明，本发明的目的和效果将变得更加明显。 

本发明中病毒协同进化的思想：将病毒个体看作具有生命力的个体，以一定的概率依次对主群体中的个体执行感染操作，如果感染后主群体个体适应度值增大，那么就用感染后的个体替换主群体中对应的个体，根据感染效果的好坏，病毒个体会执行复制或删减操作动态调整病毒个体适应度和感染率，并不断地对主群体进行感染进化，直到满足结束条件。病毒进化包括三种基本操作：感染操作、复制操作、删减操作。 

1、感染操作(Infection)：取一个病毒个体以概率p_infect来感染一个主体，感染后如果发现主体的适应度增大了，则用感染后的主个体替换感染前的主个体，操作过程如图1所示。 

2、复制操作(Copy)：随机选择一个主个体以概率p_copy代换病毒染色体中对应的基因，产生新的有效字符更多的一个病毒个体。病毒个体的初始化也是通过无感染能力的病毒个体去复制主个体部分基因操作来生一个新的病毒个体。操作过程如图2所示。 

3、删减操作(Cut)：以删减操作概率p_cut将病毒个体某位换成通配符号来代替。删减操作目的主要是在于减少这个病毒的部分活动的资源模式的数量从而来降低这个病毒的进化计算能力并得以动态调整感染率，操作过程如图3所示。 

汽油调合系统主要由组分存储罐、调合器、成品罐和分析仪四个设备组成，这些设备通过管道、流量计和阀门连接在一起，调合过程如图4所示。图中，重整油、直馏油、丁烷、催化油和烷基化油五种组分油按照一定比例调合，产生普通汽油和高级汽油两种产品。汽油产品的质量标准包括辛烷值、挥发性、硫含量、芳烃含量、雷德蒸汽压(RVP)等，但是实际生产中最重要的油品质量指标是辛烷值和蒸气压。辛烷值主要有两种形式：马达法辛烷值（Motor Octane Number，MON）和研究法辛烷值（Research Octane Number，RON），表示汽油的抗爆特性。使用ETHYL RT-70法计算汽油调合产品辛烷值，如下式所示： 

$\left(\begin{array}{c}{P}_{\mathrm{RON}}=r^{T}x+{\alpha }_1({\mathrm{cr}}^T\mathrm{diag}\left(s\right)x-\frac{\left(r^{T}x\right)\left(s^{T}x\right)}{e^{T}x})\\ +{\alpha }_2(o_s^{T}x-\frac{{\left(o^{T}x\right)}^2}{e^{T}x})+{\alpha }_3(a_s^{T}x-\frac{{\left(a^{T}x\right)}^2}{e^{T}x})\end{array}\right)---\left(1.1\right)$

$P_{\mathrm{MON}}=m^{T}x+{\alpha }_4(m^T\mathrm{diag}\left(s\right)x-\frac{\left(m^{T}x\right)\left(s^{T}x\right)}{e^{T}x})+$$\left(1.2\right)$

${\alpha }_5(o_s^{T}x-\frac{{\left(o^{T}x\right)}^2}{e^{T}x})+\frac{{\alpha }_6}{10000\left(e^{T}x\right)}{(a_s^{T}x-\frac{{\left(a^{T}x\right)}^2}{e^{T}x})}^2$

式中，x为组分流率，r为每种组分的RON，m为每种组分的MON，s=r-m，o是每种组分的蜡含量，o_s是蜡含量的平方，e为单位列向量，a为参调组分的芳烃含量，a_s是芳烃含量的平方，α_i为模型参数（其中，α₁=0.03224，α₂=0.00101，α₃=0，α₄=0.0445，α₅=0.00081，α₆=-0.00645），P_RON为产品的RON值，P_MON为产品的MON值。 

我们使用雷德蒸气压调合指数来计算雷德蒸气压（RVP），具体公式如下所示： 

$P_{\mathrm{RVP}}={\left(\underset{i=1}{\overset{n}{\Sigma }}{u}_i{\left({\mathrm{RVP}}_i\right)}^{1.25}\right)}^{0.8}---\left(1.3\right)$

式中，u_i为组分在i产品中的比例含量，n为参调组分数量，P_RVP为产品的RVP。汽油调合问题是按一定配方比例使用组分油，在满足组分油物料平衡和产品指标的条件下，使产品的利润最高，优化目标可由下式表示： 

$\max f\left(x\right)=\underset{t=1}{\overset{T}{\Sigma }}\underset{n=1}{\overset{N_p}{\Sigma }}({\mathrm{Cp}}_{n,t}·{\mathrm{Vp}}_{n,t}-\underset{m=1}{\overset{{\mathrm{Nc}}_n}{\Sigma }}{\mathrm{Cc}}_{m,t}·{\mathrm{Vc}}_{n,m,t})---\left(1.4\right)$

式中，等式右边第一项为调合产品的价值，第二项为参调组分的成本。T为调合 调度的时间尺度，Nc_n是产品n的组分种类，Np为产品种类，Vc_n，m,t为产品n中组分m的含量，Vp_n，t是产品n的调合量，Cc_m，t为组分m的成本，Cp_n，t为产品n的价格。 

汽油调合的配方优化应满足产品质量指标和产品组分两方面的约束条件，前者包括产品的辛烷指标约束和蒸汽压指标约束，后者包括产品与调合组分之间的物料平衡约束，产品的物料平衡约束、组分储量约束、产品需求量约束，下面列出具体约束指标。 

产品辛烷值指标约束A 

P_RON,n≥RON_n，min                                           （1.5） 

产品辛烷值指标约束B 

P_MON,n≥MON_n，min                                           （1.6） 

产品蒸汽压指标约束 

P_RVP，n≤RVP_n,max                                           （1.7） 

产品和调合组分之间的物料平衡约束 

$\underset{n=1}{\overset{N_p}{\Sigma }}{\mathrm{Vp}}_n=\underset{n=1}{\overset{N_p}{\Sigma }}\underset{m=1}{\overset{{\mathrm{Nc}}_n}{\Sigma }}{\mathrm{Vc}}_{n,m}---\left(1.8\right)$

调合产品的物料平衡约束 

${\mathrm{Vp}}_{n,j}={\mathrm{Vp}}_{n,j-1}+\underset{m=1}{\overset{{\mathrm{Nc}}_n}{\Sigma }}{\mathrm{Vc}}_{n,m,j}---\left(1.9\right)$

组分储量约束 

${\mathrm{Vc}}_{m,L}\le \underset{n=1}{\overset{N_p}{\Sigma }}{\mathrm{Vc}}_{n,m}\le {\mathrm{Vc}}_{m,U}---\left(1.10\right)$

产品需求量约束 

Vp_n=Mp_n        （1.11） 

上述各式中，P_RON，n为产品n的RON值，RON_n，min为产品n的RON最小值，P_MON,n为产品n的MON值，MON_n，min产品n的MON最小值，P_RVP，n为产品n的RVP值，RVP_n,max为产品n的RVP最大值，Vp_n，j为产品n第j天的储量，Vc_n，m,j为第j天参与调合产品n的组分m的体积，Vc_m，L为组分m的罐存储量下限，Vc_m，U为组分m的罐存储量上限，Mp_n为产品n的市场需求量。 

经以上分析，需要油品调和过程中成品油和组分油的参数以及需求量等相关参数如表1、表2和表3所示。 

表1：成品油参数 

  价格$/bbl RON最小值 MON最小值 RVP最大值

[0097] 

普通汽油 33.0 88.5 77.0 10.8 高级汽油 37.0 91.5 80.0 10.8

表2：组分油参数 

表3：三天产品需求量参数 

本发明基于VERNA-GA汽油调合优化调度方法，步骤如下： 

1、在高预测精度的调合模型和生产调度约束条件下，根据调合组分油、产品油和市场供求的信息，构成以利润和满足约束为目标的汽油调合调度模型，将总生产利润作为VERNA-GA寻优搜索的目标函数。 

本实施例中，根据表1数据以及公式(1.1)-(1.11)得目标函数式： 

maxf(x)=33.0x_r+37.0x_h-34.0(x_r1+x_h1)-26.0(x_r2+x_h2)-10.3(x_r3+x_k3)-31.3(x_r4+x_h4)-37.0(x_r5+x_k5)式中，x_r为普通汽油产量，x_h为高级汽油产量，x_ri为调合低级汽油过程中用到的各组分油的量，x_hi则为调合高级汽油过程中用到的各组分油的量，i＝1～5。 

物料平衡约束： 

y₁＝x₁+x₂+x₃+x₄+x₅； 

y₂＝x₆+x₇+x₈+x₉+x₁₀； 

成品油需求约束： 

g₁(x)=y₁-7000； 

g₂(x)=8000-y₁； 

g₃(x)=10000-y₁； 

g₄(x)=y₂-10000。 

组分油供应不等式约束： 

g₅(x)=12000-(x₁+x₆)； 

g₆(x)=6500-(x₂+x₇)； 

g₇(x)=3000-(x₃+x₈)； 

g₈(x)=4500-(x₄+x₉)； 

g₉(x)=4500-(x₅+x₁₀)；(x_i≥0,i=1,...,10) 

成品油指标不等式约束： 

g₁₀(x)=RON₁-88.5； 

g₁₁(x)=MON₁-77.0； 

g₁₂(x)=10.8-RVP₁； 

g₁₃(x)=RON₂-91.5； 

g₁₄(x)=MON₂-80.0； 

g₁₅(x)=10.8-RVP₂。 

RON₁、MON₁、RVP₁、RON₂、MON₂、RVP₂，分别按下式计算： 

$\left(\begin{array}{c}{P}_{\mathrm{RON}}=r^{T}x+{\alpha }_1({\mathrm{cr}}^T\mathrm{diag}\left(s\right)x-\frac{\left(r^{T}x\right)\left(s^{T}x\right)}{e^{T}x})\\ +{\alpha }_2(o_s^{T}x-\frac{{\left(o^{T}x\right)}^2}{e^{T}x})+{\alpha }_3(a_s^{T}x-\frac{{\left(a^{T}x\right)}^2}{e^{T}x})\end{array}\right)$

$P_{\mathrm{MON}}=m^{T}x+{\alpha }_4(m^T\mathrm{diag}\left(s\right)x-\frac{\left(m^{T}x\right)\left(s^{T}x\right)}{e^{T}x})+$

${\alpha }_5(o_s^{T}x-\frac{{\left(o^{T}x\right)}^2}{e^{T}x})+\frac{{\alpha }_6}{10000\left(e^{T}x\right)}{(a_s^{T}x-\frac{{\left(a^{T}x\right)}^2}{e^{T}x})}^2$

$P_{\mathrm{RVP}}={\left(\underset{i=1}{\overset{n}{\Sigma }}{u}_i{\left({\mathrm{RVP}}_i\right)}^{1.25}\right)}^{0.8}$

其中，α₁=0.03224，α₂=0.00101，α₃=0，α₄=0.0445，α₅=0.00081，α₆=-0.00645。将以上利润目标和约束代入如下目标函数计算： 

$f^{\prime }(x,\mu )=\max f\left(x\right)+\mu \underset{i=1}{\overset{m}{\Sigma }}\min {(0,g_i\left(x\right))}^2$

其中，惩罚项通常比f大10倍左右，当惩罚项约束值越小或者惩罚项代表的约束越重要，μ取值越大，在此例中μ=10000。 

2、根据产品需求和实际数据设置参数：运行的最大代数G_max，参数编码长度l，个体编码长度L为每个参数编码长度l乘以参数的个数n，种群数N，病 毒数M，自交叉概率p_adr，置换交叉概率p_c，位置自适应变异概率p_m，感染操作概率p_infect、复制操作概率p_copy、删减操作概率p_cut，衰减系数decay，以终止准则。 

本实施例中，设置运行参数：最大代数G_max为1000，个体编码长度L为20，种群大小N为60，病毒数M为5，自交叉概率p_adr＝0.5，置换交叉概率p_c＝0.5，感染操作概率p_infect=0.01，复制操作概率p_copy=0.1，删减操作概率p_cut=0.1，生命力递减率decay＝0.9，最优解连续无变化切换至SQP算法代数为g₁＝50，自适应变异概率p_m依据下式设置： 

$\left(\begin{array}{c}{p}_{\mathrm{mh}}=a_1+\frac{b_1}{1+\exp (-a(g-g_0))}\\ p_{\mathrm{ml}}=a_2+\frac{b_2}{1+\exp (-a(g-g_0))}\\ p_m=\frac{i}{L}[\alpha ,\beta ]\left(\begin{array}{c}{p}_{\mathrm{mh}}\\ p_{\mathrm{ml}}\end{array}\right)\end{array}\right);\left(\begin{array}{c}\alpha =1\\ \beta =0,\end{array}\right)1=<i=<\frac{L}{2};\left(\begin{array}{c}\alpha =0\\ \beta =1,\end{array}\right)\frac{L}{2}<i=<L$

式中，i为在个体中随机产生的一个碱基位置，L为个体长度20，g为当前进化代数，g₀为hotspot和coldspot的转折点，a为变化速率，a₁和a₂分最别为p_mh和p_ml终变异概率，b₁和b₂分别为p_mh和p_ml的变化范围，a₁=a₂=0.02,b₁=b₂=0.2，g₀=G/2，a＝20/G,G为最大进化代数1000。 

3、运行VERNA-GA，对化工动力学模型中的总生产利润进行求解，通过最大化目标函数，得到满足物料平衡和产品指标要求的一定配方比例的组分油使用量。 

该步骤通过以下子步骤来实现： 

3.1随机生成包含N个长度为L的RNA序列的初始种群，每一个RNA序列代表组分油使用量参数的可能解。通过无感染能力的病毒个体去复制主个体部分基因操作来产生M个病毒个体。组分油使用量参数均由字符集{0,1,2,3}编码为一个长度为l的RNA子序列，设未知参数个数为n，则一个RNA序列的编码产度为L=ln； 

3.2将每一个RNA序列解码为组分油使用量参数并计算其对应的目标函数值，再将目标函数值转化为对应的适应度值，适应度值最大的个体定义为这一代的最优个体。 

3.3根据适应度值优劣程度进行排序，将种群分为两类，一半为优质群体 SuG，另一半为劣质群体InG，排名第一的个体为最优个体，定义Ncnew为经交叉操作产生新个体数目，令Ncnew=0。如果最优解连续g₁代无变化，则切换到SQP算法进行寻优。 

3.4执行自交叉操作、置换交叉操作、自适应变异操作，以及病毒感染操作，如图4所示。 

3.5重复步骤3.2至步骤3.4进行迭代，直到满足终止准则； 

3.6将最终种群的最优个体解码为组分油使用量。 

所述步骤3.4中的RNA交叉操作具体为： 

3.4.1随机从劣质群体InG中选择一个RNA序列作为父代个体，以概率p_adr执行RNA自交叉操作，产生一个新个体,Ncnew=Ncnew+1； 

3.4.2随即从优质群体SuG中选择两条RNA序列作为父代个体，以概率p_c对它们置换交叉操作，产生两个子代个体，选择其中适应度高的一个个体Ncnew=Ncnew+1； 

3.4.3重复步骤3.4.1和3.4.2，直到产生的新个体数达到N/2，即Ncnew>0.5Size。 

所述的步骤3.4中的自适应变异操作为： 

根据每个碱基在个体结构中的位置信息i计算p_m，即随机产生一个位置i，将种群个体序列分成两部分，高于i位部分和低于i位部分序列里的碱基分别以概率p_mh和p_ml进行变异操作，变异概率p_mh和p_ml依进化代数g改变，总变异概率p_m依碱基位置i变化： 

$\left(\begin{array}{c}{p}_{\mathrm{mh}}=a_1+\frac{b_1}{1+\exp (-a(g-g_0))}\\ p_{\mathrm{ml}}=a_2+\frac{b_2}{1+\exp (-a(g-g_0))}\\ p_m=\frac{i}{L}[\alpha ,\beta ]\left(\begin{array}{c}{p}_{\mathrm{mh}}\\ p_{\mathrm{ml}}\end{array}\right)\end{array}\right);\left(\begin{array}{c}\alpha =1\\ \beta =0,\end{array}\right)1=<i=<\frac{L}{2};\left(\begin{array}{c}\alpha =0\\ \beta =1,\end{array}\right)\frac{L}{2}<i=<L$

式中，i为在个体中随机产生的一个碱基位置，L为个体长度，g为当前进化代数，g₀为hotspot和coldspot的转折点，a为变化速率，a₁和a₂分最别为p_mh和p_ml终变异概率，b₁和b₂分别为p_mh和p_ml的变化范围。 

所述的步骤3.4中的病毒感染操作步骤为： 

A、假设当前进化代数为t，病毒个体为i，令i＝1。 

B、对于viruspop(t)中个体i以概率p_infect对hostpop(t)中每个个体执行感染操作，hostpop(t)中被病毒i感染后的群体为U。 

C、对U中每个个体k计算主个体k被病毒i感染前的适应度为计算如果则用个体k替换hostpop(t)中相应的个体。 

D、计算如果则从U中随机选择一个主个体对病毒个体i执行概率为p_copy的复制操作；否则，对病毒个体i执行概率为p_cut的删减操作，将病毒个体i记入viruspop(t+1)。 

E、计算life_i(t+1)，若life_i(t+1)<0，则从hostpop(t)中随机选择一主个体执行概率为p_copy的复制操作，产生新的病毒个体替换viruspop(t+1)中的病毒个体i且令λ_i(t+1)=0，life_i(t+1)=0。 

F、令i＝i+1，如果i≤M，转至步骤2，否则结束子循环，并转至RNA遗传操作。 

本实施例中，对满足约束的汽油调合优化调度模型的总生产利润进行求解，通过最大化目标函数，得到满足物料平衡和产品指标要求的一定配方比例的组分油使用量如表4所示，表5列出了基于此调合配方得到的产品汽油质量指标结果及最终利润。按所得配方验证后进行生产得到的产品汽油总利润为：$75744.7122。 

表4：一天配方结果 

表5：一天配方的质量指标和利润 

三天或者更多天汽油调合配方优化，需在一天优化基础之上，根据每一天汽油需求量以及组分油供应量作进一步计算即可，具体过程为： 

1）根据第一天两种汽油需求量以及组分油供应量进行第一天的汽油调合配方优化，计算结束后输出结果并将剩余组分油加入第二天的供应量； 

2）根据第二天两种汽油需求量以及组分油供应量计算第二天汽油调合配 方，计算完毕后输出结果并将剩余组分油加入第三天的供应量； 

3）根据第三天两种汽油需求量以及组分油供应量计算第三天组分油配方，计算完毕后输出结果并将剩余组分油加入下一天的供应量。以此类推可以得到多天配方结果。本实施例中得到的三天配方结果和质量指标分别列于表6和表7。三天的总利润为：212379.4430$，并且每一天质量指标都符合要求。 

表6：三天配方结果 

表7：三天配方的质量指标和利润 

根据配方中计算的各组分油的使用量购进组分油，即可以节省库存量，实现最小库存费用。出现以下情时可以按提示进行操作： 

（1）当要实现更长周期的调合调度时，方法相同，相应扩展约束条件和惩罚项即可。 

（2）当要实现更多产品油的调合，在目标函数中增加相应的利润目标相和惩罚项即可。 

（3）当组分油或供需变化时，更改相应惩罚项即可。