CN110084355B - 大量相互作用粒子运动仿真系统的网格尺度寻优方法 - Google Patents

Abstract

本发明公开了大量相互作用粒子运动仿真系统的网格尺度寻优方法，包括如下步骤：S1、初始化场景：空间维度，边界长度；S2、初始化粒子：直径1b，质量，作用力截止距离，作用力方程，运动方程；S3、系统参数测算：A1、设置系统粒子数；A2、设定粒子占有空间的边长分别为2b、3b、4b、5b、6b五个等级；A3、对于每个密度等级，数值实验重复三次以上，寻找密度等级对应的最优格子尺寸；A4、基于最小二乘法估计回归方程参数，本发明基于理论分析，提供了一个计算方法，使得我们可以根据粒子的密度、“截止距离”以及空间维度来直接计算最优网格尺度，从而可以将获得的最优网格尺度应用于更大规模的粒子运动模拟系统中，节省计算时间。

Description

大量相互作用粒子运动仿真系统的网格尺度寻优方法

技术领域

本发明涉及网格尺度优化技术领域，具体为大量相互作用粒子运动仿真系统的网格尺度寻优方法。

背景技术

在大量粒子运动仿真系统中，邻域搜索(NS)是该系统中最基础的部分，用于粒子搜索其周围的邻近粒子，使得能够进而计算粒子之间的相互作用力，而原胞表(Cell-List)方法是加速NS进程的重要方法之一，该方法使得NS算法复杂度从最初的O(N²)变为O(N)。

有部分研究结果表明，在原胞列表(Cell-List)方法中如果网格尺度从“截止距离”(Cutoff distance)修改到某些更小的值时，系统运行速度会进一步加快，但是这些结果是基于经验的，缺乏理论分析的支撑，人们只能通过数值试算的方法来寻找到一个较优的网格尺度。

发明内容

本发明提供大量相互作用粒子运动仿真系统的网格尺度寻优方法，可以有效解决上述背景技术中提出在原胞列表(Cell-List)方法中如果网格尺度从“截止距离”(Cutoffdistance)修改到某些更小的值时，系统运行速度会进一步加快，但是这些结果是基于经验的，缺乏理论分析的支撑，人们只能通过数值试算的方法来寻找到一个较优的网格尺度的问题。

为实现上述目的，本发明提供如下技术方案：大量相互作用粒子运动仿真系统的网格尺度寻优方法，包括如下步骤：

S1、初始化场景：空间维度，边界长度；

S2、初始化粒子：直径1b，质量，作用力截止距离，作用力方程，运动方程；

S3、系统参数测算：

A1、设置系统粒子数；

A2、设定粒子占有空间的边长分别为2b、3b、4b、5b、6b五个等级；

A3、对于每个密度等级，数值实验重复三次以上，寻找密度等级对应的最优格子尺寸；

A4、基于最小二乘法估计回归方程参数。

根据上述技术方案，所述步骤S1中空间维度包括场景模型和网络模型；

所述场景模型场景是三维空间，设定其各维度相应的边界长度，从而确定场景的空间范围；

设置场景为R^D:∏L_S ^d；

其中1≤d≤D是正整数；

D表示空间维度；

L_S ^d表示场景R^D的第d 个维度的边界长度，即场景为一个立方体；

所述网格模型以场景立方体的某个顶点做为坐标原点，以长度L_C为间距，在各个维度上进行均分，得到网格化的空间，在这个网格化的空间中，每一个小立方体的边长为L_C的立方体称之为原胞，在场景的网格空间中，某个维度方向上，原胞的数量为

其中符号

表示向上取整。

根据上述技术方案，所述步骤S2中初始化粒子通过邻域搜索进行轨迹搜索；

设系统粒子数为N，为了计算这些粒子的运动轨迹，需要获得标记粒子，指当前被分配给计算资源的粒子，周围的其他粒子的信息，从而计算标记粒子和周围粒子之间的相互作用力，进而基于运动方程计算粒子运动轨迹；

