有限元逆算法（精选7篇）
有限元逆算法 篇1
0 引言
钣金零件是指利用金属的塑性变形, 通过冲压、压弯、滚弯、落压等方法, 将板材加工成需要的零件。钣金零件的制造在飞机及汽车工业中占有重要地位。据统计, 一架中型飞机上, 钣金零件的数量达六万多件, 制造工时占全机总工时的12%左右。
大多数钣金零件都是从平板原材料加工而来。出于形状匹配和加工限制等原因, 很明显的, 板材上只有特定形状部分会成为最终零件, 多余的部分会被丢弃。对于特定加工方法, 一个零件在板材上所利用的区域形状应该是固定的, 如果能够根据零件反推出对应的板材利用区域, 就可以以此调整加工方法, 提高材料利用率, 减少加工时间, 从而在降低成本的同时提高生产效率。这种反推方法统称为展开放样。
为了确保成形件的最终品质, 必须得到几何尺寸精确, 形状准确的展开图。传统钣金展开方法有:作图法, 计算法, 系数法, 计算机辅佐设计法[1]。这些方法都是根据制件的已知尺寸和几何条件, 通过解析计算, 算出绘制展开图所需的几何尺寸, 并按其绘出展开图。然而, 这些方法没有考虑到钣金加工成形过程中的力学因素, 回弹、畸变和残余应力等对钣金零件几何尺寸和形状的影响, 导致展开精度过低, 不能满足加工工作的需要。现代机械设计需要一种稳定, 高速, 高品质的钣金展开方法。同时, 还需要有可以方便快捷应用这种方法处理钣金展开的工具, 使之可以真正投入实用。
1 寻求解决方案
解决方案应包含两个方面:1) 包含力学因素考虑的钣金展开方法, 2) 采用了该种方法, 同时符合工业实际工作流模式的设计工具。
首先确定方法。从力学角度上, 钣金成形是一个包括几何流体、材料物理和边界摩擦的复杂非线性问题。由于工件形状, 材料属性等的任意性, 无法在宏观层次对任意工件给出统一的解算方案, 这对计算机软件来说是不可接受的。在必须考虑力学因素的情况下, 利用微分原理的有限元法是唯一的应对方案。有限元法将任何工件都视为最小力学单元的组合体。工件的受力变形都是这些最小单元上受力变形的集合。由于最小单元的力学特性简单, 有限元法就可以用只针对最小单元的算法, 依赖计算机高速计算能力去解决任意的工件受力, 运动与变形问题。
目前, 应用于板材分析的有限元法主要有基于流动理论的增量有限元和基于形变理论的全量有限元[2,3]。增量有限元计算时要求给定所有的冲压成形工艺参数, 计算精度高但所需计算时间过长。全量法根据已知的零件的最终几何形状和初始平板料的厚度, 通过对整个成形过程引入一系列简化假设 (基于全量理论的本构方程、等效外力等) 将零件最终构形和初始平板料进行对比就能简单而快速地求出零件最终构形上的厚度、应变和应力以及初始板料的形状。全量法在精度足够的情况下提供了远比增量法快得多的速度, 更加符合工业要求, 故采用全量法进行钣金展开分析计算。
良好的算法还需要有方便应用的软件工具才能真正发挥作用。当前已经有许多商品化的钣金成形模拟软件投入使用, 如:DYNA3D, DYNAFORM, PAMSTAMP, ABAQUS/EXPLICIT等。但是这些有限元软件只专注于计算分析, 而没有足够强大的建模功能。实际应用时, 需要使用另外的CAD软件建模, 将模型导入这些软件进行分析, 计算完的结果还需要再次导入回CAD软件中去。该流程中有工业价值的实际上仅有建模和分析两项, 其他的操作不仅耗费时间, 还增加了出错可能。另外, 作为专门的有限元软件, 其算法设定宏大复杂, 对于钣金展开这样专门的单项工作来说过于冗长, 操作难度很高, 加重工作人员的负担。直接套用现成有限元软件显然是不合适的, 最好的方法的将分析算法直接集成于CAD软件中。在实现高精度钣金展开分析目标的同时, 对现有工作流的影响也可以做到最小化。在此考虑到工作主要面向航空工业, 所以选择CATIA做为算法集成对象。
CATIA 是法国达索系统 (Dassault Systemes) 公司开发的具有界面友好、功能强大、易于操作等特点的三维CAD/CAM 软件[5], 广泛应用于汽车、航空航天、造船等领域。现根据用户设计开发的需求, 在CATIA平台上进行二次开发, 将钣金展开有限元逆算法嵌入到CATIA环境中, 结合零件设计和原始板料展开功能, 提供交互操作, 适合大批量单项工作场合, 为用户提供更大的便利。
2 有限元逆算法在CATIA环境下的实现
2.1 有限元逆算法基本方程
假定板料的弹塑性大变形是塑性变形体积不可压缩, 并且其变形过程是比例加载的。并进一步将板料和模具之间的相互作用力和接触条件简单地通过等效外力加以考虑。在已知的变形终态构形上采用虚功原理建立如下平衡方程:
其中:{u*}——虚位移;
{ε*}——由虚位移所产生的虚应变;
{σ}——柯西应力;
{f}——等效外力。
设节点虚位移为:U*, V*, W*, 则:
对于三角形单元[5,6]:〈U*n〉=〈U*iV*iW*i〉 (i=1, 2, 3) 且W*i为已知。
对于这种非线性问题, 一般采用带有收敛因子的newton—raphson (N-R) 迭代算法求解。假设整体残余应力为外力为{R (U) }, 且为外力与内力的差值:
{R (U) }=∑ ({Feext}-{Feint})
假设在第i个迭代步, 其近似的初始解{Ui}, 可以得到:
{R (Ui) }={Fext (Ui) }-{Fint (Ui) }≠{0}
根据N-R法, 将{R (Ui) }在{Ui}处泰勒展开得:
根据有限元逆算法实际情况, 选择位移准则和不平衡 (残余应力) 准则来进行收敛性判断:
undefined
式中, nd为节点总数, ε0为给定的收敛判断因子。当迭代满足收敛条件时得到的点的序列便是毛料的初始位置。
2.2 软件架构
项目依靠如下几个模块实现预定功能:模型处理、接口程序、有限元逆运算和结果后处理。各个模块通过CATIA及程序提供的特定接口, 互相传递数据进行通信, 协调合作实现在CATIA环境下的钣金展开功能。
模型处理主要功能是对钣金零件的part文件进行网格划分, 并生成划分后的网格信息。首先将要划分的零件表面作为Publication进行发布, 为后续操作提供曲面数据。由于网格划分要在Analysis模块进行, 因此新建将原part模型作为Linked Documents特征的Analysis文件。进入该模块后, 根据已发布的Publication表面添加AnalysisMeshParts属性, 并按用户要求通过修改其中的MSHPartBasicSurf, GlobalMethod, GlobalSize等参数, 最后通过Update方法完成网格划分。
接口程序提供CATIA内部数据与外部程序的通信接口。网格数据通过该接口输出至有限元计算模块, 计算完成的零件毛坯模型再通过这个接口返回CATIA。网格数据包括单元数组、节点总数、节点坐标等信息。由于在有限元计算中钣金件默认冲压方向为负z轴方向, 因此在进行运算前要调整钣金法向与z轴平行。
有限元运算模块处理从CATIA中得到的零件网格信息, 进行展开逆运算得到原始毛坯零件的形状, 然后再发送回CATIA。采用的有限元法为上文介绍的全量有限元逆运算法。
最后, 结果后处理模块将上述处理得到的原始几何结构信息重新导入回CATIA环境, 使用CreateSketch及CreatePad等接口方法在CATIA中生成Part格式的原始板料文件。并可以在此基础上完成测量板料平板上孔到边沿距离, 及查看源文件是否已被展开, 避免重复计算等后处理功能。
3 算例应用
为检验程序的可用性和实用性, 测试如图1所示的标准矩形盒零件, 经过网格划分、有限元运算等一系列操作, 可以得到如图2所示的原始板料零件。通过分析原始板料文件, 可以发现有限元逆运算结果具有较高的精度, 且运行时间短, 完全符合钣金生产实践的需要。
4 结论
通过结合有限元逆算法和基于CAA的CATIA二次开发技术, 实现了CATIA环境下嵌入钣金有限元展开功能。并通过实例应用得出以下结论:
1) 因为CATIA本身机构的层次化, 利用CAA进行二次开发具有很大的难度和复杂性, 但由于其具有直接调用CATIA底层的API等功能, 与其他二次开发方法相比CAA能实现更大功能和更深处次的开发。
2) 尽管有限元逆算法模拟精度与増量法相比略低, 但它是一种高效的数值模拟工具, 并随着算法的不断完善, 其模拟精度可进一步提高。
3) 软件实现了钣金展开模块与CATIA环境的无缝连接, 方便了钣金设计人员从钣金设计、展开到验证的一体化流程工作。
摘要：传统几何钣金展开方法缺乏对成形过程中力学因素的考虑, 导致对复杂不可展曲面展开结果精度不足。针对这种问题, 将有限元技术引入展开计算中, 提出一种基于有限元逆算法的钣金展开方案。通过对钣金成形过程中应力、应变的分析, 构建有限元控制方程并对其求解, 之后使用CAA API编程将其集成到CATIA环境中。最后通过对一不可展矩形盒零件的展开, 验证了该方法的有效性和通用性。
关键词：钣金展开,有限元逆运算,CATIA二次开发
参考文献
[1]毛昕, 黄英, 宋萃娥.钣金实用技术问答[M].北京:北京出版社, 2000.
[2]鲍益东, 杜国康, 陈文亮, 等.冲压成形模拟中有限元方程组求解算法[J].南京航空航天大学学报, 2009.41 (5) :606-610.
[3]杨晨, 谢兰生, 童国权.基于有限元逆算法的钣金展开设计[J].机械制造与研究, 2003. (3) :13-16.
[4]周仙娥, 鲁墨威, 赵海星.基于CAA的CATIA二次开发的研究[J].IT技术论坛, 2008. (36) :73-75.
[5]Kim SH, Huh H.Finite element inverse analysis for the design ofimtermediate dies in multi-stage deep-drawing processes withlarge aspect ratio[J].Journal of Materials Processing technology.2001, 113:779-785.
[6]Batoz J L, Guo YQ, Mercier F.The inverse approach with simple tri-angular shell elements for large strain predictions of sheet metalforming parts[J].Engineering Computations.1998, 15 (7) :864-892.
大型广义逆反演中求逆算法的比较 篇2
关键词：广义逆反演,共轭梯度法,不完全LU分解,不完全乔列斯基分解
0 引言
大地测量和地球物理中的反演问题都是从数据中推导模型参数,其主要目的就是在一定的理论基础上根据数据来重构模型。对于一个实际反演问题,经过一系列处理后得到一组线性方程,只要其系数矩阵满秩,就可以通过严格的数学方法来求解。但很多情况下,因为其系数矩阵非常之大而不能直接求逆,这时我们最为常见的就是用最速下降法(Steepest Descent,SD)和共轭梯度法(Conjugate Gradient,CG)等迭代方法来求解。然而,这两种方法的收敛速度却时常难以满足实际大地测量或地球物理反演问题求解的需要。为此,本文提出利用不完全LU分解(Incomplete LU Decomposition Conjugate Gradient,ILUCG)及不完全乔列斯基分解(Incomplete Cholesky Decomposition Conjugate Gradient,ICCG)的预条件共扼梯度法来优化模型求解问题。下面首先介绍线性方程的共轭梯度迭代解法,然后引出两种预条件解法并简介其特性,最后作分析比较。
1 反演问题求解模型
实际问题中各种模型都是以多元的连续函数表示的,对一个问题做反演,就得先把连续函数进行离散化(参数化),线性化(泰勒级数展开)[1],之后就可以得到如下线性方程组:
其中d为数据向量,m为模型参数,A为系数矩阵。一般的地球物理反演问题中,有时对部分参数数据有矛盾的信息,而对另一部分参数又缺少数据信息,所以(1)式一般是不能直接求解的。这时需要用最小二乘或者最小范数,再或最小阻尼等准则[1,2]来约束(1)式,从而得到能求得唯一解的方程。
当数据对参数有矛盾信息时(A阵列满秩,独立数据个数大于未知参数个数),运用最小二乘准则得到:
当有部分参数缺少必要的数据信息时(A阵行满秩,独立数据个数小于未知参数个数),运用最小范数准则得到:
当数据对部分参数有矛盾的信息,而对另一部分参数缺少数据信息(A既非行满秩也非列满秩),运用最小阻尼方法得到:
γ就是阻尼系数,它是一个经验值,当γ为零时它就是最小二乘解。一旦建立了以上的(2)、(3)、(4)式中某一式,就可以运用数学方法来求解反演问题了。但因为在地球物理或大地测量反演问题中,系数阵A的阶数往往都很大,所以要避免对右边项的直接求逆,一般都是运用迭代方法进行求解。其中下面介绍的共轭梯度法(CG)是目前最为常用的迭代求解方法。
1.1 共轭梯度法
下面以(4)式为例介绍一般的共轭梯度迭代解法。先设则有：
为了避免对B求逆,就要利用其他方法来求解(5)式、共轭梯度法就是基于变分原理的迭代方法,先建立二次函数:
括号内的运算(x,y)表示两向量的内积,经证明当B为对称正定阵时,使得二次函数F(x)取得最小值的x即是(5)式的解[3]。
显然经过处理后(5)式中的为对称正定矩阵,现在的问题就转化为怎样来求解(6)式最小值点了,下面来看共轭梯度法(CG)的做法:
先设:r(0)=y-Bx(0),它就是x(0)的负梯度方向,另外再引入一个参数q,设为:q(0)=r(0)
第1步:α(k)=(r(k-1),r(k-1))/(q(k),Bq(k))
前面两步是根据函数极值定理来求得一个α,用它来对初值x(n-1)进行改正,使得它沿x(n-1)的负梯度方向下降。这也就是最速下降法的迭代原理。但是负梯度方向只是从某直线方向上来搜索最佳收敛方向。而共轭梯度法是在一个面上来搜索最佳收敛方向,所以它具有更佳的收敛性。以下几步就是在最速下降法的基础上对q向量做改进:
第3步是计算新的负梯度方向,因为当r(k)较小时,由于舍入误差的影响,实际计算得到的r(k)会偏离最速下降的方向。所以第4、5步就是利用几何方法来求得一个改进量β,并对q做改进[3],以保证它的最速下降。
1.2 预条件的共轭梯度法
虽然共轭梯度法的收敛速度比最速下降法改进了不少,但是其收敛速度还是受到B矩阵条件数的制约[3,4]。而在大地测量或地球物理问题中的数据矩阵一般都是大型稀疏矩阵,并且其条.件数都非常大。即使经过正则化处理条件数会有所降低,但其降幅是有限的。巨大的条件数必然使得迭代算法的收敛速度变得非常慢。为了加速收敛,有必要对算法进行预条件处理。所谓预条件处理就是降低条件数。自上世纪70年代以来已经发展了多种预分解方法[4,5,6],而预条件法的优劣直接影响到迭代收敛性的快慢,目前预条件方法也被运用于大地测量及地球物理的许多领域[7,8,9]。下面我们将介绍两种预条件方法,然后再提出其中一个很少被提及而又非常重要的因素。
因为B矩阵是稀疏矩阵,可以将B进行一种近似分解—不完全LU分解[3]:B=LU,其中L为单位下三角矩阵,U为上三角矩阵。为了便于理解,先把(5)式变形为:
Bx=y (12)
这样对于(12)式,两边同乘L-1,得
L-1 BU-1 Ux=L-1y (13)
并设:,X=Ux,g=L-1 y,所以(13)式可化简为:
这样一来因为,所以F的条件数接近于最小值1。可想而知这时的迭代速度会非常之快,下面的实验也证实了这点。而其算法就是对(14)式用共轭梯度法求解,在迭代求得它的解后再对它做一个U-1求逆算子就可以得到(12)式的解。下面是ILUCG算法流程:
因为B阵是对称正定矩阵,这时还可以对其进行另一种具有更好性质的近似分解—不完全乔列斯基分解[3]。乔列斯基分解是将矩阵分解为一个上三角矩阵及一个下三角矩阵,并且两个矩阵是互为转置矩阵:。那么,F的条件数接近于1,其算法只要将上面的U改为L即可。虽然算法上这两种方法很相似,但是实际运算证明,基于不完全乔列斯基分解的预条件算子的预优性质是与不完全LU分解算子是有一定差别的。
最后我们来看看两种预条件算子的可靠性,在对矩阵进行近似分解时为了减小存储空间及加快运算速度必然会设置一个非负的阈值[6](drop tolerance)来抛弃一些元素。M.Benzi等[10]提出矩阵在任何阈值下的近似分解存在并能用于预条件算子的充分条件是矩阵是对称的H型矩阵[11]。对一般的对称正定矩阵如果抛弃阈值设置不当,对其近似分解时可能出现主元为零,而使得分解失败;或主元为负数,使得分解得到的矩阵为非正定矩阵,从而不能成为预条件算子,这也是很多情况下求逆算法不能成功收敛的原因。下面构造一个简单的非H矩阵来验证这一点,当抛弃阈值大于0.05时对它进行不完全乔列斯基分解时就会失败,
2 试验
本节利用M atlab编写了上述各种算法的程序并比较它们的特性。这里借鉴C.J.LIN和R.S-aigal的做法[12]构造B矩阵,即由M atlab中的sprandsym(n,density,rc,kind)函数生成稀疏矩阵。该函数可以控制矩阵的稀疏率(density),条件数(rc)-1及正定性(kind=1时为正定)。由于硬件条件限制,这里设计的矩阵阶数是逐渐增高,从而来推估各算法对超大型矩阵的性质。计算中迭代终止的阀值为10e-6,且为了保证矩阵能近似分解,这里设定了较小的抛弃阀值0.001。表1是几种算法比较的数据结果(注:n为迭代次数,t为解算所耗的CPU时间,N为矩阵的阶数,Cond为迭代因子的条件数。所有程序都是基于Matlab7.1,P4-
2.8 GH,内存736MB的硬件环境下进行的)。
各种算法比较试验结果表1
虽然用这几种算法最后得到的结果精度都非常好,但运算效率不尽相同。由图1可以看出,就运算时间而言,ICCG所消耗的CPU时间一直是最少的。ILUCG所耗CPU时间随着矩阵阶数的增加而迅速增加,而ICCG所耗时间的增加趋势相对较平缓。这跟它们的预分解算法及分解后的相似性有关,一般高精度的不完全LU分解会减弱原矩阵的稀疏性,需要花费较多的构造时间与存储空间。不完全乔列斯基分解产生的是互相转置两个矩阵,比起LU分解来其分解速度肯定具有优势。
(横坐标为矩阵阶数,纵坐标为迭代次数)
由图2可以发现在迭代次数上,ILUCG是非常具有诱惑力的,这样随着计算机硬件环境不断升级这种方法会更能体现其价值,所以它主要被运用于并行计算中,而ICCG的收敛迭代次数也是可以接受的,而CG法由于没有经过预条件处理,所需的迭代次数远远大于另外两者,这是在预料之中的情况。
3 结论
有限元逆算法 篇3
随着电力电子技术飞速发展,各种敏感、变频设备在电力系统中的广泛应用,电压暂降问题现已成为威胁设备安全和正常稳定运行的严重动态电能质量问题。即使短时电压暂降也可能引发敏感设备故障或停运,造成重大的经济损失[1]。本文主要采用电压暂降频次[2]作为指标来准确刻画和估计电网中电压暂降信息。
母线节点上安装电能质量监测仪可直接记录此节点发生的电压暂降频次,但出于经济性考虑,为使电网建设成本最小化,优化配置后电网中监测仪安装的数量变得有限[3,4,5,6],如何利用监测节点发生的电压暂降频次来准确估计未装监测设备母线节点的暂降频次,进而获得整个电网的电压暂降信息,是一项具有重要现实意义的研究课题[7,8,9]。
系统故障受网络拓扑、故障类型、故障点、故障阻抗、变压器接线方式等多种不确定性因素影响,增加了评估难度。现有方法主要分两类: 基于实测的统计分析法[6,8,9,10]和基于随 机建模的 随机估计法[7,8,9]。实测统计法合理的精确度需要很长的测量周期,需要安装监测设备,成本较高,且统计故障率易受天气因素、绝缘子污秽程度以及人为因素等外界影响[11]。随机估计法主要分临界距离法和随机故障点法[9]。基于此我国学者王宾较早提出电压暂降状态估计方法[4],该方法采用最小二乘法搜索故障点所在路径,进而估计非监测点电压暂降幅值特征,但该方法只局限于简单辐射性电网,不适用于环网; 文献[9]采用随机故障点法将线路故障段均一划分得到电压暂降状态估计方程,但此方法获取较高的精确度需要设置较多的故障点,计算量较大,并且对如何合理地设置故障点数目和位置至今没有统一标准。本文将解析法[12,13]和随机故障点法相结合,利用节点阻抗矩阵,根据电压暂降幅值对传输线故障距离的解析式计算得到各种短路故障情况下各电压限值对应的临界故障点,合理划分线路故障区间可大量减少计算量的同时满足评估精度的要求。
2 VSSE 的基本原理
状态估计中,通用的数学表达式为:
式中,H为量测矩阵,本文和文献[5]相似,H中每个元素表示装有量测设备的母线节点记录的电压暂降频次,暂降电压对应于一个预先设定的临界电压值; 状态变量X基于故障位置的概念,X中每个元素表示线路故障区间发生故障频次[14]; M为量测量和状态变量的关系矩阵; I为测量误差,可忽略不计[8]。式( 1) 中相关变量矩阵具体形成阐述如下。
3 VSSE 建模
首先,假设电网有N个节点,L条支路,监测仪的数量设为C,显然,本文讨论的问题中C < N。
3. 1 量测矩阵 H 的建立
H由C个监测仪记录量测的电压暂降频次数据直接形成。通常,暂降表示故障发生时母线残余电压小于预先设定的临界电压值。若是三相不对称故障,母线残余电压取其中电压幅值最小的一相。
给定一个临界电压值V,得到相对应的量测量HV。HV中每个元素表示母线残余电压低于临界电压值V时相应的监测仪记录的暂降次数。因此,HV中共包含C个元素,且都为整型。
3. 2 状态变量矩阵 X 的建立
随机故障点法的基本思想为: 预先随机设定系统中的故障点,通过计算这些故障点故障所引起的电压暂降,来具体分析其具体的特征量信息。如图1所示,图中横轴代表故障线路的百分比,t1表示预先设定的临界电压阈值。若故障点数选取过少,线路故障区间随机划为0 ~ p1,p1~ p2,p2~ 1,由于实际母线电压包含在p1~ p2区间内,这对于发生在p1~ p2区间故障点检测与估计将造成不准确,进而导致评估结果精度变差。反之一味增加故障点数则会加大计算冗余度。如何设置故障点的位置与数目对电压暂降状态估计分析计算至关重要。
为了合理划分线路故障区间,可以利用解析法分析系统短路故障,基本原理如图2所示。假定线路p-q上的节点r处发生故障,节点k为监测点。
监测点处的电压Ek根据叠加原理可以用故障前的电压Ek0和节点阻抗矩阵表示:
式中,Zkr、Zkq、Zkp、Zpq分别为阻抗矩阵Z中节点k、r、p、q间的互阻抗; Zrr、Zpp、Zqq代表节点r、p、q的自阻抗; zpq代表线路p-q段的阻抗。
基于此,本文将解析法和故障点法相结合合理确定线路故障区间。如图3所示,为简单起见,此事例中只包含一条被监测的母线m1和一条待估计母线m2,图中表示系统中沿某条线路li发生故障时在母线m1和m2求得的残余母线电压的幅值。对于预先设定的电压阈值为t1,利用上述方法推导得到系统故障时节点电压对输电线路故障距离的函数表达式,再对应不同的电压暂降临界值得出线路上的临界故障点,将线路li划分为0 ~ p1,p1~ p2,p2~ p3,p3~ 1四个区段。
用上述方法遍历系统中所有线路L则可获得全网在电压阈值t1下线路总的故障区间P1。状态变量矩阵X1中每个元素表示每个故障区间内发生的故障次数。
3. 3 关系矩阵 M 的建立
由式( 1) 可知,矩阵M1( 对应临界电压值t1) 表示状态变量矩阵X1与量测量矩阵H1之间的关系。建立M1具体步骤如下:
( 1) 由3. 2节方法求得全网在电压阈值t1下线路总的故障区间P1;
( 2) 根据电网的参数和拓扑结构,随机模拟发生P1个故障区间内短路故障。
( 3) 分别计算出故障区间P1内发生短路故障时在C个安装在母线节点监测仪记录的母线残余电压。
将计算得出的电压值与相应的临界电压值比较形成关系矩阵M1。M1为二进制矩阵,M1中每个元素对应一条被监测母线与一个故障区间发生故障的关系,其值为1或0,即:
式中,A表示故障区间p发生故障时在监测母线产生的残余电压低于临界电压值t1; B表示故障区间p发生故障时在监测母线产生的残余电压高于临界电压值t1。
通过图3,易理解M1的形成,图中假设只有一条被监测的母线m1,则M1为:
由于p1,p4所在故障区间内发生故障求得母线m1残余电压小于t1,则,同时表示在母线m1引发电压暂降;
同理m1( 1,2) = m1( 1,3) = 0,因为在p2所在故障区间内发生故障求得母线m1残余电压大于t1,表示未引发母线m1电压暂降。
考虑整个电网所有输电线路,新的列将被加入矩阵M1,每一列对应线路的一个故障区间。值得注意的是,矩阵M1仅仅取决于被考虑电网的性质,对于一个给定电网和一个给定临界电压值只需要估算一次。
需要说明的是,电网中监测仪的安装数量与位置会影响VSSE的精度,本文监测仪安装配置采用文献[6,12]的方法,保证系统中任何位置发生故障至少能被一台监测仪测量记录。如此,M1中每列至少包含一个非0元素,即M1行满秩。
3. 4 VSSE 一般化
设定一个临界电压值t,则
考虑整个电网,VSSE一般形式为:
式中,Ht1,Ht2,…,Ht T为量测量,表示在监测母线在给定临界电压值t = t1,t2,…,tT时监测仪记录相对应的电压暂降频次; Mt1,Mt2,…,Mt T为临界电压值t= t1,t2,…,tT时相对应的二进制关系矩阵; X为状态变量。
式( 3) 可简化写为:
4 基于广义逆算法求解 VSSE
4. 1 状态变量 X 的求取
传统的状态估计方程往往是超定的,然而由于量测量不足,本文中由于监测仪的数量少于母线节点数,即C < N,式( 7) 为欠定方程或者叫病态方程。针对该问题,文献[12]采用了整数线性规划( ILP)的方法,但此方法求解速度较慢,文献[15,16]分别采用神经网络算法( ANN) 和遗传算法( GA) 通过迭代寻优的方法求解欠定方程组,相比ILP算法提高了计算效率。但GA算法对于求解高维度、多局部极值复杂问题,在迭代的过程中可能导致算法陷于局部最优的情况,在收敛速度和精度上难以达到期望要求[9],且该方法在交叉和变异时参数的选择目前大部分依靠主观经验,这严重影响解的品质。基于此本文则采用广义逆算法[15,16,17],此方法本身没有算法参数约束问题,且相对于ILP、ANN、GA算法,该方法不需要迭代,不存在收敛性问题,由于此方法计算复杂度较小,计算历时明显减少,同时也适合任何复杂电网。
基于广义逆算法理论,对于式( 7) ,易证明关系矩阵M满足Moore-Penrose方程[9],本文中由于优化配置使得电网中任意位置发生故障至少能被一台监测仪所检测,因此关系矩阵M中每列元素至少有一个非零元素,即M行满秩,所以存在M的Moore-Penrose逆为[11]:
结合式( 8) 则式( 7) 解为:
定理可证明方程组MX = H为欠定方程组时,解X = M+H不但是最小二乘解,而且是具有极小范数的最小二乘解,或最佳逼近解[12,16],简记为LNLS解。利用广义逆算法的求解结果实质上是满足以下约束条件:
式中,min‖x‖_2表示状态变量的最小二次范数。
4. 2 电压暂降频次的估计
为估计全网的电压暂降频次,还需建立关系矩阵Mtnm,和3. 3节中M矩阵建立相似,不同的是Mtnm中考虑的是未安装监测仪的母线节点。Mtnm的元素形成如下:
式中,A表示故障区间p发生故障时在未被监测母线产生的残余电压低于临界电压值t; B表示故障区间p发生故障时在未被监测母线产生的残余电压高于临界电压值t。
将求得的关系矩阵Mtnm与状态变量相乘即可得带估计母线节点的电压暂降频次:
结合式( 10) ,则式( 11) 转化为:
式中,Htnm中每个元素表示电压阈值t下,未安装监测仪母线节点的电压暂降频次。
5 算例仿真
利用上述方法应用于IEEE30节点标准测试系统,系统数据由文献[18]提供,电网接线图如图4所示,由6个发电厂、30条母线、37条传输线和10台变压器组成。本文模拟故障为三相短路故障,不对称故障参看文献[15],本文不再赘述。利用GA算法和广义逆算法仿真结果对比如图5 ~ 图7所示。
图5 ~ 图7分别描绘电压阈值设置为0. 7pu、0. 9pu、0. 8pu时实际母线电压暂降频次和估计的暂降频次,误差结果见表1和表2,误差定义为:
表1和表2误差结果表明广义逆算法精度明显优于遗传算法,利用该方法估计出平均电压暂降频次和实际值非常接近,其有效性和准确性显著。表3为针对不同阶数网络两种算法在Matlab中运算时间结果对比。从表3明显可以看出,由于广义逆算法无需迭代,计算历时相比遗传算法大幅度减少,比如针对IEEE-30节点电压阈值设为0. 7pu时,运用GA算法需要1087s,但通过广 义逆算法 只需要9. 53s的短暂时间。
( 单位: s)
另外,注意到随着设定电压阈值的降低,平均误差也随之升高,这是由于电压阈值设置越低,对应监测仪的监测范围随之减小,评估结果误差会相应增加。
6 结论
有限元逆算法 篇4
1 小波网络逆变换
小波分析[5]是一种窗口大小固定、形状可变,且时间窗和频率窗均可改变的时频局域化分析方法。设待分析信号为f(t),基本小波函数为φ(t),则小波变换形式如下:
式(1)可以记成内积形式WTf(a,b)=驥f(t),φa,b(t)驦,其中,φa,b(t)为小波母函数,a为尺度因子,b为平移因子,令φa,b(t)=。小波变换的时频分布规律同自然界信号的时频特性相符合,适宜分析任意尺度的信号,而小波变换的卷积形式需满足:
式中,A和B为常系数,当A和B越接近时,逆变换误差越小。式(2)为{φj,k(t)}(j,k∈Z)构成了一个L2(R)的框架[6]。若小波变换满足式(1),则小波逆变换存在,如下式:
式中,(t)为φj,k(t)的对偶框架。此式是小波逆变换的基本式,反映了细节信号的重构。
由于信号中往往含有奇异点,造成小波多分辨率分析的频带具有不均匀性,所以基于离散小波的多分辨率分析算法的谐波测量频带分割范围不均匀,导致在不同频带内,谐波的数量不一致,容易造成信号的混叠。为此,可以用小波模系数的方法提取小波系数,避免混叠,本文采用二进小波模极大值计算方法,其模的计算公式如下:
对于高频和低频信号的奇异点,在进行小波分解时,小波系数具有模极大值,利用此特点可以有效检测特征信号。
本文结合小波变换和神经网络的优点提出一种新思路:原始信号在信号抽取时,利用小波系数模极大值原理分别对奇数和偶数进行抽取,以增加细节、减少和消除信号混叠现象,再用BP神经网络对抽取后的信号进行学习,得到每个高真度的细节和逼近信号,最后通过小波逆变换重构信号,确定基波和谐波分量。从而提高谐波检测的准确度。在这里采用BP神经网络对抽取后的信号进行训练[7],对小波逆变换的系数进行调整,设BP神经网络是一个3层BP网络,有4个输入节点,3个隐含节点,2个输出节点。系统结构如图1所示。
图1中,h0(n)、h1(n)分别为低通和高通滤波器,g0(n)、g1(n)为相应的重建滤波器,f(n)为原始信号,y(n)为输出信号。
典型环节中的二抽取是对偶数坐标位置元素的抽取不同,本算法利用小波系数模极大值原理,同时抽取数组低频、高频段的奇数和偶数坐标位置的元素,避免未抽取的部分和已抽取部分产生信号混叠。抽取后经过BP神经网络的输出为:
式(5)中yl的混叠分量o0p(-z)f(-z)g0(z)和yh中的混叠分量o1p(-z)f(-z)g0(z)被消除,最终输出结果可以达到抗混叠条件,同时兼顾信号保持与混叠抑制两方面因素,可以有效消除小波混叠误差,从而为非稳态谐波检测提供了有力手段。
2 谐波检测方法
本文提出谐波检测及补偿方法为:利用小波变换对多频电网谐波信号进行分解,将各次谐波分量分解到不同频带的子频带信号中,构成多个子空间,从中检测出含有基波分量的子频带区域,其余子频带区域均含有谐波分量。对含有各谐波分量的子频带区域的小波分解系数取负数,基波所在区域小波系数不变,利用新得到的小波系数对信号进行重构,则重构信号中除了含有基波分量的区域之外,其余各次谐波分量均己进行了反相。将重构信号和原始谐波信号相减,则得到谐波补偿信号[8]。实际应用中,通过谐波检测方法检测出电网中的谐波成分,并通过智能算法计算出谐波补偿信号,将所得到的补偿信号转变成反相PWM,再通过逆变装置注入到电网,即可实现谐波抑制。谐波检测与补偿控制的结构图如图2所示。
实际情况下,电压和电流波形产生了高次谐波,电力系统为三相交流电系统,偶次谐波基本消除,因此这里只考虑奇次谐波。另外谐波次数越高,则含量越低。电压信号用小波系数表示为:
式中,dj0(k)表示尺度函数的系数,dji(k)表示小波函数的系数(也称变换系数)。电力系统的电压有效值在某一尺度j上可以用小波变换系数模值表示为:
小波分解是将信号按尺度函数和小波函数进行划分,利用小波系数模建立模极大值的特征向量,并对特征量按照隶属函数划分。不同频率的信号根据尺度的不同被划分到不同的频段中,对各频段分别进行奇数和偶数抽取,得到信号细节a(2k)、a(2k-1)和d(2k)、d(2k-1),从而分离出各次谐波。用3层神经网络对细节信号进行逼近训练,再确定综合滤波器g0、g1,然后用小波逆变换对信号重构,得出各个采样时刻的基波值和谐波值。
3 试验
在电网中电压和电流的基波频率均为f0=50 Hz,本文选择最常见的含有3、5、7、9次谐波的情况。设单相电压信号的数学表达式为:
式中,频率fi依次为基波、3次、5次、7次、9次谐波,幅值Ai依次为1 V、0.3 V、0.15 V、0.1 V、0.08 V,其信号仿真波形如图3～图5所示。
图3为含高次谐波的电力系统单相电压波形,图4为单相电压波形的频谱图,图5为分离出的谐波成分。小波网络逆变换算法能准确地将给定信号的基波信号和谐波信号分离出来,各尺度体现的频率成分变化趋势各不相同,表明没有出现混叠和泄露现象。使用离散小波变换提取子频带内的信息,利用3层神经网络对信号进行逼近,可以精确地量化谐波信号的频率和幅值。实验数据如表1所示。
通过表1可以看出,利用该算法分解出的各次谐波频率值误差率在10-5数量级,幅值的误差率在10-3数量级,完全符合谐波分析的精度要求,从而验证了基于谐波小波分析电力系统谐波分量是可行的。
信号通过小波分解到各个尺度空间的细节信号,利用小波系数模极大值可以有效分离出基波和谐波分量,用修正的系数重构原始信号。通过小波分解系数的重构就可以测量电力系统中的各个频带内的谐波频率和幅值。通过算法可以确定出信号中的各次稳态谐波以及谐波的含量,提高分析的可靠性,满足系统对精度的要求,在电力系统谐波检测中具有较好的应用前景。
参考文献
[1]帅定新,谢运祥,王晓刚.电网谐波电流检测方法综述[J].电气传动,2008,38(8):71-76
[2]杨君,王兆安,邱关源.一种单相电路谐波及无功电流检测新方法[J].电工技术学报,1996,11(3):42-46.
[3]何英杰,邹云屏,李辉,等.用于有源滤波器的一种新型谐波检测算法[J].电力电子技术,2006,40(2):56-58.
[4]戴君.基于小波变换的电力系统谐波分析与检测方法研究[D].无锡:江南大学,2008.
[5]飞思科技产品研发中心.小波分析理论与Matlab7实现[M].北京:电子工业出版社,2005.
[6]HAMID E Y,KAWASAKI Z I.Wavelet packet transformfor RMS and power measurement[J].IEEE Power Engineer-ing Review,2001,21(9):49-51.
[7]雷锡社,刘敏.基于神经网络的高精度电力系统谐波分析[J].电测与仪表,2007,44(3):17-21.
有限元逆算法 篇5
归纳与演绎的学习方法是机器学习与推理的策略,在人工智能领域得到广泛的研究与应用。Jevons在研究归纳问题时,认为归纳问题比演绎问题困难复杂,提出了逆演绎的思想。文献[1,2]将逻辑与程序统一为逻辑程序,在此基础上文献[3]将机器学习与逻辑程序结合,提出了归纳逻辑程序。归纳逻辑程序的研究主要包括理论研究与系统实现,最具代表性的是文献[4]给出了归纳逻辑程序的标准语义,并建立了基于逆蕴涵与精化算子的Progol系统。在国内,相关专家在归纳逻辑程序基础上提出了一些新的有关归纳逻辑程序的改进的设计方法,如文献[5]将遗传算法与逻辑程序结合,提出了遗传归纳逻辑程序设计方法;文献[6]向约束方向扩展归纳逻辑程序,提出了约束归纳逻辑程序设计方法等。本文在归纳逻辑程序相关研究启发下,针对传统归纳学习的困难,提出了一种新的逆演绎的学习算法,并将之应用到智能通信网学习推理系统中,来完成智能通信网的学习任务。实验结果表明,在逆演绎规则约束下,逆演绎算法能够高效地归纳出逻辑程序,学习能力强,能够完成智能通信网的学习任务。
1 逆演绎学习算法
1.1 逆演绎学习算法的相关概念
归纳是从特殊样例E推导出一般的理论T(特殊到一般);演绎则是从先验理论T推出样例E(一般到特殊)。虽然归纳是演绎的逆过程,但演绎具有可靠的推理规则,而归纳具有不可靠的假想成分,需要训练样本的统计假设支持。训练样本构成的集合为样例E,部分先验理论称为背景B,训练样本的统计称为假设H,从演绎角度,它们具有关系B∧H|=E。在此基础上,可以给出逆演绎学习的定义。
定义1 逆演绎的学习 一个操作因子q(B,E)=H使用样例E={<x,f(x)>}和背景B作为输入,输出一个满足(∀<x,f(x)>∈E)(B∧H∧x)|=f(x)的假设H,定义为逆演绎的学习。其中x为样例E中的一个样本,f(x)为x的目标。
该定义通过假设H和背景B来解释样例E={<x,f(x)>},学习过程就是发现假设H,使得样本x的目标f(x)从H、B中演绎出来。下面对进行学习算法过程中所涉及到逻辑程序、模式声明和霍恩状态函子作如下定义。
定义2 霍恩子句与逻辑程序 霍恩子句是带有最多一个肯定文字的子句。有且仅有一个肯定文字的子句叫明确子句,没有任何肯定文字的子句叫目标子句。现行逻辑程序的基本语句属于一阶谓词演算的一个子集,即霍恩子句集。霍恩子句的一般形式为:A1,…,A□←B1,B2,…,B□,其中A□(1≤□≤n)、B□(1≤□≤n)都是原子公式,分别代表结论和前提的形式。前提部分是各原子的合取式,构成子句体,结论部分最多只有一个原子,称为子句头。由此可将霍恩子句分成两个基本类型:①有头霍恩子句(用来代表一条规则),例如,grandmother(z,x)←mother(z,y),mother(y,x)代表:z是y的母亲且y是x的母亲,则z是x的祖母。有头无体的霍恩子句是一断言(用来代表一个事实),例如,father(A,B)代表:A是B的父亲;②无头霍恩子句,称为目标语句(用来代表结论的否定式),例如,←grandfather(A,C)代表:A不是C的祖父。
定义3 模式声明 模式声明M是由模式头mod eh(z,p)和模式体mod eb(z,p)构成的集合,其中z≥1为返回值,p为原子,p的项分为常规型和标记型2种,常规型有常量符或函数符构成,标记型有输入类符+type、输出类符-type或常量类符#type构成,例如,对于动物分类的模式声明,模式头为mod eh(1,class(+animal,#class)),模式体为mod eb(1,hasEggs(+animal))等;由文字和模式M确定的子句为C=h:-b1,…,bn,若子句属于模式语言即C∈L(M),当且仅当满足①h是M的模式头mod eh(z,p)的原子且+type和-type由变量替代,#type由项替代;②bi是M的模式体mod eb(z,p)的原子且+type和-type由变量替代,#type由项替代。
定义4 霍恩状态函子 B是霍恩子句集的背景,M是模式声明,a←b1,…,bn是定义子句;C⊥是模式语言L(M)中的最特殊子句,C⊥的基为自然数m满足1≤m≤n;C是L(M)中满足θ包容Cθ⊆C⊥的子句;若<C′,θ′,m′>∈ρ(<C,θ,m>),当且仅当C′=C,m′=m+1,θ′=θ,m<n,则ρ(<C,θ,m>)定义为霍恩状态函子。在霍恩状态函子定义的状态空间,学习算法进行子句搜索,可利用定义如下激励信息和调整函数来提高效率。
定义5 激励信息和调整函数 对霍恩子句表示的样例E、背景B、模式M、状态s=<C,θ,m>、最大子句体长度c、集合S的基
1.2 逆演绎学习算法
根据逆演绎的相关定义,逆演绎的学习任务就是:在背景B引导下,利用E为训练样例,寻找一个假设空间H,使得H对样例E是充分完备的。逻辑程序的核心问题是:如何从背景知识中择优选择谓词来构造满足约束的假设H。针对学习任务的算法如下。
输入:背景B,样例E,模式M。
输出:假设H{H1,…,Hn}。
Step 1 初始化H={ }。
Step 2 反复计算H,if E=Φ,返回H,else e=first(E)。
Step 3 构造实例e=first(E)的最特殊子句C⊥,取第一个头模式声明,建立哈希表,将项变量与哈希表作出映射,初始C⊥为空,头部项为空,搜索深度m为0,通过搜索深度m和C⊥匹配返回;利用模式M的头部项构建C⊥的头部和体部,输出最特殊子句C⊥。
Step 4 最大压缩量的子句Cm的搜索,顶结点为CΤ、底结点为C⊥、子句结点为CΤ≤H≤C⊥,激励信息为fs=ns-(cs+hs),有效样例ns,句体长度cs,文字长度hs;利用霍恩逻辑函子ρ(s)构造状态空间s=<C,θ,k>,利用Open和Close表进行待检测状态和已检测状态的存储,利用调整函数best(Open),prune(s)和ter min ated(Open,Close)对表操作,在状态空间搜索覆盖有效样例且压缩量为
Step 5 消去Cm体中的所有等式H:=unflattering(Cm)。
Step 6 覆盖规则集计算H:=H∪{H}。
Step 7 返回Step2继续操作。
Step 8 输出假设H={H1,…,Hn}。
2 智能通信网学习系统的任务和模式声明
2.1 智能通信网的学习任务声明
智能通信网中包含有各种子通信网:电话通信网、移动通信网、卫星通信网、支撑网、用户接入网、数据通信网、计算机通信网和综合业务数字网。智能通信网的这些子通信网络又有各种具体的网络,例如,电话通信网中的市话网、长途网和国际网;移动通信网中的模拟和数字移动通信网;卫星通信网中的FDMA卫星通信网、TDMA卫星通信网、SDMA卫星通信网;支撑网中的数字同步网、电信管理网和信令网;用户接入网中的光接入网、光纤/同轴电缆接入网和无线接入网;数据通信网中的分组交换数据网、帧中继网和数字数据网;计算机通信网中的局域网、都市网和广域网;综合业务数字网中的窄带ISDN、宽带ISDN和ATM交换网。我们把智能通信网的这些子通信网以及这些子通信网的各种具体网络都看作结点,这些子网之间的关系看作边,所有的这些结点和边就构成了智能通信网的知识库。对智能通信网的学习任务就是在智能通信网知识库的基础上,利用某个通信资源样例,来寻找一个规则空间,使得规则空间对于该资源样例是充分完备的,并通过知识获取来扩展丰富学习内容,完善智能通信网的知识库。
2.2 智能通信网模式和背景的声明
在智能通信网分类关系的学习中,本文通过优化参数来设定激励信息[f、n、h],其中f为压缩函数、n为覆盖有效样例、h为输入输出关系文字长度。需要学习的通信网分类关系用目标谓词表示为class(X,Y),通信网的特征用谓词hasPhone(电话业务)、hasComputer(计算机通信业务)、hasSatellite(有卫星天线支持)、hasUser(有用户接入到核心网)、hasTopology(拓扑结构)、hasPart(组成结构)、workspace(工作空间)表示,则模式声明M如下(限于篇幅仅举知识库中的几个样例):
3 逆演绎学习算法在智能通信网中的实例测试与分析
3.1 智能通信网的学习算法测试
针对逆演绎的学习算法,利用Visual-Prolog编写一个智能通信网的学习推理系统来进行测试。系统读取具有模式、背景和样例的规格说明文件,通过规格说明告知系统的解释器,指出所有模式声明并增加背景、样例到子句库;模式语言是对系统的解释器的查询,采用Prolog[7]的查询语句,诸如:-modeb(1,hasPhone(+communicationnetwork)) 。模式声明的语句可用作系统的假设H子句,如样例语句class(local-call-network,telephone-network)(市话网属电话通信网);背景语句是Prolog任意代码段,如hasTopology(local-call-network,star)(市话网的拓扑是星型);系统生成与模式头modeh相关的谓词,如语句:-modeh(1,class(+communicationnetwork,#class) 生成class(X,#):-…形式,其中X是类型communicationnetwork/1输入变量,#是类型class/1常量,…将通过模式体modeb声明:-modeb(1,hasPhone(+communicationnetwork)) 来构造。从样例的格结构的顶结点class(local-call-network,telephone)出发,通过搜索算法可学出如class(A,telephone-network):-hasPhone(A)的关系。
对其中一个有效样例结点class(local-call-network,telephone-network),学习算法测试的运行结果如下:
3.2 测试结果分析
逆演绎算法在PC机上测试,CPU为Pentium 4,主频为2.0GHz,内存为1GB,编译系统为Visual-Prolog。算法测试结果表明:逆演绎的学习算法包括构造最特殊子句、搜索状态空间和计算规则集;利用模式M和最特殊子句构造算法,以最泛化子句为顶结点、最特殊子句为底结点、其他为中间结点,通过霍恩状态函子构建状态空间,搜索得到激励信息的压缩最大的有效样例结点;再利用规则集算法不断覆盖有效样例,计算出最优子句,重复计算出具有最大压缩率的假设H,这就是搜索的结果,也是学习的结果。
本文将逆演绎算法应用于智能通信网学习的知识获取,并与传统的归纳学习算法—AQ算法相比。如图1(a)所示,采用逆演绎学习算法的学习系统的学习效率明显提高,在样例容量较大的情况下还能保持较高的搜索精准度,结果的改进是明显的。这是由于采用逆演绎学习算法的学习系统充分利用了谓词的模式,这使得算法的搜索空间大大减小,另一方面,该学习系统使用相关的背景知识对搜索空间加以适当约束,也缩小了搜索范围,如图1(b)所示。从图1显示表明,在学习系统样例容量较大时,采用逆演绎学习算法的搜索空间比传统归纳学习算法缩减了将近一半,搜索精准度则提高了一倍,从而提高了搜索效率。
4 结 语
本文在对智能通信网的学习过程中,针对传统的归纳学习方法的困难,提出一种逆演绎的学习算法,实现了逆演绎的学习算法在智能通信网学习系统中的成功应用。利用逆演绎的学习算法,我们能对智能通信网的各个子网结点的相关分类等信息进行学习,从而织就出一张有关智能通信网的知识网络。通过知识获取,不断往知识网中增加新的结点,并将之融入进已学习的网络中,扩展智能通信网中的结点信息,设计实现并完善智能通信网的学习推理系统。
参考文献
[1]Kowalski R.Algorithm=logic+control[J].CommunicationsoftheACM,1979,22(1):424-436.
[2]Colmerauer A.AnintroductiontoprologⅢ[J].CommunicationsoftheACM,1990,33(7):69-90.
[3]Muggleton S.Inductivelogicprogramming[J].NewGenerationComputing,1991,8(4):295-318.
[4]Muggleton S.Inverseentailmentandprogol[J].NewGenerationComputing,1995,13(5):245-268.
[5]杨新武,刘椿年.遗传归纳逻辑程序设计中规则的位串表示法[J].北京大学学报,2001,27(3):297-302.
[6]郑磊,刘椿年.约束归纳逻辑程序方法的研究[J].计算机工程与应用,2003,39(10):63-66.
有限元逆算法 篇6
椭圆密码算法ECC (Elliptic Curve Cryptography) [1]具有每bit最高的安全性[2], ECC算法的主要优点有[3]:安全性能高;计算量小且运算速度快;占用存储空间小;对数据运算带宽要求低。目前常用的定义椭圆曲线的有限域有两种:素数域GF (p) 和二进制域GF (2m) 。在二进制域上元素可以直接用二进制位向量表示, 这样元素的加、减、乘等操作在计算机中可以并行实现, 所以二进制域的椭圆曲线加密通常被用于硬件实现的椭圆曲线加密系统中。GF (2m) 上的非奇异椭圆曲线是Weierstrass公式定义的:y2+xy=x3+ax2+b。其中a, b是GF (2m) 中的元素, 且b≠0。曲线E同时还包括一个无穷远点, 通常用O表示, 为曲线上点的加法单位元。有限域上的模逆运算是椭圆曲线密码系统中的关键性计算, 但是其运算效率却是最低的, 据M.Brown, D.Hankerson等人粗略估计[4], 求逆运算与乘法运算 (带有快速取模运算) 的资源占用比是80∶1, 因此通过改进加速求逆运算减少其计算时间可以在很大程度上提高椭圆加密的运算效率和吞吐量。
GF (2m) 上的求逆算法有两类:基于费尔玛小定理使用乘法的算法、扩展欧几里德算法及其变体。Scliroeppel等人在文献[5]中提出了Almost求逆算法。如果椭圆加密算法使用仿射坐标, 那么求逆运算将是制约标量乘法性能提高的主要瓶颈。扩展欧几里德算法及其变体适合加密算法的软件实现;基于扩展欧几里德算法求逆电路会对核心电路和控制器并行性要求较高, 所以硬件实现的研究还较少。
本文基于二元Euclidean算法提出了改进, 根据资源和需求的不同, 通过修改算法执行结构, 将原算法中制约效能的步骤进行改进, 并引入并行执行方式, 实现了求逆运算的加速。
1求逆运算的原理
域元素的求逆运算无论在软件和硬件实现上都有很大复杂度, 相比于乘法运算其计算效率要低很多。单纯的求逆算法有两种:基于扩展Euclidean方法的算法、基于Fermat小定理的算法。基于Fermat小定理的算法是将求逆运算转化为升幂运算, 算法中的模乘和模平方运算决定其运算速度, 算法的评估软件得到结果是完成一次求逆运算需要「log2 (m-l) 」+w (m-1) -1次模乘和m+1次模平方, 这种方法不但会增加软硬件设计的复杂度, 算法中的运算需要反复执行, 性能受到很大影响[7]。而扩展Euclidean方法便于硬件实现, 没有大量的乘法, 只涉及移位、判断和加法, 可以方便并行化设计, 因此, 本文通过研究三种基于扩展Euclidean方法的算法, 得到新的算法:对二元Euclidean方法进行并行优化和改进。
Euclidean算法本质是计算多项式的最大公因式。在多项式a1 (x) 和a2 (x) 的最大公因式gcd ( (a1 (x) , a2 (x) ) 时, 用a1 (x) 除a2 (x) 得到商p (x) 和余数q (x) , 满足a2 (x) =a1 (x) p (x) +q (x) , 其中q (x) 的次数低于a1 (x) 的次数, 由Euclidean算法得gcd (a1 (x) , a2 (x) ) =gcd (a1 (x) , q (x) ) 。
定理1设a1 (x) , a2 (x) 和a3 (x) 是GF (2m) 上的任意多项式, 求a1 (x) 和a2 (x) 的最大公因式有:gcd (a1 (x) , a2 (x) ) =gcd (a1 (x) , a2 (x) -a3 (x) a1 (x) ) 。对于二元扩域上的逆元也可以使用该算法。二元扩域GF (2m) 上的元素a, 其多项式表示为a (x) , 那么a的逆元a-1可用a (x) 的扩展Euclidean算法求出。
使用传统的Euclidean算法计算正整数a, b的gcd原理就是当b≥a时, 用a去除b得到商q和余数r, 满足:b=qa+r, 且0≤r<a, 根据定理知gcd (a (x) , b (x) ) =gcd (a (x) , r (x) ) , 所以求gcd (a (x) , b (x) ) 的问题就转化为gcd (a (x) , r (x) ) , 求其中元素 (r, a) 小于 (a, b) 。上述过程重复执行, 直到 (r, a) 中的一个元素为0, 此时得到gcd (0, d) =d, 由于非负余数严格递减, 此算法肯定收敛, 再者, 如果a的位宽为k位, 则除法步数至多为2k, 所以该算法效率还是较高的。Euclidean算法可以扩展到寻找整数x和y满足ax+by=d, 其中d=gcd (a, b) 。
算法1求二进制多项式的Euclidean算法
INPUT:非零多项式a和b同时deg (a) ≤deg (b) 。
OUTPUT:d=gcd (a, b) 以及二进制多项式g, y满足ag+bh=d。
假定f为一个m次二进制不可约多项式, 非零多项式a的次数不超过m-1, 这时gcd (a, f) =1, 在算法中, 如果输入是a和f, 则在步骤3中, 最后的非零的输入u=1, 那么在步骤4中有ag1+fh1=1。所以ag1≡1 (mod f) and so a-1=g1。说明:在求取g1的过程中, 不需要求出h1和h2, 基于此导出如下的求GF (2m) 上的求逆的扩展Euclidean算法。
算法2扩展Euclidean算法求逆[6]
算法3二元Euclidean算法[3]
另外还有一种改进的二进制算法求逆方法叫殆逆算法, 运算过程是先计算一个多项式h (x) 以及一个正整数k, 满足a (x) h (x) ≡xkmod f (x) , 然后得到a-1 (x) =h (x) x-kmod f (x) 。如果多项式是某些特定的结构, 其约减速度殆逆算法会比二元Euclidean算法快, 如果多项式为任意结构, 其约减速度比较的话二元Euclidean算法适应性更强。
2求逆算法的改进及其硬件实现
2.1改进求逆算法
在实现二元Euclidean算法时, 总共需要四个m+1位的寄存器u, v, g1, g2, 其中deg (u) 和deg (v) 分别表示u和v中域元素的多项式次数, 在第2步中判断u和v是否能整除x, 只要判断u和v的最低比特位是否为0, 在整除时只要对应的寄存器进行逻辑右移, 第3步是域元素的加法, 因此制约二元Euclidean算法的执行效率的因素是对循环的判断, 以及比较deg (u) 和deg (v) 。文献[8]给出了一种方案, 这种算法能在一个时钟周期内比较deg (u) 和deg (v) 的函数, 其硬件实现结构为优先级译码器, 这种设计利用硬件资源换取了速度的提高, 但是当选取域GF (2m) 的m值较大的时候, 其延迟也很大。因此, 提出一种改进算法有效加速deg (u) 和deg (v) 的比较, 能在很大程度上提高求逆速度。本文就是基于这一思想, 增加了两个寄存器Count_u和Count_v, 用于存储并检测u和v中元素次数的变化, 设计了一种并行的硬件结构, 同时执行循环中的判断和移位、加法操作。本文提出的改进的算法二元Euclidean算法描述如下:
算法4改进的二元Euclidean算法
利用设计的独立判断和移位元件并行执行2.1和2.2
本文增加的两个寄存器Count_u和Count_v只需在初始化的时候赋予二进制项的系数, 然后随着移位的进行检测u和v的次数的变化, 也就是在步骤2.1和2.2中根据相应的跳转条件减1计数, 在硬件中比较两个 (log2 (m+1) ) 位的寄存器Count_u和Count_v是容易实现的。程序状态图如图1所示。
在循环的判断、移位、域加法的操作执行时增加的两个并行部件, 这两个部件同时执行, 而且他们的出口都是连接到并行判断元件上去, 一旦跳出条件满足, 两个部件同时终止, 这就避免了串行循环带来的大延迟, 以增加了很小的硬件资源的代价换来了求逆速度的极大的提升。
2.2改进求逆算法的结构图
根据2.1节的改进算法, 运算流程图如图2所示。
参数初始化之后进入并行执行模块, 在并行执行模块内部, 首先判断执行条件, 当条件满足时, 在两个运算部件中分别判断Count_u和Count_v以及u或者v次数大于1的条件满足时, 分别单独执行。在并行执行的同时, 并行判断部件也在运行, 当Count_u或Count_v的跳转条件满足时, Count_u>Count_v或者Count_u<Count_v就进入并行执行部件, 当Count_u=1或Count_v=1时程序运行结束。这样经过至多2m-1个时钟周期后, 就可以得到程序的输出结果。
由于改进的求逆算法在计算没有太多的I/O限制, 而且计算密集, 判决条件具有并行性, 运算步骤可以流水处理和并行处理, 根据上述的运算流程, 本文基于FPGA设计了高效的执行结构, 利用FPGA内部丰富的逻辑布线资源和算法功能块实现结构图如图3所示。
2.3硬件仿真与实现
根据2.2节改进原理以及实现结构图, 本文选取有限域GF (2192) , 编写相应的Verilog代码, Xilinx Virtex-5XC5VFX70T的硬件开发板上对改进算法进行板级验证。本文设计的求逆代码仿真结果如图4所示。
针对一次求逆所需的时间, 与现有的实现效果进行了对比, 比较如表1所示, 数据表明:在使用资源规模相当的情况下, 本文提出的改进算法运算效率提高约20.9%。
3结语
椭圆曲线密码进行实现时, 有限域上的求逆运算占有很大比例。由于求逆运算是标量乘法的重要组成部分, 提高求逆算法效率可以在很大程度上提高点乘的运算速度, 进而提高整个椭圆加密的吞吐率, 本文提出的改进二元Euclidean算法, 简化了其中的耗时运算, 并且引入了并行化操作机制, 在此基础上设计了一套适合FPGA实现的求逆运算器结构。采用Verilog语言, 将此求逆算法进行并行优化, 并且编程实现, 使用Xilinx ISE 12.4对整个算法仿真、综合, 在Xilinx Virtex-5XC5VFX70T的开发板上实现了域内标量乘法的验证, 比较结果表明该改进算法有效地提高了椭圆曲线加密算法的硬件运算速度。
摘要：针对二进制域上现有求逆算法计算量大、并行度小、速度慢的缺点进行改进, 基于二元Euclidean算法提出了改进, 设计了相应的乘法器硬件结构, 并且分析了其运算效能和资源占用情况。将此求逆计算器的并行改进算法使用Verilog语言编程实现, 利用Xilinx ISE 12.4对整个求逆算法综合仿真 (行为级) , 在Xilinx Virtex-5 XC5VFX70T的硬件平台上验证求逆算法的运算效率, 结果表明对求逆算法的改进有效地提高了求逆运算的速度。
关键词：椭圆加密,二进制域,求逆,扩展欧几里得算法
参考文献
[1]W.DIFFIE H M.New direction in cryptography[J].IEEE Transactions on information Theory, 1976, 22 (6) :644-654.
[2]王峰.GF (2163) 上椭圆曲线密码体制的FPGA实现[D].广州:广州大学, 2006.
[3]HANKERSON Darrel, MENEZES Alfred, VANSTONE Scott., Guide to Elliptic Curve Cryptography[M].张焕国, 译.北京:电子工业出版社, 2005.
[4]BROWN M, HANKERSON D, LOPEZ J, et al.Software implementation of the NIST elliptic curves over prime fields[C]//Topics in Cryptology-CT-RSA.[S.l.]:Springer, 2002:250-265.
[5]SEHROEPPEL R, ORMAN H, O′MALLEY S, et al.Fast key exchange with elliptic curve systems[C]//Advances in Cryptology-CRYPTO.[S.l.]:[s.n.], 1995:43-56.
[6]IEEE.IEEEP1363-2000 IEEE standard specifications for public-key cryptography[S].USA:IEEE, 2000.
[7]ZHANG Yu, CHEN Dong-dong, CHOI Younhee, et al.A high performance ECC hardware implementation with instruction-level parallelism over GF (2163) [J].Microprocessors and Microsystems, 2010, 34:228-236.
[8]高献伟, 欧海文, 董秀则, 等.基于FPGA的GF (2) 域求逆算法的设计研究[J].计算机工程与应用, 2005 (9) :135-137.
[9]秦帆, 戴紫衫.有限素域上椭圆曲线模逆运算的设计与实现[J].计算机工程与应用, 2008, 44 (23) :117-119.
有限元逆算法 篇7
关键词：SCARA机器人,运动学,逆解
1 引言
工业机器人是一种能够自动定位控制、可重复编程的、多功能的、多自由度的设备。现在机器人已经成功应用于多个领域。工业机器人可以完成材料搬运、分拣、组装等作业。可以用于危险、不适宜人工作的环境代替人进行作业, 从而减少对人力的需要, 减少对人体安全的损害。提高生成效率。SCARA机器人是专门为工业要求而开发的机器人系统, 适合在平面范围内实现对物体的快速取放或者装配等。
2 SCARA机器人简介
SCARA机器人是一种四轴机械手, SCARA (Selective Compliance Assembly Robot Arm, 中文译名:选择顺应性装配机器手臂) 是一种圆柱坐标型的特殊类型的工业机器人。SCARA机器人很类似人的手臂的运动, 它包含由肩、肘和腕关节连接的臂部、肘部和手部。可以实现平面内的定位。它具有四个关节, 三个转动关节和一个移动关节。它可以实现沿三个坐标轴移动和沿Z轴转动。
3 SCARA机器人的逆解算法设计
机器人由连杆以及衔接连杆的关节组成, 关节一般装有驱动装置, 用于驱动连杆运动, 从而实现希望的功能。连杆确定了其连接的关节之间的位置关系。
SCARA机器人是一种串联机器人。串联机器人的运动学分析通常采用D-H方法。根据D-H方法, 首先确定机器人每个关节的结构参数和运动参数, 即连杆的扭角、连杆的杆长、偏距和关节转角。对于每个关节结构参数是常量, 运动参数是变量。下面是采用D-H法建立SCARA型搬运机器人的坐标系:
由坐标系示意图得知如下的D-H参数表:
首先需要知道坐标系{i}相对于坐标系{i-1}的连杆变换。坐标系{i-1}经过如下4次有序的相对变换可得到坐标系{i}:
连杆变换矩阵:
将D-H表中的参数带入变换矩阵, 可得到相邻两坐标系的位姿变换矩阵如下, 公式中si=sinθi, ci=cosθi:
根据创建好的连杆坐标系和整理的机器人D-H参数, 计算SCARA机器人末端执行器参考坐标系相对于基坐标系的齐次变换矩阵。
逆运动学解法分为三类:代数法、几何法和数值解法。代数法和几何法的步骤和最终表达式, 对于不同的机器人不同。数值解法计算量大, 计算比较耗时, 目前在实际应用中往往不能满足控制要求。
对应SCARA机器人的运动学逆解采用代数法, 应用矩阵变换来求解。先计算SCARA机器人的运动学正解的各变换矩阵。然后计算各矩阵的逆矩阵, 最后解代数方程, 计算出用末端执行器空间位姿表示的各关节变量表达式。
假设已知第四关节臂终端位置[px, py, pz]T, 求解各关节变量[θ1, θ2, θ3, θ4]T
上面SCARA的动力学方程的第四列为第四节关节臂终端的位置矢量, 因此有:
再联合解上面的方程组, 可得:
式中代表:
由动力学方程可得:
4轨迹规划
轨迹规划根据实际的工作要求, 根据作业过程中选取的一些关键点的位置、速度、加速度, 为机器人末端执行器从一个关键点运动到另一个关键点选择满足某些要求的机器人运动轨迹, 从而控制机器人的运动。轨迹规划是运动学反解的具体应用。
轨迹规划方法包括:关节空间轨迹规划和直角坐标空间轨迹规划。关节空间轨迹规划只要求机器人手部经过要求的设定点, 至于设定点间的路径基机器人手部的位置及姿态未加要求。该方法的优点是计算简单, 缺点是运动期间各连杆位置不易判定。直角坐标空间轨迹规划以轴位置 (转角) 来描述轨迹, 所以无论起点、终点还是途径点都必须采用逆运动学将直角坐标下机器人末端执行器的位姿转变成各关节的转角才可以进行轴坐标控制。每根轴随时间的转动都可以用一个函数表示。代表各轴运动的函数可以不同。直角坐标空间轨迹规划的优点是1) 机器人手部运动平滑, 2) 手部的位姿在整个路径上均有定义, 3) 适用于有障碍物的工作空间。缺点是1) 要将空间路径转换成轴关节转动, 计算复制、耗时, 2) 计算中有奇异点产生。
本文应用的是基于关节空间, 以下是关节空间路径规划的一些基本步骤:
1、通过示教获得轨迹的起点和终点。
2、计算起点和终点对应的运动学逆解。
3、进行轨迹规划。
4、对规划的轨迹进行插值。
轨迹规划时, 为了使电机的运动平稳。在各关节的运动过程中, 在过渡点我们希望电机的速度变化平滑。也就是在起点、终点或方向改变的点希望速度变化平稳, 其它的点速度可以比较快一点, 以提高效率。为了实现这个目的, 我们才用如下的轨迹规划方法:每个关节的运动分成从起始点到提升点, 从提升点到下放点, 从下方点到终点三段。其中第一段使用四次多项式规划, 第二段采用立方规划, 第三段采用四次多项式规划。
并且为了计算方便, 采用归一化时间变量t∈[0, 1], 这样我们可以对不同的关节都采用同样的处理方法。
应该满足的边界条件如下:
根据这些边界条件可以根据上面几段式子进行求解。有下面结果:
其中:
5 示教规划设计及实现
通过示教盒规划作业流程及机器人动作, 编制机器人控制程序, 然后机器人控制器根据机器人编好的机器人控制程序, 进行机器人运动轨迹规划。机器人运动控制器把规划好的数据通过串口发送给多轴运动控制卡, 运动控制卡控制各个电机沿规划好的轨迹运动。
6 结束语