并行化仿真

2024-08-07

并行化仿真（精选七篇）

并行化仿真篇1

大规模的实时人群仿真在计算机游戏、电影、虚拟现实等方面有着非常重要的应用, 实时人群仿真是其中一个重要的研究方向[1,2]。基于Agent的人群仿真模型, 比较符合真实世界中人群运动的特点, 并且具有非常灵活的可控性, 比如可以控制每个Agent的每个参数。因此基于Agent的人群仿真模型被使用得非常广泛。但是对于大规模的实时人群仿真, 由于每个Agent都要根据其周围的环境状况来更新自身状况, 并且每个Agent相对其他Agents来说都是一个动态障碍物, 系统计算量随着仿真人数的增多而快速增加[3], 时间复杂度为Agent数目n的平方。如此一来, 当n比较大时, 就会严重影响计算每一帧所需的时间, 导致无法实时更新每个Agent的状态。

近年来, 图形处理单元GPU由于其强大的并行计算能力及可编程性, 被广泛地应用在各种通用计算领域。传统意义上的GPU, 是一个专用于各种图形处理的器件, 主要的计算单元包括vertex shader和pixel shader, 在图形流水线中分别对vertex和pixel执行并行操作, 技术的不断进步, 使得他们的可编程性以及单个芯片上可集成的计算单元数越来越高, 慢慢地演变成为通用计算单元, 使其能更容易地应用到各种非图形学领域的计算, 这些计算单元一般采用单指令多数据并行计算模式。单从计算能力上来说, 现在的GPU已经要远胜于CPU了。在人群仿真中, 我们发现单个Agent的状态更新只需要考虑它周围某个范围内的世界状况, 基于上述观察, 我们设计了一种并行化的人群仿真方法, 采用部分并行化更新策略, 通过把场景均匀分割成很多bin (每个bin为一个正方形区域, Agents分布于各个bins之中) , 从而确定每个Agent所需查询的范围。对于那些非相邻的bins, 它们之间的Agents可采用并行化更新策略, 共分四个批次的并行更新, 这样就利用了GPU并行计算能力。实验结果表明这种策略可以大大提高更新的效率。

在现实生活中, 每个人走在路上都会根据他自己周围的环境状况作出相应决策, 比如前面有一辆车过来, 他就会根据车的行驶方向、行驶速度, 及距离的远近, 作出相应的判断, 或者保持现状或者改变自己的行走方向及速度。行人之间也会尽量避免发生碰撞。文献[4]中提出的Agent模型可以很好地模拟这些情况。本文采用了类似但相对简单的Agent模型, 该Agent受到三个力的作用[5]。首先是躲避障碍物的作用力, 其次是避免与其他Agent相撞的作用力, 最后是路径跟随的作用力。三种作用力的优先级依次降低。计算每个Agent受到这三个力的合力, 根据这个合力来更新行走速度。

本文提出的基于GPU的并行化人群实时仿真算法, 它充分挖掘了人群行为中的并行性, 使得大规模人群实时仿真成为可能。并介绍本文在近邻查询和人群更新方面所采用的并行化算法。

本文采用的GPU编程工具为NVIDIA公司开发的CUDA (Compute Unified Device Architecture) [6]语言, 其设计思想来源于斯坦福大学设计的Brook架构[7,8]。CUDA语言在C语言基础上进行了一些扩展。CUDA程序主要分成两个部分, HOST和DEVICE, HOST运行在CPU上, DEVICE运行在GPU上, HOST程序调用DEVICE程序来完成并行计算功能, 两者的通信主要通过CUDA的RUNTIME API来完成。

1并行化的近邻查询

在基于Agent的人群仿真中, 每个Agent都有自己的可感知范围。在每次更新时, 每个Agent必须查询自己的可感知范围, 以获取该范围内的其他Agents或者障碍物的信息, 然后根据查询到的周围状况计算得到将作用于该Agent的力, 调整自己的速度大小及方向, 并更新位置。一般基于CPU的近邻查询采用串行化查询, 即轮到当前Agent更新时才执行近邻查询, 然后是下一个Agent。在本文算法中, 我们采用并行化的近邻查询, 查询的主体是bin, 而不是每个Agent, 通过这样的近邻查询, 就可以得到该bin中Agents需要考虑的所有潜在邻居Agents。执行并行化近邻查询前, 对整个场景进行均匀分割, 得到所有bin的信息。

1.1场景分割

在进行并行化近邻查询前, 本文对整个场景进行分割, 分割过程在初始化时完成。我们整个场景的大小为1200×1200, 设定每个bin为边长是4的正方形, 这样整个场景就被分割成300×300个bins。采用边长为4的正方形, 是因为本文设置Agent的可感知范围是一个查询半径为4的圆形区域。其计算公式如下:

queryRadius=predTime×maxVelocity×2

predTime为前向预测时间, 设定为1s, maxVelocity为Agent的最大行走速度, 为2m/s, 在查询过程中, 我们只需要考虑以该bin为中心的8个bin即可。

如图1所示, 场景被均匀分割成很多个bin, 场景中的Agents分布在各个bin中, 既有含有多个Agents的bin, 又有一个Agent也不包含的bin。

1.2近邻查询

完成场景分割并统计完Agents在各个bin之中的分布情况之后。我们以bin为主体, 查询每个bin的近邻Agents, 在场景分割步骤中, 由于我们设定每个Agent的观察范围是半径为4的圆形区域, 分割出来的bin边长也为4, 所以我们只需要考察该bin周围的8个bin及其自身即可。如图2所示, 编号为49的绿色bin, 它的观察范围为虚线圆圈, 需要查询的bin编号分别为37, 38, 39, 48, 49, 50, 59, 60, 61。由于bin之间的近邻查询只是对数据的读取, 不存在写操作, 所以我们采用并行化的查询方案, 查询结果保存在相应的数据结构里。

在具体的计算过程中, 每个bin用8个threads计算得到与它的近邻Agents, 由于CUDA的kernel函数在GPU上执行是以32个threads组成一个wrap, 运行时以wrap为一个单元, 为了充分利用资源, 我们把4个bin总共32个查询threads组成一个block, 比如bin0, bin1, bin11, bin12, 我们把计算结果存储到一个数据结构里—binNeighberAgent, 详细的 CUDA kernel伪代码见图3, 在伪代码中首先getBinId函数根据当前的blockIdx与threadIdx得到currentBinIdex, 然后用8个线程并行查询其周围的8个bin, 如果其不为空, 则把该bin内的Agents暂时放入tempAgentNeighber数组里。等所有的相邻bin同步查询完成后, 在threadIdx.x为0的线程里把tempAgentNeighber存放的Agents及currentBin里的Agents一起存到binNeighberAgent里。

2并行化的人群更新

人群速度和位置的更新是在一个Agent的近邻查询完之后执行的, 然后再是下一个Agent的更新。在前一步并行化的近邻查询中, 我们已经查询完所有bin的近邻, 但是由于相邻bin内的Agents会相互作用, 这一步的更新操作无法完全并行化。比如bin 48和bin 49, 如果两者并行更新, 那就会产生读写矛盾。因为bin48内的Agents在更新状态时要考虑bin49内的Agents, 读到的数据可能是过时的bin49内的Agent数据, 反之亦然。所以本文采用部分并行化更新策略, 具体步骤如下所述。

2.1Agent模型

参照现实世界中行人的行为规则, 及以往文献[5,9,10,11]中的Agent模型资料, 我们设计了本文所采用的Agent模型。本文的Agent模型受3类力的作用:1) 障碍物躲避力, 2) Agents之间的避让力, 3) 路径跟随力。每类力的形象展示如图4所示。

三种作用力的优先级依次下降。障碍物的躲避行为具有最高的优先级, 比如当Agent k检测到他将会碰撞到obstacle d, 但同时他又有可能与Agent g相撞, 那么他就会只考虑躲避obstacle d。因为障碍物是固定不动的, 如不及时调整自己的速度, 就一定会撞上, 而Agent g是会移动的, 就算Agent k没有针对其调整, Agent g本身也会针对Agent k调整速度, 躲避力的计算方法如下:Agent的向前预测距离与其观察半径组成一个矩形区域, 在该区域内找到会发生碰撞的且距离最近的obstacle, 图中为obstacle d, 躲避力 (Force k) 为当前Agent的前向拉力与障碍物排斥力的合力方向。另一种情况是当Agent j向前预测到j′位置, 并且预测附近的Agent i将会同时到达i′位置, 当j′与i′的中心点之间距离小于它们的观察半径之和时, 就认为会发生碰撞, 那么他们就不会受路径跟随力的作用, 而是受到Agent相互之间避让力的影响, 该避让力为一个侧向的拉力 (Force j) , 使其偏离当前行使方向。在图4中路径有一个路宽 (与路径平行的虚线和路径之间的距离) , 当前待更新Agent会向前预测一段距离, 把预测的位置P投影到路径上P′点, 如果P′与P之间的距离大于路宽, 则说明Agent的行走方向偏离了路径方向, 否则没有偏离, 如图Agent h, 他会受路径跟随力的作用 (Force h) , 为投影点P′与Agent的中心连线方向, 使其慢慢靠拢到路径方向。

2.2人群更新

本文分四次并行更新Agent速度。如图2所示, 总共有四种颜色的bin, 每次并行更新颜色相同的bin, 直到所有颜色的bin都更新完后, 才完成一次总体更新。同一个bin内的Agents相互之间会有影响, 无法并行更新的, 所以只能采用串行更新。在具体执行过程中, 每个block代表一个bin, 我们遵循上述Agent模型的同时, 根据各类力的优先级, 进行了适当的调整。具体速度更新程序伪代码如图5所示, 首先在threadId.x为0的线程中进行一些初始化工作, 得到binAgentNum, AgentNeighberNum, 以及AgentVelocity数组。在计算作用力的过程中, 由于障碍物避让力优先级最高, 并且计算该力对每个Agent来说都是独立的, 所以我们采用并行计算, 一个thread代表一个bin内Agent, IsIntersectWithObstacle函数根据Agents的速度及minTimeToCollision, 计算得到最近的相撞Obstacle, 如果存在这样的Obstalce, 则根据模型中的介绍计算得到避让力, 然后更新Agent速度。接下来就串行更新每个Agent, 遍历所有该bin内的Agents, 如果没有被更新过, 则先利用函数IscollisionWithNeighber计算得到与所有neighber相撞需要的时间与距离, 从中找到最有可能会相撞的neighber, 如果没有找到, 则计算路径跟随力, 相应的作用力计算方法参考上节, 然后更新Agent的速度。最后把更新的结果写到Agent的数据库里。

上述人群并行化更新算法伪代码如图 5所示。

在CPU主程序中, 调用四次updateAgentVelocity (offset) kernel, 更新所有Agents的速度, 然后根据更新后的速度得到新的position, 重新分配Agents到所有bin中。这样就完成了一次整体更新。

3实验结果

本文的实验平台为一台PC电脑, CPU为4核 (Q9550) , 每个核主频都为2.83GHZ, 内存4G, 显卡为NVIDIA GTX280, 含有30个multiprocessors, 每个multiprocessors8个核, 共240核。

如图6所示为使用GPU并行人群仿真算法和CPU串行人群仿真算法的帧速率比较。

从图6中可以看出, 基于GPU并行化计算的人群仿真速度明显要高于基于CPU的人群仿真, 随着Agent数目的增多, 基于CPU的FPS下降得非常快, 当人数达到5000以上时, 已经无法实现实时仿真的效果, 而GPU程序在Agent数目为10000以上时, FPS还能达到24帧以上, 完全满足了实时性, 在Agent个数为1000时, GPU的FPS为60。

从图中曲线还可以看出, 虽然基于GPU算法的FPS明显高于基于CPU算法的FPS, 但仍然没有达到十倍甚至数十倍的倍率。这主要是由于在Agent速度更新时, 存在着大量的条件性判断语句, 这类语句的存在降低了基于GPU算法的并行化计算效率。

4结论

通过基于GPU的并行化计算算法, 本文实现了大规模人群的实时仿真。利用场景的均匀分割, 将每个Agent根据他自身位置分配到相应的bin里, 将原本必须串行化的计算转化为一种并行化计算。每个bin作为一个独立的更新单元, 充分利用了GPU在并行化计算方面的能力。实验结果表明, 本文方法在更新效果上与基于CPU的方法一致, 但是在效率上有了很大的提高。