寻找到与标记粒子距离小于截止距离的其他粒子的计算行为称之为邻域搜索；

用公式表示邻域NA和邻居NP，

NA_i(t)＝{x|||x-x_i(t)||＜R_C},x∈R^D (1)；

NP_i(t)＝{j|||x_i(t)-x_j(t)||＜R_C,j≠i} (2)；

其中，符号NA_i(t)表示编号为i的标记粒子在t时刻的球体邻域；

半径为截止距离R_C；

x_i(t)表示i粒子的空间坐标位置；

i是粒子编号且为正整数1≤i≤N；

符号NP_i(t)表示由球体邻域NA_i(t)内的其他粒子的编号所构成的集合。

根据上述技术方案，所述步骤S3中A1粒子通过元胞列表法搜索；

所述原胞列表法为了寻找到粒子i的邻居NP_i(t)，通常需要遍历系统中其它(N-1)个粒子，计算它们之间的距离，并使用(2)式作为判别标准，获得邻居NP_i(t)，这种方法的复杂度为O(N²)，其基本方法是将空间用边长为L_C的原胞网格化，原胞存储有占据其中的粒子编号的列表，标记粒子只需要在其周围的摩尔邻域内的原胞中寻找其原胞列表，将这些列表所对应的粒子构成的集合NP′_i(t)做为遍历对象，进而基于(2)式计算出NA_i(t)，在原胞列表法中，原胞边长设置为L_C，可以保证不会有比截止距离小的相邻粒子被漏掉，但是摩尔邻域显然是大于邻域NA_i(t)的，这将导致有些大于截止距离的相邻粒子进入到遍历集合NP′_i(t)中。

根据上述技术方案，所述摩尔邻域由粒子集合NP′_i(t)尽量接近球体邻域内的粒子集合NP_i(t)，选择更小的原胞边长，来构造网格空间；

所述摩尔邻域用数学公式来描述扩展邻域的定义如下：

其中，符号

表示向下取整，

表示向上取整；

d 表示D维空间中的第d 个维度，L_C表示原胞的边长。

根据上述技术方案，所述扩展邻域的期望原胞数和期望粒子数为设定比率l＝L_C/R_C为原胞边长和截止距离的比值，令比率l＝n+m，其中整数

小数m＝R_C/L_C-n，0≤m＜1；

设粒子在空间任意一个位置出现概率相等，则扩展邻域NE_i(t)内所包含原胞数量LC_i(t)的期望值可以表示为，

E(LC_i(t))＝P₁·(2(n+1))+P₂·(2n+1),m∈[0,0.5)；

E(LC_i(t))＝P′₁·(2(n+1))+P′₂·(2n+3),m∈[0.5,1)；

P₁＝2m,P₂＝1-2m,P′₁＝2(1-m),P′₂＝2m-1,

E(LC_i(t))＝2l+1 (4)；

进一步设空间中粒子均匀分布，即密度ρ为空间位置x∈R^D的常数，则可知道扩展邻域NE_i(t)内所包含的粒子数量LP′_i(t)，

E(LP′_i(t))＝ρL_C ^D(2l+1)^D＝ρR_C ^D(2+l^-1)^D (5)。

根据上述技术方案，所述粒子进行一次遍历并更新位置的运算中，邻域搜索运算和计算距离和作用力是单个粒子的基本运算；

设循环模块对应时间参数分别为σ₁和σ₂；

则一个粒子对应的循环模块运算量为(σ₁LC_i(t)+σ₂LP′_i(t)+τ_c)，其中τ_c是一个时间常量，由式(4)和(5)得到循环模块运算量的期望值为，

T_CE(L_C)＝σ₁(2l+1)^D+σ₂ρR_C ^D(2+l^-1)^D+τ_c (6)；

式(6)是一个关于原胞边长L_C的函数，我们可以进一步寻找满足下式时，

