CN106875001A - 具有冯·诺依曼结构的随机漂移粒子群优化方法 - Google Patents

Abstract

本发明公开了一种具有冯·诺依曼拓扑结构的随机漂移粒子群优化方法。针对生物网络如合成基因回路建模，需要建立生物系统的微分方程模型并根据观测值推理模型参数。参数计算过程首先需要定义合适的目标函数和约束条件，同时需要得到优化算法的支持。本发明将微分方程的参数估计等价为带有微分‑代数约束的非线性规划问题，通过寻找合适参数集来拟合量测数据。针对参数估计过程的寻优计算易陷入局部最优的局限性，该发明在随机漂移粒子群优化算法中引入冯·诺依曼结构，以增强全局搜索能力。该改进算法针对全连接的全局拓扑结构的局限性进行改进，采用局部拓扑结构作为粒子间的信息共享方式，增强算法全局搜索能力，能够给出可信度高的模型参数。

Description

具有冯·诺依曼结构的随机漂移粒子群优化方法

技术领域

本发明涉及生物网络辨识与参数估计领域，尤其涉及一种具有冯·诺依曼结构的随机漂移粒子群优化方法。

背景技术

生物技术的快速发展为食品医药、环境治理、新能源等领域带来持续动力，各国政府对以合成生物学等为代表的生物技术关注度和投资力度逐步增加。合成生物学的目的在于建立人工生物系统，并让它们按照特定规律发挥作用，从基因片段、DNA分子、基因调控网络到细胞的人工设计与合成。波士顿大学的柯林斯等人开发一种“套环开关(ToggleSwitch)”的装置，所选择的细胞功能可随意开关。加州大学的埃洛维茨等人开发的合成基因振荡回路，当某种特殊蛋白质含量发生变化时，细胞能在发光和非发光状态之间切换，起到有机振荡器的作用。合成生物装置或系统的合理有效设计有赖于有效的数学模型，因此可靠的数据源和参数推理方法将是生物系统建模和设计中的有效工具。

但是由于生物系统的复杂程度及问题求解规模，用于参数估计的推理算法能力之间的矛盾日益加剧，因此迫切需要性能经过提升的推理算法来解决具有较多变量的生物系统，如基因调控网络、代谢网络等。对生物网络进行建模既为探索未知生物对象提供有力工具，同时能够辅助合成生物系统的构建，在合成生物学建模领域，相当部分的工作围绕着工程化基因回路建模展开。通常生物系统含有数百个甚至上千个变量，比如DREAM 5平台公布的关于E coli的微阵列表达数据，涉及4511个基因。但是当前的参数估计方法能够应对的系统一般不超过十个未知参数，在未知参数个数增加时求解难度和计算所需资源呈现快速上升趋势，其中用于参数估计的最优化方法的计算能力属于制约因素之一。除计算速度外，估计参数的可靠性也是网络推理过程必须考虑的问题，所得到的模型参数应最大程度上接近真实值，这就要求寻优算法具有较好的全局搜索能力。考虑到参数估计过程涉及到目标函数优化计算，开发高效、全局搜索能力强的优化算法在一定程度能够缓解这种矛盾。

粒子群算法(Particle Swarm Optimization,PSO)及其改进方法凭借其寻优能力和适应性在优化计算领域持续受到关注和应用，作为一种基于种群的优化算法，粒子群优化算法通过迭代最小化候选解的目标函数，这些候选解在算法中用粒子等价表示。粒子在解空间移动的过程既受到自身经验的影响，同时也受到当前种群最优位置的影响，在基本PSO方法的速度更新公式中，每个粒子的位置都会根据个体历史最优位置和种群全局历史最优位置不断进行调整，从而收敛到最优解。基本的PSO算法能够有效处理低维度的优化问题，但在处理高维度的复杂优化问题时易陷入局部最优。由于通过改变粒子速度的更新模式能够在一定程度提高其搜索能力，Jun等人提出随机飘移粒子群(RDPSO)算法，模拟自由电子向具有最小势能位置的运动，使种群中粒子不断调整自身位置。本发明在此基础上，通过对算法的拓扑结构进行调整，将全连接方式用部分拓扑连接即冯·诺依曼结构代替。

发明内容