目前, 本文所使用的Agent模型相对较为简单, 如何使用更复杂的模型来实现更加逼真的实时人群仿真, 如何减少基于GPU算法框架中逻辑判断型计算的出现次数, 如何进一提高算法的并行性, 这些都是下一步工作的重点。

摘要：针对大规模实时人群仿真往往需要在极短时间内完成大量的计算, 而一般的串行算法无法完成该类计算任务的现状, 提出了一种基于GPU (Graphic processing unit) 的并行化人群实时仿真算法。该算法利用人群仿真中人群个体的独立性, 将个体近邻查询, 个体位置更新两个计算过程并行化, 从而实现了大规模人群的实时仿真。

关键词：并行计算,并行近邻查询,人群仿真,图形处理单元

参考文献

[1]Farenc N, Musse S R, Schweiss E Kallmann.AParadigmfor ControllingVirtual Humans in Urban Environment Simulations[J].Applied Artifi-cial Intelligence, 2000, 14 (1) :69 91.

[2] Tecchia F, Loscos C, Conroy R.Agent Behaviour Simulator (ABS) : A Platform for Urban Behaviour Development[C]//Proceedings of Games Technology Conference 2001.

[3] Lavalle S M.Planning Algorithms[M].Cambridge University Press, 2006.

[4]Reynolds C W.Flocks, herds, and schools:A distributed behavioralmodel[C]//Computer Graphics (SIGGRAPH’87 Proceedings) , MCStone.1978, 21:25 34.

[5] Craig W. Reynolds.steering behaviors for Autonomous Characters[C].Sony Computer Entertainment America 919 East Hillsdale Boulevard Foster City, California 94404.

[6]Nvidia Corporation:Cuda Compute Unified Device Architecture Pro-gramming Guide[G].2007.

[7]Buck I, Foley T, Horn D, et al.Brook for GPUs:Stream Computing onGraphics Hardware[J].ACM Transactions of Graphics, 2004:777786.

[8] Buck I, Hanrahan P.Data Parallel Computation on Graphics Hardware [R]. Tech. Report 2003-03, Stanford University Computer Science Department, 2003.

[9]Sallach David, Charles Macal.The simulation of social Agents:an in-troduction[J].Special Issue of Social Science Computer Review, 2001, 19 (3) :245 248.

[10] Charles M Macal, Michael J North.Tutorial on Agent-based modeling and simulation part 2: how to model with Agents[C]//Association for Computing Machinery, Winter Simulation Conference, 2006.

并行再生制动控制策略仿真研究篇2

随着能源危机的加剧,混合动力汽车和电动汽车已经成为新一代汽车的发展方向,而再生制动技术作为混合动力汽车和电动汽车一项关键的节能技术,具有突出优点,已经得到越来越大的重视。所谓再生制动是指汽车在制动或下坡时将储存于车身上的势能和动能,通过电机转化为电能,并储存于电储能装置中的过程[1]。电池中,实现能量回收,同时产生车辆所需的部分制动力,既实现了车辆的减速和制动,又有效地降低了整车的燃油消耗和污染物排放,还减少了制动器摩擦片的磨损。

国内外对再生制动的研究已经开展了10多年时间,国外在混合动力汽车和电动汽车再生制动系统的研究上取得了很大的进展,特别是各大汽车公司,已经在量产的混合动力汽车和电动汽车上普遍采用该系统。国内各汽车厂商、科研院都在这一方面进行了研究,也取得了一定成果,但由于混合动力汽车和电动汽车的研究尚处于起步阶段,因此再生制动这一关键技术也显得相对薄弱,目前国内再生制动技术主要面临如何提高汽车制动稳定性和能量回收有限等一些问题。随着混合动力汽车和电动汽车技术的发展,再生制动技术也会日趋成熟,国产的混合动力汽车和电动汽车上也会越来越多地采用再生制动系统。

再生制动对电动汽车的能量利用率的提高和混合动力汽车的燃油经济性、排放性和行驶安全性有直接影响,是电动汽车和混合动力汽车的关键技术之一。

1 再生制动系统结构与工作原理

再生制动系统的结构如图1所示,由驱动轮、主减速器、变速器、电机、AC/DC转换器、DC/DC转换器、能量储存系统以及控制器组成[2]。

汽车在制动或滑行过程中,根据驾驶员的制动意图,由制动控制器计算得到汽车需要的总制动力,再根据一定的制动力分配控制策略得到电机应该提供的电机再生制动力,电机控制器计算需要的电机电枢中的制动电流,通过一定的控制方法使电机跟踪需要的制动电流,从而较准确地提供再生制动力矩,电机的电枢中产生的电流经AC/DC整流再经DC/DC控制器反充到电化学蓄能装置中保存起来。

2 并行制动力分配控制策略

汽车制动时前、后轮同时抱死,对附着条件的利用、制动时汽车的方向稳定性均有利。对轿车来说,只用比例阀就可以满足制动稳定性和附着系数利用率高的要求。结合混合动力车辆的特点,提出一种分配控制策略——并行制动。所谓并行制动是指再生制动与机械制动以固定的关系分享驱动轮制动力。也就是说,驱动轮制动力等于再生制动力与机械制动力总和。并行制动的控制策略见图2所示。

当需要的制动力较小时,由再生制动力单独作用,其中包括对发动机制动的模拟;当需要的减速度增大时,再生制动力所占的比例逐渐减小,机械制动力开始起作用;当总制动力大于一定值(图例为0.7G)时意味着这是一个紧急制动,再生制动力减小到零,机械制动提供所有的制动力;当所需的制动减速度在两者之间时,再生制动与机械制动共同作用。

3 并行制动力分配控制策略实现

ADVISOR中控制算法的实现是通过修改制动力分配控制模型来实现的。修改后的并行制动力分配控制模型如图3所示。

控制算法的输入量为总制动力,由制动踏板力传感器得到。再生制动力由存储在芯片中的再生制动力控制曲线得到。前后轮总摩擦制动力等于总制动力减去再生制动力,前后轮摩擦制动力通过制动回路中的液压比例阀实现分配。控制算法的输出量为前轮(驱动轮)再生制动力、前轮摩擦制动力和后轮摩擦制动力。其中前轮(驱动轮)再生制动力由电机控制单元实现控制。

4 仿真结果及分析

汽车正常行驶时,频繁的启动、停车、加速、制动,据统计我国市区车辆的行驶工况符合图4的速度时间历程。为了研究并行制动控制策略的有效性,有必要对循环工况进行仿真分析。

图5是电池荷电状态(SOC)情况,可见电池的SOC值在制动过程中有上升的趋势。

表2给出了并行制动力分配控制策略在十五工况下回收的总能量和回收效率。从表中也可以看出并行制动控制策略可以回收50.09%的制动能量,有效能量回收效率为16.7%。

5 结论

本文提出了一种再生制动控制策略,在ADVISOR中建立了仿真模型,并且嵌入到Saturn SL1车型并联混合动力汽车。在循环工况下进行了仿真,结果表明并行制动控制策略可以回收50.09%的制动能量,有效能量回收效率为16.7%。

摘要：再生制动系统是混合动力汽车和电动汽车特有的系统。该系统可以将汽车制动过程中消耗的汽车动能和势能通过电机发电的方式储存到电池中,在启动和加速过程中加以利用。装备有再生制动系统的汽车在制动时,主要由前轴摩擦制动力、后轴摩擦制动力和再生制动力三个力作用。提出一种控制策略来解决这三个力的分配问题,结合Saturn SL1车型,在市区十五工况下进行了仿真,结果表明,并行控制策略既可以达到制动安全性的要求又可以达到最大限度的回收制动能量的目的。

关键词：混合动力汽车,再生制动,控制策略

参考文献

[1]陈全世,仇斌,谢起成.燃料电池电动汽车[M].北京:清华大学出版社,2005.

[2]Sasaki Y,Mtomo A,Kawahata F,et al.Toyota Braking System for Hybrid Vehicle with Regenerativ System[C]//Proceedings of the14th International Electric Vehicle Symposium(EVS-14),Florida,1997.

并行化仿真篇3

关键词：SQP,并行优化,多领域仿真,Modelica

0 引言

Modelica语言对不同领域物理系统的模型进行统一表述, 实现了统一建模, 并支持层次化结构建模[1]。文献[2,3]基于Modelica 语言的复杂机械系统建模与仿真平台, 实现了工程实际中复杂机械系统多领域统一可视化建模、仿真和优化。目前国内外已经推出了基于Modelica的复杂机械系统的建模与仿真软件。

在MWorks中, 优化模块是其核心之一, 目前集成了约束变尺度法、拟牛顿法、Hook-Jeeves算法和遗传算法 (genetic algorithm, GA) 。约束变尺度法又被称为序列二次规划 (sequential quadratic programming, SQP) 算法, 具有收敛快、效率高、可靠性与整体收敛性好、适应能力强等一系列优点[4], 并被成功应用于工程优化问题的求解。序列二次规划算法通过求解一系列的二次规划子问题 (quadratic programming, QP) 来获得原问题的最优解, 而其二次规划子问题的精度依赖于原问题的目标函数和约束函数梯度值的精度, 并对问题的收敛性有较大影响。对于基于多学科仿真的多实例 (即多工况或一个系列产品) 多目标仿真优化问题, 当采用差分法计算梯度时, 由于每次迭代需要多次进行仿真函数评估 (评估次数等于设计变量维数的两倍乘以实例数+1) , 优化求解过程很费时。因此研究SQP并行优化算法很有必要。

1 SQP算法的并行化

SQP算法并行化研究, 国内外已取得了一定的成果。High等[5]从全局优化、函数评估、梯度评估与Hessian矩阵的有限差分计算等方面讨论了SQP算法并行化, 利用SCHEDULE软件在Alliant FX/8并行机上进行了大量的算例运算。张帆等[6]讨论了大规模化工过程系统中, 优化算法SQP的并行化和机群计算效率的问题, 并通过算例证明其合理性。陈忠[7]叙述了关于非线性规划问题的求解有几种主要并行思想, 简要地综述了值优化的进展。

总之, SQP算法并行化研究主要集中在传统优化的基础上并结合具体的问题来展开, 而基于多学科仿真多实例多目标的SQP优化算法并行化研究甚少, 且绝大多数的研究是针对专用并行机的。因此, 本文从基于多学科仿真的多实例多目标一般性优化出发, 研究SQP梯度计算并行化和并行仿真评估。

在优化算法中, 梯度计算是最耗时间的步骤, 同时也是设计优化的瓶颈[8], SQP法也不例外。然而, 因为各个变量的灵敏度计算是互相独立的, 因此并行化很适合SQP法。

SQP并行优化算法步骤如下:

(1) 给定初始值X (0) 、λ (0) , 其中, X为设计变量, λ为Lagrange乘子。

(2) 并行仿真评估、并行计算梯度值[6]。应用集中式动态负载平衡技术[9]实现并行计算梯度。工作池技术能很好地实现上述过程, 如图1所示。

(3) 构造二次规划子问题。

(4) 求解二次规划子问题, 确定新的乘子向量λ (k) 和搜索方向d (k) 。

(5) 沿d (k) 方向进行一维搜索, 确定步长αk, 得到新的约束极小值点。

(6) 满足[6]:

‖∇xL[X (k+1) , λ (k+1) ]‖2≤ε1 (1)

式中, L为Lagrange函数;ε1为收敛精度。

或者满足:

$\begin{array}{l} c_{i} [X^{(k + 1)}] \leq ε_{2} i \in E \cup Ι \\ \frac{| f [X^{(k + 1)}] | - f [X^{(k)}] |}{f [X^{(k)}]} \leq ε_{3} \end{array}} (2)$

式中, f为优化的目标函数;c为优化约束;E为等式约束集合;I为不等式约束集合;ε2、ε3均为收敛精度。

则停止计算。否则转步骤 (2) 。

(7) 求解Hessian近似阵[6], 令k←k+1, 返回步骤 (3) 。

2 并行调度算法

2.1 问题描述

在分布式环境下, 由于各个处理机的性能和功能不尽相同, 一个任务需要一个或多个处理机同时执行才能完成, 这样的任务被称为多处理机任务。一般地, k个处理机的多处理机系统P={p1, p2, …, pk}, 多处理机任务集合为J={j1, j2, …, jn}, 其中每个任务ji (i=1, 2, …, n) 具有二元属性 (si, ti) , si⊆P为执行任务需要拥有的一组处理机, 通常称之为处理机模式[10], ti为这组处理机同时处理该任务所需要的时间量, 称为时间长度。