最优原胞边长

式子(7)即为我们的优化目标函数；

对函数(6)求导数，

可以得到：

或：

为待求的最优原胞边长，即最优网格尺度；

由方程式(9)，可以得到，

进而可以将(10)式化为一元线性回归方程形式，

其中，δ是回归方程的参数。

与现有技术相比，本发明的有益效果：本发明基于理论分析，提供了一个计算方法，使得我们可以根据粒子的密度、“截止距离”以及空间维度来直接计算最优网格尺度，从而可以将获得的最优网格尺度应用于更大规模的粒子运动模拟系统中，节省计算时间。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：

图1是本发明的步骤示意图；

图2是本发明的摩尔邻域示意图；

图3是本发明的扩展邻域示意图；

图4是本发明的包含五个原胞的一维扩展邻域示意图；

图5是本发明的二维扩展邻域示意图；

图6是本发明的m∈[0,0.5)时，蓝色原胞的三个区域划分示意图；

图7是本发明的m∈[0.5,1)时，蓝色原胞的三个区域结构示意图；

图8是本发明的通道行人粒子交通仿真示意图；

图9是本发明的截至距离为1.2m条件下的回归方程示意图；

图10是本发明的扩展邻域中的原胞个数随LC减小而增加示意图；

图11是本发明的与数组L^a对应的原胞边长LC的150个值示意图；

图12是本发明的扩展邻域中的邻居个数示意图；

图13是本发明的扩展邻域的面积示意图；

图14是本发明的单个粒子随LC减小的平均时间耗费示意图；

图15是本发明的单个粒子随LC减小的平均总时间耗费(T2+T3) 的面积示意图。

具体实施方式

以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

实施例：如图1所示，本发明提供技术方案，大量相互作用粒子运动仿真系统的网格尺度寻优方法，包括如下步骤：

S1、初始化场景：空间维度，边界长度；

S2、初始化粒子：直径1b，质量，作用力截止距离，作用力方程，运动方程；

S3、系统参数测算：

A1、设置系统粒子数；

A2、设定粒子占有空间的边长分别为2b、3b、4b、5b、6b五个等级；

A3、对于每个密度等级，数值实验重复三次以上，寻找密度等级对应的最优格子尺寸；

A4、基于最小二乘法估计回归方程参数。

根据上述技术方案，步骤S1中空间维度包括场景模型和网络模型；

场景模型场景是三维空间，设定其各维度相应的边界长度，从而确定场景的空间范围；

设置场景为R^D:∏L_S ^d；

其中1≤d≤D是正整数；

D表示空间维度；

L_S ^d表示场景R^D的第d 个维度的边界长度，即场景为一个立方体；

网格模型以场景立方体的某个顶点做为坐标原点，以长度L_C为间距，在各个维度上进行均分，得到网格化的空间，在这个网格化的空间中，每一个小立方体的边长为L_C的立方体称之为原胞，在场景的网格空间中，某个维度方向上，原胞的数量为

其中符号

表示向上取整。

根据上述技术方案，步骤S2中初始化粒子通过邻域搜索进行轨迹搜索；

设系统粒子数为N，为了计算这些粒子的运动轨迹，需要获得标记粒子，指当前被分配给计算资源的粒子，周围的其他粒子的信息，从而计算标记粒子和周围粒子之间的相互作用力，进而基于运动方程计算粒子运动轨迹；

寻找到与标记粒子距离小于截止距离的其他粒子的计算行为称之为邻域搜索；

用公式表示邻域NA和邻居NP，

NA_i(t)＝{x|||x-x_i(t)||＜R_C},x∈R^D(1)；

NP_i(t)＝{j|||x_i(t)-x_j(t)||＜R_C,j≠i} (2)；

其中，符号NA_i(t)表示编号为i的标记粒子在t时刻的球体邻域；

半径为截止距离R_C；

