技术领域
[0001] 本
发明属于大规模群体建模与仿真技术领域，特别涉及一种大规模群体仿真中的高速并行计算方法，本发明群体指的是人群、
鸟群、车辆群体与
航天器群体等。
背景技术
[0002] 群体仿真在计算机游戏、电影、
虚拟现实等方面有着非常重要的应用，实时群体仿真是其中一个重要的研究方向。传统的群体仿真方法大多是针对单个Agent寻找路径，这类
算法在群体规模较少时非常有效，但是计算量随着仿真人数的增多而快速增加，因此不适用于大规模群体仿真。
[0003] 全局的路径规划是当前群体仿真的重要
基础和前提。路径规划帮助Agent选择从当前
位置前往目标位置的合理平滑路径。路径规划通常要在不同路径之间做出权衡，一般来说，Agent希望能够减少路径的长度、运动时间等。为了事先对全局路径进行计算，图论的某些方法常常被广泛运用在群体仿真的路径规划当中。当下使用较多的路径规划算法是Dijkstra算法和A*算法，它们都是在加权图上寻找最
短路径的著名算法，然而这些算法都是从起始位置开始搜索，其运算规模与群体的规模成正比，因而运行效率比较低，这表示它们并不适用于大规模群体的路径规划。
发明内容
[0004] 为了克服上述
现有技术的缺点，本发明的目的在于提供一种大规模群体仿真中的高速并行计算方法，基于OpenCL使用GPU并行计算方法计算Agent在任意位置前往目标地点的最优方向，可以有效地提高运算效率，突破大规模群体的计算
瓶颈。
[0005] 为了实现上述目的，本发明采用的技术方案是：
[0006] 一种大规模群体仿真中的高速并行计算方法，对大规模群体仿真构建一种面向全局的离散
势能场，采用数值函数来描述仿真区域的几何结构，通过搜索势能场的负梯度方向获取Agent在任意位置前往目标地点的最优运动方向，势能场模型的计算复杂度与群体规模无关，只与地图的精细度及复杂度相关，从而减少了计算量。为了提高运算速度，对快速
迭代算法进行改进，基于OpenCL来进行GPU并行计算，以较高的效率实现大规模群体运动仿真。算法包括如下步骤：
[0007] 引入分区运算规则，使用并行化的快速迭代算法求解势能场模型，包括：
[0008] 1.将仿真区域的所有单元格划分为预定义大小的区
块，在GPU运算术语中称为tile；
[0009] 2.设定运算列表，将包含目标区域单元格的tile移入运算列表；
[0010] 3.对仿真区域进行循环计算：
[0011] 1)将tile映射为OpenCL中的工作组，程序将单元格的势能值存储在局部内存(local memory)中以便提高运算速度。每个tile中的单元格势能值的计算由一个工作项负责。对运算列表中的tile进行更新运算，对tile中所有的单元格计算其势能值，每计算完一次后比较当前结果与前一次的结果，如果差值小于
阈值，则将此单元格记为收敛，循环计算，直到该tile中的所有单元格都趋于收敛。
[0012] 2)对tile检查其四向邻tile，若其中有单元格势能值发生变化，将该tile重新加入计算列表中，若没有变化，将该tile从运算列表中删除。
[0013] 3)重复1)和2)直至运算列表为空。此处设定tile永远处于三种状态之一：更新(Update)、休眠(Sleep)、收敛(Converged)。Update状态表示该tile必须在下一个外循环迭代运算中被更新计算，Sleep状态表示该tile应该挂起而不进行任何更新，Converged状态表示该tile经过上一次运算已经收敛。
[0014] 根据目标区域单元格来初始化各tile的初始状态，即对于含有目标位置单元格的tile，其状态为Update，其余tile则设定为Sleep。
[0015] 当tile处于Update状态时，计算单个单元格此次与前次势能值的差值，如果该差值小于阈值，则该单元格为收敛。遍历tile中的所有单元格，如果均为收敛，则该tile收敛。则状态由Update转变为Converged。
[0016] 当tile处于Sleep状态时，如果该tile的四向邻tile均处于Sleep状态，则该tile仍然保持Sleep状态，因为该tile挂起时周围没有发生信息变化，否则对该tile进行循环一次计算其所有单元格的势能值；如果tile内所有单元格满足此次与前次势能值的差值小于阈值，则表示邻tile的信息虽然发生变化，但并没有对该tile发生影响，该tile的状态仍然保持Sleep；如果tile内任意单元格满足此次与前次势能值的差值大于阈值，说明其邻tile的信息更新将导致本tile内的信息更新，则对该tile进行循环多次计算其所有单元格的势能值，如收敛则状态变为Converged，否则变为Update。
[0017] 当tile处于Converged状态时，如果该tile的四向邻tile均处于Sleep状态，则该tile仍然保持Converged状态，因为该tile收敛时周围没有发生信息变化，否则对该tile进行循环一次计算其所有单元格的势能值；如果tile内所有单元格满足此次与前次势能值的差值小于阈值，则表示邻tile的信息虽然发生变化，但并没有对该tile发生影响，该tile的状态仍然保持Converged；如果tile内任意单元格满足此次与前次势能值的差值大于阈值，说明其邻tile的信息更新将导致本tile内的信息更新，则对该tile进行循环多次计算其所有单元格的势能值，如收敛则状态变为Converged，否则变为Update。
[0018] 与现有技术相比，本发明的有益效果是：
[0019] 1.本发明方法利用快速迭代算法的思想，对其局限性进行改进和优化，引入分区运算规则，将原有的整体计算区域划分为预定义大小的区块，对每个区块的求解是并行的，提高了算法的效率。
[0020] 2.在此基础上实现在GPU的多个处理单元并行计算群体仿真，从而大幅度地提高了大规模群体仿真的计算速度。
[0021] 3.与传统的基于CPU算法相比，充分利用了GPU强大的并行计算能
力，降低了大规模群体建模仿真的经费投入。
附图说明
[0022] 图1是本发明离散势能场示意图。
[0023] 图2是本发明仿真区域离散化分割图。
[0024] 图3是本发明tile示意图。
[0025] 图4是本发明实验仿真区域示意图。
[0026] 图5是本发明势能场梯度全局图。
[0027] 图6是本发明势能场梯度局部放大图。
[0028] 图7是本发明最优路径示意图。
具体实施方式
[0029] 下面结合附图和
实施例详细说明本发明的实施方式。
[0030] 本发明采用势能场模型，采用一个数值函数来描述仿真区域的几何结构，通过搜索势能场的负梯度方向获取最优运动方向，如图1所示。势能场模型的计算复杂度与群体规模无关，只与地图的精细度及复杂度相关，从而减少了计算量。为了提高运算速度，采用GPU来
加速势能场模型的计算，实现了大规模的群体运动仿真。
[0031] 人在选择前往目标的路径时，总是尽量使得某些代价最小。通常情况下，人们倾向于选择路程最短的路径，这是最常见的想法。人们有时更希望能够以自己最理想的速度运动，寻找耗时最短的路径。人们有时会有自己更偏爱的区域，也有自己不希望去到的区域，例如，有些区域群体聚集，即使此时路径较短，人们还是会选择其他人少的区域通行，本发明用不适度指数来表示这一特性。
[0032] 代价函数与路径长度、耗费时间、不舒适度三个指标有关，耗费时间与路径长度均和Agent的运动速度有关，由此推导代价函数的计算公式：
[0033]
[0034]
[0035] 式中：l——Agent的运动路径；El——路径l的代价函数；wl——路径长度所占权重；wt——耗费时间所占权重；wc——不适度值所占权重；c——不舒适度函数；v——Agent的运动速度；e——路径的单位代价函数。
[0036] 在离散化的仿真区域内，路径可以看成一系列的有序单元格序列，Agent从某个单元格沿着某方向运动可以进入其上、下、左、右的相邻单元格，因此，单条路径的代价的积分运算可以用从单元格进入上、下、左、右的相邻单元格的单位代价值之和进行近似，因此公式(1)可以转化为：
[0037] El＝∑d∈{W,N,E,s}ei,j,d (2)
[0038]
[0039] 式中：ei,j,d——从单元格(i，j)沿着方向d进入相邻单元格的代价值； ——单元格内Agent沿着方向d运动的可达速度；ci′j′——对应方向的相邻格内的不舒适度指数，如超出边界则为正无穷。
[0040] 公式(1)定义了仿真区域内从当前位置到目标位置的一条路径的代价值，则平面空间内从起始位置S到目标区域G的最优路径C(s)的搜索过程可以转化为如下的约束优化问题：
[0041] min∫e(C(s))ds
[0042] s.t.C(s)＝{x(s),y(s)},C(0)＝S,C(L)＝G, (4)
[0043] 式中：e——路径的单位代价函数；L——路径的长度；S——初始位置；G——目标区域。
[0044] 可以证明，平面上A点到B点的最优路径，即满足公式(4)的曲线为标量场p的梯度下降曲线，其中势能场p满足方程：
[0045]
[0046] 仿真区域中，目标位置往往是单一且明确的，因此势能场的计算通常从目标位置开始进行回溯计算，由于我们关心的是势能场的梯度而非其绝对数值，因此为简化计算，设定目标区域的势能值为0。综上，势能场的构建函数可用如下公式表示：
[0047]
[0048] 式中：G——表示目标区域，p表示势能场，e为路径的单位代价函数，p(G)＝表示目标区域的势能值， 表示势能场梯度的模值。
[0049] 对该函数进行求解可以得到区域内任意位置的势能值。
[0050] 在群体仿真系统中，Agent的数目与计算复杂度成正比，为了实时完成大规模群体的势能场运算，本方法采用OpenCL编程模型来实现并行化快速迭代算法求解公式(6)，从而使用GPU进行并行运算。
[0051] 引入分区运算规则，将仿真区域离散化处理，如图2所示，再将所有单元格划分为预定义大小的区块，在GPU运算术语中称为tile，如图3所示。其中每个tile中的单元格的势能值的求解是并行进行的，从而实现GPU上的并行运算。
[0052] 使用并行化的快速迭代算法求解势能场模型，包括如下步骤：
[0053] 1)定义：
[0054] (1)设定t表示将仿真区域进行等尺寸划分的单个tile；
[0055] (2)设定运算列表L，用于存储需要更新的tile；
[0056] 2)初始化：
[0057] (1)将所有单元格均匀划分为相同大小的区块tile；
[0058] (2)将所有包含目标区域单元格的tile移入运算列表L；
[0059] 3)循环计算：
[0060] (1)对运算列表L中的每一个区块t进行更新运算：
[0061] a)循环n次计算t中的所有单元格的势能值(此循环运算简称为内循环运算，下同)；每计算完一次后比较当前的结果与前一次的结果，如果差值小于阈值，则该网格标记为收敛。
[0062] b)遍历t中的所有单元格，如果所有单元格都被标记为收敛，则t标记为收敛，否则继续进行上一步中的循环。
[0063] (2)检查区块t的四向邻tile：
[0064] a)对上一步中的已收敛的所有区块t，对其所有相邻区块进行一次势能值计算，检查t的邻区块中是否有单元格势能值发生变化；
[0065] b)将上步中发生变化的区块t重新加入运算列表L，势能值没有发生变化的t从L中删除；
[0066] (3)对L包含的所有未收敛的tile，重复步骤(1)和(2)直至L为空(此循环运算称为外循环运算，下同)。
[0067] 基于OpenCL，将并行快速迭代算法在GPU上实现，包括如下步骤：
[0068] 1)内循环更新实现方法
[0069] 根据OpenCL的工作原理，势能场的kernel根据二维索引空间处理仿真区域中的单元格值来计算，上述并行快速迭代算法中的tile即可映射为OpenCL中的工作组，程序将单元格的势能值存储在局部内存(local memory)中以便提高运算速度。每个tile中的单元格势能值的计算由一个工作项负责。势能场模型按照目标不同将Agent分为不同的组，设定Agent最多可以分为四组，即每个单元格中存储四个分别属于不同组的势能值，每个单元格的四组势能值(P0，P1，P2，P3)将由OpenCL的向量化数据类型float4来表述。
[0070] 势能场计算的kernel程序将势能值载入local memory以进行内循环运算。单元格的势能值计算需要其上、下、左、右邻单元格的势能值参与。
[0071] 根据程序定义的内循环次数n，kernel对所有local memory中存储的单元格进行势能值的计算。
[0072] 2)外循环运算实现方法
[0073] 外循环运算是针对tile的循环运算，设定两个整形变量StateBefore和StateAfter分别存储各tile经过处理之前和之后的状态，设定tile永远处于三种状态之一：更新(Update)、休眠(Sleep)、收敛(Converged)。Update状态表示该tile必须在下一个外循环迭代运算中被更新计算，Sleep状态表示该tile应该挂起而不进行任何更新，Converged状态表示该tile经过上一次运算已经收敛。
[0074] 算法根据目标区域单元格来初始化各tile的初始状态，即对于含有目标位置单元格的tile，其状态为Update，其余tile则设定为Sleep。
[0075] 在每一次外循环更新中，kernel读取StateBefore的值来决定对tile的处理方式。实现时，每个工作组中处理该tile的第一个工作项读取该tile的状态信息，存入local memory。该工作组中的其余工作项则依据该状态信息处理其余单元格。该tile经过处理后，其状态发生转变，其新状态仍由首个工作项写入StateAfter，这个状态转换过程记作StateBefore→StateAfter。
[0076] 针对单个tile的更新，定义两种操作：MultiStep和SingleStep。MultiStep操作表示对单个tile连续进行多次内循环运算，SingleStep操作表示对其进行单次内循环运算。
[0077] 当tile处于Update状态时，进行一次MultiStep操作，计算单个单元格此次与前次势能值的差值，如果该差值小于阈值，则该单元格为收敛。遍历tile中的所有单元格，如果均为收敛，则该tile收敛。则状态由Update转变为Converged。
[0078] 当tile处于Sleep状态时，如果该tile的四向邻tile均处于Sleep状态，则该tile仍然保持Sleep状态，因为该tile挂起时周围没有发生信息变化，否则对该tile进行SingleStep操作；如果tile内所有单元格满足此次与前次势能值的差值小于阈值，则表示邻tile的信息虽然发生变化，但并没有对该tile发生影响，该tile的状态仍然保持Sleep；如果tile内任意单元格满足此次与前次势能值的差值大于阈值，说明其邻tile的信息更新将导致本tile内的信息更新，则对该tile进行一次MultiStep操作，如收敛则状态变为Converged，否则变为Update。
[0079] 当tile处于Converged状态时，如果该tile的四向邻tile均处于Sleep状态，则该tile仍然保持Converged状态，因为该tile收敛时周围没有发生信息变化，否则对该tile进行循环一次计算其所有单元格的势能值；如果tile内所有单元格满足此次与前次势能值的差值小于阈值，则表示邻tile的信息虽然发生变化，但并没有对该tile发生影响，该tile的状态仍然保持Converged；如果tile内任意单元格满足此次与前次势能值的差值大于阈值，说明其邻tile的信息更新将导致本tile内的信息更新，则对该tile进行循环多次计算其所有单元格的势能值，如收敛则状态变为Converged，否则变为Update。
[0080] 建立如图4的多障碍物迷宫环境，图中三
角形表示出发位置，五角星表示目标位置，多边形代表障碍物。
[0081] 设定环境中没有其余Agent，仿真区域被划分为256*256的离散网格，设定OpenCL的工作组大小为16*16。根据势能场的计算模型，此处采用快速迭代算法计算势能场。
[0082] 实验程序设定将各单元格的梯度方向以箭头表示，如图5所示，其中图6为图5中阴影矩形框区域的放大图。
[0083] 任意位置前往目标位置的最优方向为势能场的负梯度方向，因而从初始位置的单元格沿负梯度方向连线，即可得到最优路径，如图7中折线所示。