并行调度就是为提交给系统的每一个任务分配一定的资源, 并指派占用这些资源的起止时间, 目的是在总体上达到时间跨度最小, 这是典型的多处理机任务调度问题, 记为Sk|fix|Cmax[11], Sk为k个独立的处理机, fix为调度实例的约束, 包括处理机的分配方式、任务之间的优先性、资源需求等方面的约束;Cmax表示要使最后完成时间达到最小, 提高系统的总吞吐量。

在多实例多目标并行优化中, 调度模型Sk|fix|Cmax有些特殊, 其中fix表示每一个任务只需要一个处理机。如3个变量、3个实例的多实例、多目标优化问题, 每迭代一次有21个并行任务, 可以描述为J={j1, j2, …, j21}, 其中, J中每一个任务都只需要一个处理机。

现在我们研究N个处理机环境下的并行优化问题。假设优化问题有n个设计变量 (调节参数) 、m个模型实例, 那么多实例、多目标并行优化并行任务可以描述为

Jiter={ji= (si, ti) |i=1, 2, …, m (2n+1) ;

si=pi} (3)

其中, iter为优化迭代次数。由于优化的每一次迭代任务的模式都相同, 只是其中的数值不同, 因此我们只分析其中的一次迭代就够了。执行时间为titer, 且

$t_{i t e r} = \max {\sum_{k = 1}^{q_{i}} t_{k}, i = 1, 2, \dots, Ν$ ;qi<m (2n+1) } (4)

显然, 任意一个调度的时间跨度不小于N个性能最好的处理机执行的总时间, 我们得到一个下界, 令

B=max{tk|k=1, 2, …, N}

则Opt (J) ≥B, 其中, Opt (J) 表示优化函数。

2.2 调度算法

下面讨论多实例、多目标并行优化S5|fix|Cmax问题的一个调度算法。

输入:任务实例J={j1, j2, …, j21}, 其中, ji= (pi, tk) ;i=1, 2, …, m (n+1) ;k=1, 2, …, N。

输出:实例J的一个调度。

调度算法步骤如下:

(1) 得到所有连接的客户端, 评估其性能, 从高到低排序, 生成处理机性能队列。

(2) 以单位性能任务数和处理机性能数为依据分配任务。将单位性能处理机看作有两个线程的计算节点, 那么其他处理机以其性能数为依据, 决定其启动的线程数。如5处理机系统, 其性能数分别为1.0、1.5、2.0、1.0、2.5, 那么它们同时启动的线程数分别为2、3、4、2、5。每一个节点处理机可以同时启动两个或者两个以上的线程, 实现二级并行计算。由于任务数大, 要分散操作多次才能完成所有任务。每一次分散操作的总任务数为所有处理机可以同时启动的线程数之和。如5处理机系统, 可以同时启动的线程数为16, 那么一次分散操作只能完成16个任务。这样我们的任务分配可以描述为

J={{j1, j2, …, jL}1, {jL+1, jL+2, …, j2L}2, …,

{jx, jx+L, …, jm (2n+1) }i, m (2n+1) >N} (5)

式中, L为所有处理机可同时启动的线程数。

(3) 分散操作。分散任务并执行, 即将任务队列中的每一组任务按其线程数分别发送给调度算法步骤 (1) 中建立的处理机队列对应的处理机;最后一次分散操作比较特殊, 可能所有处理机的总线程数大于所剩任务数, 如果我们按性能好的处理机优先进行分散操作。如5处理机系统, 其性能数分别为1.0、1.5、2.0、1.0、2.5, 那么他们同时启动的线程数分别为2、3、4、2、5, 问题规模为21个任务, 应该进行两次分散操作, 第一次每一个处理机分配和同时启动线程数相等的任务数, 共分配了16个任务, 第二次5个处理机16个线程只有5个任务, 那么这时就有两种策略:①每一个处理机分配一个任务;②只分配给性能好的部分处理机, 显然这两种策略所花费总时间是不相等的。此问题一般性的描述如下:假设有k个处理机, 按其性能从高到低排序p1, p2, …, pk, 它们分别可以同时启动的线程数为L1, L2, …, Lk, 分别启动一个线程计算所花时间为t1, t2, …, tk, 任务数为J, 当 $J ＜ \sum_{i = 1}^{k} L_{i}$ 时, 问题就转化为

$\begin{array}{l} \min \max {t_{i} x_{i}} \\ s . t . \sum_{i = 1}^{k} x_{i} = J \\ x_{i} \in S_{i}; S_{i} \in S \end{array}} (6)$

其中, xi表示第i个处理机启动的线程数, 且只能取整数; $J ＜ \sum_{i = 1}^{k} L_{i}$ 、ti均为给定的常数。

假设 $f (X) = \max {t_{i} x_{i}}, G (X) = \sum_{i = 1}^{k} x_{i}$ , 那么模型变为

$\begin{array}{l} \min f (X) \\ s . t . G (X) = J \\ X \in S_{i}; S_{i} \in S \end{array}} (7)$

这是等式约束离散变量优化问题。应用方向差商法[12]可求得最优解, 即最优分配策略。

如有8台处理机, 它们可同时启动的线程分别为2、3、4、2、5、6、3、5, 任务是3个设计变量, 3个模型实例的多实例、多目标优化问题。每一次迭代的任务数为21, 而8台计算机可以同时启动的线程总数为30, 假设每个处理机启动一个线程计算一个任务需要的时间分别是5s、4s、3s、5s、2s、1s、4s、2s。那么数学模型为

$\begin{array}{l} \min \max {5 x_{1}, 4 x_{2}, 3 x_{3}, 5 x_{4}, 2 x_{5}, x_{6}, 4 x_{7}, 2 x_{8}} \\ s . t . \sum_{i = 1}^{8} x_{i} = 21 \\ x_{1} \in {0, 1, 2} \\ x_{2} \in {0, 1, 2, 3} \\ x_{3} \in {0, 1, 2, 3, 4} \\ x_{4} \in {0, 1, 2} \\ x_{5} \in {0, 1, 2, 3, 4, 5} \\ x_{6} \in {0, 1, 2, 3, 4, 5, 6} \\ x_{7} \in {0, 1, 2, 3} \\ x_{8} \in {0, 1, 2, 3, 4, 5} \end{array}} (8)$

通过方向差商法求得数据如表1所示。

在可行域中, 找到最小值为8的有8组。

(4) 集中操作, 结束。当都执行完分散的任务时, 将每一个客户处理机的结果收集到服务器。

3 算法分析[9]

SQP每次迭代都构造一个QP问题, 且每次迭代的时间复杂度一样, 所以我们只考虑一次迭代就可以。

顺序执行时间为

ts=m (2n+1) (9)

顺序时间复杂性为O (n m) 。

对于并行优化而言, 设处理机总数为p, 那么从处理机有p-1个。

(1) 通信。首先, 向从处理机广播仿真评估模型, 并返回收到模型确认, 即

tconmm1= (p-1) (tstartup+tdata1) +tstartup+tdata0 (10)

式中, tstartup为形成消息和启动传递所需的时间;tdata0 、tdata1为传递不同数据所需的时间。

再执行分散操作, 即

tconmm2= (p-1) (tstartup+tdata2) (11)

式中, tdata2为分散操作需要的平均时间。

(2) 计算。然后所有处理机并行地完成它们的仿真评估任务, 即

$t_{c o m p} = \frac{2 n m}{p} = Ο (\frac{n m}{p}) (12)$

(3) 通信。执行集中操作, 即

tconmm3=tstartup+tdata3 (13)

式中, tdata3为集中操作需要的平均时间。

(4) 总执行时间。总执行时间可表示为

$t_{p} = \frac{2 n m}{p} + 2 p t_{s t a r t u p} + (p - 1) (t_{d a t a 1} + t_{d a t a 2}) + t_{d a t a 3} + t_{d a t a 0} (14)$

(5) 加速系数。加速系数可表示为

$\begin{array}{l} a_{i n d e x} = \frac{t_{s}}{t_{p}} = \\ \frac{2 n m}{\frac{2 n m}{p} + 2 p t_{s t a r t u p} + (p - 1) (t_{d a t a 1} + t_{d a t a 2}) + t_{d a t a 3} + t_{d a t a 0}} (15) \end{array}$

如果问题规模足够大, 加速系数可能接近p。

(6) 计算/通信比。计算时间与通信时间比可表示为

$c = \frac{2 n m}{p [2 p t_{s t a r t u p} + (p - 1) (t_{d a t a 1} + t_{d a t a 2}) + t_{d a t a 3} + t_{d a t a 0}]} (16)$

(7) 并行效率。并行效率可表示为

$\begin{array}{l} η = \frac{t_{s}}{t_{p} p} = \\ \frac{2 n m}{2 n m + p [2 p t_{s t a r t u p} + (p - 1) (t_{d a t a 1} + t_{d a t a 2}) + t_{d a t a 3} + t_{d a t a 0}]} (17) \end{array}$

如果问题规模足够大, 并行效率可接近100%。

从算法分析中, 我们可以看出通信开销占有相当的比例, tdata0、tdata1、tdata3的大小和问题规模成正比增加。从加速系数可知, 当问题有一定规模时, 通信时间相对减少, 加速系数增大。因此, 并行优化算法是很有价值的。

4 SQP并行优化过程

并行执行流程如图2所示, 图2中, X为设计变量, C为约束函数, O为目标函数, W为加权系数。并行执行流程中涉及如下若干关键过程:

(1) 客户机和服务器转换角色。用户以其中一个MWorks平台为操作对象, 将其转换成 (或增加启动本机作) 服务器, 如果原服务器空闲, 将所有现有的客户端与现在的服务器连接, 包括原服务器上添加一个执行器, 共有执行器数等于客户端个数与原服务器端个数之和。如原服务器忙, 则当前客户机执行单机运算。

(2) 生成仿真模型与参数文件。编译仿真模型, 生成仿真参数mat文件 (MATLAB数据格式文件) , 用于保存所有参数在各仿真步的参数值;将仿真模型生成执行文件Simulator.dll, 初始化时从仿真参数mat文件中读出参数值, 这样每次迭代执行仿真不需要重新生成仿真执行文件。

(3) 客户端创建任务。若角色转换成功, 首先创建并行任务存入队列 (即向量容器vector〈SimulateData〉, 其中SimulateData是结构体) , 如图1所示的工作池的队列, 然后分配任务创建执行器 (即map容器, map〈Socket *, vector〈SimulateData〉〉mExecutors, 其中SimulateData是结构体) , 使客户机和任务对应, 便于管理。

(4) 启动服务器任务分配与负载平衡程序, 进行并行计算, 计算结束通知客户端。

(5) 并行优化任务。并行优化任务基类Task (如C++伪代码) , 其用于存储和设置每一个任务共有的信息。其他任务均从此基类派生。

class Task {//任务

protected:

int m_uType; //命令类型

string m_sClientIP; //客户机名称

int m_nPort; //端口号

bool m_bFinished; //任务完成标志

char * m_sCommandLine; //命令流, 去报头

string m_sAppendMsg; //显示文本信

string m_sModifListStr; //移除的列表信息

static CMySocket * m_pSocket;//套接字

public:

Task () ;//构造函数

～Task () ;//析构函数

int GetType () ;//获得命令类型

string GetClientIP () ;//获得客户机名称

string GetAppendMsg () ;//获得消息

void GetListModStr (StringArray &a) ;//获得移除信息

virtual bool Parse () =0;//解析命令, 将字符流转化成对象内容

virtual string GenerateStream () =0;//将命令对象变成字符流

virtual bool Execute (CPtrList &pSockets) =0;//执行命令

void SetSocket (CMySocket* pSock) ;//设置套接字

void InitInfo (const string* sClient, int nPort, char* pBuf) ;//初始化

void SetInfo (const string* sClient, int nPort, char* pBuf) ;//设置基本数据

……

};

几个关键任务设计, 包括初始化仿真器任务 (InitSimulateDLLTask) , 初始化参数文件任务 (InitParaMatTask) , 执行仿真任务 (SimulateRunTask) , 和仿真结果任务 (SimulateResultTask) , 其中InitSimulateDLLTask用于存储初始化模型动态库文件信息;InitParaMatTask用于存储初始化参数文件信息;SimulateRunTask用于执行仿真任务, 读取初始化参数文件, 仿真评估结束后将结果写入结果数据;SimulateResultTask用于存储结果数据。