为了提高生物网络参数估计的可靠性，解决现有粒子群寻优算法所面临的局部最优问题，本发明提供一种具有冯·诺依曼结构的随机漂移粒子群优化方法，能够作为参数推理算法较好处理基因网络等生物系统建模问题，在参数寻优过程提升了全局搜索能力。

本发明解决其技术问题采用的技术方案是：一种具有冯·诺依曼结构的随机漂移粒子群优化方法，该方法包括以下步骤：

(1)初始化粒子群，包括粒子群规模m和待解决问题的维度D，粒子群中每个粒子都包含3个矢量表征自身特性，即当前位置X_i、当前速度矢量V_i和个体历史最优位置P_besti，为所有粒子随机初始化当前位置矢量i表示种群大小为m的粒子群中第i个粒子，i＝1,2,…,m，设置步数k＝0；

(2)根据目标函数f(X_i(k))来计算各粒子的初始代价函数值，令初始位置为个体历史最优位置P_besti(k)，并计算种群的初始全局历史最优位置G_gbest(k)和冯·诺依曼邻域中所有粒子的平均历史最优位置C_mbesti(k)；

(3)判断循环终止条件，当代价函数值误差满足设定的精确度，或者迭代次数达到最大值，停止搜索并输出得到的最优解，否则执行步骤(4)；

(4)设置迭代步数k＝k+1，更新种群中各粒子位置；

(5)更新粒子i的个体历史最优位置P_besti(k)，重新计算种群中各粒子的代价函数值f(X_i(k))，若当前f(X_i(k))小于代价函数f(P_besti(k-1))，则将f(X_i(k))置为粒子i的历史最优位置，即P_besti(k)＝X_i(k)，否则P_besti(k)＝P_besti(k-1)；在粒子群算法中P_besti(k)和G_gbest(k-1)的更新方程如下所示：

(6)更新种群全局历史最优位置，将每个粒子对应代价函数值f(X_i(k))与当前全局历史最优位置的代价函数值f(G_gbest(k-1))相比较，若满足条件f(X(k))＜f(G_gbest(k-1))，则最优位置更新G_gbest(k)＝X(k)；否则G_gbest(k)＝G_gbest(k-1)；

(7)计算每个粒子的冯·诺依曼邻域中所有粒子的平均历史最优位置C_mbesti(k)；

(8)返回步骤(3)。

进一步地，在第k次迭代时，每个粒子i都有热运动和向局部吸引因子p_i(k)的定向运动，这两个运动的速率在第d维度上分别表示为和粒子i在第d维度上速度表示为：

假设服从双指数分布，利用随机模拟方法将表示为：

其中，

其中，为分布的标准差，s和是在区间(0,1)上服从均匀分布的两个不同的随机数。

进一步地，漂移运动的速度采取线性形式，即当时，随着k趋于无穷，将趋于因此保证粒子在整体上向运动；β为漂移系数，反应粒子向定向运动的能力。

进一步地，粒子的速度和当前位置更新采用如下形式：

其中，α为热系数，其大小体现了算法的全局搜索能力，其值越大则全局搜索能力越强。热运动部分的表示冯·诺依曼邻域中所有粒子个体历史最优位置平均值。

进一步地，热系数α设为S型函数，表达式如下：

α＝0.3+0.6/(1+(4*k/n)³)

其中，k为迭代次数，n为最大迭代次数。

本发明的有益效果是：(1)粒子间信息共享方式采用高效冯·诺依曼结构形式；(2)提高多极值寻优问题和生物网络参数辨识中算法的全局寻优能力。若粒子与其他粒子采用全连接结构，将限制算法全局搜索能力，采用局部拓扑结构进行粒子间信息交互。为了使具有冯·诺依曼邻域的RDPSO算法中粒子向全局最优值运动，考虑从和分布形式和计算入手，本发明根据迭代次数的增加相应改变热系数α，在开始阶段进行较强的全局搜索，随着计算进行，全局搜索能力逐步下调而局部搜索能力提升。

附图说明

图1为具有冯·诺依曼结构的RDPSO方法流程图；

图2为环形拓扑结构图；

图3为星形拓扑结构图；

图4为冯·诺依曼拓扑结构图；

图5为热系数α与迭代次数关系；

