[发明专利]基于脉冲投放捕食‑被食模型的多目标组合调度优化方法

摘要

本发明公开的是一种基于脉冲投放捕食‑被食模型的多目标组合调度优化方法，假设在一个生态系统中种植了一种贵重药用植物，以该药用植物为食的害虫种群有若干种；这些害虫种群存在同一种天敌种群，天敌种群与害虫种群构成捕食‑被食关系；害虫种群数量的减少会导致天敌种群数量的减少；而天敌种群数量的减少会导致害虫种群数量的剧增；定期人工放养天敌种群会使天敌种群数量瞬间突增，突然增加的天敌种群会逐步抑制害虫种群数量的剧增；天敌种群数量的脉冲增加而导致的害虫种群数量的突然变化使得搜索能够快速跳出局部最优解陷阱；利用害虫种群的生长变化和脉冲投放捕食‑被食模型能够快速找到多目标组合调度优化问题的全局最优解决方案。

主权项

一种基于脉冲投放捕食‑被食模型的多目标组合调度优化方法，即MOSLO_CIPPSD方法，其特征在于：设要解决的多目标组合调度优化模型的一般形式如下：$\begin{array}{c}\min \{O_1{f}_1\left(X\right),O_2{f}_2\left(X\right),...,O_M{f}_M\left(X\right)\}\\ \begin{array}{cc}s.t.& \left\{\begin{array}{c}{g}_{i_a}\left(X\right)\&GreaterEqual;0,i_a\&Element;I\\ h_{i_b}\left(X\right)=0,i_b\&Element;E\\ X\&Element;H\&Subset;R^n\end{array}\right.\end{array}\end{array}---\left(1\right)$式中：(1)Rⁿ是n维欧氏空间，n为该优化模型所包含的变量总数；(2)X＝(x₁，x₂，…，x_m，x_m+1，…，x_n)是一个n维决策向量，其中，前m个变量x₁，x₂，…，x_m是连续实数型变量，用来表示模型中涉及到的资源型参数；后n‑m个变量x_m+1，…，x_n是0、1整数型变量，一个0、1整数型变量又称为一个要素，即对于任意x_j∈{x_m+1，…，x_n}，若x_j＝1，则表示第j个要素被选中为该最优组合调度中的一个要素，若x_j＝0，则表示第j个要素未被选中；(3)f₁(X)，f₂(X)，…，f_M(X)为M个目标函数，用来表示选择组合调度策略时的M个控制目标要求；(4)O₁，O₂，…，O_M为M个目标函数的优先级，优先级次序要求满足O₁>O₂>…>O_M，即目标函数f₁(X)首先要求达到最小，其次是f₂(X)，再其次是f₃(X)，依次类推，最后要求达到最小的是目标函数f_M(X)；(5)表示要素选择时所需满足的第i_a个不等式约束条件；I为不等式约束条件编号的集合；(6)表示要素选择时所需满足的第i_b个等式约束条件；E为等式约束条件编号的集合；(7){f_i(X)，i＝1，2，…，M}、的数学表达式没有限制条件；(8)H为搜索空间，又称解空间；(9)计算时，决策向量X也称为试探解；若试探解X不满足约束条件，则令f(X)＝+∞；将多目标组合调度优化模型式(1)转换成如下单目标组合调度优化模型：$\begin{array}{c}\min \{F\left(X\right)=\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{O}_k{f}_k\left(X\right)\}\\ \begin{array}{cc}s.t.& \left\{\begin{array}{c}{g}_{i_a}\left(X\right)\&GreaterEqual;0,i_a\&Element;I\\ h_{i_b}\left(X\right)=0,i_b\&Element;E\\ X\&Element;H\&Subset;R^n\end{array}\right.\end{array}\end{array}---\left(2\right)$式中，O_k＝10^M‑k；k为目标函数的编号；所述MOSLO_CIPPSD方法采用常数脉冲投放多种群捕食‑被食系统动力学理论，假设在一个生态系统中种植了一种贵重药用植物，以该药用植物为食的害虫种群有若干种；为了确保该贵重药用植物天然品质，人为依靠自然力量来控制害虫种群对该贵重药用植物的危害；这些害虫种群存在同一种天敌种群，该天敌种群以这些害虫种群为食，即天敌种群与害虫种群构成捕食‑被食关系；害虫种群和天敌种群的自我繁衍过程是一个连续的过程，害虫种群数量的减少会导致天敌种群数量的减少；由于害虫种群数很多，天敌种群数量的减少会导致害虫种群数量的剧增；为了有效控制害虫种群，必须定期人工放养天敌种群，使天敌种群数量瞬间突增，突然增加的天敌种群会逐步抑制害虫种群的数量剧增；害虫种群的生长变化相当于搜索空间的试探解从一个位置转移到另外一个位置；天敌种群数量的脉冲增加会导致害虫种群数量的突然变化，此相当于搜索空间的试探解从一个位置猛烈跳到另外一个位置，这种性质有利于使搜索跳出局部最优解陷阱；时期t，害虫种群P_i的生长能力强弱用种群生长指数PGI来表示，害虫种群P_i的PGI指数计算方法为：式中，X_i(t)为时期t害虫种群P_i所对应的试探解；N为害虫种群数；i表示害虫种群P_i的编号；所述MOSLO_CIPPSD方法包括如下步骤：(S1)初始化：a)令时期t＝0；按表1初始化本方法中涉及到的所有参数；b)随机确定N个害虫种群的初始密度y₁(0)，y₂(0)，…，y_N(0)；c)随机确定天敌种群的初始密度z(0)；d)随机确定N个试探解X₁(0)，X₂(0)，…，X_N(0)；表1 参数的取值方法(S2)执行下列操作：(S3)令时期t从0到G，循环执行下述步骤(S4)～步骤(S23)，其中G为演化时期数；(S4)计算：c^t＝Rand(c₀，c₁)，d^t＝Rand(d₀，d₁)，Q^t＝Rand(Q₀，Q₁)；式中：c^t，d^t，Q^t分别为参数b_i，η_i，α_i，c，d，Q在时期t的取值；b_i表示害虫种群P_i的增长率，b_i>0，b₀和b₁表示取值的下限和上限，且满足b₀>0，b₁>0，b₀≤b₁；η_i表示害虫种群P_i的减少率，η_i>0，η₀和η₁表示取值的下限和上限，且满足η₀>0，η₁>0，η₀≤η₁；α_i表示害虫种群之间的竞争参数，α_i>0，α₀和α₁表示取值的下限和上限，且满足α₀>0，α₁>0，α₀≤α₁；d表示天敌种群的转化率，d>0，d₀和d₁表示d^t取值的下限和上限，且满足d₀>0，d₁>0，d₀≤d₁；c表示天敌种群的增长率，c>0，c₀和c₁表示c^t取值的下限和上限，且满足c₀>0，c₁>0，c₀≤c₁；Q表示天敌种群的投放量，Q>0，Q₀和Q₁表示Q^t取值的下限和上限，且满足Q₀>0，Q₁>0，Q₀≤Q₁；Rand(A，B)表示在[A，B]区间产生一个均匀分布随机数，A和B为给定的常数，要求A≤B；(S5)按式(5)计算r_i(t)：$r_i\left(t\right)=\frac{y_i\left(t\right)}{\underset{s=1}{\overset{N}{\mathrm{\&Sigma;}}}{y}_s\left(t\right)},i=1,2,...,N---\left(5\right)$式中，r_i(t)为时期t害虫种群P_i在所有害虫种群中所占的比例，r_i(t)又称为害虫种群P_i的占比；y_i(t)，y_s(t)分别表示时期t害虫种群P_i和P_s的密度，y_i(t)≥0，y_s(t)≥0；(S6)令i从1到N，循环执行下述步骤(S7)～步骤(S20)；(S7)生成特征种群集合AS、PM、SM；时期t，对于当前害虫种群为P_i，特征种群集合AS、PM、SM的生成方法如下：a)产生高密度害虫种群集合AS：从N个害虫种群中随机挑选出L个害虫种群，其编号形成集合使得对于所有s∈{s₁，s₂，…，s_L}，满足r_s(t)＞r_i(t)；L又称为施加影响的害虫种群数；b)产生优势害虫种群集合PM：先从N个害虫种群中随机挑选出L个种群，这些害虫种群的PGI指数比当前害虫种群P_i的PGI指数高，形成集合其中g₁，g₂，…，g_L是这些害虫种群的编号；c)产生强势害虫种群集合SM：从N个害虫种群中随机挑选出L个种群，这些害虫种群的PGI指数和占比要比当前种群P_i的PGI指数和占比高，形成强势种群集合其中h₁，h₂，…，h_L是这些害虫种群的编号；即对于所有s∈{h₁，h₂，…，h_L}，有PGI(X_s(t))>PGI(X_i(t))，且占比r_s(t)>r_i(t)；(S8)若t不能被T整除，则按式(6)计算害虫种群P_i的密度y_i(t+1)和天敌种群的密度z(t+1)；否则，若t能被T整除，则按式(7)计算y_i(t+1)和z(t+1)；$\left\{\begin{array}{cc}{y}_i(t+1)=y_i\left(t\right)+y_i\left(t\right)(b_i^t-y_i(t)-\underset{s=1,s\&NotEqual;i}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_s^t{y}_s(t)-{\mathrm{\&eta;}}_i^{t}z(t\left)\right)& i=1,2,...,N\\ z(t+1)=z\left(t\right)+z\left(t\right)(-c^t+d^t\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&eta;}}_i^t{y}_i(t\left)\right)& \end{array}\right.---\left(6\right)$$\left\{\begin{array}{cc}{y}_i(t+1)=y_i\left(t\right)& i=1,2,...,N\\ z(t+1)=z\left(t\right)+Q^t& \end{array}\right.---\left(7\right)$式中，T表示天敌种群的投放周期；z(t)表示时期t天敌种群的密度，z(t)≥0；式(6)、式(7)来自于式(4)所描述的带有常数脉冲投放的多种群捕食‑被食系统动力学模型：$\left\{\begin{array}{cc}\left.\begin{array}{cc}\frac{{\mathrm{dy}}_i\left(t\right)}{dt}=y_i\left(t\right)(b_i-y_i(t)-\underset{s=1,s\&NotEqual;i}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_s{y}_s(t)-{\mathrm{\&eta;}}_{i}z(t\left)\right)& i=1,2,...,N\\ \frac{dz\left(t\right)}{dt}=z\left(t\right)(-c+d\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&eta;}}_i{y}_i(t\left)\right)& \end{array}\right\}& t\&NotEqual;kT\\ \left.\begin{array}{cc}{\mathrm{\&Delta;y}}_i\left(t\right)=0& i=1,2,...,N\\ \mathrm{\&Delta;}z\left(t\right)=Q& \end{array}\right\}& t=kT\end{array}\right.---\left(4\right)$(S9)令j从1到n，循环执行下述步骤(S10)～步骤(S18)；(S10)计算：p＝Rand(0,1)，其中p为害虫种群P_i被天敌种群捕食和害虫种群P_i与其它同类种群相互竞争时，其生长特征受到影响的实际概率；(S11)若p≤E₀，则执行步骤(S12)～(S16)，其中E₀为害虫种群因被天敌种群捕食和同类种群相互竞争时，其生长特征受到影响的最大概率；否则，转步骤(S17)；(S12)计算：q₀＝Rand(0,1)，其中q₀为食饵算子、天敌算子、优势算子、强势算子被执行的实际概率；(S13)若q₀≤1/4，则按式(8)执行食饵算子，得到v_i,j(t+1)；$v_{i,j}(t+1)=\left\{\begin{array}{cc}{r}_{i_1}\left(t\right)x_{i_1,j}\left(t\right)+r_{i_2}\left(t\right)x_{i_2,j}\left(t\right)-r_{i_3}\left(t\right)x_{i_3,j}\left(t\right)& j\&le;m\\ Great(AS,i,j)& j\&gt;m\end{array}\right.---\left(8\right)$式中：v_i,j(t+1)为时期t+1当前害虫种群P_i的特征j状态值；和分别为时期t害虫种群和的特征j的状态值；i₁，i₂，i₃是从{s₁，s₂，…，s_L}中随机选取出来的，且满足i₁≠i₂≠i₃；Great(W，j)的含义是：k₁，k₂，k₃是从集合W中随机选取的三个不同种群的编号，即满足k₁≠k₂≠k₃；(S14)若1/4<q₀≤1/2，则按式(9)执行天敌算子，得到v_i,j(t+1)；$v_{i,j}(t+1)=\left\{\begin{array}{cc}{r}_i\left(t\right)x_{i,j}\left(t\right)+\underset{s\&Element;AS}{\max }\left\{r_s\left(t\right)x_{s,j}\left(t\right)\right\}-\underset{s\&Element;AS}{\max }\left\{r_s\left(t\right)x_{s,j}\left(t\right)\right\}& j\&le;m\\ Great(AS,i,j)& j\&gt;m\end{array}\right.---\left(9\right)$(S15)若1/2<q₀≤3/4，则当j≤m时按式(10)执行优势算子，得到v_i,j(t+1)；当j>m时按式(11)执行优势算子，得到v_i,j(t+1)；$v_{i,j}(t+1)=\left\{\begin{array}{cc}{r}_{g_a}\left(t\right)x_{g_a,j}\left(t\right)+r_{g_b}\left(t\right)x_{g_b,j}\left(t\right)-r_{g_c}\left(t\right)x_{g_c,j}\left(t\right)& \left|PM\right|\&GreaterEqual;3\\ r_i\left(t\right)x_{i,j}\left(t\right)+r_{g_a}\left(t\right)x_{g_a,j}\left(t\right)-r_{g_b}\left(t\right)x_{g_b,j}\left(t\right)& \left|PM\right|=2\\ r_{g_a}\left(t\right)x_{g_a,j}\left(t\right)& \left|PM\right|=1\\ x_{i,j}\left(t\right)& \left|PM\right|=0\end{array}\right.---\left(10\right)$$v_{i,j}(t+1)=\left\{\begin{array}{cc}Great(PM,i,j)& \left|PM\right|\&GreaterEqual;1\\ x_{i,j}\left(t\right)& \left|PM\right|=0\end{array}\right.---\left(11\right)$式中，g_a、g_b、g_c在{g₁，g₂，…，g_L}中随机选择，要求g_a≠g_b≠g_c；(S16)若3/4<q₀≤1，则当j≤m时按式(12)强势算子，得到v_i,j(t+1)；当j>m时按式(13)执行强势算子，得到v_i,j(t+1)；$v_{i,j}(t+1)=\left\{\begin{array}{cc}{r}_{h_a}\left(t\right)x_{h_a,j}\left(t\right)+r_{h_b}\left(t\right)x_{h_b,j}\left(t\right)-r_{h_c}\left(t\right)x_{h_c,j}\left(t\right)& \left|SM\right|\&GreaterEqual;3\\ r_i\left(t\right)x_{i,j}\left(t\right)+r_{h_a}\left(t\right)x_{h_a,j}\left(t\right)-r_{h_b}\left(t\right)x_{h_b,j}\left(t\right)& \left|SM\right|=2\\ r_{h_a}\left(t\right)x_{h_a,j}\left(t\right)& \left|SM\right|=1\\ x_{i,j}\left(t\right)& \left|SM\right|=0\end{array}\right.---\left(12\right)$$v_{i,j}(t+1)=\left\{\begin{array}{cc}Great(SM,i,j)& \left|SM\right|\&GreaterEqual;1\\ x_{i,j}\left(t\right)& \left|SM\right|=0\end{array}\right.---\left(13\right)$式中，h_a、h_b、h_c在{h₁，h₂，…，h_L}中随机选择，要求h_a≠h_b≠h_c；(S17)若p>E₀，则令v_i,j(t+1)＝x_i,j(t)；(S18)令j＝j+1，若j≤n，则转步骤(S10)，否则转步骤(S19)；(S19)按式(14)执行生长算子，得到X_i(t+1)；式中：X_i(t)＝(x_i,1(t),x_i,2(t),…,x_i,n(t))；V_i(t+1)＝(v_i,1(t+1),v_i,2(t+1),…,v_i,n(t+1))；(S20)令i＝i+1，若i≤N，则转步骤(S7)，否则转步骤(S21)；(S21)若新得到的全局最优解X^*t+1与最近一次获得的全局最优解之间的误差满足最低要求ε，则转步骤(S24)，否则转步骤(S22)；(S22)保存新得到的全局最优解X^*t+1；(S23)令t＝t+1，若t≤G，则转上述步骤(S4)，否则转步骤(S24)；(S24)结束。