(6) 集中操作。服务器收集客户机的仿真结果, 并保存在仿真结果mat文件中。由仿真结果mat文件来计算目标或约束变量值, 对于未激活的不考虑, 变量值需要根据其取值类型 (如最终值等) , 通过仿真时间变化曲线来确定。

(7) 目标聚合。对于多实例、多目标问题, 设q为所有的归一化目标值 $q_{i j} = \frac{c_{i j}}{d_{i j}}, i, j \in S_{\min}$ 组成的向量, 其中, Smin表示所有实例中的最小化的实例目标集合, qi j, ci j, di j分别表示第j个实例的第i个目标的归一值、目标值和期望值。采用向量最大范式, 优化的总目标就是获得一个q中最大值最小的归一向量所对应的调节参数集合的值。

令

α∶=max qi j

这样如果α1<α2则认为第一个设计比第二个设计更满意。如果α≤1, 则该设计是满足约束的, 因为此时各个目标值都小于限定的约束权值。

综上所述, 优化问题可认为是寻找最小的α, 即

α*=min α

由此, 该多实例优化问题的数学模型描述如下:

$\begin{array}{l} \min_{t} \max_{i j \in S_{\min}} \frac{c_{i j} (t)}{d_{i j}} \\ s . t . c_{i j} (t) \leq d_{i j} i, j \in S_{i n e q u a l} \\ c_{i j} (t) = d_{i j} i, j \in S_{e q u a l} \\ t_{\min} \leq t \leq t_{\max} \end{array}} (18)$

其中, t为时间, tmin和tmax分别为仿真时间上下限, Sinequal为不等式实例约束, Sequal为等式实例约束。

5 实例分析

下面以MWorks系统建立的F14战机简易模型的控制系统参数优化为例, 说明多实例多目标问题仿真优化并行求解。该模型用于对战机模型的纵向运动响应进行仿真与分析, 优化设计的目的是确定合适控制参数, 使该响应符合操控要求。MWorks提供了建立优化模型的统一界面。用户只需从编译仿真模型后得到的参数变量结构树中选择并设置优化模型信息, 亦可读取保存在仿真模型优化信息, 在此基础上构建新的优化模型。

5.1 设计参数

选择3个控制器增益作为设计参数, 并设置初始值和边界值, 即Kf∈[-10, -0.5], 取初始值为-2;Ki∈[-10, -0.5], 取初始值为-2;Kq∈[0.1, 10]取初始值为0.5。其中, Kf为攻击角的比例增益;Ki为攻击角的积分增益;Kq为俯仰速率的比例增益。

5.2 模型实例

选取战机模型的5个空气动力学参数作为实例参数, 并定义了如表2的两个模型实例, 其中实例2的实例参数值是实例1参数值变动+10%后得到的。该多实例优化的目的是在战机模型参数扰动的情况下仍能获得满意的优化结果。模型实例见表2。表中, Ma为俯仰角对攻击角的变化系数;Md为俯仰扭矩对单位攻击角偏差的变化系数;Mq为俯仰扭矩对俯仰速率的变化系数;Za为单位攻击角的提升系数;Zd为单位攻击角偏差的提升系数。

5.3 实例目标和约束

选择攻击角的超调量f1、上升时间f2和调节时间f3作为实例目标, 以攻击角最大偏差c为约束, 并设定各实例期望值见表3。

5.4 优化模型数学表达式

建立以上优化模型后可以得到如下优化模型的数学表达式:

$\begin{array}{l} f i n d X = [Κ_{f}, Κ_{i}, Κ_{q}] \\ \min \sum | \frac{f_{i j}}{d_{i j}} | i = {1, 2, 3}, j = {1, 2} \\ s . t . c_{1 j} \leq 2 \\ - 10 \leq Κ_{f} \leq - 0.5 \end{array}} (19)$

5.5 计算机性能参数

利用2节点机群系统求解。单机及机群系统的机器配置为:Pentium (R) D CPU 3.0GHz, 1.00GB RAM。网络环境为100M局域网。

分别采用序列二次规划法串行和并行求解问题式 (19) , 经过15次迭代后收敛, 优化结果为X=[Kf, Ki, Kq]=[-3.0432, -3.3956, 0.7723]。优化求解后得到的目标和约束值见表4。从表4可以看出所有的目标和约束均满足表3设定的实例期望。单机耗时339s, 机群耗时249s, 可以看出2节点机群计算时间有明显下降。我们以3、4、5节点机群系统分别求解, 得出如图3所示, 计算节点增多, 计算时间缩短, 加速比增大, 但是并行效率有所下降。

在两节点机群系统中, 求解F14战机算例的加速比和并行效率。

加速系数为

$a_{i n d e x} = \frac{t_{s}}{t_{p}} = \frac{339}{249} = 1.36$

并行效率为

$η = \frac{t_{s}}{t_{p} p} = \frac{339}{249 \times 2} = 68.07 %$

根据Amdahl定律, 加速极限为

$a_{\lim i t} = \frac{1}{\frac{t_{s} - t_{p}}{t_{s}}} = \frac{1}{\frac{339 - 249}{339}} = 3.77$

说明加速比的极限远没有达到, 可以通过增加处理器的数量来使加速比进一步提高。仅单项仿真评估任务并行处理和梯度计算并行化, 就可得到比较好的效率。若提高机群性能, 网络环境等硬件设备, 或者改进算法, 将进一步提高计算效率。

6 结论和展望

本文从SQP优化算法的并行化以及实现等方面研究了SQP并行参数优化设计方法。研究表明SQP方法中梯度计算和模型仿真评估并行化效果较好, 明显加速了优化设计的进度, 在基于多学科仿真的多实例多目标优化中体现的更加明显;SQP方法成功与否依赖于梯度计算的准确性和有效性;利用机群系统计算可用较低的成本获取较大的计算能力, 解决仿真优化中高耗时的问题。系统可以很容易地引入那些可买到的最新的处理器, 扩充机群系统。SQP并行算法的进一步研究, 应考虑全局并行化和线性搜索并行化, 以及全局、函数评估、梯度计算与线性搜索混合并行化和伴随梯度法[5,6]。

参考文献

[1]吴义忠, 吴民峰, 陈立平.基于Modelica语言的复杂机械系统统一建模平台研究[J].中国机械工程, 2006, 17 (22) :2391-2396.

[2]赵建军, 丁建完, 周凡利, 等.Modelica语言及其多领域统一建模与仿真机理[J].系统仿真学报, 2006, 18 (2) :570-573.

[3]丁建完, 陈立平, 周凡利.复杂陈述式仿真模型的指标分析[J].系统仿真学报, 2006, 18 (8) :2092-2096.

[4]余俊, 周济.优化方法程序库OPB-2—原理及应用[M].武汉:华中理工大学出版社, 1997.

[5]High K A, LaRoche R D.Parallel NonlinearOpti mization Techniques for Chemical ProcessDesign Problems[J].Computers&ChemicalEngineering, 1995, 19 (6-7) :807-825.

[6]张帆, 邵之江, 仲卫涛.化工过程系统优化的分布式并行计算[J].化工学报, 2001, 52 (5) :396-400.

[7]陈忠.关于非线性规划问题的并行算法[J].长江大学学报 (自科) , 2006, 4 (12) :1-8.

[8]Ling J, Lorenz T, Biegler V, et al.Design andOpti mization of Pressure Swing Adsorption Systemswith Parallel I mplementation[J].Computers andChemical Engineering, 2005, 29:393-399.

[9]Wilkinson B, Allen M.并行程序设计[M].陈鑫达, 译.北京:机械工业出版社, 2005.

[10]黄金贵, 陈建二, 陈松乔.并行环境下基于多处理机任务的调度模型与调度算法[J].计算机科学, 2002, 29 (4) :1-3.

[11]Hochbaum D S.Approxi mation Algorithms for NP-hard Problems[M].London:PWS PublishingCompany, 1997.

并行化仿真篇4

随着新能源和电力电子技术的发展和广泛应用, 以及用户对电能质量及供电可靠性等要求的不断提高, 现有交流配电网将面临分布式新能源接入、负荷及用电需求多样化, 以及高供电可靠性等方面的巨大挑战。国内外研究资料表明, 基于直流的配电网在输送容量、经济性、分布式能源接入以及提高供电质量等方面具有比交流更好的性能, 可以有效地解决现有交流配网面临的一些问题[1]。

直流配电系统是包含中压配电网和用户侧配电网的公共配电网络, 为了保证直流配电系统的安全、稳定运行, 发展直流配电网仿真技术, 研究其复杂的暂态过程和动态响应就变得非常重要。

在直流配电网中, 各器件都反应迅速, 需要精细仿真, 因此对仿真步长要求很高;同时, 在网络中控制系统逻辑复杂, 需要采用仿真步长更短、基于元件线性化的电磁暂态仿真[2,3]。

在直流配电网中, 主要有两类元件:一次系统中的各类电气元件和二次系统中的各类控制元件。在进行电磁暂态仿真的时候需要分别考虑, 各自求解。

电气网络的电磁暂态仿真利用了仿真步长短的特点, 将单个元件进行线性化, 再对线性化后的系统列写KVL方程, 通过求解线性方程YU=I来求解节点电压[4,5,6]。

电磁暂态仿真中, 控制系统仿真算法有两类:第一种是将各控制元件单独建模, 根据信号流动方向, 依次对每个元件进行仿真;第二种是将各元件离散化, 建立整个控制网络的网络方程求解[7]。第一种单元件方法系统所有元件串行处理, 无法进行并行计算, 速度很难得到提高[8]。而第二种网络求解方法实现了系统的网络化求解;但是因为非线性元件无法通过常规方法离散化, 必须单独计算, 而且随着网络规模的扩大, 计算量大规模增加, 需要采取并行计算方法, 将系统降阶, 提高运算速度[9]。

为了提高电力系统电磁暂态仿真的效率, 很多学者对分网并行进行了研究。目前分网并行算法主要有3种方式[10,11,12,13]:长距离输电线路解耦法, MATE (MultiArea Thevenin Equivalent) 方法, 节点分裂法。但是这3种算法都是针对电力系统一次侧电气网络而提出的, 并不适用于二次侧控制网络仿真。

针对上述问题, 本研究在传统的控制网络电磁暂态仿真算法 (Transient Analysis of Control Systems, TACS) 的基础上, 提出控制系统的时延环节解耦并行仿真算法, 使得在直流配电网仿真中, 控制系统并行仿真得到实现;在本研究中, 电气网络和控制网络二者分别采用MATE算法和时延环节解耦算法进行并行计算。

1 直流配电网络元件建模

目前, 国内外对直流配网系统的研究还处于起步阶段, 现有文献中也没有见到包含各种分布式电源、各类型负荷的完整直流配网系统的仿真模型。因此, 为了深入研究直流配网的特性, 有必要建立与外部交流电网相连、包含分布式电源、储能、各类型负荷的直流配网系统的仿真模型。

目前已在电磁暂态仿真程序中初步建立了电压源换流器 (VSC) 、直流变压器、锂电池储能单元、光伏电池等系统元件模型。

1.1 电压源换流器 (VSC) 模型

电压源型换流器 (VSC) 的主电路拓扑结构如图1所示, 由6个全控型的功率开关管 (IGBT) 构成, 每个功率开关管都反并联一个续流二级管。

R—变流器功率损耗的等效电阻, L—变流器等效电感, C—直流侧滤波电容

本研究对电压源换流器的控制采用直接电流控制, 具有双闭环的控制结构:外环控制根据变流器待实现的控制功能, 如定功率控制或定电压控制等, 确定变流器电网侧输出电流分量的目标值;内环控制用于实现变流器输出调制电压的控制。控制过程较为复杂, 包括:ABC/DQ变换, 双闭环控制, DQ/ABC变换, SP-WM控制[14,15,16]。

1.2 直流变压器模型

双向直流变压器的拓扑结构如图2所示[17]。直流变压器的原边和副边都由4个可控的开关管和4个与可控开关管反并联的二极管构成, 原边和副边通过一个高频变压器相连。开关管S1和S2、S3和S4、S5和S6、S7和S8同时开关, 每个开关管以近似0.5的恒定占空比工作。

1.3 锂电池储能单元模型

单个锂电池的仿真模型可以根据锂电池输出电压表达式用一个可控的电压源和一个电阻相串联来建立。

锂电池的输出电压表达式为:

式中:Vbatt—锂电池出口电压;i—锂电池出口电流;E0—锂电池开路电压;R—锂电池内阻;K—极化常数;Q—锂电池容量;it—锂电池实际电荷量;A, B—锂电池特性曲线中指数区电压幅度和时间常数的逆。

1.4 光伏电池模型

光伏电池是利用半导体材料PN结的电子特性把太阳光能直接转换成电能的一种固态器件。光伏电池的等效电路如图3所示, 由文献[18]可得到光伏电池的模型。

光伏电池的输出电流表达式如下式所示:

式中:Iph—单晶硅的光电流强度;Ios—电池反向饱和电流;Rs—光伏电池的等效串联电阻;Rsh—光伏电池的等效并联电阻;q—电子电荷, 且q=1.602×10-19C;a—二极管特性因子, 通常取值为1~1.5;k—玻耳兹曼常数。

1.5 交、直流负荷模型

交流负荷模型采用的是有刀闸的定功率负荷, 示意图如图4所示。

开关BRK在0时刻的闭合的, 此时交流负荷输出功率为7 MW;在2 s时刻断开, 交流负荷输出功率变为1 MW。

直流负荷被简化为1Ω的电阻。

分析以上模型可以看出, 这些模型的电气网络并不复杂, 但二次侧控制网络却都比较复杂。例如VSC整流器中, 电气网络有节点12个, 但控制6个电力电子开关的控制网络却很复杂, 包括32个节点, 27个线性的控制元件和一个较复杂的非线性元件, 即SPWM (Sinusoidal Pulse Width Modulation) 脉宽调制元件。

2 直流配电网控制系统并行仿真算法

本研究采用网络化求解方法进行控制系统求解, 并提出了控制系统的时延解耦仿真算法。

2.1 控制元件电磁暂态建模方法

对于有可传递函数的函数块, X (s) =G (s) U (s) , 都能够利用和电气网络类似的梯形差分法来进行差分化。

令:

得到X (s) ·D (s) =K·N (s) ·U (s) , 转化成时域的方程为:

经过线性化后可以得到:

在完成每一时步的求解后, 需要更新n个历史项, 以得到下一时步的H。因为计算比较复杂, 一般采用递归方法。

2.2 控制系统电磁暂态分析算法

将元件进行离散化后, 每个元件能够得到的一阶线性方程如下式所示:

式中:c, K, d—常数;hist (t) —历史项[19,20]。

包含4个传递函数以及两个求和函数的系统如图5所示, 其中, 系统输入信号u=[u1u2], 输出信号xout=[x4]。

下面以图5的算例来说明, 当元件离散化完成之后, 是如何进行系统的网络求解的。

根据离散化后的元件格式, 可以得到4个方程组成的方程组如下:

可简化成:

在进行下一步计算之前, 首先要对[Axx]和[Axu]进行三角分解, 然后按照如下的步骤求解未知数[x]:

(1) 计算右端项[hist];

(2) 从上向下作前代运算;

(3) 回带求出[x];

(4) 更新每一块的历史项。

2.3 控制系统电磁暂态仿真并行算法

2.2节中介绍了单个控制系统网络电磁暂态仿真的分析算法。在实际的电力系统中, 由于电气元件的隔离, 会出现各个子控制网络相互解耦的情况, 如果采用传统的分析方法, 公式 (7) 中的系数矩阵规模将会很大, 不利于计算。另外在控制网络中, 存在较多的时延环节, 当仿真步长小于等于时延时, 时延环节两端的网络是解耦的。为了减小公式 (7) 中系数矩阵的维数, 提高整个仿真程序的并行效率, 本研究提出了控制网络的并行仿真算法。

控制网络分网示意图如图6所示。子网络ABC为2.2中计算的3个线性控制子网络, 根据上述算法, 对于3个子网络分别有:

各参数含义如图6和2.2节说明所示。

可得整体控制网络的方程如下所示:

式中:PAB—相对应输入/输出的关联矩阵;例如PAB—子网络B的输入节点与子网络A的内部节点的关联矩阵。

当各个子网络之间无相互关联时, 各个P矩阵都为0, 很容易可以看出, 这个时候各个子网络相互解耦, 每个时步仿真只需直接分网计算各个子网络内部数据即可。

当各个子网络之间以时延环节相互联系时, 公式 (9) 会发生如下变化:

式中:xAhist, xBhist, xChist—时延时刻之前的xA, xB, xC值。

在T时刻时是已知的, 与当前状态下需要求的xA, xB, xC相互独立。因此可得调整后的方程如下, 可得此时3个子网络仍然相互解耦。计算时首先根据子网络之间的数据交换来计算方程右侧历史项, 然后各个子网络分网并行计算。

3 算例分析

本研究采用的测试算例系统为如图7所示的双端直流多电压等级配网系统, 包含光伏电池、锂电池储能系统和交直流负载。系统两端都与无限大电网相连, 左端VSC整流器实现的是定电压控制, 使得整流器输出电压稳定为15 k V, 右端VSC整流器进行定功率控制, 使得整流器交流侧稳定输出功率为10 MW。保证在负荷、光伏、锂电池出力变化的时候, 系统都能够稳定运行, 并且直流母线上各个节点的电压都保持在15 k V左右。直流母线上通过直流变压器连接光伏电池和锂电池向系统输送电能, 以及通过逆变器和直流变压器分别于交、直流负载相连接。系统中所有的元件都基于EMTP程序进行建模, 模型如第1小节所述。为了保证仿真的准确性, 仿真步长为5μs, 仿真总时长为1 s。

如图7所示, 上述算例共包括3个VSC整流器, 3个DC/DC直流变换器, 光伏电池, 锂电池, 直流负载和交流负载。在自然情况下, 上文中已经提到, 在这些模块中, VSC整流器电气网络和控制系统都较其他模块为复杂, 尤其是控制系统, 包括32个节点, 27个线性的控制元件和一个较复杂的非线性元件, 即SPWM脉宽调制元件。因为非线性元件不能按照一般的离散化方法来建模, 需要单独计算。如果将整个VSC模块作为一个子分区, 那么分区内部, 线性部分和SPWM的仿真仍是串行的, 极大降低了仿真效率, 也导致了系统中分区的不平衡, 降低了系统的并行度。因此要将VSC模块拆分成线性部分和SPWM两个子分区。在VSC中, 为了模拟真实的系统运行状况, 在SPWM之前, 本研究增加了一个时步的时延。利用这个时延, 就可以将线性部分和SPWM进行解耦, 作为两个分区并行计算。

3.1 程序准确性验证

首先验证仿真程序的正确性。笔者选取代表性的线路选取直流线路节点电压Udc1, 将电磁暂态仿真串、并行程序的仿真数据与PSCAD的计算结果进行比较。串、并行程序对比图如图8所示, 线路电压由0上升, 经过0.25 s左右的时间, 达到稳定点15 k V, 并保持稳定。在图8中3条曲线是完全重合的, 而PSCAD与并行仿真程序的误差不超过10-8, PSCAD与并行程序误差分析如图9所示。该数据表明电磁暂态仿真串、并行程序的仿真结果是正确的。

3.2 并行仿真程序效率分析

为了检验时延解耦并行仿真算法能否对程序的效率有提高作用, 本研究采用了3种不同的分区方案, 对系统进行仿真。以下3个方案都采用10个计算机核心进行仿真:

方案1:不采用时延解耦并行仿真算法, 电气网络仿真采用传统的MATE算法, 每一个模块作为一个子分区, 共分为10个小的电气网络子分区, 所有模块的控制网络作为一个大的分区进行串行计算;

方案2:采用延时解耦并行仿真算法, 每一个模块作为一个子分区, 共分为10个小的电气网络子分区, 以及10个控制网络子分区;为了最大效率地提高并行度;

方案3:采用延时解耦并行仿真算法, 并对VSC整流器模块进行解耦拆分, 拆成VSC线性部分和SPWM两个子分区。这样共分为10个小的电气网络子分区, 以及13个控制网络子分区包括3个VSC线性部分和3个SPWM子分区, 相应地采用13个核心进行并行计算。

仿真时长为1 s时, 各个方案的仿真时间如表1所示。

从以上结果可以得知, 与不采用时延解耦算法的方案1相比, 采用了并行计算的方案2和方案3都在保证程序的正确性的基础上, 提升了程序的计算效率;而方案3中, 对VSC控制器的解耦使得各个分区的规模大小变得更为相似, 减弱了系统各个分区规模大小的不平衡程度, 进一步提升了计算效率。

为了进一步研究控制系统并行仿真策略对控制系统仿真的影响, 单步仿真的耗时统计如表2所示。

由表2可以看出, 应用延时解耦并行算法之后, 控制网络的单步耗时减少了58.14%, 效率提升了一倍以上, 充分证明了本研究提出的算法是有效提高了仿真效率。算例中, 控制系统包含13个分区, 采用了13个计算核心对系统进行并行计算, 但耗时的减少并未达到相应的大规模减小。造成该问题的主要原因可能是在实际的电力系统中, 并行计算分区很难实现理想的平衡状态。时延解耦并行算法中, 只能在子分区之间没有直接相连或者通过时延环节相连时, 才能进行分区间的解耦, 分网方式还不够灵活, 需要进一步研究。

4 结束语

本研究开发了直流配电网电磁暂态仿真程序, 另外搭建了光伏电池、锂电池、VSC变流器等系统的电磁暂态仿真模型;针对直流配电网的特点, 提出了适用于控制网络的时延解耦分网并行仿真算法, 并在仿真程序中得到了实现, 验证了算法的正确性, 并且利用该算法有效提高了直流配电网的仿真效率。

本研究搭建了包含光伏电池和储能的双端直流配电网算例, 仿真结果证明该程序能够进行一定规模的直流配电网电磁暂态仿真。而本研究提出的网络并行计算方法也是正确而高效的。

本文引用格式:

张达, 江全元.直流配电网电磁暂态仿真并行算法研究[J].机电工程, 2014, 31 (9) :1185-1190.

动态碰撞网格的DSMC及其并行化篇5

1.1 DSMC算法的背景

通常气体动力学中对气体的处理采用连续介质假设, 但当气体密度很低时, 气体的粒子效应会变得十分显著, 从而破坏连续介质假设。这种情况下需要采用稀薄气体动力学的方法才能得到正确的结果。在导弹、飞船和航天器的设计中, 稀薄气体动力学有着重要的应用。

稀薄气体的运动规律由波尔兹曼方程描述, 但这是一个带积分项的偏微分方程, 除了某些特殊情况, 极难得到解析解。并且由于带有积分项, 因此偏微分方程数值算法也不能很好地适用, 这使通过数学模型方法来研究稀薄气体问题陷入困境。

DSMC方法是由Bird[1,2,3]提出。该方法没有从波尔兹曼方程出发, 而是直接通过模拟分子的运动和分子间相互碰撞, 得到气体的宏观情况。与传统的模型方法相比, 它简化了模型, 易于引入各种物理和化学模型, 且得到了实验数据的很好验证。现在DSMC已经成为了求解稀薄气体问题的一种全能算法。

1.2 DSMC算法简介

DSMC方法用大量模拟分子来模拟真实气体中分子的运动过程。由于气体分子运动的确定性计算量非常大, DSMC通过对分子运动和分子间碰撞进行解耦, 且在计算中采用了Monte Carlo方法, 减小了计算量, 使大数据量的分子模拟成为可能。

图1为DSMC的算法流程图, 其中的具体步骤说明如下:

1) 初始设置根据流场的初始条件 (温度、密度、初始速度等) , 以及模拟分子的条件 (模拟比例、能量模型等) , 设置模拟分子的位置、速度、能量。并根据流场和固壁的条件划分网格, 使网格中的分子数量及网格的大小在一个合理的范围之内。

2) 分子移动计算分子在其当前速度下, 经过一个时间t之后的位置, 分子的移动采用确定性的算法。同时, 在这一步要考虑分子与固壁的碰撞和能量交换, 分子与固壁的碰撞计算采用Monte Carlo算法, 保证碰撞过程前后动量和能量的期望守恒。

3) 分子标定根据分子移动后的新位置, 标记其所属的网格。

4) 分子碰撞对流场中的每一个网格, 根据网格的信息计算碰撞次数N。然后随机选取网格中的N对模拟分子, 按照碰撞模型重新计算分子对的速度和能量分布。

5) 判断结束重复以上2) -4) 步直到满足结束条件。

6) 输出宏观量采样网格中的分子信息, 计算得到整个流场的宏观性质。