图6为PSO、RDPSO和Von-RDPSO三种方法对Griewank函数的寻优收敛曲线；

图7为PSO、RDPSO和Von-RDPSO三种方法对Ackley函数的寻优收敛曲线

图8为参数推理中收敛速度对比；

图9为输出y₃结果对比；

图10为输出y₄结果对比；

图11为输出y₅结果对比；

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细说明。

如图1所示，本发明提供的一种具有冯·诺依曼结构的随机漂移粒子群优化方法，该方法包括以下步骤：

(1)初始化粒子群，包括粒子群规模m和待解决问题的维度D，粒子群中每个粒子都包含3个矢量表征自身特性，即当前位置X_i、当前速度矢量V_i和个体历史最优位置P_besti，为所有粒子随机初始化当前位置矢量i表示种群大小为m的粒子群中第i个粒子，i＝1,2,…,m，设置步数k＝0；

(2)根据目标函数f(X_i(k))来计算各粒子的初始代价函数值，令初始位置为个体历史最优位置P_besti(k)，并计算种群的初始全局历史最优位置G_gbest(k)和冯·诺依曼邻域中所有粒子的平均历史最优位置C_mbesti(k)；

(3)判断循环终止条件，当代价函数值误差满足设定的精确度，或者迭代次数达到最大值，停止搜索并输出得到的最优解，否则执行步骤(4)；

(4)设置迭代步数k＝k+1，更新种群中各粒子位置；

(5)更新粒子i的个体历史最优位置P_besti(k)，重新计算种群中各粒子的代价函数值f(X_i(k))，若当前f(X_i(k))小于代价函数f(P_besti(k-1))，则将f(X_i(k))置为粒子i的历史最优位置，即P_besti(k)＝X_i(k)，否则P_besti(k)＝P_besti(k-1)；在粒子群算法中P_besti(k)和G_gbest(k-1)的更新方程如下所示：

(6)更新种群全局历史最优位置，将每个粒子对应代价函数值f(X_i(k))与当前全局历史最优位置的代价函数值f(G_gbest(k-1))相比较，若满足条件f(X(k))＜f(G_gbest(k-1))，则最优位置更新G_gbest(k)＝X(k)；否则G_gbest(k)＝G_gbest(k-1)；

(7)计算每个粒子的冯·诺依曼邻域中所有粒子的平均历史最优位置C_mbesti(k)；

(8)返回步骤(3)。

进一步地，在第k次迭代时，每个粒子i都有热运动和向局部吸引因子p_i(k)的定向运动，这两个运动的速率在第d维度上分别表示为和粒子i在第d维度上速度表示为：

假设服从双指数分布，利用随机模拟方法将表示为：

其中，

其中，为分布的标准差，s和是在区间(0,1)上服从均匀分布的两个不同的随机数。

进一步地，漂移运动的速度采取线性形式，即当时，随着k趋于无穷，将趋于因此保证粒子在整体上向运动；β为漂移系数，反应粒子向定向运动的能力。

进一步地，粒子的速度和当前位置更新采用如下形式：

其中，α为热系数，其大小体现了算法的全局搜索能力，其值越大则全局搜索能力越强。热运动部分的表示冯·诺依曼邻域中所有粒子个体历史最优位置平均值。

进一步地，热系数α设为S型函数，表达式如下：

α＝0.3+0.6/(1+(4*k/n)³)

其中，k为迭代次数，n为最大迭代次数。

实施例

本实施例提供一种具有冯·诺依曼结构的随机漂移粒子群优化方法，包括下述步骤：

步骤1，参数估计问题在数学上等价为一个带有微分-代数约束条件的非线性规划问题，其目标是找到未知参数集θ使目标函数J取得最小值，寻优过程需要满足给定之约束条件；

步骤2，利用优化算法对带有未知参数的系统进行参数辨识时，实际上通过优化算法最小化目标函数，同时不断搜索、更新参数的值。每一组估计的模型参数都会被代入到描述生物系统的微分方程组中，然后采用四阶的龙格-库塔方法求解该微分方程组，得到一组模型的预测值y_pred。利用优化算法不断搜索新的可行参数，选择目标函数值最小时的模型参数集。目标函数值越小，表示得到的参数估计值和参数的标准值之间匹配度越高；