x_i(t)表示i粒子的空间坐标位置；

i是粒子编号且为正整数1≤i≤N；

符号NP_i(t)表示由球体邻域NA_i(t)内的其他粒子的编号所构成的集合。

根据上述技术方案，步骤S3中A1粒子通过原胞列表法搜索；

原胞列表法为了寻找到粒子i的邻居NP_i(t)，通常需要遍历系统中其它(N-1)个粒子，计算它们之间的距离，并使用(2)式作为判别标准，获得邻居NP_i(t)，这种方法的复杂度为O(N²)，其基本方法是将空间用边长为L_C的原胞网格化，原胞存储有占据其中的粒子编号的列表，标记粒子只需要在其周围的摩尔邻域内的原胞中寻找其原胞列表，将这些列表所对应的粒子构成的集合NP′_i(t)做为遍历对象，进而基于 (2)式计算出NA_i(t)，在原胞列表法中，原胞边长设置为L_C，可以保证不会有比截止距离小的相邻粒子被漏掉，但是摩尔邻域显然是大于邻域NA_i(t)的，这将导致有些大于截止距离的相邻粒子进入到遍历集合NP′_i(t)中。

根据上述技术方案，如图2-3所示，摩尔邻域由粒子集合NP′_i(t)尽量接近球体邻域内的粒子集合NP_i(t)，选择更小的原胞边长，来构造网格空间；

摩尔邻域用数学公式来描述扩展邻域的定义如下：

其中，符号

表示向下取整，

表示向上取整；

d 表示D维空间中的第d 个维度，L_C表示原胞的边长；

如图4所示；其中，

表示右侧扩展邻域的半径；

如图5所示，标记粒子位置为蓝色原胞中的红色点，其中

分别表示为标记粒子在X和Y轴正、反方向上的扩展邻域半径。

根据上述技术方案，扩展邻域的期望原胞数和期望粒子数为设定比率l＝L_C/R_C为原胞边长和截止距离的比值，令比率l＝n+m，其中整数

小数m＝R_C/L_C-n，0≤m＜1；

如图6所示，当m∈[0,0.5)时，蓝色原胞的三个区域划分；

如图7所示，当m∈[0.5,1)时，蓝色原胞的三个区域划分。

设粒子在空间任意一个位置出现概率相等，则扩展邻域NE_i(t)内所包含原胞数量LC_i(t)的期望值可以表示为，

E(LC_i(t))＝P₁·(2(n+1))+P₂·(2n+1),m∈[0,0.5)；

E(LC_i(t))＝P′₁·(2(n+1))+P′₂·(2n+3),m∈[0.5,1)；

P₁＝2m,P₂＝1-2m,P′₁＝2(1-m),P′₂＝2m-1,

E(LC_i(t))＝2l+1 (4)；

进一步设空间中粒子均匀分布，即密度ρ为空间位置x∈R^D的常数，则可知道扩展邻域NE_i(t)内所包含的粒子数量LP′_i(t)，

E(LP′_i(t))＝ρL_C ^D(2l+1)^D＝ρR_C ^D(2+l^-1)^D (5)。

根据上述技术方案，粒子进行一次遍历并更新位置的运算中，邻域搜索运算和计算距离和作用力是单个粒子的基本运算；

设循环模块对应时间参数分别为σ₁和σ₂；

则一个粒子对应的循环模块运算量为(σ₁LC_i(t)+σ₂LP′_i(t)+τ_c)，其中τ_c是一个时间常量，由式(4)和(5)得到循环模块运算量的期望值为，

T_CE(L_C)＝σ₁(2l+1)^D+σ₂ρR_C ^D(2+l^-1)^D+τ_c (6)；

式(6)是一个关于原胞边长L_C的函数，我们可以进一步寻找满足下式时，

最优原胞边长

式子(7)即为我们的优化目标函数；

对函数(6)求导数，

可以得到：

或：

为待求的最优原胞边长，即最优网格尺度；