1.3 DSMC算法的不足

尽管DSMC算法的优点很多, 但其缺点也是很明显的, 主要有两点:

1) 网格处理复杂 DSMC方法中的网格是固定的。为了保证计算结果的准确性, 对网格的大小及网格中的分子数都有一定的要求。但实际中流场的情况千变万化, 实际计算中需要根据算例的情况采用复杂的网格生成算法甚至是手工生成网格。在流场发生变化时, 还需要不断调整并重新划分网格, 十分不方便。

2) 计算量大除了网格问题, 实际使用中DSMC方法最大的问题就是计算量大、计算时间长, 以至于到了现有单机无法承受的地步。一个典型的DSMC算例需要计算几万个时间步长, 计算几十万的网格和几百万的模拟分子, 这样的计算量在单个CPU上所花费的计算时间是无法承受的。

1.4 本文对DSMC算法的改进

针对DSMC方法存在的不足, 本文在两方面对DSMC算法进行了改进。 (1) 将原先DSMC算法中的固定碰撞网格部分改为动态生成网格, 使DSMC算法可以做到对复杂流场和流场变化自适应。 (2) 对动态网格的DSMC进行了共享内存模型的并行化, 提高了其在单机上的计算效率。

2 DSMC网格和并行化相关工作

根据网格的不同, 现有的DSMC方法大致可以分为两类, 一类是非结构的适体网格程序, 另一类是以位置元方案为主的直角结构网格[3]。这两种方法各有优缺点。非结构网格可以适应复杂的流场条件, 但是追踪分子的算法复杂且效率低下。位置元方法是直角坐标网格的进一步改进, 但存在分子与表面元的碰撞判断复杂低率的问题, 且影响了DSMC的精度和准确性。

在DSMC的动态网格方面, 就笔者所知, 仅有Olson等人研究了无网格 (动态网格) 的DSMC方法[7]。该文使用八叉树 (Octree) 作为动态空间划分的结构, 但其对网格的处理较为繁复, 也未对无网格算法的性能进行比较实验。

在DSMC算法的并行方面, 由于实际应用中对DSMC算法的需要以及DSMC本身所要求的大量计算, 对DSMC算法的并行化研究很早就展开了, 如文献[4,5,6]。随着计算机硬件的发展, 多核CPU成为发展的趋势。如何使DSMC算法充分利用现有多核架构的计算能力已经成为新的研究课题。

3 DSMC的动态网格算法

3.1 DSMC中网格的分析

DSMC方法的重要实质是将分子的运动与分子间的碰撞解耦。而碰撞网格的作用就是划分模拟分子。根据文献[1], 理想的DSMC方法中的网格需要满足如下条件:

1) 很好的计算效率;

2) 能适合复杂的流场和固壁条件;

3) 网格可以自适应流场局部的梯度及分子自由程;

4) 当流场发生变化时, 网格能自适应地调整;

5) 可以有效地降低碰撞分子对的平均距离;

6) 可以高效地确定分子所属的网格。

传统的DSMC算法中使用固定的碰撞网格, 显然难以满足以上要求。而对固定网格算法的改进也集中在对已有网格进行局部修改, 当流场变化剧烈时无法适用。

3.2 DSMC的动态网格算法

通过对DSMC算法的分析可以发现, 碰撞网格和采样网格可以分开, 且每一轮的计算中只涉及到碰撞网格。碰撞网格的作用是为了有效地计算分子碰撞数以及选取分子碰撞对, 因此碰撞网格从概念上可以作为碰撞过程中的一部分。在DSMC算法的每一轮的分子碰撞过程中, 先根据分子的分布情况动态构建网格, 标记分子所属网格, 再遍历所有的动态网格计算分子碰撞, 即得到动态网格算法。整个算法流程图如图2所示。

动态划分算法采用kd-tree算法。具体的划分算法如下:

算法1 划分网格算法

输入: 网格根节root, 分子数组moles。

输出: kd-tree, 叶子节点为碰撞网格及其所属分子。

步骤:

1) 如果moles的分布和数量以及root的属性满足构成一个碰撞网格的条件, 则将root的左右子节点设为空, 设置root中指向moles的数据, 返回root。

2) 根据空间划分策略, 计算空间划分的轴和划分位置。

3) 根据第2步的结果, 将moles分成两部分:moles-left和moles-right。

4) 构造两个新的kd-tree节点node-left和node-right, 并设置相应属性。

5) 在node-left、node-left和node-right、 node-right上递归调用划分网格算法, 并设置node的左右子节点为node-left和node-right, 返回root。

该算法中仍有两点需要确定, (1) 如何决定输入的节点和模拟分子是否可以构成一个碰撞网格, 即终止条件。 (2) 如何有效地进行空间划分。本文采用的终止条件为分子数小于30个分子;空间划分策略为按网格最大维度的几何中心为划分平面。

根据分子动力学, 单个网格在Δtm时间进行分子碰撞数为:

$Ν_{t} = \frac{1}{2} Ν_{m} n \bar{σ}_{Τ} \bar{g} Δ t_{m}$

其中Nm是网格中的分子数, $\bar{σ}$ T是分子平均碰撞截面积, $\bar{g}$ 是分子平均相对速度。但计算Nt 需要计算每对分子的相对速度, 计算量过大。为了解决这个问题, Bird提出了非时间计数器方法[2], 即先从网格中选择NNTC对候选分子, 再按概率 $(\bar{σ}_{Τ} \bar{g}) / (\bar{σ}_{Τ} \bar{g})_{\max}$ 决定该对分子是否碰撞。NNTC按下式计算:

$Ν_{Ν Τ C} = \frac{Ν \bar{Ν}}{2} F_{Ν} (\bar{σ}_{Τ} \bar{g})_{\max} Δ t_{m}$

其中 $\bar{Ν}$ 是碰撞网格的平均分子数, FN是模拟分子所代表的真实分子数。但是在每轮碰撞动态划分网格的情况下, 无法得到有意义的 $\bar{Ν}$ 值。 Bird[8]证明了, 在碰撞网格中的分子数量N满足Poisson分布的条件下, 可以用N (N-1) /2代替上式中的 $Ν \bar{Ν} / 2$ , 不会影响正确的分子碰撞数计算。实验结果也验证了这种处理的正确性。

3.3 算法的复杂度分析

根据上述的动态网格划分算法, 设该算法的复杂度为T (N) , N为总的分子数, 则复杂度的递推式为:

T (N) =P1 (N) +P2 (N) +T (N1) +T (N2) N1+N2=N

其中P1 (N) 为计算划分维度和划分位置的开销, P2 (N) 为划分模拟分子的开销, T (N1) 和T (N2) 则为分别递归划分两部分分子的开销。根据算法, 采用网格最大维度上的几何中心作为划分位置, 则P1 (N) =O (1) , 线性划分分子的开销为P2 (N) =O (N) , 则有:

T (N) =O (N) +T (N1) +T (N-N1)

上式的解为T (N) =O (Nlog2N) (平均情况下) , 这个复杂度是可以接受的。同时注意到, 在简单直角网格的情况下, 标定分子的时间复杂度为O (N) , 因此动态划分网格的复杂度显然大于固定网格算法。但通过良好的程序实现, 可以使得动态划分网格的DSMC算法的速度在可接受的范围之内。

4 DSMC的多核并行

DSMC算法是一个非常适用于共享内存并行的算法。在分子移动和分子碰撞的计算中, 不同的分子和不同的网格之间不存在相关性, 因此可以简单地分别对分子和网格的规模进行分解, 然后分别计算, 达到很好的并行效果。具体算法这里略去。

在动态网格划分的算法中, 由于划分模拟分子需要对内存进行写入, 因此针对每一趟划分进行并行比较困难, 但显然一旦划分完毕, 则划分得到的两部分分子可以分别进行再划分, 互相之间不存在相关性。以一百万分子、四核并行为例, 设t为单步划分所需时间, 易得串行程序的运行时间为20t, 而并行程序理论上的运行时间为6t, 加速比3.33, 并行效率为83%, 并行效率是比较高的。

5 实验结果及讨论

实验的算例为超音速横掠平板流动问题[1], 具体设置如下:

超音速氮气流过一热平板, 计算区域尺寸L=1m, H=0.6m, 底部的平板长度L=0.9m, 来流气体速度为4.0Ma, 温度300K, 平板的温度为500K。划分采样网格数为100×60, 模拟比例为0.175×1016, 流场稳定时模拟分子数约为3.7×104, 每一轮计算时间约为4.0×10-6秒, 迭代计算10000次。

5.1 DSMC动态网格算法的正确性验证

将DSMC动态网格算法与原有DSMC固定网格算法所得到的结果进行对比如图3-图8所示。可以看出, 动态网格的实验结果与固定网格的实验结果吻合得很好。

5.2 DSMC动态网格算法的单机性能

将动态网格的DSMC算法 (C++) 与Bird的DSMC程序 (FORTRAN) 进行对比。在Intel Core平台上进行实验结果如表1所示。

从表1可以看出, 静态网格比动态网格要快75%左右, 这一差异是可以接受的。通过对动态网格程序实现进行优化 (更好的数学库, 更好的随机数生成器以及优化内存存取) , 这一差异还可以进一步缩小。

5.3 DSMC动态网格算法的并行性能

在Windows平台上使用Windows Thread API进行DSMC动态网格算法的并行化, 在Intel Quad Core平台上进行实验, 得到的运行时间以及DSMC算法的各步骤运行时间如表2所示。

从表2可以看出, 网格构建和分子碰撞这两个过程的加速性能明显好于分子移动。这是因为分子移动的计算量较小, 此时瓶颈在于内存带宽, 并行加速有限。另外, 通过对该程序进行性能剖析表明, 有大量CPU运算在随机数上, 而要实现线程安全的随机数生成器, 需要频繁地进行同步, 影响了性能。这就是理论上加速比较低的网格构建反而有较高的加速比的原因 (网格构建不需要随机数) 。通过换用更好的线程安全的随机数生成器可使DSMC算法的并行性能可以得到进一步的提高。

6 结论及未来工作

本文分析了传统的固定网格算法的DSMC方法的特点和不足, 提出了一种新的动态构建网格的DSMC算法, 改进了固定网格算法无法应付复杂流场的缺点。本文通过与固定网格算法的对比实验, 进一步验证了该算法的有效性和正确性, 通过实验肯定了动态网格的速度性能。最后, 为了进一步提高该算法的性能, 本文对该算法在共享内存模型下进行了并行化, 得到了比较满意的结果。

尽管多核CPU可以有效地提高DSMC算法的计算效率, 但是实际的DSMC算例的计算量还是大大超过了单机可以承受的范围。通过多机并行来提高DSMC算法的效率已是必然要求。但多机并行DSMC算法面临着通信量大、通信频率高和负载同步等一系列问题。本文提出的动态网格的DSMC算法, 在多机并行的环境下可以省去同步网格的开销, 有效地降低了多机并行中的通信量, 并且降低了负载平衡的难度。相信动态网格的DSMC算法可以为多机并行DSMC效率的提高提供新的思路。

摘要：DSMC (Direct Simulate Monte Carlo) 方法是处理稀薄气体问题的一种有效的方法, 但现有的DSMC存在着网格生成及处理复杂、算例需进行大量人工调整、计算量大、耗时长等缺点。分析DSMC方法中对网格的要求以及网格在整个DSMC方法中所起的作用, 提出了动态划分碰撞网格的DSMC算法, 有效地解决了复杂流场条件下网格自适应的问题, 并通过实验验证了该算法的正确性。同时, 针对DSMC算法计算量大的特点, 利用共享内存的并行模型对动态网格的DSMC算法进行了并行化, 得到了较好的结果。

关键词：DSMC,动态网格,并行化

参考文献

[1]BIRD G A.Molecular Gas Dynamics[M].Oxford:Clarendon Press, 1976.

[2]BIRD G A.Molecular Gas Dynamics and the direct simulation of gasflow[M].Oxford:Clarendon Press, 1994.

[3]BIRD G A.Application of the DSMC method to the full shuttle geome-try[R].AIAA-90-1692.

[4]LAUX M.Optimization and parallelization of the DSMC method on un-structured grids[R].AIAA-97-2512.

[5]KIMMG, KIMHS.A parallel cell based DSMC method with dynamicload balancing using unstructured Adaptive meshes[R].AIAA2003-1033.