步骤3，算法评价指标，采用算法运行时间、算法收敛速度、迭代结束时的目标函数值以及每次实验的最佳参数估计值与标准值的比较和最佳参数估计值下对象的状态变量的名义输出值和最佳预测值的比较作为评价算法的主要指标。

本实施例涉及对象采用微分方程组进行描述，该对象表示具有反馈通路的代谢通路，包含5个基因，对本发明方法的求解效率进行验证。该对象可用非线性微分方程描述

其中p_i(i＝1,…,5)表示待估计参数，物理意义为基因之间的动态转化常数；而y_i(i＝1,…,5)表示5个基因的浓度，该代谢通路模型有关参数如下表所示：

参数	取值
			5
	10
			10
	10
			10

针对此非线性对象，选择的目标函数为

需要满足的约束条件可描述为

x(t₀)＝x₀

h(x,y,θ,v)＝0

g(x,y,θ,v)≤0

其中θ表示待估计的模型参数矢量，y_msd(t)和y(θ,t)分别表示输出状态变量的测量值和模型预测值，而函数h和g则分别表示对系统性能的额外要求。为便于寻优计算，代价函数采用离散化形式，即实验测量值和模型预测值差值的平方和

其中l表示实验中单个状态的采样数据点数量，n表示状态变量数目，y_exp和y_pred分别为表示系统状态变量的实验值矢量和模型预测矢量。将优化算法得到的估计参数代入对应对象的微分方程组中，求解微分方程组时，本发明使用四阶龙格-库塔方法。附图2是环形拓扑结构图。m为算法中粒子的数目，实线为粒子间直接相连接，虚线表示中间有省略的粒子。

附图3是星形拓扑结构：其中一个粒子作为整个拓扑结构的中心，与其他所有的粒子相连，但其他粒子间不连接。此结构加快了粒子间的信息传播速度，使得算法的收敛速度较快，多用于多模函数优化，但也较易陷于局部最优。

附图4是冯·诺依曼拓扑结构图。这是由许多方形组成的立体结构，每个个体都与上下左右的个体邻居相连接。附图3左侧是冯诺依曼的立体结构，是一个四方网格，顶点相连形成的环面。附图3右侧是冯诺依曼结构展开的形状，是一个由上下左右相邻粒子组成的二维矩阵形式。

附图5是PSO、RDPSO和Von-RDPSO三种方法对Griewank函数的寻优收敛曲线,其中“Von-RDPSO”表示冯·诺依曼拓扑结构的RDPSO算法，横坐标为算法迭代次数，纵坐标为目标函数值，即当前函数值与函数的全局最小值之间的误差

附图6是PSO、RDPSO和Von-RDPSO三种方法对Ackley函数的寻优收敛曲线，其中“Von-RDPSO”表示冯·诺依曼拓扑结构的RDPSO算法。PSO算法在对两个函数进行寻优时，其收敛速度均最快，RDPSO算法次之，冯·诺依曼拓扑结构的RDPSO算法则最慢，这也满足设计这些算法时的初衷，即增强算法的全局搜索能力，但是在得到的寻优结果方面，RDPSO算法并没有取得很好的效果，即它有可能陷入局部最优，而冯·诺依曼拓扑结构的RDPSO算法在此方面表现较好。在三种算法中，冯·诺依曼拓扑结构的RDPSO算法对两个函数的寻优结果最好，由此也可以看出冯·诺依曼拓扑结构的RDPSO算法较好的性能。

附图7是参数推理中收敛速度对比，横坐标表示迭代次数，纵坐标代表目标函数值。alpha0表示热系数α为定值，alpha1采用线性函数，表达式为α₁＝0.9-0.6*(k/n)，其中k为迭代次数，n为算法中设置的最大迭代次数；alpha2采用S型函数表示，表达式为α₃＝0.3+0.6/(1+(4*k/n)³)。在拓扑结构相同的情况下，热系数α为alpha1时收敛速度最快，热系数α为alpha0和alpha2时全局搜索能力基本接近，前期alpha0收敛较快，后期alpha2对应曲线快速下降。在热系数相同的情况(均设为alpha2)下，具有全连接拓扑结构的算法收敛速度最快，冯·诺依曼拓扑结构下次之，环形拓扑结构下收敛最慢，但全连接结构的算法较不稳定。