由方程式(9)，可以得到，

进而可以将(10)式化为一元线性回归方程形式，

其中，δ是回归方程的参数。

如图8-13所示，以重复三次为例，我们可以得到回归方程式(10)

两个变量所构成的数组

的三十组数据，其中维度D ＝2，密度数组ρ^a＝[6^-2,8^-2,10^-2,12^-2,14^-2,16^-2,18^-2,20^-2,24^-2]， K＝length(ρ^a)＝10。

运用最小二乘法，可对于这三十组实验数据进行参数估计，即令： X＝ρ^-1，

则有方程参数，

如图9所示，水平坐标表示密度数组ρ^a＝[6^-2,8^-2,10^-2,12^-2,14^-2,16^-2,18^-2,20^-2,24^-2]中的元素的出现顺序，纵坐标方向表示变量

的值，刻度为平方分米；图中黑点表示结果来自方程(4)和(5)的理论计算值，绿十字星表示的数据，来自于扩展邻域内的原胞个数及其邻居个数的实际测量值。

与现有技术相比，本发明的有益效果：本发明基于理论分析，提供了一个计算方法，使得我们可以根据粒子的密度、“截止距离”以及空间维度来直接计算最优网格尺度，从而可以将获得的最优网格尺度应用于更大规模的粒子运动模拟系统中，节省计算时间。

最后应说明的是：以上所述仅为本发明的优选实例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.大量相互作用粒子运动仿真系统的网格尺度寻优方法，其特征在于：包括如下步骤：

S1、初始化场景：空间维度，边界长度；

S2、初始化粒子：直径1b，质量，作用力截止距离，作用力方程，运动方程；

S3、系统参数测算：

A1、设置系统粒子数；

A2、设定粒子占有空间的边长分别为2b、3b、4b、5b、6b五个等级；

A3、对于每个密度等级，数值实验重复三次以上，寻找密度等级对应的最优格子尺寸；

A4、基于最小二乘法估计回归方程参数。

2.根据权利要求1所述的大量相互作用粒子运动仿真系统的网格尺度寻优方法，其特征在于，所述步骤S1中空间维度包括场景模型和网络模型；

所述场景模型场景是三维空间，设定其各维度相应的边界长度，从而确定场景的空间范围；

设置场景为R^D:∏L_S ^d；

其中1≤d≤D是正整数；

D表示空间维度；

L_S ^d表示场景R^D的第d 个维度的边界长度，即场景为一个立方体；

所述网格模型以场景立方体的某个顶点做为坐标原点，以长度L_C为间距，在各个维度上进行均分，得到网格化的空间，在这个网格化的空间中，每一个小立方体的边长为L_C的立方体称之为原胞，在场景的网格空间中，某个维度方向上，原胞的数量为

其中符号

表示向上取整。

3.根据权利要求1所述的大量相互作用粒子运动仿真系统的网格尺度寻优方法，其特征在于，所述步骤S2中初始化粒子通过邻域搜索进行轨迹搜索；

设系统粒子数为N，为了计算这些粒子的运动轨迹，需要获得标记粒子，指当前被分配给计算资源的粒子，周围的其他粒子的信息，从而计算标记粒子和周围粒子之间的相互作用力，进而基于运动方程计算粒子运动轨迹；

寻找到与标记粒子距离小于截止距离的其他粒子的计算行为称之为邻域搜索；

用公式表示邻域NA和邻居NP，

NA_i(t)＝{x|||x-x_i(t)||＜R_C},x∈R^D (1)；

NP_i(t)＝{j|||x_i(t)-x_j(t)||＜R_C,j≠i} (2)；

其中，符号NA_i(t)表示编号为i的标记粒子在t时刻的球体邻域；

半径为截止距离R_C；

x_i(t)表示i粒子的空间坐标位置；

i是粒子编号且为正整数1≤i≤N；