[6]姜恺, 黄良大.直接模拟蒙特卡洛 (DSMC) 方法的并行化[J].高性能计算发展与应用, 2005 (2) .

[7]Olson S E, Christlieb AJ.Gridless DSMC[J].Journal of Computation-al Physics, 2008, 227 (17) .

基于聚类算法的并行化研究篇6

由于数据挖掘是从海量数据中提取有用信息,处理效率问题成了对海量数据处理的瓶颈之一,传统的单机串行算法效率较低;由于部分聚类算法中蕴涵并行性,所以为了解决处理效率问题,将并行化的程序设计思想(并行处理)引入聚类算法,同时降低算法的复杂度,使用机群系统进行并行计算,从而有效的缩短聚类的时间。

1 K-means算法

1.1传统K-means聚类算法

K-means算法以k为输入参数,把包含n个对象的集合分为k个簇,使得结果簇内的相似度高,而簇间的相似度低。簇的相似度是关于簇中对象的均值度量,可以看做簇的质心或重心[2]。

传统K-means算法的处理流程如下:

输入:k:簇的数目

D:包含n个对象的数据集

输出:k个簇的集合

方法:

1)从D中任意选择k个对象作为初始簇重心

2)Do

3)根据簇中对象的均值,将每个对象(再)指派到最相似的簇

4)更新簇均值,即计算每个簇中对象的均值

5)while数据集中所有对象的平方误差和E不再发生变化

通常,采用平方误差准则,其定义如下:

其中,E是数据集中所有对象的平方误差和,p是空间中的点,即给定对象,mi是簇Ci的均值(p和mi都是多维的)。换言之,对于每个簇中的每个对象,求对象到簇中心距离的平方再求和。这个准则试图使得生成的k个结果簇尽可能的紧凑和独立。

1.2并行化K-means改进算法

随着并行处理技术的快速发展,越来越多的研究人员尝试将并行处理方法应用于提高聚类算法的效率,通过研究发现K-means算法具有很大的并行性。首先,可将待挖掘的数据集N划分为t个数据子集,t为并行处理环境中处理机的数目;然后将划分后t个数据子集分别发送到t台处理机进行数据聚类处理;最后主机将收到的节点机的聚类结果计算平方误差准则函数E的值,并将前后两次结果做差,如果差的绝对值小于阈值10-6,则处理结束,否则继续循环处理。并行K-means算法的流程如图1所示。

1.3实验结果与分析

我们搭建工作站机群系统,通过以太网卡等连接5台PC机(Intel P4.17GHz、256MB RAM,安装LINUX redhat OS),采用Master/Slave模式的数据并行策略,建立基于消息传递的工作站机群系统,用MPI进行算法编程验证实验。

本实验的主要目的是验证并行化后的K-means算法的执行时间和效率,所以为了简单起见,本实验中的数据是通过计算机随机产生的整型数据。同时,我们将并行与串行算法的实验结果相比较,当进行算法比较时,把程序运行10次并取平均值进行作图比较(如图2)。

从图2中我们可以看出并行K-means在数据集较大时表现出比串行K-means更好的执行效率,而当数据集较小时,主要由于并行计算中PC间通信时耗较大,所以单机串行算法表现出相对更高的执行效率。实验可以证明K-means算法在并行机群上具有了良好的并行性和可扩展性。

2 PAM算法

2.1 PAM聚类算法

PAM是k中心点(k-medoid)算法之一,它试图确定n个对象的k个划分。在随机选择k个初始代表对象之后,该算法反复试图选择簇的更好的代表对象。分析所有可能的对象对,每对中的一个对象看作是代表对象,另一个看做非代表对象。对于每个这样的组合,计算结果聚类的质量。对象oj被那个可以使误差值减少最多的对象所取代。再一次迭代中产生的每个簇中最好的对象集合成为下次迭代的代表对象。最终集合中的代表对象便是簇的代表中心。PAM算法的处理流程如下[2]:

输入:k:结果簇的个数

D:包含n个对象的数据集合

输出:k个簇的集合

方法:

1)从D中任意选取k个对象作为初始的对象或种子

2)repeat

3)将每个剩余对象指派到最近的代表对象所代表的簇

4)随机地选取一个非代表对象Orandom

5)计算用Orandom交换代表对象Oj的总代价S

6)if S<0 then用Orandom替换Oj,形成新的k个代表对象的集合

7)until不在发生变化

2.2并行化PAM改进算法

为了使问题简单化,首先我们选择任意的当前k个对象作为节点{Ol,…,Ok}。对于PAM算法,当每一步结束时,一种情况是找到一个代价最小的相邻节点,另一种情况是算法结束(当前节点代价最小)[3]。如果我们需要从当前节点移动到一个新的节点,我们必须交换一个已选对象和一个未选对象。为了保证已选对象在前k位,我们交换他们的下标。这样{Ol,…,Ok}会一直作为当前节点,而且不会受到当前节点移动的影响。

PAM的主要任务是检查当前节点的所有相邻节点,而且必须在划分的同时检查[3]。假设在p个进程(记为p1,p2,…,pp)上运行PAM算法。算法描述为:

1)将所有相邻节点写在列表中并按下标(升序)排序;

2)前[k(n-k)/p]个相邻节点指派给p1,接着的[k(n-k)/p]个相邻节点指派给p2,…,最后的[k(n-k)/p]个相邻节点指派给process p;

3)p个进程并行,并且报告各自相邻节点n1,…,np;

4)如果没有相邻节点被报告,算法结束(当前节点的代价最小);

5)从n1,…,np中选择代价最小的节点,将此节点改为当前节点,重复第一步。

下面举一个例子简单说明该算法,给定一个对象集{1,2,3,4,5,6,7},假设k=4,“1234”相邻节点为(用上述方法得到):1235,1236,1237(i=4);1245,1246,1247(i=3);1345,1346,1347(i=2);2345,2346,2347(i=1)。

每个进程被指派任务后,各自查找代价最小的节点,最后所有的进程(除了p1)将得到的节点报告给p1,由p1作比较工作。

2.3实验结果与分析

利用2.3中搭建的工作站机群系统,此时用3台PC机,进行PAM算法的执行效率验证,并对比串行和并行PAM的执行时间(如图3),由于PAM算法不适用于大量数据集的处理,所以实验n取1000以内的数值。

从图3中我们可以看出并行PAM的执行时间比串行PAM的执行时间长,并没有提高算法的执行效率,由此我们可知K-means算法有比PAM更好的并行性和可扩展性。

3具有并行性的其他聚类算法

聚类算法中除了上述K-means、PAM算法具有潜在的并行性和可扩展性外,还有一些算法可以进行并行化处理例如:并行硬聚类算法中的K-mediods,面向大规模数据库系统的BIRCH算法,处理非数值属性聚类的CACTUS算法,子空间聚类算法ENCLUS等[4],以及模糊聚类算法中的FCM等算法,理论上也具有在并行机群上的加速性。

4进一步研究方向与展望

近年来诞生了聚类算法中的一个崭新分支和研究热点—谱聚类算法,谱聚类算法建立在谱图理论之上,其实质是将聚类问题转化为图的最优划分问题,相对于传统的聚类算法有许多优势,并在实践中取得了很好的效果。由于谱聚类算法一般可以归纳总结为三个步骤[5]:

步骤一:构造数据集表示矩阵Z;

步骤二:计算Z的前k个特征值和特征向量,构造特征值的向量空间;

步骤三:利用K-means或其它传统聚类算法对特征向量进行聚类。

由于谱聚类算法研究中可以运用K-means算法等具有并行性的聚类算法进行特征向量的聚类,所以本文对K-means算法并行化的研究也可以运用于谱聚类的并行化,提高谱聚类算法的执行效率,是很有前景的研究问题。

参考文献

[1]Tan P N,Steinbach M,Kumar V.Introduction to Data Mining[M].Beijing:POSTS&TELECOM PRESS,2006.

[2]Jia W H,Micheline Kamber.Data Mining Concepts and Techniques[M].Beijing:China Machine Press,2006.

[3]Yan Z P,Shao X L,Zhang F.Parallel Clustering Methods for Data Mining[J].Acta Scientiarum Naturalium Universitatis Nankaiensis,2008(4):106-111.

并行化仿真篇7

视频压缩编码是图像处理领域的热门方向。为得到较高的压缩率使得图像数据易于存取和传输, 需要大量复杂的运算, 使得在CPU处理器上难以进行大尺寸图像的实时编码。

为此, 出现了许多实现高速运算的解决方案。当前的主流方案平台是基于ASIC、DSP, 并且已有基于这些平台的编码器相关实例。如TI公司的DM642 DSP[1], 可实现40fps的CIF格式视频图像编码;富瀚微公司的多通道视音频编码芯片FH8735, 具有1路1080p, 2路720p或者8路标清的实时编码能力。但DSP运算能力有限, 无法实现高清画质视频的实时编码;ASIC虽然运算能力较强, 但开发周期长、前期投入大、风险高、灵活性差, 无法适应市场多样性的需求。

GPU是目前发展迅速的高速处理器, 主流GPU的单精度浮点处理能力已达到933GFlops, 为同期CPU的10倍左右;其外部存储器带宽141MB/s, 为同期CPU的5倍左右。2007年NVIDIA公司推出了CUDA平台, 利用CUDA平台, DCT变换可实现8倍速度的性能提升[2];Canny边缘检测可实现22倍速度的性能提升[3];而H.264编码中的ME模块运算速度可提高12倍[4]。CUDA平台开发成为大数据量计算、存储的最优方案, 在性能和开发时间上都具有巨大的优势, 资源的利用率大大提高。

本文基于CUDA平台实现了一种高度并行化MPEG-2高清视频编码器, 无论从系统架构上, 还是模块算法与内存管理上, 都进行了最优于并行化的处理, 并在CUDA平台上予以实现, 从大数据量的运算和传输两个方面上达到了实时要求。

1 并行化编码系统

本文方案选用MPEG-2视频标准作为并行化编码的对象, 其主要原因在于:第一, MPEG-2相对MPEG-4, H.264以及AVS来说, 具有相对较低的复杂度[5], 而对于现在发展迅速的网络带宽以及海量存储, 完全可以适应码率较高的MPEG-2[6]码流。第二, 在MPEG-2中, 宏块间独立性很强, 便于并行化, 非常适合CUDA平台的架构。第三, 由于MPEG-2标准已被ATSC, DVB和ISDB等多个HDTV标准采纳, 作为HDTV节目的信源压缩算法, 在工业界已被广泛的使用。因此, 对其研究具有实际的意义和价值。本文方案综合MPEG-2编码系统以及CUDA平台并行资源的特点, 从系统架构设计层面以及模块算法层面均做了并行化设计, 以满足大数据量实时运算的要求。

1.1 分级并行编码

CUDA平台具有不同级别的并行运算单元。其最小并行运算单元为Thread;若干个Thread组成一个Block, 而Block与Block间同样为并行执行的运算单元;若干个Block组成一个Grid, 每个Grid顺序执行各自的kernel编码模块主函数。对于可并行化的每个编码主模块, 将一帧图像映射到CUDA的一个Grid中。从而在GPU内部实现了两个级别的并行, 即Block级和Thread级。在MPEG-2编码器的系统架构中, 由于不存在帧内预测, 除了VLC模块, 其余编码模块宏块间无相关性, 可做并行化设计。且在各编码模块内部, 算法同样具有可并行化之处。因此, 将待编码图像的不同运算单元映射到CUDA平台的相应并行运算单元中去, 便可实现最优的并行化设计。

本文方案分别针对编码器系统采用了3个级别的并行化处理, 最大限度地利用可并行资源。

1.1.1 处理器级并行

根据CPU与GPU两处理器间的异构并行特性, 本文方案采用两者并行处理系统架构。CPU负责逻辑性强的事务处理和以及读写数据的串行计算, 而GPU负责高度并行化的数据处理。两处理器之间的通信主要通过以下两方面工作完成:一是CPU调用在GPU端执行的kernel函数, 实现指令的调度和执行;二是CPU与GPU间相互的内存复制, 实现数据的交换。

图1为并行化编码系统架构设计图, 虚线左侧为CPU执行部分, 负责系统模块的调度, 以及视频图像序列帧数据的读取、输出码流的写回等工作。虚线右侧为GPU执行部分, 依次执行几个编码主模块的并行运算工作。而数据交换的过程和CPU函数以及GPU的kernel函数的执行是相互异步的, 即用GPU的kernel函数的执行时间来掩盖另外两个过程的执行时间。