附图8为输出y₃结果对比；附图9为输出y₄结果对比；附图10为输出y₅结果对比，对于输出y₃，y₄和y₅，存在四条曲线，从图中可以看出：在全连接拓扑结构下，无论热系数α为alpha0，alpha1或alpha2时，得到的输出结果与标准输出均出现了较大偏差，即未找到满足条件的参数值，而环形拓扑结构和冯·诺依曼拓扑结构下的输出结果与标准输出几乎相同。全连接的拓扑结构下，取得的参数p₃和p₄与其标准值存在较大差异，而输出y₃和y₄分别与参数p₃和p₄之间存在着直接联系，而y₅虽然未被其直接影响。输出y₄和y₅之间存在相互影响，即可能是y₄对y₅的间接影响导致y₅的估计值在这三种算法下与标准值产生了较大偏差。

除使用冯·诺依曼拓扑结构的RDPSO算法外，还将使用以下四种算法与其进行对比，分别是：①原始RDPSO算法，②将原始RDPSO算法中的热系数α设为线性函数得到的算法，③将原始RDPSO算法中的热系数α设为S型函数得到的算法，④将原始RDPSO算法中的热系数α设为S型函数，同时拓扑结构采用环形拓扑结构得到的算法。可以看出，这五种算法中，②和③是保持原始RDPSO算法中的全连接的拓扑结构不变，④和冯·诺依曼拓扑结构的RDPSO算法是采用③中的热系数，改变算法中的拓扑结构。其中S型函数所采用的线性函数的表达式如下所示

α₁＝0.9-0.6*(k/n)

为说明具有冯·诺依曼拓扑结构的RDPSO算法的有效性，将其与以下四种方法进行对比，这四种算法分别是：①原始RDPSO算法，②将原始RDPSO算法中的热系数α设为线性函数得到的算法，③将原始RDPSO算法中的热系数α设为S型函数得到的算法，④将原始RDPSO算法中的热系数α设为S型函数，同时拓扑结构采用环形拓扑结构得到的算法。将五种算法分别对基因回路进行参数估计实验，运行20次取得的最优目标函数值进行汇总得到表1，

表1 20次最优目标函数值分析表

评价指标	alpha0(Global)	alpha1(Global)	alpha2(Global)	alpha2(Ring)	alpha2(Von)
						平均值	31.98	24.87	11.83	6.81	8.06
均方差	25.62	26.39	18.56	15.39	15.33
						最好值	0.21	0.26	0.0064	0.0065	0.15
最差值	53.91	55.72	52.62	52.62	54.01

拓扑结构相同时，热系数α为alpha2时得到的最优目标函数值最小，由此可以看出对于此对象，α做S型变化较适合，热系数α相同时，环形拓扑结构下得到的最优目标函数值最小，由此可以看出环形拓扑结构下算法的全局搜索能力较强，冯·诺依曼拓扑结构的RDPSO算法次之，两者均比全连接拓扑结构下得到的最优目标函数值小；从均方差来看，冯·诺依曼拓扑结构的RDPSO算法最小，环形拓扑结构下热系数α为alpha2时次之，由此可知，跟其他情况相比，冯·诺依曼拓扑结构的RDPSO算法的稳定性较好。

与此同时，表2给出五种算法20次运行试验取得的最佳参数估计值与真值的对比情况

表2 20次实验最佳参数估计值分析表

参数	alpha0(Global)	alpha1(Global)	alpha2(Global)	alpha2(Ring)	alpha2(Von)	标准值
							p1	5.121	5.16	4.96	5.02	5.08	5.00
p2	10.06	10.11	9.99	10.01	10.06	10.00
							p3	9.58	9.50	9.60	9.94	9.73	10.00
p4	9.50	9.37	9.60	10.01	10.00	10.00
							p5	9.98	10.04	9.99	9.97	9.99	10.00

表2中的最佳参数估计值与参数标准值的比较中可以看到，环形拓扑结构和冯·诺依曼拓扑结构下热系数为alpha2时得到的最佳参数估计值与标准值较接近，体现出算法的高精度参数估计能力；另外三种情况下得到的最佳估计参数值与标准值相差大一些，尤其是参数p3和p4与其标准值相差较大。

Claims (5)