符号NP_i(t)表示由球体邻域NA_i(t)内的其他粒子的编号所构成的集合。

4.根据权利要求3所述的大量相互作用粒子运动仿真系统的网格尺度寻优方法，其特征在于，所述步骤S3中A1粒子通过原胞列表法搜索；

所述原胞列表法为了寻找到粒子i的邻居NP_i(t)，通常需要遍历系统中其它(N-1)个粒子，计算它们之间的距离，并使用(2)式作为判别标准，获得邻居NP_i(t)，这种方法的复杂度为O(N²)，其基本方法是将空间用边长为L_C的原胞网格化，原胞存储有占据其中的粒子编号的列表，标记粒子只需要在其周围的摩尔邻域内的原胞中寻找其原胞列表，将这些列表所对应的粒子构成的集合NP′_i(t)做为遍历对象，进而基于(2)式计算出NA_i(t)，在原胞列表法中，原胞边长设置为L_C，可以保证不会有比截止距离小的相邻粒子被漏掉，但是摩尔邻域显然是大于邻域NA_i(t)的，这将导致有些大于截止距离的相邻粒子进入到遍历集合NP′_i(t)中。

5.根据权利要求4所述的大量相互作用粒子运动仿真系统的网格尺度寻优方法，其特征在于，所述摩尔邻域由粒子集合NP′_i(t)尽量接近球体邻域内的粒子集合NP_i(t)，选择更小的原胞边长，来构造网格空间；

所述摩尔邻域用数学公式来描述扩展邻域的定义如下：

其中，符号

表示向下取整，

表示向上取整；

d 表示D维空间中的第d 个维度，L_C表示原胞的边长。

6.根据权利要求5所述的大量相互作用粒子运动仿真系统的网格尺度寻优方法，其特征在于，所述扩展邻域的期望原胞数和期望粒子数为设定比率l＝L_C/R_C为原胞边长和截止距离的比值，令比率l＝n+m，其中整数

小数m＝R_C/L_C-n，0≤m＜1；

设粒子在空间任意一个位置出现概率相等，则扩展邻域NE_i(t)内所包含原胞数量LC_i(t)的期望值可以表示为，

E(LC_i(t))＝P₁·(2(n+1))+P₂·(2n+1),m∈[0,0.5)；

E(LC_i(t))＝P′₁·(2(n+1))+P′₂·(2n+3),m∈[0.5,1)；

P₁＝2m,P₂＝1-2m,P′₁＝2(1-m),P′₂＝2m-1,

E(LC_i(t))＝2l+1 (4)；

进一步设空间中粒子均匀分布，即密度ρ为空间位置x∈R^D的常数，则可知道扩展邻域NE_i(t)内所包含的粒子数量LP′_i(t)，

E(LP′_i(t))＝ρL_C ^D(2l+1)^D＝ρR_C ^D(2+l^-1)^D (5)。

7.根据权利要求4所述的大量相互作用粒子运动仿真系统的网格尺度寻优方法，其特征在于，所述粒子进行一次遍历并更新位置的运算中，邻域搜索运算和计算距离和作用力是单个粒子的基本运算；

设循环模块对应时间参数分别为σ₁和σ₂；

则一个粒子对应的循环模块运算量为(σ₁LC_i(t)+σ₂LP′_i(t)+τ_c)，其中τ_c是一个时间常量，由式(4)和(5)得到循环模块运算量的期望值为，

T_CE(L_C)＝σ₁(2l+1)^D+σ₂ρR_C ^D(2+l^-1)^D+τ_c (6)；

式(6)是一个关于原胞边长L_C的函数，我们可以进一步寻找满足下式时，

最优原胞边长

式子(7)即为我们的优化目标函数；

对函数(6)求导数，

可以得到：

或：

为待求的最优原胞边长，即最优网格尺度；

由方程式(9)，可以得到，

进而可以将(10)式化为一元线性回归方程形式，

其中，δ是回归方程的参数。

Images