非线性处理（精选8篇）
非线性处理 篇1
0 引 言
科研人员在解决科技或工程中的一些实际问题时,要做大量试验,得到大批数据,然后还要对这些数据进行分析处理,以求找到一些变量间的因果关系。为了得到变量之间关系,有时要探求依据,进行理论推导,得到所要求的方程或公式;有些时候则是基于试验数据,利用拟合方法来确定一些变量间的因果关系,在此过程中,经常会遇到试验数据的非线性拟合。拟合因变量与自变量的关系曲线并导出其计算公式比较复杂,需要利用专门针对数据处理和制图的Origin或MATLAB等专业软件进行和完成。这里以三种线形为例,简要介绍利用Origin 6.0[1]进行非线性拟合的推演过程。三种线形分别选择三种函数,进行非线性拟合,得到了拟合后的公式。一些论文中常会出现这三种线形,所以提出的非线性拟合方法在实际应用中具有很广阔的适用范围。同时,本文也介绍了拟合时的注意事项及拓展应用等问题,以求与同行们彼此借鉴,共同探讨非线性拟合方法,推进非线性拟合的研究进程。
1 用自然指数(Exponential)曲线拟合
硅光电池是不需外加电源而可直接将太阳辐射能转换为电能的器件,研究其伏安特性,对了解和使用各种光电器件具有十分重要的意义。设硅光电池的输出电流为I,输出电压为U1,有一组全暗情况下硅光电池的伏安特性的实验数据[2],如表1所示。
现在,假设要寻找因变量I受自变量U1的影响的变化规律,就可以对这组数据进行非线性拟合,找到I受U1的影响的变化规律。过程如下:
打开Origin 6.0,在工作表的A[X]列输入U1/V数据,在B[Y]列输入I/μA数据,然后,在菜单栏的Analysis(分析)中选Non-Llinear Curve Fit(非线性拟合),就会出现非线性拟合菜单(如果打开的是基本模式,单击“More”按钮,切换到高级模式),这个菜单上部有可见的菜单栏和工具栏。拟合可按下列步骤进行。
(1)选择函数。
首先打开Select Functions(选择函数)对话框,在此对话框左部的列表框Categories(类型)中选择Exponential(自然指数),而右部列表框Functions(函数)中选取Exp3P1,在函数三种预览方式(Equation 公式、Sample Curve 示样曲线和Function File 函数文件)中可预览选中的函数,公式的形式如下:
(2)编辑函数。
在同上的对话框菜单栏Function中选择Edit,出现Function Edit对话框,按“Save”按钮;
(3)指定变量。
在同上的对话框菜单栏Action中选择Dataset,出现Select Dataset对话框,在上部列表框单击Y变量,在中部列表框单击Data1_b(X变量对应于Data1_a数列已自动指定),然后单击“Assign”按钮,此时会出现一幅曲线图。由出现的曲线图可知,在选择函数时是因何选取Exponential(自然指数);
(4)曲线模拟。
在同上的对话框菜单栏Action中选择Simulate,出现Simulate Curves对话框,将其上的a、 b和c(即拟合时所要寻求的系数值)文本框内的1、1和0.5分别改写为0. 5、-25和-5,单击“Create Curve”按钮,在曲线图上会出现一条摸拟线;再将a、 b和c分别改写为0.1、-20和-4,单击“Create Curve”按钮,在曲线图上又会出现一条摸拟线;
(5)实施拟合。
在同上的对话框菜单栏Action中选择Fit,出现Fitting Session对话框,单击下部“Chi-Sqr” 按钮,再单击“10Iter” 按钮(单击一次进行10次迭代),可再单击数次“10Iter”,直至其中的a、 b和c的数值不变为止;
(6)显示结果。
在同一对话框的菜单栏Action中选择Result,出现Genarate Result 对话框,单击“Param. Worksheet”按钮,生成Parameters工作表,关闭对话框,会出现Result Log窗口;
(7)查看结果。
在Parameters工作表中可以看到拟合后的参数a、b与c的值。若未显示出数值,可将列加宽,并双击列名(Value)出现Worksheet Column Format对话框,在其下部调整Column Width数值,若要求给出科学记数形式,可在Format下拉菜单中选定 Scientific:1E3,在Numeric Display下拉菜单中选取Set Decimal Places=,同时将其后面的数字定为5,单击“OK”按钮,Parameters工作表中的a、b与c的值就会转换成符合要求的位数。在Result Log窗口和Graph1图上的文本框中都有参数a、b与c的值。至此,得到拟合的参数值,a=0.067 65,b=-19.145 04,c=-3.587 85。此外有些结果是供分析拟合误差之用,如有必要可以查看,此处仅只讨论拟合方法。
为了使拟合曲线更加突出,要删掉模拟线,可在下部窗口中直接删除exp3p11与exp3p12(exp3p1是函数名,后面的1和2是顺序编号)。初始设定a、b与c值的两条模拟线删除后,留下原始数据曲线和拟合曲线。另外,Graph1图上还有显示拟合参数的文本框,可留也可删。此时,下部窗口文件夹中有数据表、拟合曲线图、拟合数据表和参数表四份文件。
遵照上述步骤进行拟合,拟合前的I-U1曲线图,如图1所示。 绘制有拟合曲线的I-U1曲线图,结果如图2所示。
这里要说明的是,自然指数曲线拟合成功与否,需要视拟合曲线与原始数据曲线符合的程度,符合程度越高,拟合越成功,这与初始参数值的设定有着极其密切的关系,而初始参数值的设置则是通过试验最终确定的。如果使用Simulate Curves对话框给定的值a=b=1、c=0.5,拟合会失败;如果只用a=0.5、b=-25、c=-5拟合,单击次数会太多;而给出a=0.1、b=-20、c=-4,就是为了拟合得更好。因为拟合是由最后给出的模拟数据启动的,若先给出a=0.1、b=-20、c=-4,后给出a=0.5、b=-25、c=-5,先前的数据将不起作用;若能直接给出a=0.1、b=-20、c=-4,当然最好。
确定a、b与c的值之后,就得到了I-U1的关系计算公式,如之前选择Exponential(自然指数)那样,拟合曲线的公式为:
经数学变换,等号两边取自然对数,并加以整理,也可以表示为:
具有此种线形的实例很多,但不一定采用自然对数拟合,也并非都需要拟合。这种线形在很多论文都可以看到,比如拉弯联合载荷下弹塑性J积分估算方法研究,J积分与载荷比的关系[3],有多条曲线均属这一类型,水面上的饱和蒸气压和冰面上的饱和蒸气压与温度的关系[4],也是这种线形。
2 用幂函数(Power)曲线拟合
现有照度变化率随焦距变化(光照面Φ=5mm时)曲线[5] 如图3所示。有拟合曲线的照度焦距曲线如图4所示。
为此,假设要寻找最大照度变化率(Y/%)随焦距变化(f/mm)的相应规律,就可以对这条曲线进行非线性拟合,得到Y受f影响的变化规律。
拟合步骤与曲线拟合1大致相同,但在Categories中选择了Power(乘幂),在Functions中选取了Pow2P2,其公式形式是:
y=a(1+x)b (4)
在曲线拟合时,只需确定a与b两个系数。
在曲线模拟时,将a与b的1、1第一次变为600和﹣1.6,第二次变为1 000和﹣1.8,绘制有拟合曲线的Y -f曲线图,如图4所示。
在查看结果时,得到a=2 728.900 84、b=﹣2.138 31。由此,得到Y-f关系式为:
Y=2 728.900 84(1+f)-2.138 31 (5)
具有如此线形的实例也很多,但不一定采用Power拟合,也并非都需要拟合。这种线形在相当多的论文中也都可以看到,比如中厚板辊式淬火机淬火过程的温度场分析[6],讨论不同换热系数下钢板心部温度与时间的变化曲线,许多曲线都是这种类型曲线;再比如,在自旋对量子中束缚磁极化子性质的影响[7],则有更多的同种类型曲线。
3 用玻耳兹曼函数(Boltzmann)曲线拟合
现有催化剂用量对罗丹明B溶液降解率的影响曲线图[8],如图5所示。图5的拟合曲线如图6所示。
在此,假设要寻找因变量降解率(Y/%)受自变量催化剂含量(X/g)的影响的变化规律,就可以对这条曲线进行拟合,得到Y受X的影响的变化规律。
拟合步骤与曲线拟合1大致相同,但在Categories中选择了Origin Basic Functions,而在Functions中选取了tzmann(玻耳兹曼),公式形式是:
在曲线拟合时,确定A1、A2、x0和dx四个系数。
在曲线模拟时,将A1、A2、x0和dx的0、1、0和1第一次变为0、90、0和0.15,第二次变为0、100、0和0.20,绘制有拟合曲线的Y﹣X曲线图,如图6所示。
在查看结果时,得到A1=–18.418 41、A2=100.395 54、x0=0.084 61和dx=0.179 21。如此得到Y-X关系式为:
具有同种线形的实例也是很多的,但也不一定采用Boltzmann(玻耳兹曼)拟合,也并非都需要拟合。这种线形在许多论文中也都可以看到,比如在基于现场试验的超长桩端阻力承载性状研究[9],调研了很多建筑物的超长桩端阻力,对端阻力-桩端沉降的关系,采用Boltzmann(玻耳兹曼)拟合进行了归一化处理,并给出了九种建筑物(六座桥梁三座大楼)端阻力-桩端沉降曲线图及其Boltzmann(玻耳兹曼)拟合结果。
4 结束语
对试验数据进行非线性拟合时,需要选择函数,有的函数存在于多种类别之中,如Lorents就分别出现在Origin Basic Functions、Peak Fuctions和Spectroscopy中。但文中所关注的只是函数,而并不考虑其归属类别。此外,有的函数可以拟合多种线形,如Hyperbl和Exp3P1Md就是这样的函数,因其系数变化,就会产生不同的线形。但文中关注的则为其是否满足将欲拟合的线形要求。
试验数据非线性拟合的成败,首先与选用的函数有关。一条试验数据曲线,可以选择几种函数进行拟合,但不一定均获成功。拟合曲线与原始数据曲线符合的程度越高,即越成功,而在拟合成功的函数中还要选择更为优异的,这就需要进行拟合的误差讨论。在不进行误差讨论的情况下,则需凭直接观察做出判断,拟合程度最高的就是最好的选用函数。
试验数据非线性拟合的成败,其次还与初始参数值的设定有着密切的关系,而初始参数值的设定是通过多次试验选取得到的。拟合是从最后给出的模拟数据开始的,若使用Simulate Curves对话框提供的系数值,或者给出的模拟数据与试验数据曲线偏离较远,拟合都会失败;若给出的模拟数据与试验数据曲线非常接近,则拟合时点击次数会很少,会很快得到成功的结果,否则拟合时点击次数太多甚至会失败。
文中讨论的是试验数据的非线性拟合,实际拟合还有Linear(线性拟合)、Polynomial(多项式拟合)、Sigmoidal (S拟合)等,这些拟合可以综合运用,比如一个因变量有两个自变量,分别进行拟合,得到两个公式,最后将其综合在一起,成为一个关系式,这在科学试验和工程中也有机会遇到。
摘要：科技人员经常要对大量试验数据进行分析处理,以求找到一些变量间的因果关系。在进行数据处理时,会经常遇到非线性拟合问题。用三种线形为例,分别选择Exponential(自然指数)、Power(乘幂)和Boltzmann(玻耳兹曼)三种函数,以简单的操作步骤,进行非线性拟合,找出了变量间的因果关系,并举出应用实例。同时也介绍了拟合时应注意的事项及引申应用等问题。
关键词：数据处理,非线性拟合,变量间关系,应用实例
参考文献
[1]郝红伟,施光凯.Origin6.0实例教程[M].北京:中国电力出版社,2000.
[2]张玮,杨景发,闫其庚.硅光电池特性的实验研究[J].实验技术与管理,2009(9):42-46.
[3]白永强,汪彤,等.拉弯联合载荷下弹塑性J积分估算方法研究[J].工程力学,2010(3):6-9.
[4]李英干,范金鹏.湿度测量[M].北京:气象出版社,1990:6.
[5]梁世安,胡友旺,段吉安.阵列波导器件封装中机器视觉系统的光源分析与设计[J].光电子技术,2009(1):55-59.
[6]袁国,王国栋,王黎筠,等.中厚板辊式淬火机淬火过程的温度场分析[J].东北大学学报:自然科学版,2010(4):527-530.
[7]李志新,肖景林.在自旋对量子中束缚磁极化子性质的影响[J].固体电子学研究与进展,2009(1):10-13.
[8]周云龙,胡志彪,等.TiO2/竹炭复合材料研究(Ⅱ)光催化降解性能[J].功能材料,2010(2):197-199.
[9]蒋建平,章杨松,高广运.基于现场试验的超长桩端阻力承载性状研究[J].工程力学,2010(2):149-160.
非线性处理 篇2
摘 要：BLAS库是基本线性代数子程序库，是许多大型科学与工程计算的核心计算程序，FitenBLAS库是在多核多线FT1000微处理器上开发的基本线性代数库，其研制对FT1000微处理器在科学与工程计算中的应用具有重要意义.根据多级存储结构和寄存器的数目，设计了向量与向量、矩阵与向量和矩阵与矩阵运算的多级循环展开方法，采用指令调度、数据预取等通用优化技术，优化BLAS库串行程序.对于BLAS3子程序，设计了矩阵乘无冗余数据拷贝分块算法，采用指令重排、访存与计算的重叠、分块等技术优化矩阵乘子程序，基于矩阵乘子程序实现了其他BLAS3子程序.研制了汇编线性代数程库FitenBLAS，其核心子程序矩阵乘的双精度计算性能达到6.91Gflops，是峰值性能的86.4%.
关键词：FT1000微处理器；BLAS库；性能优化
中图分类号：TP332.2 文献标识码：A
基本线性代数子程序BLAS（Basic Linear Algebra Subroutines）库，提供最基本的线性代数函数接口[1]，分为三级：BLAS 1（Level 1）包括向量与向量操作子程序，如点积、向量相加等.BLAS 2（Level 2）包括矩阵与向量操作子程序，如矩阵向量相乘等.BLAS 3（Level 3）包括矩阵与矩阵操作子程序，如矩阵与矩阵相乘等.
BLAS库是每款微处理器要移植和优化的数学库，是许多大型科学与工程计算的核心计算模块，同时BLAS库子程序可以反映许多应用程序的计算特点，如果BLAS库可以在微处理器上获得高性能，同样的应用程序也可以获得好的性能，BLAS库程序可以验证微处理器的功能和计算性能.因此各个厂家在新型号微处理器推出时，都会配套针对微处理器特点研制、优化和推出高性能BLAS库，BLAS库已经成为微处理器的必备数学库之一.
湖南大学学报（自然科学版）2015年
第4期迟利华等：FitenBLAS：面向FT1000微处理器的高性能线性代数库
Intel公司针对通用CPU开发了MKL基本数学运算库[2]，包含采用多线程进行并行计算的函数库，可以在结点内实现高性能；AMD公司针对通用CPU开发了ACML基本数学运算库[3]，具有MKL类似的特点；IBM公司开发了ESSL基本数学运算库[4].上述厂商开发的基本数学运算库，均包含了使用汇编语言手工优化的高性能BLAS库.GotoBLAS[5-6]是开源的针对不同类型的微处理器开发的BLAS库，提供了OpenMP多线程并行版本，是性能最好的BLAS库之一，2008后不再更新，不支持最新推出的微处理器；ATLAS[7-8]针对不同的微处理器提供可以自动优化的BLAS库.OpenBLAS是针对龙芯等微处理器开发的高性能BLAS库[9-10].
FT1000是由国防科学技术大学研制的单芯片多线程（CMT）处理器，是天河1A计算结点采用的微处理器之一[11].BLAS库在FT1000上获得高性能是需要研究的重要问题，目前没有可以在FT1000上运行的高性能BLAS库，本文结合FT1000多核多线微处理器特点，设计了并行计算方法和数据结构，研制了手工汇编子程序，进行了针对性的性能优化，研制了高性能线性代数库FitenBLAS.
1 FT1000微处理器
FT1000微处理器包含8个处理器核，每核包含8套硬件现场，支持8个线程.每个线程有1个完整的寄存器文件，大部分ASI，ASR和特权寄存器都是每个线程1份.
每个核包含2条整数流水线、1条浮点流水线和1条存储流水线.8个线程共享浮点流水线和存储流水线.8个线程分成2组，每组4个线程，共享1条整数流水线.虽然8个线程同时运行，但是在任意时刻，最多两个线程是活跃的，这两个线程发射的指令只可能是下面的组合：1对整数操作、1个整数和1个浮点、1个整数和1个存储、1个浮点和1个存储.每个组内的线程按照LRU算法每个周期进行切换.
每个核内有独立的一级数据cache和一级指令cache.L1I Cache大小为 16 kB，8路组相联，块大小32字节.L1D Cache大小为8 kB，4路组相联，块大小32字节.指令TLB为64项全相联，数据TLB为128项全相联.8个线程共享L1I，L1D和TLB，通过TLB中的自动释放机制使多个线程更新TLB.
FT1000微处理器采用SPARC V9指令集，现有的BLAS线性代数库不能直接运行，需要重新设计.
2 程序优化方法
2.1 循环展开方法
循环展开可以减少循环体的循环次数，减少分支执行的时间，为流水线提供更多的并行机会，是一种通用的程序优化方法，是手工汇编优化BLAS子程序的主要方法.在循环展开中，展开因子的选择是核心问题，目前的编译器一般只对循环体很小的循环进行循环展开，且使用的展开因子是很小的常数，这可能损失一些计算性能，特别是对于多层嵌套循环，编译器的优化效果不明显，只能采用手工的循环展开优化.
BLAS 1（Level 1）包括向量与向量操作子程序，包含SUM，DOT，AXPY，COPY，SCAL等子程序.当跨步为1时，BLAS1程序中的计算操作针对一维数据进行连续访问，具有良好的空间数据局部性，可以采用手工循环展开来优化性能.采用汇编编程时寄存器的数目限制了展开次数，为此，本文根据寄存器的数目进行多级展开，展开因子可以随意调节，在展开循环体内，循环使用寄存器.以AXPY为例来说明多级展开方法，图1给出了多级展开循环体，双精度浮点寄存器的数目为n，涉及2个向量运算，平均使用向量寄存器，每个向量使用的寄存器最大数目为n/2，循环展开因子可以取为m“*”n/2.ldd，fmuld，faddd和std是汇编指令，分别表示取数、相乘、相加和存数.为了表示方便，图1中使用了循环控制变量，在实际展开的循环体中，没有循环控制变量，而将图1中的循环体重复m“*”n/2次.
AXPY使用寄存器的多级展开循环体：
根据寄存器数目，图1中给出的是AXPY多级循环展开的一般方法，BLAS1其它子程序中可以根据图1中给出的展开方法，进行多级循换展开.针对FT1000微处理器，双精度寄存器数目为32，对于AXPY子程序，每个向量至多使用16个寄存器.使用16个寄存器时，循环展开因子就是16的倍数，同样地循环展开因子可以是12，13，14和15等数的倍数，可以根据实际测试情况，进行调整，以找到最佳的循环展开因子.
BLAS2包括矩阵与向量操作子程序，如矩阵向量乘、rank1和rank2矩阵校正等，涉及二维数组A和一维数组x和y间的计算操作，计算结果为一个一维数组.以矩阵向量乘GEMV子程序为例进行说明，为了复用一维数组x，对数组x进行分段，对数组A进行分块.在分配寄存器时，二维数组A，一维数组x和y以及保存临时变量的一维数组z给分配寄存器总数的1/4.针对FT1000微处理器，双精度寄存器数目为32，对于GEMV子程序，数组A至多使用8个寄存器，一维数组x和y分别使用8个寄存器.使用8个寄存器时，循环展开因子就是8的倍数，同样地循环展开因子可以是4，5，6和7等数的倍数，可以根据实际测试情况，进行调整，以找到最佳的循环展开因子.图2给出了BLAS2子程序GEMV多级展开循环体.
BLAS 3（Level 3）包括矩阵与矩阵操作子程序，包含GEMM，SYMM，SYRK，SYR2K，TRMM和TRSM等子程序.矩阵乘GEMM子程序是BLAS3的核心，其它BLAS3子程序可以基于GEMM来实现.下面以GEMM子程序为例来进行说明.所要求解的矩阵乘形式如下：
C=βC+αA×B
其中A∈Rm×k，B∈Rk×n，C∈Rm×n，α和β是实数.
矩阵乘在多级存储结构上获得高性能的基本方法是分块算法，将A，B和C划分成如下子矩阵：
原矩阵乘算法转化为多个子矩阵相乘，根据Cache和TLB的大小来选择子矩阵的大小.具体实现时选取Aip的大小为L2 Cache大小的一半，同时保证不发生TLB失效，并驻留在二级Cache中，直到不再使用为止，即完成如下操作：
Ci，0…Ci，N-1=Ai，p×
Bp，0，…，Bp，N-1.
计算时，子矩阵Bp，j（j=1，…，N）以流水的方式进入L1D Cache.Ci，j的数据不重用，不必长时间保存在L1D和L2 Cache中.
把寄存器总数的一半分配给矩阵C使用，另一半寄存器被矩阵A，矩阵B和临时变量平均使用.针对FT1000微处理器，双精度寄存器数目为32，对于GEMM子程序，数组C至多使用16个寄存器，数组A和B分别使用4个寄存器，剩下的作为临时变量寄存器使用.图3给出了BLAS3子程序GEMM按照4*4*4展开的循环体.
通过指令的合理调度、数据的预取和寄存器的合理使用，当Ai，p，Bp，j和Ci，j在Cache中时，可以发挥CPU的峰值性能.
2.2 数据预取
矩阵和向量是存放于内存中的，计算过程中，通过多级存储结构，将数据取到寄存器中开始运算，数据的预取可以很好地实现数据传输和计算的重叠，是提高数据空间局部性的性能优化方法.开始执行计算时，参与运算的矩阵和向量存放在内存中，开始的取数不会命中Cache，FT1000上从内存中取数到二级Cache的延迟大于100拍（即100个CPU时钟周期），Cache块大小为32字节，每次内存访问，同时会把内存中连续的32字节分别存放到二级Cache和一级数据Cache中.
数据预取就是将后面要用的数据提前取到L2 Cache中，将访存的数据传输操作和乘加计算操作重叠起来，如果计算时间大于数据传输时间，那么整个计算就可以得到CPU的峰值性能，如矩阵相乘.通过参数可以调整数据预取的间隔，获得最佳性能.
BLAS1，BLAS2和BLAS3均用到数据预取，本文仅以双精度AXPY为例来说明数据预取方法，图4给出了带数据预取的展开循环体.
AXPY带数据预取的展开循环体：
令k=8，PREFECHSIZE=8，SIZE=8
循环展开因子为m*8
数组a（*）、b（*）和c（*）表示不同的寄存器
for i ← 0 to m*8 step 8 do
for j ← 0 to 4 do
ldd x[m*i+j]， a[j]；
fmuld a[j]， alpha， a[j]；
ldd y[m*i+j]， b[j]；
faddd a[j]， b[j]， b[j]；
std b[j]， y[m*i+j]；
endfor
prefetch [&x[m*i]+PREFETCHSIZE * SIZE]， 0
prefetch [&y[m*i]+PREFETCHSIZE * SIZE]， 0
for j ← 5 to 8 do
ldd x[m*i+j]， a[j]；
fmuld a[j]， alpha， a[j]；
ldd y[m*i+j]， b[j]；
faddd a[j]， b[j]， b[j]；
std b[j]， y[m*i+j]；
endfor
prefetch [&x[m*i+4]+PREFETCHSIZE * SIZE]， 0
prefetch [&y[m*i+4]+PREFETCHSIZE * SIZE]， 0
endfor
图4 BLAS1子程序AXPY带数据预取的展开循环体
Fig.4 Multilevel loop unrolling with the prefetching
data for subroutine AXPY of BLAS1
FT1000通过预取指令对内存数据进行操作，SIZE表示一个数据所占的内存字节数，PREFECHSIZE表示提前预取数据间隔，可以根据需要改变大小.图4中，预取的数据提供给下一次循环体使用，将数据从内存中，取到二级Cache中，预取数的时间和图4中的两个小循环进行时间重叠，用计算重叠数据传输.
3 多线程并行计算方法
对BLAS1采用向量一维平均划分的方式组织并行计算.对BLAS2中涉及的二维矩阵采用按行或按列一维划分.对BLAS3中涉及的二维矩阵采用按行和按列二维划分.
下面重点对BLAS3中的矩阵乘并行计算方法展开讨论.
P 个线程按 pr×pc 划分成二维拓扑结构， 满足 P=pr×pc.假设 P（r， c） （0≤r 多线程矩阵乘并行算法： 全局数组 局部数组 r，c， 每一个线程P（r， c），其中 ，执行： for i→r×mr to r×mr+mr do for l←0 to k-1 do 第r行线程中的任一线程将子矩阵 拷贝到 for j←c×nc to c×nc+nc do 将ij拷贝到 Call GEMM_kernel（bm，bn，bk，，，ij， LDC） end for end for end for 图5 避免重复拷贝A子矩阵的多线程矩阵乘并行算法 Fig.5 Multithread parallel matrix multiplication algorithm avoiding the redundant packing of A BLAS3中的SYMM，SYR2K，SYRK，TRMM和TRSM子程序并行算法均基于GEMM实现，参考文献＼[12＼]. 4 性能测试与分析 测试环境的硬件平台为8核FT1000微处理器，主频为1 GHz，双精度浮点性能8Gflops.操作系统为银河麒麟操作系统，编译器为gcc-4.5.1. 图6给出了BLAS1的双精度性能测试结果，固定向量的长度n不变，统一取为256 000 000.从图6可以看出，从1线程到8线程各BLAS1程序具有明显的加速效果，各程序在16或32线程时，达到最高性能，64线程时性能略有下降.造成32或64线程性能下降的主要原因是：1）BLAS1程序的计算访存比是1或2，受限于访存；2）BLAS1程序访存有空间局部性，没有时间局部性，也就是不存在数据的复用；3）FT1000微处理器的访存带宽在16或32线程时达到最大. No of threads 图6 BLAS1不同线程数性能测试结果（n=256 000 000） Fig.6 Computation performance for BLAS1 on difference number of threads（n=256 000 000） 图7给出了BLAS2的双精度性能测试结果，固定矩阵的长度n不变，统一取为16 000.从图7可以看出，浮点性能明显高于BLAS1的性能，性能最好的SYMV在64线程时达到最高性能，3.11Gflops/s，是峰值性能的38.8%，SYMV的计算访存比是6，存在空间和时间局部性.BLAS2中的其他子程序的计算访存比是3，只有空间局部性，不存在时间局部性. No of threads 图7 BLAS2不同线程数性能测试结果（n=16 000） Fig.7 Computation performance for BLAS2 on difference number of threads（n=16 000） BLAS3双精度浮点性能测试结果由图8给出，由于BLAS3的计算访存比大，对于阶为N的方阵，计算量为2N3，访存量为3N2，在64线程时获得最高性能，对不同的计算规模展开测试.从图8中可以看出，矩阵乘GEMM的最高性能达到6.91Gflops，是峰值性能的86.4%.SYMM，SYR2K，SYRK，TRMM和TRSM子程序的最高性能分别达到6.75，6.73，6.74，6.75和6.69Gflops，和矩阵乘GEMM的性能接近. 影响矩阵乘GEMM性能的主要因素包括：1）数据拷贝开销.为了充分发挥多级Cache的性能，采用了分块算法，分块子矩阵在内存中的存储是不连续的，为了减少对TLB冲突，需要将矩阵A和B的数据预先拷贝到一个连续的内存空间.2）数据延迟时间.每个子矩阵块相乘前，每个线程需要进行数据的填充，数据需要从内存中传输到L2 Cache，从L2 Cache到L1 Cache，从L1 Cache到寄存器，每个数据大概需要150拍，150拍过后，每拍可以取一个双精度数据.3）数据对存储带宽的竞争.64线程并行计算时，需要同时从内存中取数，数据带宽成为影响性能的因素.4）计算不能全覆盖数据的存取.FT1000微处理器中每核启动8线程，8线程共享一套硬件计算资源，靠多线程的切换来屏蔽访存延迟.这种情况主要存在每个循环展开启动阶段，8线程需要同时取数，此外循环展开计算完成后，需要把计算结果保存到矩阵C中，此时存取数据量和计算量处于相同量级.
M=N=K
图8 BLAS3不同计算规模性能测试结果（64线程）
Fig.8 Computation performance
of difference scale of matrixes on 64 threads
相对于其他X86微处理器，FT1000是多核多线体系结构，每8线程共享一个计算核，通过线程切换来屏蔽访存延迟，相对而言更难发挥其峰值性能，我们认为矩阵乘GEMM双精度能发挥峰值性能的86.4%已经是能达到的最好浮点性能.
5 结论与展望
提炼共性基础数值算法，研制高性能计算库，统一编程接口，是用户充分发挥微处理器峰值性能的重要手段.BLAS库是基础线性代数库，需要根据CPU的具体特点，设计高性能的算法和数据结构，经过高度的手工汇编优化，其中的BLAS3可得到接近峰值的浮点性能，满足用户的性能要求.
本文基于国防科学技术大学研制的多核多线FT1000微处理器，研制了基本线性代数汇编子程序库FitenBLAS，根据寄存器的数目，设计了向量与向量、矩阵与向量、矩阵与矩阵运算的多级循环展开方法，采用计算指令流水线调度、数据预取等通用优化技术，优化BLAS库串行程序.对于BLAS3子程序，设计了矩阵乘无冗余数据拷贝分块算法，采用指令重排、访存与计算的重叠、分块等技术优化矩阵乘子程序，基于矩阵乘子程序实现了其他BLAS3子程序.其核心子程序矩阵GEMM乘的双精度计算性能达到6.91Gflops，是峰值性能的86.4%.
下一步，面向短向量飞腾微处理器，研制支持向量优化运算的BLAS库，吸收并行编程框架底层调用的矩阵向量运算和稀疏线性数值算法，完善FitenBLAS库，支持更加广泛的数值模拟运算.
参考文献
[1] DONGARRA J. Basic linear algebra subprograms technical forum standard [J]. International Journal of High Performance Applications and Supercomputing， 2002， 16（1）： 1-111.
[2] Intel MKL homepage. http：//software.intel.com/zhcn/intelmkl/
[3 AMD ACML homepage. http：//developer.amd.com/cpu/ libraries/acml/
[4] IBM ESSL homepage. http：//www03.ibm.com/systems/software/essl/
[5] GotoBLAS homepage. http：//www.tacc.utexas.edu/taccprojects/gotoblas2
[6] GOTO K， VAN DE GEIJN R. Anatomy of highperformance matrix multiplication [J]. ACM Transactions on Mathematical Software， 2008， 34（3）： 1-25.
[7] ATLAS homepage. http：//mathatlas.sourceforge.net/
[8] WHALEY R， PETITET A， DONGARRA J. Automated empirical optimizations of software and the ATLAS project [J]. Parallel Computing， 2001， 27（1）： 3-35.
[9] OpenBLAS homepage. http：//xianyi.github.com/ OpenBLAS
[10]张先轶，王茜，张云泉. OpenBLAS：龙芯3A CPU的高性能BLAS库 [J]. 软件学报， 2011， 22（zk2）： 208-216.
ZHANG Xianyi， WANG Qian， ZHANG Yunquan. OpenBLAS： a high performance BLAS library on Loongson 3A CPU [J]. Journal of Software， 2011， 22（zk2）： 208-216. （In Chinese）
[11]About FT1000 processor， http：//www.nscctj.gov.cn/ resources/resource_1.asp， 2010-11-17.
非线性处理 篇3
脉冲压缩理论始于二战初期, 随着脉冲压缩技术的发展以及元器件性能的进一步提高, 目前, 脉冲压缩技术已经比较成熟, 并在现代雷达中得到了广泛的应用[1]。现在较常用的信号形式是线性调频和相位编码信号, 随着技术的进一步发展非线性调频信号在脉冲压缩处理领域也发挥着越来越重要的作用。
1 脉冲压缩处理技术
脉冲压缩处理就是通过信号处理技术来获得相应的脉压输出信号, 这样就能在不提高雷达发射功率的条件下, 得到较高的输出信号信噪比。按其实现方法可分为时域脉压和频域脉压处理法。
时域脉压处理系统其本质就是采用FIR滤波器来实现, 时域处理是直接对雷达回波信号进行卷积运算处理, 用公式表示为[2]:
式中, x (n) 为模拟回波信号经过A/D采样后得到的离散信号, N为采样点数;h (n) 为匹配滤波器的冲激响应;y (n) 为脉压输出信号。时域脉压实现的原理图如图1所示。
时域脉压对小的时宽信号能充分体现其优越性, 但随着脉冲宽度的增加, 信号处理部分的运算量变得很大。因此, 对大时宽信号的处理需要用频域脉冲压缩处理方法去完成。
频域脉冲压缩处理首先对接收到的经过采样得到的数字信号s (n) 做FFT变换, 获得其频谱函数S (ω) , 并对匹配滤波器的冲击响应h (n) 做FFT变换, 从而得到其频谱函数H (ω) , 将S (ω) 与匹配滤波器函数H (ω) 相乘, 然后再做IFFT变换, 最后得到脉压输出信号y (n) , 其表达式为:
频域处理法原理图如图2所示, 它在完成大时宽信号处理时其设备量并不会增加很大。这是由于采用了高效快速傅立叶变换 (FFT) 器件, 当采用FFT算法时, 能够使N点DFT的乘法的计算量由N2次降为 (N/2) log2N次[3]。因此, 对于大时宽脉压信号的处理, 采用频域脉压处理有其更大的优越性。随着DSP和FPGA技术的快速发展及其在快速傅立叶中的广泛应用, 采用频域脉冲压缩处理的运算能力将更强, 处理速度更快以及实时性将更高。
2 脉冲压缩信号形式
当前, 在工程中常用的脉冲压缩信号有线性调频信号和相位编码信号, 非线性调频信号在工程上应用较少。下面分别对这三种信号进行分析。
2.1 线性调频信号
线性调频信号作为雷达一种常用的发射波形, 当具有大的时宽带宽积时, 脉冲压缩系数D=BT, 线性调频信号的时宽T越大, 则雷达发射的总能量越大, 在干扰环境下, 雷达的作用距离就越远。同时, 线性调频信号的频宽B越大, 则压缩后的脉冲宽度越窄, 雷达距离分辨力也越高。因此, 选择线性调频信号作为一种雷达发射波形在脉压系数一定的条件下, 应该在雷达作用距离和分辨率之间进行折衷选择, 从而满足设计的需求。
线性调频信号矩形脉冲信号的复数表达式为[4]:
式中, T为脉冲宽度, 当时, rect (t) =1, 当t为其它值时, rect (t) =0。f0为载波频率, k为调频斜率。因此, 调制信号f (t) 的瞬时频率可以表示为:
当线性调频信号的时宽带宽积较小时, 就会导致其带内的菲涅尔纹波较大, 即使通过加权去平滑其矩形谱的边缘部分, 也不能使带内的纹波得到抑制。基于以上的原因, 小时宽带宽积信号加权结果与理论值相比有较大差距。所以, 在今后设计中可以通过谱修正技术来抑制带内纹波。
2.2 相位编码信号
相位编码信号就是把宽脉冲分成多个等宽度的窄脉冲, 它和线性调频信号最大的区别是线性调频信号是“连续型”信号, 而相位编码信号属于“离散型”信号[5]。同时, 相位编码信号对载波信号进行调制。
当前广泛应用的是二进制编码, 即二项码 (2PSK) , 它由“0”和“1”或“+1”和“-1”组成, 发射信号的相位在0°和180°之间交替变化, 其变化的规律由各脉冲是“1”还是“0”来决定。二项码的表达式为:
式 (5) 中, φ是0°或180°由数字信息是“1”或“0”决定。
对于小时宽带宽积信号, 相位编码信号旁瓣抑制比很高, 脉冲压缩处理效果较好, 同时, 信号波形变化灵活, 能够很好的实现雷达的低截获性能, 其不足之处是对多普勒很敏感, 对于回波信号存在多普勒频移时, 信号的脉冲压缩性能影响很大, 因此, 相位编码信号多用于多普勒频率范围变化较小的领域。
2.3 非线性调频信号
非线性调频信号是利用加权窗函数通过反设计得到的, 因而它的谱本身就具有加权窗的频谱特性。其最大的优点是不需要加权, 匹配滤波器输出就能得到很高的旁瓣抑制比, 这样就能避免失配引起信噪比的损失。非线性调频信号的缺点就是对多普勒频移很敏感, 但对于海上的低速目标非线性调频信号可以作为一种重要的雷达发射信号。
设非线性调频信号s (t) =a (t) exp[jθ (t) ]的频谱及匹配滤波器传输函数分别为S (ω) 和S* (ω) , 那么脉压输出信号y (t) 的频谱为[6]:
在上式中, 当把窗函数作为脉压输出信号y (t) 的频谱时, 就可以得到脉压输出信号, 由此可以推出非线性调频信号的调频函数为:f (t) =T-1 (f) (7)
式 (7) 中, T-1 (f) 为其群延迟函数T (f) 的反函数, 同时, 群延迟函数T (f) 和窗函数W (x) 又存在如下关系为:
式 (8) 中, m为常数。因此, 我们可以由信号的功率谱S (f) 2和选定的窗函数W (f) 来设计非线性调频信号。这样得到的非线性调频信号, 不但能避免失配损失, 同时, 还能获得很低的距离旁瓣。
3 线性调频信号脉冲压缩处理仿真分析
当前非线性调频和相位编码脉冲压缩处理已经有了较深入的研究, 线性调频信号作为一种较常用的脉冲压缩处理信号, 本文主要通过对线性调频信号脉冲压缩处理进行仿真分析。
通过对回波信号进行数字下变频以及脉冲压缩处理后的结果进行仿真实验, 这里时宽分别为T=5us, 20us和80us, 带宽B=4MHz, 信号中心频率f0=35MHz。采样频率fs=80MHz。固定作用距离为12公里的点目标进行仿真。不同时宽回波信号脉冲压缩处理结果的仿真图如下图3、图4和图5所示。
从雷达回波信号脉冲压缩处理仿真图可以看出, 随着时宽带宽积和脉压处理点数的增加, 经过处理后的回波信号的主瓣得到提高, 通过对仿真图的进一步分析和计算, 可以得出T=5us时的旁瓣抑制比为13.8297d B;T=20us时的旁瓣抑制比为13.4732d B;T=80us时的旁瓣抑制比为13.4056d B。因此, 对于固定距离目标点来说, 随着时宽带宽积的增加, 回波信号的分辨率有所降低, 这也是导致中远量程脉冲压缩雷达杂波大的一个重要原因。
随着发射信号时宽带宽积的增大, 其脉压处理增益分别为:T=5us时的处理增益为13.0103d B;T=20us时的处理增益为19.0309d B;T=80us时的处理增益为25.0515d B。同时, 处理增益还要考虑雷达的重复频率、天线转速和积累方式等因素。所以, 对于处理增益的提高应是以上多种因素的折中考虑。
4 总结
脉冲压缩处理技术作为一种重要的雷达信号处理方法, 已经获得了广泛的应用。本文通过对线性调频、非线性调频和相位编码信号进行详细的分析, 在总结其各自优缺点的基础上, 对常用的线性调频信号进行了仿真研究, 并对三种不同时宽带宽积信号进行了深入的对比研究, 得出了一些有价值的结论。下一步的工作就是将在实际工程中, 对该方法不断的优化, 从而获得更好的效果。
参考文献
[1]承德保.现代雷达反对抗技术[M].北京:航空工业出版社, 2002.
[2]吴太亮, 刘铮.基于FPGA的时域脉冲压缩器研究[J].制导与引信, 2007, (4) :46-50.
[3]胡广书.数字信号处理理论、算法与实现[M].北京:清华大学出版社, 2001.
[4]罗军辉, 罗勇江, 白义臣等.MATLAB7.0在数字信号处理中的应用[M].北京:机械工业出版社, 2005.
[5]杨静.相位编码信号的脉冲压缩及旁瓣抑制.南京:南京理工大学硕士学位论文, 2007.
非线性处理 篇4
在工业废水和生活废水的处理中, 生化法愈来愈受到人们的重视。生化处理是通过微生物的作用将有机污染物转化为无毒的气体 (如CO2、NO2、N2) 、液体 (如H2O) 及固体产物 (如生物污泥) 。微生物生长过程中的一个有效控制环节就是溶解氧 (DO) , 对DO的有效控制不仅影响到出水的水质而且还涉及到能耗。如何选择合适的空气流量, 使得水中的DO不仅能满足微生物生长的需要, 还可使其浓度保持在一定值 (一般为2mg/L) , 是DO控制的关键。DO在水中的传递模型表现出很强的非线性, 传统的PID控制对这种非线性、耦合性的模型难以实现有效的控制。本文针对DO的非线性特性提出了一种非线性的设计方法??反馈线性化, 验证表明它具有理想的控制效果。
2 适于控制系统的DO动态建模
2.1 活性污泥法污水处理流程
典型的污水处理过程见图1[1], 主要由生化反应池和沉淀池组成。其中:Q、Qr、Qw分别为进水、污泥回流和污泥排放的流量 (m3/h) ;X, Xr分别为混合液悬浮固体浓度 (MLSS) 和回流污泥浓度 (mg/L) ;S0, S, Sr分别为:进水、生化池和回流污泥的底物浓度 (mg/L) ;O, Or分别为生化池和回流污泥的溶解氧浓度 (mg/L) ;V为生化池有效容积 (m3) 。
2.2 DO动态建模
国际水质协会 (WA) 推出的活性污泥法数学模型 (ASM) 系列虽然可用于工艺流程的选择、建筑尺寸及运行参数的确定, 但是很难运用于控制系统的设计。为了得到简化模型, 将各种物质归并为微生物和有机物两类, 即溶解性物质 (有机底物) 和颗粒性物质 (微生物) , 并作如下假设: (1) 进水中的微生物浓度、DO浓度为零; (2) 沉淀池中没有微生物的代谢作用。这样溶解氧的传递模型[2]为
其中, Osat─与温度有关的溶解氧饱和浓度。
KLa─氧的传质系数, 一般可近似为空气流量QA的线性函数。
ROU─氧气吸收率, 与进水的底物浓度S0有关[3]
当空气流量比较大时, KLa是氧饱和度Ka的函数[4], 即KLa=KaQA。
将KLa代入 (1) 式, 有
其中x=O氧为状态变量, QA=U为控制变量, y为输出。从 (2) 式可以看出, 这是一个一阶非线性系统。
3 反馈线性化原理
反馈线性化是一种非线性控制设计的方法, 其核心思想就是把一个非线性系统转化为一个 (全部或部分) 线性系统, 以便可以应用线性系统的技巧。反馈线性化和普通 (雅可比) 线性化的根本区别在于, 它不是通过系统的线性逼近而是通过状态变换和反馈得到的[5]。对于单输入-单输出非线性系统, 有
进行输入输出线性化需要对y微分, 直至有如下形式:
而p是使LgLf h (x) ≠0的最小整数, 称为系统的相对阶。很显然, 由状态反馈控制律
可以使输入─输出映射简化为输出y和输入v之间的线性关系, 即
对于本文中的DO控制而言, 因系统是一单输入单输出非线性系统, 如果令e=yref-y, 其中yref为设定输入值, 则可以直接设计控制律:
u=1LgLpf-1h (.x) [-Lpfh (x) -k (yref-y) ] (7)
k值可以通过e=ke的特征方程确定。
4 DO非线性控制器的设计
由于控制目标是使生化池的DO控制在设定值, 在 (2) 式中取反馈量为O氧, 则按反馈线性化理论对输出y求导, 注意到y=x, 有
乘积项, 即消除非线性项, 取
从 (9) 式可以看出, 为了得到控制变量, 必须约定Osat≠x。事实上, 由于Osat为饱和DO的浓度, 而x为生化池溶解氧浓度, 二者不可能相等。控制时为使输出能跟随yref, 也就是使得鼓风机的能量消耗最低, 我们选取yref=2mg/L, 而对于k, 假定在稳态时有yref=y, 于是可以通过确定与设定值yref相关得特征方程的系数来确定变量x的变化即.e=ke。其中 (e=yref-y) , 移项有e.-ke=0, 选择常数k使得λ-k=0的根位于左半平面, 以保证系统稳定。这样, 利用原控制输入u来镇定原非线性系统 (2) 的问题, 已经转化为利用新控制律-k (yref-y) 来镇定系统的问题。如选取k=-3, 计算后可以得到控制方程为
如此可以得到非线性控制结构如图2所示。
5 系统仿真及结论
为了验证此方案的正确性, 采用Matlab语言对控制系统进行了仿真。根据实际经验可以取ROU=1.021kg/L?d, Ka=0.08/m3, Osat=6.25mg/L, QA=3m3/s, 这样得到的仿真曲线如图3所示。
图3中, 曲线1是经过反馈线性化处理的系统, 即将式 (10) 代入 (2) 仿真所得的动态响应, 它能够在很短的时间内跟踪输入且在达到系统所要求的溶解氧 (2mg/L) 后保持稳定;而曲线2是没有经过反馈线性化的系统, 即直接由式 (2) 仿真得到的动态曲线, 显见系统响应呈指数规律上升, 调整时间很长, 且存在稳态误差。由此不难得到结论:基于状态反馈线性化控制方案使得控制系统具有比传统控制更优越的动态性能。
参考文献
[1]ANDREWS J F.Dynamic models and control strategies forwastewater treatment processes[J].Water Research, 1974, 8 (3) :261-289
[2]杜树新.污水生化处理过程建模与控制[J].控制理论与应用, 2002, 19 (5) :660-666
[3]MARSILI—LIBEUi S.Adaptive estimation of bioactivities inthe activated sludge process[J].IEE Proc-D.1990, 137 (6) :349-356
[4][瑞]G.乌尔松, [澳]B.纽厄尔.污水处理系统的建模、诊断和控制[M].北京:化学工业出版社, 2004:33
非线性处理 篇5
1 非线性处理器设计
非线性自适应处理器的原型是自适应波形预测器, 它使用两组单频正弦信号作为线性组合器的输入信号, 不断地调整权向量来合成一个模拟的MSK信号。正弦信号的频率是载波的中心频, 两组信号的相位差为π2, 每一个采样时刻的输出信号都是它们的线性组合, 但是需要调整权值来使输出信号逼近MSK信号。通过非线性处理器对接收信号进行自适应处理, 其目的是将脉冲噪声从接收信号中去除, 并尽可能地还原MSK信号的波形, 减小信号失真。图1为本文非线性处理器的结构图。
1.1 自适应非线性组合器
自适应非线性组合器是信号处理器中的重要原件, 它可以将多个输入映射为若干个输出, 本文采用4个信号输入映射为2个输出的结构, 如图2所示。
在输入端由x1, x2, x3, x4组成输入向量, 可调权向量为w1, w2, 以及求和输出为y1, y2。用X和W的形式可表示为:
式中:T表示向量转置;k为抽样时间, 向量的元素均取自第k个抽样时间。由图所示映射关系得到输出如下:
随着抽样时间的变化, 权向量也不断进行变化:
式中:ΔWk为Wk+1与Wk权向量的改变量, 由自适应预估值与接收机接收值的误差通过LMS运算决定ΔWk取值, 使预估值更准确地接近MSK信号, ΔWk取值如下:
由ΔWk设置输出量yk的判决条件, 最后得出自适应预估值:
1.2 误差处理
由于接收机接收到的为混有脉冲噪声的MSK信号, 与自适应非线性组合器输出的预估值之间存在误差ε, 在高斯白噪声的影响下运用LMS计算可以使信号误差ε趋近于0。但当大幅度的脉冲噪声影响时, 误差ε会变得非常大, 达到MSK信号峰值的几倍甚至几十倍, 会造成误差收敛速度缓慢, 使预估的MSK信号产生失真。因此当遇到脉冲干扰时需对信号进行限幅滤波, 判断信号是否受到脉冲噪声干扰还要设置判决门限。
由式 (3) 和式 (4) 得:
式中:ε (k) =s (k-1) -y (k-1) 。设置一个脉冲判决门限B, 输出为n, 有:
用n与误差信号相乘得到ε′:
经过误差处理后信号消除了脉冲, 波形较为平滑, 相位依然连续, 如图3所示。
2 系统仿真及结果分析
仿真模拟在脉冲噪声干扰条件下通信系统误码率随信噪比改变的变化情况, 改进后的接收机在解调前对接收信号进行非线性处理, 消除脉冲噪声的影响。其结构如图4所示。
通信系统中脉冲噪声模块可设置参数SNR和Vd改变脉冲噪声的强度, SNR为信噪比, Vd为电压偏差, 定义电压偏差Vd为:
对Vd的认识给出了sn2的值, 而Vd的值在国际无线电咨询委员会 (CCIR) 报告中被估计出。
在仿真中分别设置Vd=2, 5, 10, 对比在经过信号处理器和未经过信号处理器通信系统的误码率。
图5 (a) 中脉冲噪声的干扰强度较弱, 不采用非线性处理器对接收机的误码率没有影响。当干扰增大, 如图5 (b) 所示, 采用非线性处理器的接收机误码率降低的速度明显快于不使用非线性处理器的接收机。而当脉冲噪声非常严重时, 如图5 (c) 所示, 没有非线性处理器的接收机无法有效降低误码率, 有非线性处理器的接收机可以有效降低误码率。
3 结语
改进后的通信接收机通过对接收信号进行非线性处理, 在接收机解调信号前用削波法削除脉冲噪声大幅度干扰, 采用自适应算法尽可能恢复通信信号。通过对系统的仿真及结果分析表明, 对于通信有严重干扰的脉冲噪声, 采用非线性处理器处理脉冲噪声的干扰能够有效提高系统的接收性能。
摘要：脉冲噪声由于其冲击能量和幅度较大, 会对通信系统接收机产生严重影响, 为了提高通信系统接收机的性能, 必须对混有脉冲噪声的接收信号进行处理, 需要改进接收机。以MSK信号为例, 通过对改进后的接收机性能仿真, 验证了采用非线性自适应处理器改进通信接收机的可行性。
关键词：通信接收机,脉冲噪声,非线性处理,误差处理
参考文献
[1]付贞, 温东, 孙晓磊.削波后的甚低频大气噪声中的通信接收性能[J].电讯技术, 2010 (1) :57-60.
[2]GOGOI A K, MALLIK RK, MAHANTA A, et al.Performanceof a coherent BPSK receiver in an impulsive noise environment[J].IEEE Transactions on Communication System, 1999, 15 (1) :471-475.
[3]沃特.无线电工程[M].北京:国防工业出版社, 1973.
[4]刘玉昌, 韩恩权, 陈妍.基于Matlab的VFL/LF通信系统仿真[J].电讯技术, 2004 (5) :140-143.
[5]肖明波, 杨松光, 许芳.通信系统仿真原理与无线电应用[M].北京:机械工业出版社, 2005.
非线性处理 篇6
物理现象中物理量之间存在着相互制约关系, 根据这种关系可总结验证物理定律、定理及物理公式.有些物理量之间存在着线性关系, 有些物理量之间存在着非线性关系.通常图象法是解决上述问题的一个比较直观的途径, 但图象法在解题过程中遇到非线性变化关系时, 就只能定性研究而难以定量计算.实际上在处理这样的一些非线性数据或问题时, 我们可以采用转化的思想, 可以将非线性变化关系转化为线性变化关系, 将曲线变化规律转化为直线变化规律, 这样有助于分析问题和定量求解.
二、非线性转化为线性的处理策略
在物理实验和解题时, 常常会遇到一些非线性的两个物理量, 若直接以这两个物理量为坐标, 所得的图象为曲线, 此时我们无法找到两个物理量之间的定量关系.这时我们就需要去寻找呈线性关系的自变量和因变量.再以这两个变量为坐标, 将图象由曲线转变为直线.那么如何去转化?其策略是:从物理规律和实验原理出发, 找出非线性的两个物理量的关系式, 并通过数学移项变形, 找到合适的线性关系, 然后利用其斜率和截距来求解.比如力学实验中通常需要测量的数据是位移s和时间t, 而在很多情景中, 这两个物理量成非线性关系, 下面我们就在具体的问题中来探讨其转化的策略.
三、非线性转化为线性策略的应用
1. 转化为图象
例1利用如图1所示的装置可测量滑块在斜面上运动的加速度.一斜面上安装有两个光电门, 其中光电门乙固定在斜面上靠近底端处, 光电门甲的位置可移动, 当一带有遮光片的滑块自斜面上滑下时, 与两个光电门都相连的计时器可以显示出遮光片从光电门甲至乙所用的时间t.改变光电门甲的位置进行多次测量, 每次都使滑块从同一点由静止开始下滑, 并用米尺测量甲、乙之间的距离s, 记下相应的t值;所得数据如下表所示.
完成下列填空和作图:
(1) 若滑块所受摩擦力为一常量, 滑块加速度的大小a、滑块经过光电门乙时的瞬时速度vt、测量值s和t四个物理量之间所满足的关系式是;
(2) 根据表中给出的数据, 在图2给出的坐标纸上画出s/ t-t图线;
(3) 由所画出的s/t-t图线, 得出滑块加速度的大小a= m/s2 (保留2位有效数字) .
解析:本题是利用光电门的装置在斜面上测物体运动的加速度, 实验中, 用米尺测量甲、乙之间的距离s, 并用计时器显示出遮光片从光电门甲至乙所用的时间t, 设加速度为a, 光电门乙时的瞬时速度vt, 则由运动学公式可得s=vtt-2/1at2, 由公式可知s和t成二次关系, 作s-t图象是曲线, 为了研究问题, 处理的策略是:作s/t-t, 即表达式:, 作出图象是直线, 成线性关系, 其中图象的斜率k=-21a, 截距b=vt.
2. 转化为s-t2图象
例2图3是“研究匀变速直线运动”实验中获得的一条纸带, O、A、B、C、D和E为纸带上六个计数点, 加速度大小用a表示.
(1) OD间的距离为cm.
(2) 图4是根据实验数据绘出的s-t2图线 (s为各计数点至同一起点的距离) , 斜率表示, 其大小为m/s2 (保留三位有效数字) .
解析:本题是利用打点计时器研究匀变速直线运动, 利用毫米刻度尺测出A、B、C、D、E到O点的距离s, 再根据打点计时器打出的计数点计算出距O点的时间t, 则由运动学公式可得s=1/2at2, 可知s和t成二次关系, 作s-t图象是曲线, 为了研究问题, 处理的策略是:作s-t2图象是直线.其中图象的斜率k=1/2a.
3. 转化为图象
例3现要通过实验验证机械能守恒定律.实验装置如图5所示:水平桌面上固定一倾斜的气垫导轨;导轨上A点处有一带长方形遮光片的滑块, 其总质量为M, 左端由跨过轻质光滑定滑轮的细绳与一质量为m的砝码相连;遮光片两条长边与导轨垂直;导轨上B点有一光电门, 可以测量遮光片经过光电门时的挡光时间t, 用d表示A点到导轨底端C点的距离, h表示A与C的高度差, b表示遮光片的宽度, s表示A、B两点间的距离, 将遮光片通过光电门的平均速度看做滑块通过B点时的瞬时速度.用g表示重力加速度.完成下列填空和作图.
(1) 若将滑块自A点由静止释放, 则在滑块从A运动至B的过程中, 滑块、遮光片与砝码组成的系统重力势能的减小量可表示为.动能的增加量可表示为.若在运动过程中机械能守恒, 与s的关系式为=.
(2) 多次改变光电门的位置, 每次均令滑块自同一点 (A点) 下滑, 测量相应的s与t值, 结果如下表所示:
以s为横坐标, 1/t2为纵坐标, 在答题卡上对应图6位置的坐标纸中描出第1和第5个数据点;根据5个数据点作直线, 求得该直线的斜率k=×104m-1·s-2 (保留3位有效数字) .由测得的h、d、b、M和m数值可以计算出1/t2-s直线的斜率k0, 将k和k0进行比较, 若其差值在试验允许的范围内, 则可认为此实验验证了机械能守恒定律.
解析:本题是利用光电门和气垫导轨的装置来验证机械能守恒定律.实验中s表示A、B两点间的距离导轨上, B点有一光电门, 可以测量遮光片经过光电门时的挡光时间t.多次改变光电门的位置, 测量相应的s与t值.由机械能守恒定律可得, 由公式可知:s和t成二次幂函数关系, 作s-t图象是曲线, 为了研究问题, 处理的策略是:, 作图象是直线, 成线性关系, 其中图象的斜率.
四、结束语
非线性处理 篇7
不论单频信号还是多分量的单频信号, 都是假定信号的频率是不随时间变化的, 这样的信号称作时不变信号。信号的傅立叶变换与时间无关, 因此, 傅立叶变换只适合于时不变的信号。然而, 现实物理世界中的绝大部分信号的频率都是随时间变化的。频率随时间变化的信号又称时变信号, 亦称非平稳信号。这类信号的傅立叶变换反映不出信号频率随时间变化的行为, 因此, 傅立叶变换只适合分析平稳信号, 而对频率随时间变化的非平稳信号, 它只能给出一个总的平均效果。时频分析的主要思想是将时间和频率结合起来对信号进行分析和处理, 它弥补了傅立叶变换缺少定位功能的缺点。时-频分析法就是用时间和频率的联合函数来描述信号在不同时间和频率的能量密度或强度。
线性调频信号 (LFM, linearfrequencymodulation) 是一种有着典型特征的基本的非平稳信号[1], 在雷达对目标进行探测、追踪和成像中有着广泛的应用。由于它的频率在脉冲宽度内按线性规律变化, 所以它常被用作衡量一种时频分析方法是否有效的试验信号。而且, 利用有效的时频分析方法准确的评估LFM信号, 对于己方探测敌方雷达辐射信号也具有重要意义。
瞬时频率是相对复信号或解析信号而言的, 它定义为信号相位随时间的变化率。瞬时频率估计方法只适用于单分量信号, 对于多分量信号, 它得到的是只是“平均频率”。本文在阐述几种不同时频分布基本原理的基础上, 针对LFM的时频 (时间—瞬时频率) 特征展开研究, 利用线性时频表示方法和二次型时频分布来分析和处理LFM信号, 并得出了一些有益的结论。
2 时频分析方法的理论基础
2.1 线性时频分析方法
对非平稳信号的第一类时频描述方法是线性时频描述, 也称为核函数分解法 (atomicdecomposition) 。它们将信号分解为基本成分之和的形式, 这些基本成分 (即核) 在时域和频域内具有很好的局域性质。
短时傅立叶变换使用一个很窄的窗函数取出信号, 并求其傅立叶变换。令g (t) 为一个时间宽度很短的窗函数, 它沿时间轴滑动。于是, 信号x (t) 的短时傅立叶变换定义为
式中*代表复数共轭。窗函数宽度越小, 则时域分辨率越好, 频率分辨率越差;窗函数越宽, 频率分辨率越好, 时间分辨率越差。如果取无穷长的窗函数g (t) =1, 则短时傅立叶变换退化为传统Fourier变换。
Gabor在1946年提出了信号时频展开的思想, 即Gabor展开:
式中g (t) 是窗函数, Cm, n是展开系数, m代表时域序号, n代表频域序号。这实际是用时频平面离散栅格上的点来表示一个一维的信号。
线性时频分析技术如STFT和Gabor扩展的主要缺陷在于所谓的“窗效应”, 它们使用固定的分析和综合窗函数, 其时间和频率分辨率受到窗长的约束, 不能得到同时优化, 短的时间持续窗函数能使时间分辨率得到提高, 但却导致了频率分辨率的降低;相反, 长的时间持续窗函数能使频率分辨率得到提高, 而时间分辨率却降低了。这是由于不确定性原理排除了在时域和频域同时具有高分辨率的可能。因而在需要同时具有高的时间和频率分辨率的情况下, 使用STFT和Gabor展开很难得到令人满意的结果。
2.2 二次型时频分布
当我们要用时频表示来描述或分析LFM等非平稳信号时, 二次型时频方法是一种更加直观和合理的方法。Wigner于1932年首先提出了Wigner分布的概念, 并把它用于量子力学领域。直到1948年, 首次由Ville把它应用到信号分析领域[6]。
1980年, Clas s e n对WVD的定义、性质等作了全面的讨论。信号的WVD就是它的时频能量密度, 其定义为
显然, Wx (t, w) 是关于t和w的二维函数。Wx (t, w) 有着时间边缘特性和频率边缘特性等一系列好的性质, 在一定程度上解决了短时傅立叶变换存在的问题, 而且, 具有明确的物理意义, 可被看作信号能量在时域和频域中的分布。它是应用甚为广泛的一种信号时频分析方法。1970年, Mark提出Wigner-Ville分布中最主要的缺陷——交叉干扰项的存在。令x (t) =x1 (t) +x2 (t) ,
与线性时频表示不同, 二次型时频分布不再服从线性叠加原理, 任何p0二次型时频分布都满足所谓的“二次叠加原理”:对于p分量信号z (t) 的二次时频分布pz (tw) , 每个信号分量都有一个自分量, 每一对信号分量都有一个对应的互分量, 即有p个信号项以及Cp2个两两组合的交叉项。信号分量越多, 交叉项就越严重。虽然交叉项是振荡的, 对于信号的能量并无贡献, 信号的能量完全包含在自项中, 但却在时频谱中引入了模糊, 这是WVD的一个重大缺陷。
对信号的WVD施加窗函数h (t) 后关于t的傅立叶变换称为伪Wigne r-Ville分布 (PWVD, ps e udo Wigne r-Ville dis tribution)
式中, H (w) 是窗函数h (t) 的傅立叶变换。加窗的结果等于未加窗时的WVD和窗函数的频谱在频率方向上的卷积, 使x (t) 的WVD在频率方向上得到平滑, 可在某种程度上降低交叉项干扰, 但同时也会降低WVD的频率分辨率。
减少交叉项的另一种方法是直接运用H (θ) w× (t, w-θ) dθ平滑函数对WVD分布进行平滑操作, 得到所谓的平滑Wigner-Ville分布 (SWVD, smoothed Wigner-Ville distribution)
其Wg (t, w) 中是一个平滑滤波器, tw**表示对时间和频率的二维卷积。当平滑滤波器Wg (t, w) 为函数g (t) 的Wigne r-Ville分布时, 该平滑WVD等效为信号x (t) 加窗函数g (t) 做短时傅立叶变换得到的谱 (下转第214页) 图。
实际上, 往往是对PWVD分布完成频域平滑的基础上进行时域平滑, 得到的分布称为平滑的伪Wigner-Ville分布 (SPWVD, s m oothe d ps e udo Wigne r-Ville dis tribution) , 即可分离平滑
式中, g (t) 和h (t) 是两个实的偶函数, 且g (0) =h (0) =1。
3 结论
当前, 随着高分辨成像雷达的发展, 线性调频信号的作用越来越重要。对线性调频信号的识别和分析处理成为了一个重要的关键问题。本文在阐述了几种不同的时频分析方法基本原理的基础上, 利用线性时频表示方法和二次型时频分布来分析LFM的时频特征, 得出的结论是二次型时频分布的时频分辨率比起时频表示方法的高。
参考文献
[1]Ville J.Theorie et applications de la notion de signal analytique[J].Cables et Transmission.1948.
[2]Classen T A C M, Mecklenbrauker.The Wigner distribution-part[J].Philips ResJ, 1980.
[3]于凤芹, 曹家麟.多Chirp成分信号双线性时频分布的交叉项分析[J].江南大学学报, 2004.
非线性处理 篇8
现代电网运行工况复杂、难以预测。为提高计算速度,单纯地采用简化模型或改进算法无法取得突破性进展[1],为此,并行计算成为了解决此类问题的有效途径之一。
大部分电力系统并行计算都是针对稀疏线性方程组进行的。对于粗粒度并行结构而言,将线性方程组的系 数矩阵划 分为分块 对角加边 形式(BBDF),是提高并行计算效率的有效方法。文献[2-3]先把电力系统分割为若干联系松散的区域,再对各区域进行并行计算。文献[4-5]研究了电力系统潮流及暂态稳定计算的并行化方法。文献[6-7]在个人计算机(PC)机群下利用上述算法实现了暂态稳定的并行计算,主要是对方程组的求解过程实现了多机的并行化,并获得了显著的加速效果。近年来,多核CPU架构成为 总体的发 展趋势,文献[8]对此进行了研究,并在4核CPU机器上获得了3倍左右的加速效果。
并行体系架构的革新,使得把多个核心集成于同一块芯片上、以数量来弥补速度的不足成为可能。多核CPU架构是其 中的一种,众核图形 处理器(GPU)架构是另 一种。统一 计算设备 架构(computeunifieddevicearchitecture,CUDA)的问世,使得GPU具备了通用计算的能力,为大型电力系统并行计 算提供了 新的途径。 文献 [9-10]在CPU-GPU混合架构下,用稳定的双共轭梯度法实现了大型电力系统暂态稳定的超实时仿真计算。文献 [11]利用GPU以及广义 最小残差 方法(GMRES)实现了电力系统暂态稳定仿真计算,且在系统规模较大时取得了一定的成效。文献[12]利用GPU实现了基于隐式梯形积分的机电暂态稳定仿真计算。
对于大规模稀疏线性方程组并行求解问题,若采用BBDF划分方法,则在众核GPU架构下,BBDF划分是否仍能够胜任细粒度的并行结构,需要进一步研究与探讨;若采用迭代法,虽然非常适合于众核GPU架构,但预处理的选择问题却很难解决;若采用直接法,则虽然不存在难收敛的问题,但在CPU-GPU混合架构下,如何处理方程组的稀疏性,又是一个值得研究的问题。
为了解决这2个问题,本文首先论证了先排序BBDF划分方法对于细粒度并行计算的可行性,其次,利用GPU的线程块(block)和线程(thread)并行特性,对具有粗/细粒度两层并行结构的线性方程组,提出了面向稀疏线性方程组的直接求解方法。
1众核 GPU 并行环境介绍
NVIDIA公司在2007年推出了CUDA,使GPU能够用于 图像处理 以外的一 般计算。其 与CPU主要的区别在于,CPU拥有较少的核心数目,但功能强大;而GPU每个核心的处理能力不强,但一个芯片上能汇集众多的核心,如图1所示。
CUDA技术的程序执行模型如图2所示[11],其从硬件上支持两层的并行结构:线程和线程块。各个线程块间可以并行地执行,但相互之间不能高效地进行通信,可以形成粗粒度的并行。而同一线程块内的不同线程可共享内存,实现线程间同步,形成细粒度的并行。本文针对电力系统仿真计算的特点,研究基于线程块和线程上的并行,形成线程块间的粗粒度并行和线程间的细粒度并行。图2中,线程块和线程均以2维予以展示,分别用Block(i,j)(i=0,1,2;j=0,1)和Thread(i,j)(i=0,1,2,3;j=0,1,2)表示。但实际上采用1维、2维均可。
2电力系统仿真计算并行方案
考虑到电力系统有其特殊的稀疏结构,因此,直接采用针对满阵的求解方法[13],抑或是基于稠密分块的多波前方法[14],都是不合适的。
本节以暂态稳定计算为例,对电力系统仿真计算的并行方案进行简要介绍后,给出了基于先排序的BBDF划分方法能够应用于粗/细粒度并行结构的依据,然后介绍 了所提方 法的原理 和技术实现细节。
2.1电力系统仿真计算的并行方案概述
图3为用于暂态稳定计算的数学模型框架[15]。
系统中的绝大部分动态元件都是相互独立的,仅通过电力网络连接形成统一的整体。系统的数学模型可以描述 为如下一 般形式的 微分—代数方 程组:
式中:x为微分方程组中描述系统动态特性的状态变量;y为系统的代数变量。
常规解法有交替求解法和联立求解法。本文采用的是交替迭代求解法。但无论采取何种策略,在求解过程中都会形成大规模稀疏线性方程组,研究重点都是一致的(具体分析见附录A)。
2.2并行方案分析
对于粗粒度并行计算,首先要将方程组的系数矩阵A划分为形如式(3)的BBDF[2,3,4,6,8]。
式(3)的物理意义在于:将原始节点集合划分为若干互不关联的子集合和一个规模较小的联系节点集合(见图4)。其中,联系节点集合对应BBDF中的边界块Acc,各子集合对应BBDF中的对角块Aii,i=1,2,…,k。
一般观点认为,这种方法不适合于众核GPU。其原因在于,为了能够充分利用GPU的众多核心,需要将节点集合分为数量足够多的子集合,这会导致联系节点集合的规模急剧增大。然而,对于先排序的BBDF划分,事实却并非如此,具体分析如下。
假设已采用了最小度最小深度(MDML)法[16]对所有节点进行了重排序。由于所需要的子集合较多,故采用文献[8]提及的分区方法进行集合划分。在多核CPU环境下,子集合的个数Q一般取为2至4。但对于GPU来说,所需要的Q较大。以3872条母线系统为例,测试不同分区数对联系节点集合的影响,如图5所示,具体数字参考附录B。
曲线符合对数特征。对该曲线进行对 数拟合后,得到以下结果:
式中:Q为分区数;Y为联系节点集合所占百分比。
式(4)对图5曲线拟合程度较好,联系节点集合规模增加趋向平缓。当分区数为64时,联系节点集合仅占2%。因此,对大规模系统使用先排序BBDF方法将节点集合划分为几十个子集,是可行的。
2.3暂态稳定并行计算的具体设计
根据前文的研究,并将其应用到电力系统的暂态稳定分析中,可确定如下的计算方案:由CPU进行逻辑复杂的前期处理,在形成大规模稀疏方程组后,交由GPU进行方程的求解。以交替求解法为例,算法流程如图6所示,联立求解法类似。
图6中解算网 络方程的 具体流程 见图7。图7(a)为CPU并行处理 的流程,作为参考。图7(b)为本文的研究,不同于CPU,其额外拥有两步“数据”处理的操作,其原因为CPU与GPU拥有不同的内存空间,需要对其数据进行同步处理。
对于非解方程部分,其并行化则相对较为简单。由于各元件对系统的影响都较为独立,可以简单地分配到众多的核心中进行求解,具体参考 附录C。因此,仿真中并行化的难点在于前文重点分析的大规模稀疏线性方程组的求解。本文在进行因子分解时,根据GPU的特点,制定出以下两层并行结构。
1)线程块层面的并行结构
根据式(3),每个的因子分解相互独立,其可以交由GPU的多个线程块进行处理。然而与文献[8]类似,此时由于每个线程块之间不能简单地交互数据,故暂时不能修正由每个Aii因子分解时对Acc的影响,而只修正了Aci和Aic。此过程利用了GPU的线程块特性,实现了对方程组求解的粗粒度并行。
2)线程层面的并行结构
在每个线程块的因子化过程中,可派生多个线程对其进行并行化,实现细粒度的并行。为简单起见,以下面5阶的矩阵a为例说明。
在第i步,包含归一化和秩一修正2个过程。
归一化过程:i=1时,需要对a12,a14,a15进行归一化,其相互是没有关联性的。在实际操作过程当中,可以开辟32条线程。当需要归一化的元素个数少于32个时,某些线程实质上不需要进行操作。当元素个数多于32个时,某些线程会串行地对不止一个元素进行归一化。下文开辟多线程也类似处理。
秩一修正过程:i=1时,实际上需要进行处理的元素是a22,a24,a25,a42,a44,a45,a52,a54,a55。这些元素处于归一化时a12,a14,a15所在的列(第2,4,5列)与第2,4,5行的相交位置。而其他位置的元素无论其是否非零元,均不会受到影响。而需要修正的这些元素之间并不存在依赖关系,也可以实现线程级别的并行。
当i=1处理完毕时,以同样的策略处理i=2,3,…,直到处理完毕所有的行。此过程利用了GPU的线程特性,实现了对方程组求解的细粒度并行。
上述过程实现了对式(3)中各个子块的因子化。当每个分块处理完毕后,需要对Acc进行修正:
Acc 中的每个元素的修正均是独立的,可以同样地通过开辟 多个线程 进行处理。Acc修正完毕后,需要对Acc,new进行因子化,过程与之前每个分块的方法一致,这里不再重复。
当因子化过程完成以后,进入前代回 代过程。与因子化过程一样,前代也可以形成两层的并行。
1)线程块层面的并行:式(3)的每个对角块交由GPU的一个线程块处理。
2)线程层面的并行:同样以式(5)为例,假设右端项为[b1,b2,b3,b4,b5]T。当i=1时,首先对b1进行归一化处理,b2,b4,b5可以并行地进行修正。当i=2时,对b2进行归一化处理,而后b4,b5可以并行地进行修正。需要注意的是,此过程的修正仅针对Aii,i=1,2,…,k;对于边界系统,需要额外处理。
上述过程完成后,需要研究[Ac1,Ac2,…,Ack]对右端项的影响。由于此矩阵中的第i行会且只会影响到右端项对应的bi,因此可为受到影响的每个bi 开辟一个线程,进行并行处理。回代过程与之类似。因子化过 程和前代 回代过程 分别见附 录D图D1和图D2。
为了能够充分利用GPU的性能,实现过程中应尽量保证CUDA核心处于满负荷计算状态。在图D1所示的过程中多线程块多线程区域,共开辟了32×k条线程同时执行。在2.2节中已经予以说明,k可以取为较大的数值,且保持联系节点集合的规模在可以 接受的范 围内。若k>30,则32k>960,已经超过了CUDA核心的数目(本文测试环境GTX570为480个核心),足以给各个核心分配一定的计算量。这也是即便要付出增加联系系统规模的代价下,也要划分众多子网的原因。对于具体的某个CUDA核心,在某个分区不需要它的时候,可以由另外的分区 进行调用,以消除等 待时间。由于GPU可以在不同线程间几乎零延时地进行切换(即一个CUDA核心可以在不同线程间来回切换,而不需要额外的时钟周期),因此较多的线程可以在一定程度上均衡CUDA核心的计算量。
综上所述,本文结合GPU运行环境实现了两层并行结构的方法,实现了大规模稀疏线性方程组在GPU中的求解,并应用到电力系统暂态稳定计算中。
3算例测试与分析
3.1仿真环境及技术指标
本文使用VisualStudio2012作为开发 环境,nVidiaCUDAToolkit5.5作为GPU开发工具,采用C++及CUDA语言分别编写了运行于CPU上串行暂态计算程序和基于GPU的并行暂态计算程序。两种程序采用相同算法框架,并且都应用稀疏存储技术,区别在于线性方程组求解过程是在CPU还是GPU上进行。硬 件平台如 下:CPU为intelCorei52320/3 GHz,内存8 GB;GPU为nVidiaGeforceGTX570,运行环境MicrosoftWindows8.1x64。
主要测试算例参数见表1,其中7800案例为利用IEEE300系统倍增26倍形成,7744系统为利用3872系统倍增形 成。初始潮流 收敛且系 统稳定。微分方程采用改进欧拉法求解,代数方程采用高斯消元法直接求解,整体采用交替求解法。故障设定为某母线出线端发生三相接地短路,延时切除故障。计算步长为0.005s,数据采用双精度。
本文主要考察程序进行并行化后的效率提升,并用式(7)所示的加速比技术指标来表示。
3.2计算结果与分析
3.2.1计算参数设置及结果
由于在GPU上可调参数较多,例如分区个数、GPU上开辟的线程块数目、每个线程块内的线程数等均可调整,因此本文首先通过一系列的测试,确定参数的设置,具体设置见附录D。对于多线程块多线程区域每个线程块开辟32条线程、单线程块多线程区域开辟1024条线程能 够获得最 优的加速 效果。在此基础上,表2显示了不同分区数对结果的影响。
由于本文的设定为一个分区(即一个对角块)由GPU的一个线程块进行处理,因此表中的结果也即为线程块个数的设定。对多个案例测试的经验结果表明,分区个数(N为方程的总维数)时能够获得最优的效果。此时,每个分区的维数约为在上述设置下,联系系统的规模与各个子网的规模相差不大,前者约为单个后者的80%~90%。测试中发现,联系系统的运算时间(包括各个子网对联系系统的修正过程)大约为全部子网运算时间的50% ~60%,比重较大,因此再继续增大联系系统的规模就会影响并行效率。
根据以上结果,本文对测试案例进行分区时主要采取的设置参数如表3所示。
注:指导阈值为文献[8]分区方法的参数设置,即前驱节点数超过此值的节点,直接纳入到联系系统。
本文使用微软的VC编译器和英伟达的NVCC编译器分别进行了测试,结果如表4和表5所示。表4中,“GPU”指图7中由GPU进行的部分内容,其中包括了数据同步部分(即表5中的通信部分);“CPU”表示图7(a)部分内容,形成了与GPU部分的对照程序;“共同”则表示图6中为除网络方程求解外的其他内容,GPU程序和CPU对照程序均使用此部分内容,在计算加速比时将扣除此共同时间后再予以比较。
注:1四个案例仿真的总时长不一样,不可横向比较;2 共同及CPU 部分由 VC编译,GPU 部分由 NVCC编译。
注:GPU 部分扣除了通信时间。
从表4的结果可以看出,随着系统规模的增大,加速比的效率也随之增加,目前规模下可达3.84;在同等矩阵规模下,已经明显优于使用广义最小残差迭代法(GMRES)结合GPU所带来的加速比[11]。表5给出了各个案例下数据同步(即CPU、GPU内存间通信)所占据时间的比例,一般在50%~60%左右。若扣除此部分时间,只考察GPU和CPU在解稀疏方程组的加速效果,加速比则能够达到8以上。但在应用过程中,同步时间是必要的,表4的数据更具参考价值。
3.2.2计算结果分析
结合理论方法以及测试结果,具体分析如下所示。
1)经过预处理后,矩阵形成如式 (3)所示的形式;为每个对角块分配一个线程块,同时为每个线程块分配了32个线程,形成了两层的并行结构。结合稀疏技术后,每个CUDA核心仍可 以满负荷 地运行,充分发挥了GPU的性能,并获得了 良好的效果。
而在联系节点集合的处理上,由于其相对比较稠密,根据文献 [13],也能够有 效发挥GPU的作用。
2)本文最大的 加速比为3.84,小于文献 [1213]中的加速比。文献[12]使用判断某处是否为非零元素进行处理,相当于系数矩阵为满阵,本文采用此方法进行测试,结果如表6所示。表中时间为不同运行环境下,执行一次暂态稳定计算所耗时长。
从表6可以看出,结合满阵技术,在矩阵规模较小时已经能够有较高的加速比,但计算时间过长,远劣于应用稀疏技术时的结果,已经失去了加速的意义。此外,当使用满阵时,用于CPU和GPU之间的数据同步时间将会大大地增加,表6及文献[12]结果均表明,超过50%的解方程组时间将用于数据同步。当系统规模较大时,传输时间将增大到不可忍受的程度,且GPU将不足以一次性存下所有的数据。满阵的研究成果不可直接用于电力系统计算分析,因此也造成了暂态稳定并行计算的加速比比满阵研究结果稍差。
3)在硬件环境上,本文的GTX570对双精度的支持并不好。表7给出了CPU和GPU硬件特性比较[11]。
对于CPU,从单精度到双精度,浮点运算速度只下降了一半,但对于GTX570来说,下降为原来的1/8。如果GTX570的双精度性能能够达到与CPU一样,为单精度的1/2,理论上本文的加速比可在现有基础上再提升2倍以上,达到7~8倍的加速比。事实上,已经有高端的显卡Quadro6000能够做到这一点。
4)在理论上,直接法求解方程组的方法也制约了GPU的性能,在许多步 骤上出现 了“三角形 问题”,例如:假设需要对式(8)进行前代。
完成对a11的归一化后,可开辟4个线程对第2至5行进行秩一修正;但对a44归一化后,仅能开辟1个线程对第5行进行修正,即某些时刻只有少数的线程进行计算。由于单个CUDA核心的性能是远差于单个的CPU核心,因此这将会造成加速比的严重下降,拖累了整体的加速效果。此问题覆盖了因子分解和前代回代的大部分步骤,造成加速比低于理论上限。
本文的研究结果表明,尽管存在一些瓶颈,但利用了GPU实现求解电力系统中线性方程组的直接法后,在电力系统暂态稳定计算中依然能够有较好的加速效果,目前的测试环境下能够达到3~4倍的加速。当技术发展到目前高端水平的显卡能够普及时,加速比可望达到7~8倍。
4结语
本文在分析电力系统计算中大规模线性方程组的稀疏特点后,提出了稀疏线性方程组直接求解法在GPU上的并行化实现方法,形成了粗粒度和细粒度的两层并行结构,分别利用GPU的线程块应用到电力系统计算中。通过算例分析验证了方法的有效性。分析结果表明,目前在一般配置的电脑上已有3~4倍的加速效果,优于GMRES等迭代法在GPU上的效率提升,且远低于理论的上限值,存在着较大的改进空间;当矩阵规模扩大时,加速比提升更为显著,具有良好的发展前景。
向上海交通大学国家能源智能电网(上海)研发中心提供良好的科研条件表示感谢。
附录见本 刊网络版 (http://www.aeps-info.com/aeps/ch/index.aspx)。
摘要：针对电力系统大规模线性方程组的稀疏特点,提出了基于图形处理器(GPU)的直接求解方法。该方法首先利用基于先排序的分块对角加边形式(BBDF)划分方法对方程组系数矩阵进行分割,形成具有粗粒度和细粒度两层并行结构的线性方程组,然后利用GPU的线程块和线程并行特性对其分别予以求解。将上述方法应用到电力系统暂态稳定计算中,并对其加速效果进行了测试。测试结果表明,在目前普及的设备上,所提方法可获得3~4倍的加速比;在高端设备上,能够获得7~8倍的加速比。