1.一种具有冯·诺依曼结构的随机漂移粒子群优化方法，其特征在于，该方法包括以下步骤：

(1)初始化粒子群，包括粒子群规模m和待解决问题的维度D，粒子群中每个粒子都包含3个矢量表征自身特性，即当前位置X_i、当前速度矢量V_i和个体历史最优位置P_besti，为所有粒子随机初始化当前位置矢量i表示种群大小为m的粒子群中第i个粒子，i＝1,2,…,m，设置步数k＝0；

(2)根据目标函数f(X_i(k))来计算各粒子的初始代价函数值，令初始位置为个体历史最优位置p_besti(k)，并计算种群的初始全局历史最优位置G_gbest(K)和冯·诺依曼邻域中所有粒子的平均历史最优位置C_mbesti(k)；

(3)判断循环终止条件，当代价函数值误差满足设定的精确度，或者迭代次数达到最大值，停止搜索并输出得到的最优解，否则执行步骤(4)；

(4)设置迭代步数k＝k+1，更新种群中各粒子位置；

(5)更新粒子i的个体历史最优位置P_besti(k)，重新计算种群中各粒子的代价函数值f(X_i(k))，若当前f(X_i(k))小于代价函数f(P_besti(k-1))，则将f(X_i(k))置为粒子i的历史最优位置，即P_besti(k)＝X_i(k)，否则P_besti(k)＝P_besti(k-1)；在粒子群算法中P_besti(k)和G_gbest(k-1)的更新方程如下所示：

$P_{besti}\left(k\right)=\left\{\begin{array}{c}{X}_i\left(k\right)iff\left(X_i\left(k\right)\right)<f\left(P_{besti}\left(k-1\right)\right)\\ P_{besti}(k-1)iff\left(X_i\left(k\right)\right)\&GreaterEqual;f\left(P_{besti}\left(k-1\right)\right)\end{array}\right.$

$G_{gbest}\left(k\right)=P_{bestg}\left(k\right),g=\underset{0\&le;i\&le;m}{\arg min}\&lsqb;f\left(P_{besti}\right(k\left)\right)\&rsqb;$

(6)更新种群全局历史最优位置，将每个粒子对应代价函数值f(X_i(k))与当前全局历史最优位置的代价函数值f(G_gbest(k-1))相比较，若满足条件f(X(k))＜f(G_gbest(k-1))，则最优位置更新G_gbest(k)＝X(k)；否则G_gbest(k)＝G_gbest(k-1)；

(7)计算每个粒子的冯·诺依曼邻域中所有粒子的平均历史最优位置C_mbesti(k)；

(8)返回步骤(3)。

2.根据权利要求1所述的具有冯·诺依曼结构的随机漂移粒子群优化方法，其特征在于，在第k次迭代时，每个粒子i都有热运动和向局部吸引因子p_i(k)的定向运动，这两个运动的速率在第d维度上分别表示为和粒子i在第d维度上速度表示为：

$V_i^d\left(k\right)=V{1}_i^d\left(k\right)+V{2}_i^d\left(k\right)$

假设服从双指数分布，利用随机模拟方法将表示为：

其中，

其中，为分布的标准差，s和是在区间(0,1)上服从均匀分布的两个不同的随机数。

3.根据权利要求2所述的具有冯·诺依曼结构的随机漂移粒子群优化方法，其特征在于，漂移运动的速度采取线性形式，即当 时，随着k趋于无穷，将趋于因此保证粒子在整体上向运动；β为漂移系数，反应粒子向定向运动的能力。

4.根据权利要求3所述的具有冯·诺依曼结构的随机漂移粒子群优化方法，其特征在于，粒子的速度和当前位置更新采用如下形式：

$X_i^d(k+1)=X_i^d\left(k\right)+V_i^d(k+1)$

其中，α为热系数，其大小体现了算法的全局搜索能力，其值越大则全局搜索能力越强。热运动部分的表示冯·诺依曼邻域中所有粒子个体历史最优位置平均值。

5.根据权利要求4所述的具有冯·诺依曼结构的随机漂移粒子群优化方法，其特征在于，热系数α设为S型函数，表达式如下：

α＝0.3+0.6/(1+(4*k/n)³)

其中，k为迭代次数，n为最大迭代次数。