将视频图像序列的数据读入CPU内存的时间和VLC输出码流写回的过程是必须在CPU端完成的事务, 若按照原编码流程, 将有很大一部分时间GPU处于闲置状态, 等待CPU端完成上述工作。因而, 本文方案根据CPU与GPU各自的优点, 采用协同工作的方法, 利用两者的异步特性来掩盖部分CPU执行时间, 提高整体编码速度。

由于ME、MC以及VLC模块是在GPU上运行的相对耗时部分, 因而本文方案使用当前帧的VLC执行时间来掩盖下一帧的图像数据读入时间;且使用下一帧的ME执行时间来掩盖当前帧输出码流的写回时间。采用这样的系统设计架构, 使得CPU与GPU能够获得最大限度的占用率, 减少各自的闲置时间, 实现了处理器级的并行化。

1.1.2 Block级并行

MPEG-2编码器的主模块在编码过程中, 同一帧的宏块间大多无相关性, 从而为并行运算的实现提供了可能。因此, 本文方案去掉了需要双向预测的B帧以及场图, 以高并行度的高速运算来弥补压缩率的少量损失。对于ME、MC、DCT/IDCT、Q/IQ模块, 将一帧图像中每一个宏块分别映射到CUDA平台一个Block中, 从而使得图像中所有宏块的编码运算并行执行。而对于VLC模块, 由于其包含了差分编码, 当前宏块与前一宏块有一定的相关性, 相邻宏块间独立性不佳。针对VLC模块, 本文方案将图像每一行宏块数据映射到一个Block中, 使得行与行之间并行地进行运算。从而实现了Block级的并行化。

1.1.3 Thread级并行

CUDA的一个Block可有多个并行执行的Thread, 根据不同编码模块的算法特点, 合理地设置Block维度以及分配使用Thread资源, 是提高GPU资源占用率以及计算速度的关键。在CUDA中, SM的指令是以warp为单位发射的, 即一个warp中的Thread是同步执行的。因此, 设计过程中, 应尽量使得Block中Thread数为warp内最大Thread数的整数倍, 以获得更好的资源占用率。针对此特性, 本文方案采用如下的Thread级并行算法。

如表1所示, 针对不同的模块, 为其分配合适的Block维度。对于ME模块, 将一个搜索窗中各个参考宏块映射到当前Block的各个Thread中去, 在各Thread中分别进行其与当前宏块的匹配计算, 从而使得每个参考宏块与当前宏块匹配的计算并行执行。而对于MC模块, 为进一步提高并行度, 将一个宏块的Y分量划为32个4×2块, 然后映射到一个Block的32个Thread中去。

DCT/IDCT的运算, 相当于原始图像矩阵与变换核矩阵的相乘运算;而对于Q/IQ的运算可转化为高度并行化的矩阵运算。对于码率控制部分, 本文方案沿用复杂度较低, 易于并行化的TM5参考算法。因此, 本文方案将每一个像素点的计算都映射到一个Thread中去 (一个宏块划为4个8×8块, 即4×8×8个像素点) 。每个Thread分别读取A与B两矩阵的一行与一列, 计算出结果矩阵C的一个点, 从而实现高度并行化的矩阵相乘运算。

熵编码VLC模块完成在GPU上的查表工作, 对于每一行图像, 采用80个Thread使系统能为该行内全部宏块并行地查询储存在常数存储器中的表, 以达到行内的并行, 并利用GPU中cache资源提高读取速度。

采用上述一系列的设计方案, 从处理器级、Block级以及Thread级三方面依次对编码器进行并行化处理, 使得系统资源利用率和编码速度大大提高。

1.2 ME并行快速搜索算法

MPEG-2编码标准中采用了ME全搜索算法, 在预先定义的搜索区域内, 把当前宏块与所有参考宏块进行匹配运算, 寻找最佳匹配参考宏块。由于CUDA平台Block内可并行Thread最多仅能有32个, 全搜索法的覆盖范围远超过32个点, 这使得每个参考宏块的匹配过程需要多次并行才能完成, 该算法在得到最大搜索精度的同时也带来了较长的运算时间。而目前速度最快的Four Step search (FSS) 四步搜索法[7]每个搜索窗内仅有25个 (5×5模板) 或9个 (3×3模板) 参考宏块, 均小于CUDA平台上一个warp内的Thread数32, 这使得部分Thread被闲置, 并行度不高。且FSS搜索算法最大搜索覆盖范围为垂直和水平方向 (-7, +7) , 运算精度剧烈下降。因此, 本文方案提出一种基于CUDA平台的快速ME并行搜索算法, 最大程度利用可并行资源, 搜索模板如图2所示, 包含粗搜索用大模板以及细搜索用小模板。两种模板的参考点均为32个, 保证了CUDA平台warp内32个Thread全部并行地进行运算, 充分利用了资源, 并且提高了搜索精度。

1.2.1 算法步骤

本文方案提出的算法搜索最佳匹配参考宏块过程分为以下几步完成。

Step 1:以当前帧当前宏块左上顶点像素点位置为初始搜索起点。

Step 2:应用大模板在前一帧进行第一次搜索, 此为粗搜索, 在较大的覆盖范围内搜索最接近当前宏块的参考宏块可能位置区域, 得到当前宏块的最佳匹配参考宏块。

Step 3:以Step 2中得到的最佳匹配参考宏块左上顶点像素点为第二次搜索的初始搜索起点, 并应用小模板进行搜索。经过以上步骤, 最终得到最匹配参考宏块, 搜索结束。

本文方案提出的搜索算法可覆盖区域较大, 搜索范围达到了垂直方向 (-17, +17) , 水平方向 (-17, +18) , 相比FSS搜索法其搜索精度大大提高;且整个过程仅需要两步搜索, 每一步的匹配运算都是完全并行的, 从而本文方案搜索算法相比全搜索法, FSS搜索法, 速度有了大幅度提高。

1.2.2 实验结果

通过在CUDA平台上对几种ME搜索算法的实现, 对于不同运动程度的测试序列, 所得到的ME搜索所需时间以及相应PSNR如表2所示。

实验结果表明, 从搜索速度上来说, 本文方案所提出的快速ME并行搜索算法相比全搜索算法快至少6倍以上, 而相比FSS搜索算法快至少3.4倍以上。从编码质量上来说, 采用全搜索算法的编码图像质量最好, 搜索精度最高;采用本文方案算法的编码图像质量稍有下降, 但仍优于FSS搜索算法。即本文方案提出的并行化搜索算法充分利用了CUDA平台的并行资源, 搜索速度快, 搜索范围广, 得到了较高性能。

2 内存管理

为满足高清视频图像数据的实时传输存取, 无论在CPU与GPU间的数据传输上, 或是GPU内部的内存分配管理上, 都需要足够的带宽。对于720p视频图像 (灰度) , 其每帧所产生的数据量为1280×720×8/8=900kB, 每秒所产生的数据量为25×900KB=22.5MB。GPU与CPU通过PCI-E接口进行通信, 实验用PCI-E通道提供了上下行各4GB/s的带宽, 满足了两处理器之间大数据量高速传输的条件。而GPU的GDDR内存带宽为141MB/s, 使得数据可以高速地在GPU的内存间进行传输。

表3与描述了再编码过程中CPU与GPU的内存通信管理的过程。本文方案使用了Pinned memory, 在尽可能小地影响CPU运行性能的情况下, 得到了较高的CPU与GPU之间的内存复制速度。且将编码中所用到的固定参数以及熵编码查找表存入GPU端只读的Constant memory, 相比Global memory得到更高的GPU内部内存复制速度。每个SM运算单元通过存取Constant memory和Global memory中自己所需数据, 来完成块级编码, 最终通过Global memory输出码流至CPU端。

3 实验结果

为验证本文所提出的并行编码系统的性能, 选取了7个高清格式 (720p) 标准测试序列进行了测试。序列帧数为200帧, 帧格式为IPPP, GOP (Group of Pictures) 长度为15帧 (1个I帧和14个P帧) , 码率8Mbps。实验环境为: (1) CPU运行平台AMD Sempron 2800+, 主频为1.60GHz, 内存为512M的DDR; (2) GPU采用NVIDIA GeForce GTX285显卡; (3) Microsoft Windows XP sp2; (4) Microsoft Visual Studio 2005; (5) CUDA Toolkit and SDK 2.2。

表4给出了本文方案与单CPU方案的性能对比。可以看到, 单CPU方案编码速度不超过0.13fps, 距离视频序列的实时编码存在巨大差距。而本文方案的编码速率相比单CPU方案有了大幅度提升, 对于运动剧烈程度不同的测试序列, 编码速率达到26fps～35fps, 实现了高清视频序列的实时编码。

与此同时, 本文方案得到了较好的编码图像质量, 如图3所示, 对于实验中7个测试序列, 图像的平均PSNR值均大于33.5dB。相对单CPU方案, 图像质量略有下降, 究其原因, 主要是由于, 第一, 本文方案提出的利于并行运算的ME搜索算法, 相比原MPEG-2编码器的全搜索算法搜索精度下降;第二, 由于是高度并行的编码, 对于量化步长的确定是基于帧级别的, 而在原始CPU端进行的量化是串行基于宏块级的。这都对计算精度产生了一定的影响, 但其主观画质评价并未有明显下降, 失真属于人眼可接受程度。

4 结束语

本文首先分析了MPEG-2编码器的可并行性, 并基于CUDA平台的软硬件体系, 利用GPU的高存储带宽以及并行运算能力, 提出了一种具有高并行度的高速视频编码系统设计方案;其次, 对各编码主模块在CUDA平台上做了并行化优化设计, 并提出了一种具有高并行度、高搜索覆盖率的快速ME搜索算法;最后, 优化的系统模块调度以及资源分配, 使得CPU与GPU之间具有高度的并行性, 并且各编码模块内部的运算也具有不同级别的并行性。此外, 合理的内存管理, 使得大数据量的实时传输以及运算速度满足需求。最终, 测试结果显示, 基于CUDA平台的ME并行算法无论是在运算速度上还是搜索精度上都取得了优异的性能。且在测试环境下的CPU+GPU方案编码速度远高于单CPU方案, 对于720p格式视频编码速度达到25fps以上;同时保持了编码图像质量;在极小的硬件代价和开发时间的情况下, 便实现了高清视频实时编码的目标。

摘要：HD视频的编码技术具有广阔的应用前景, 为解决其大数据量实时处理的瓶颈问题, 提出一种基于CUDA平台的并行编码系统架构。根据CUDA平台软硬件结构特性, 采用三级并行机制;并提出一种高并行化快速ME搜索算法;同时合理分配内存空间, 实现大数据量的实时运算与存取。实验结果表明, 方案具有高并行度, 高编码速率的特点, 对HD视频可达到实时编码要求。

关键词：CUDA,并行计算,视频编码,ME

参考文献

[1] Lin Daw-Tung, Yang Chung-Yu. H.264/AVC Video Encoder Rea-lization and Acceleration on TI DM642 DSP[J]. PSIVT 2009, LNCS 5414, 2009:910-920.

[2] Yang Zhiyi, Zhu Yating, Pu Yong. Parallel Image Processing Based on CUDA[C]. Proceedings-International Conference on Computer Science and Software Engineering (CSSE 2008) , 2008, 3:198-201.

[3]Seung In Park, Sean P Ponce, Jing Huang, et al.Francis Quek.Low-cost, high-speed computer vision using NVIDIA’s CUDA ar-chitecture[J].Proceedings-Applied Imagery Pattern RecognitionWorkshop, 2008.

[4]Chen Wei-Nien, Hang Hsueh-Ming.H.264/AVC motion estimationimplementation on compute unified device architecture (CUDA) [C].International Conference on Multimedia and Expo 2008 (IC-ME2008) , 2008:697-700.

[5] Ishfaq Ahmad, Yong He, Ming L Liou. Video compression with parallel processing[J]. Parallel Computing.2002, 28 (7-8) :1039-1078.

[6] Le Gall, Didier. MPEG: a video compression standard for multimedia application[J]. Communications of the ACM, 1991, 34 (4) :46-58.

本文来自 360文秘网(www.360wenmi.com)，转载请保留网址和出处

【并行化仿真】相关文章：

并行化分析05-22

浅议“四个能力建设”与“户籍化管理”并行对消防管理产生的影响01-20

课程仿真化教学05-30

并行检测05-01

双轨并行05-04

并行配置05-17

并行开发05-22

并行教学07-15