自动微分范文10篇

自动微分范文篇1

关键词自动微分切线性模式数据相关分析统计准确率

1.引言

计算微分大致经历了从商微分，符号微分，手写代码到自动微分几个阶段。与其它几种微分方法相比，自动微分具有代码简练、计算精度高及投入人力少等优点。自动微分实现的基本出发点是：一个数据相对独立的程序对象（模式、过程、程序段、数值语句乃至数值表达式），无论多么复杂，总可以分解为一系列有限数目的基本函数（如sin、exp、log）和基本运算操作（加、减、乘、除、乘方）的有序复合；对所有这些基本函数及基本运算操作，重复使用链式求导法则，将得到的中间结果自上而下地做正向积分就可以建立起对应的切线性模式，而自下而上地做反向积分就可以建立起对应的伴随模式[1]。基于自动微分方法得到的切线性模式和伴随模式，在变分资料同化[2]、系统建模与参数辨识[3]、参数的敏感性分析[4]、非线性最优化以及数值模式的可预测性分析[5]等问题中有着十分广泛的应用。

迄今为止，已有数十所大学和研究所各自开发了能够用于求解切线性模式的自动微分系统，比较典型的有TAMC系统[6]、ADJIFOR系统[7]和ODYSSEE系统[8]。在一些特定的运用中，它们都是比较成功的，但在通用性和复杂问题的处理效率上还存在许多不足。通常，自动生成切线性模式的关键难题在于对象自身的强相关性，这给系统全局分析（如数据IO相关分析和数据依赖相关分析）和微分代码的整体优化都带来了很多困难。同时，对于程序对象不可导处的准确识别和微分处理，至今仍还没有一个统一而有效的算法。另外，最优或有效求解稀疏雅可比矩阵一直是衡量一个自动微分系统有效性的重要尺度。

统计准确率被我们视为评价一类自动微分工具及其微分模式代码可靠性与有效性的重要尺度。其基本假设是：如果对于定义域空间内随机抽样获得的至多有限个n维初始场（或网格点），微分模式输出的差分和微分逼近是成功的；那么对于定义域空间内所有可能初始场（或网格点），微分模式输出的差分和微分逼近都是成功的。微分模式统计准确率评价的具体方法是：在所有随机抽样得到的初始场（或网格点）附近，当输入扰动逐渐趋向于机器有效精度所能表示的最小正值时，模式输出的差分和微分之间应该有足够精度有效位数上的逼近。

DFT系统具有许多优点，它能够完全接受用FORTRAN77语言编写的源代码，微分代码结构清晰，其微分处理能力与问题和对象的规模及复杂性无关。它基于YACC实现，具有很强的可扩展性。DFT系统具有四个重要特色。它通过对象全局依赖相关分析，准确求解雅可比矩阵的稀疏结构，自动计算有效初始输入矩阵，从而可以用较小的代价求得整个雅可比矩阵。同时，它可以自动生成客观评价微分模式效率与可靠性的测试程序，对奇异函数做等价微分处理，并采用二元归约的方法，在语句级层次上实现微分代码优化。

2.系统概况

DFT系统主要由两部分组成：微分代码转换和微分代码评价，图2.1。微分代码转换部分接受用户输入指令并自动分析对象模式，生成切线性模式代码及其相关测试代码，后者直接构成微分代码评价系统的主体。微分代码评价是DFT系统的一个重要特色。DFT系统的开发小组认为，一个微分模式如果在可靠性、时间和存储效率上没有得到充分的验证，至少对实际应用而言，它将是毫无意义的。

原模式切线性模式

统计评价结果

图2.1DFT系统结构简图

2.1微分代码转换

DFT系统是基于YACC在UNIX环境下开发的，其结构图2.2所示。通过DFT系统产生的切线性模式代码成对出现，并在语句级程度上做了简化，可读性很强，如图2.4。

切线性模式

评价函数集

图2.2微分代码转换

微分代码转换部分从功能上分为四个部分：词法分析，语义分析，对象复杂性及数据相关分析和微分代码转换。对于一组具有复杂数据相关的程序模式对象，通常需要系统运行两遍才能得到有效而可靠的微分代码。这主要有两方面的考虑：其一，根据对象的复杂性（如最大语句长度、最大变量维数、子过程或函数数目、子过程或函数内最大变量数目等对象特征）选择合适的系统参数以求最优的运行代价；其二，模式内各子过程或函数之间以及一个子过程或函数内往往具有很强的数据相关性，需要事先保存对象的相关信息并且在考虑当前对象的属性之前必须做上下文相关分析。

图2.3PERIGEE源程序代码图2.4DFT系统生成的切线性代码

2.2微分代码评价

通常，评价一个编译系统的性能有很多方面，如处理速度、结果代码可靠性及质量、出错诊断、可扩展和可维护性等。对于一类自动微分系统来说，由于软件开发人力的局限以及对象模式的复杂多样性,通过自动转换得到的微分模式并非常常是有效而可靠的（即无论是在数学意义上还是在程序逻辑上应与期待的理想结果一致），因而在微分模式被投入实际应用前，往往需要投入一定的人力来对其做严格的分析测试。

对切线性模式做统计评价测试的主要内容可以简单叙述为：在网格化的模式定义域空间内，选择所有可能的网格点形成微分模式计算的初始场；在不同的网格点附近，随机选取至少个线性无关的初始扰动，对每个扰动输入分别进行网格点逼近，统计考察模式输出差分和微分在有效位数上的逼近程度。图2.5描述了整个测试过程，它包含网格点数据随机采样（1）和网格点数据逼近（2）两级循环。

图2.5切线性模式代码的测试过程

3.系统主要特色

DFT系统并不是一个完整的FORTRAN编译器，但它几乎可以接受和处理所有FORTRAN77编写的源模式代码，并且可以很方便地扩展并接受FORTRAN90编写的源模式代码。本节将着重介绍DFT系统（版本3.0）的以下几个重要特色。

3.1结构化的微分实现

DFT系统采用标准化的代码实现，切线性模式的扰动变量和基态值变量、微分计算语句和基态值计算语句总是成对出现，并具有清晰的程序结构。微分代码保持了原模式本身的结构和风格（如并行和向量特性、数据精度等），即语句到语句、结构到结构的微分实现。在奇异点或不可导处，DFT系统对微分扰动采取简单的清零处理，实践证明这对抑制扰动计算溢出具有重要意义，但并不影响评价测试结果。

3.2全局数据相关分析

DFT系统具有较强的数据相关分析能力，它包括全局数据IO相关分析、全局数据依赖相关分析、全局过程相关分析以及数据迭代相关分析几个不同方面。数据依赖相关与数据IO相关关系密切，但又存在根本不同。前者强调每个变量在数学关系上的依赖性；而后者描述了一个对象的输入输出特性，且具有相对性，即任何一个变量参数，无论它是独立变量还是依赖变量，在数学意义上都可等价为一个既是输入又是输出的参数来处理。

DFT系统记录所有过程参数的IO属性表，通过深度递归相关计算，准确计算每个过程参数的最终IO属性。DFT系统通过对数据相关矩阵做模二和及自乘迭代计算（An+1=An⊕An2）来完成数据的依赖相关分析，这种算法具有很好的对数收敛特性。DFT系统通过全局过程相关分析的结果，自动生成模式的局部或整体相关引用树结构（如图3.1），这对用户分析复杂数值模式和微分评价测试都具有很好的指导作用。DFT系统还具有分析局部数据迭代相关和函数迭代相关的能力，这两种形式的数据迭代相关是自动微分实现颇具挑战的难题之一。

图3.1GPSRayshooting模式的相关树结构片段

3.3自动生成测试程序

基于IO相关分析的结果，DFT系统自动生成微分测试代码，分别对切线性模式的可靠性和运行代价做统计评价测试。特别地，DFT系统还可将任何模式参数都视为输入输出参数，生成在数学意义上等价的测试代码，这样处理的不利之处在于往往需要极高的存储开销。

3.4基于语句级的代码优化

目前，DFT系统仅仅具备局地优化能力。在语句级微分实现上采用二元归约的方法对微分代码进行优化是DFT系统的一个重要特色。根据右端表达式的乘法复杂性及含变元数目的不同，DFT系统采取不同的分解策略。二元归约的方法避免了微分计算中的许多冗余计算，在一些复杂的非线性表达式的微分计算中具有最小的计算代价，同时也非常适合于微分系统的软件实现。同时，对于某些特殊的运算操作（除法、乘方）和特殊函数（如sqrt、exp），DFT系统较好地利用了基态值计算得到的中间结果，避免了微分实现中的冗余计算。

4.系统应用

运用自动微分工具得到的切线性模式，可以在无截断误差意义下求解函数的数值微分和导数、稀疏雅可比矩阵。同时这些结果在数值参数敏感性分析、非线性最优化以及其它数值理论分析中有着非常重要的应用。这里简单介绍切线性模式的几个基本应用。

4.1符号导数和微分

如果输入为数学关系式，DFT系统可以自动生成对应的微分表达式和梯度，而与数学关系式的复杂程度无关。例如我们输入关系式：

，（1）

DFT系统将自动生成其符号微分形式及其梯度形式分别为

，（2）

4.2数值导数和微分

切线性模式最基本的应用就是在一定扰动输入下求解输出变量的扰动（响应）。表4.1给出了DFT系统在对IAP9L模式、GPSRayshooting模式和GPSRaytrace模式三个数值模式做切线性化的具体应用中，一些不同计算粒度、不同引用深度和不同程序风格的核心子过程，以及它们的切线性模式在SGI2000上运行的统计评价测试结果，其中切线性模式的可靠性指标都准确到六个有效数字以上，在运行时间、存储开销和代码复杂性方面分别是原模式的两倍左右，比较接近于理想的微分代价结果（1.5倍）。除了IAP9L模式由于过于复杂仅做粗略统计外，其余模式都用非注释语句行数来表示各自的代码复杂性。

表4.1DFT系统在三个数值模式中的统计评价测试结果

性能指标

对象模式运行时间（10-3秒）存储开销（字节数）代码复杂性

原模式切线性

模式

原模式切线性

模式

原模式切线性

模式

Xyz2g2.5306.1605524110485589

IntCIRA1.5602.750133426614165

Dabel0.0350.072601202749

LSS8.30017.50669133879143

RP42.4085.10360572102238

Vgrad10.1000.21218564368282454

RefGr43.0086.0071865414373083578

LL2JK0.6261.350262252442232

RayFind462.70

×103125.4

×103985618212111179

EPSIMP1.76011.50445589101327

Hlimits0.8301.8802425774842543774

Int3sL26.9051.2082002916394584690

MAKE

NCEP1340392072292514458504584

Curvcent0.0130.038527542754

DYFRM3.800

×1037.250

×1035000*9500*161279

PHYSIC2.750

×1035.385

×10330005000*1399*

(含注释行)2826*

(含注释行)

适当设置输入扰动的初值，运用切线性模式可以简单求解输出变量对输入的偏导数。例如，对于一个含有个输入参数的实型函数

（3）

这里设，。运用DFT系统，可以得到对应的切线性模式

（4）

其中，为切线性模式的扰动输入参数。可以通过以下办法来求得偏导数：

（5）

其中。如果对于某个既是输入参数又是输出参数，可以类似以下过程引用的办法来处理。对于过程引用的情形，例如一个含有个输入参数的子过程

（6）

其中，为输入参数；，为输出参数；，既为输入参数又为输出参数。运用DFT系统，可以得到对应的切线性模式为

（7）

其中，，，分别为切线性模式的微分扰动输入、输出和输入输出参数。可以通过以下输入扰动设置并引用切线性模式（7）来求得偏导数：

a)设置；（，）；（）可以同时求得（）和（），其中。

b)设置（）；；（，）可以同时求得（）和（），其中。

4.3稀疏雅可比矩阵

运用上节讨论的方法来求解稀疏雅可比矩阵,具有极高的计算代价。例如，一个含个独立和个依赖参数的子过程，为求解整个雅可比矩阵就需要反复调用次切线性模式，当相当大时，这对许多实际的数值计算问题是不能接受的。事实上，如果雅可比矩阵的任意两列（行）相互正交，那么可以通过适当设置扰动输入值，这两列（行）的元素就可以通过一次引用切线性模式（伴随模式）完全得到。设和分别为雅可比矩阵的行宽度和列宽度，即各行和各列非零元素数目的最大值，显然有，。这里介绍几种常用的求解方法。

正向积分当时，通常采用切线性模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择右乘初始输入矩阵，可以获得接近的计算时间代价。DFT系统采用一种逐列（行）求解的方法，来有效求解右（左）乘初始输入矩阵。其基本思路是：按照某种列次序考察雅可比矩阵的各列；考察当前列中所有非零元素，并对这些非零元素所在行的行向量做类似模二和累加运算（即将非零元素视为逻辑“1”，零元素视为逻辑“0”），从而得到一个描述当前列与各行存在“某种”相关的标志向量（其元素都是“1”或“0”）；依据此标志向量，就很容易得到一个与之正交的列初始向量，其中与当前列序号对应的元素设置为“1”，而与标志向量中非零元素序号对应的元素设置为“0”，与标志向量中非零元素序号对应的元素设置为“-1”，显然，该列初始向量是唯一的，并且对应着当前右乘初始输入矩阵的最后一列；逐一考察已求解得到的列初始向量，如果某列初始向量与当前求解得到的列初始向量按下面定义的乘法（见过程4）正交，那么这两列就可以合并，即将当前列初始向量中非“-1”的元素按照对应关系分别赋值给该初始向量，并从记录中删除当前列初始向量；重复以上过程，继续按照给定列次序考察雅可比矩阵的“下一列”。不难说明，按照不同列次序求解得到的右乘初始输入矩阵可能不同。其中逐列求解右乘初始输入矩阵的过程可以简单叙述为：

1）将右乘初始输入矩阵所有元素的初值均设置为，，。。

2）如果，转6）。否则，如果雅可比矩阵的第列中的所有元素均为，，重复2）的判断。否则转3）。

3）计算标志向量。令，做如下计算：

，;

4）设为的列向量。在上定义乘法，对任意的，我们有：a）；b)如果，必有和。然后，做如下计算：

,;

，6);

2);

5）令，并做如下计算：

，;

令,。如果，转6）；否则，重复2）的判断。

6）对，，如果，则。取的前列，这样，我们就得到了一个维右乘初始输入矩阵。

这里需要说明的是，运用上面的方法求得的右乘初始输入矩阵不仅与求解雅可比矩阵的列序有关，而且与过程4）中的合并顺序也有关系。至于如何最优求解右乘初始输入矩阵，目前还很难讨论清楚。但是，大量模拟试验结果表明，运用上面自然次序求得的右乘初始输入矩阵宽度已经非常接近于其下界值。

反向积分当和时，通常采用伴随模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择左乘初始输入矩阵，可以获得接近的计算时间代价。其中左乘初始输入矩阵的求解过程完全可以按照上面的方法进行，但是在处理前必须先将雅可比矩阵转置，最后还需将得到的初始输入矩阵转置才能最终得到左乘初始输入矩阵。同时，其行宽度也已经非常接近于其下界值。

混合积分如果将切线性模式和伴随模式相结合，往往可以避免梯度向量运算中的诸多冗余计算。例如，ADJIFOR系统在求解雅可比矩阵时，在语句级微分实现中首先用伴随方法求得所有偏导数，然后做梯度向量积分；其计算时间代价与和模式的语句数目有关，而其存储代价为。具体讨论可参考文献[7]。

5．结论

切线性模式在无截断误差意义上计算函数的方向导数、梯度或雅可比矩阵，以及在模式的可预测性及参数敏感性分析、伴随模式构造等相关问题中有着广泛应用。DFT系统主要用于求解FORTRAN77语言编写的切线性模式，具有很强的全局数据相关分析能力。此外，DFT系统还具有其它几个重要特色，如结构化的微分实现、自动生成微分测试程序以及基于语句级的微分代码优化。本文简单给出了DFT系统在求解数值和符号导数和微分、稀疏雅可比矩阵中的应用。为评价一类自动微分系统，本文初步提出了统计准确率的概念。

参考文献

[1]AndreasGriewank.OnAutomaticDifferentiation.InM.IriandK.Tanabe,editors,MathematicalProgramming:

RecentDevelopmentsandApplications.KluwerAcademicPublishers,1989

[2]LeDimet,F.XandO.Talagrand,Variationalalgorithmsforanalysisandassimilationofmeteorological

observations:theoreticalaspects,Tellus,1986,38A,97-110

[3]P.Werbos,Applicationsofadvancesinnonlinearsensitivityanalysis,InsystemsModeling

andOptimization,NewYork,1982,SpringerVerlag,762-777

[4]ChristianBischof,GordonPusch,andRalfKnoesel."SensitivityAnalysisoftheMM5WeatherModelusing

AutomaticDifferentiation,"ComputersinPhysics,0:605-612,1996

[5]MuMu,etal，Thepredictabilityproblemofweatherandclimateprediction,ProgressinNatureScience,accepted.

[6]GieringR.etal.RecipesforAdjointCodeConstruction.ACMTrans.OnMath.Software.1998，24（4）：

437-474.

[7]C.Bischof,A.Carle,P.Khademi,andG.Pusch."AutomaticDifferentiation:ObtainingFastandReliable

Derivatives--Fast"inControlProblemsinIndustry,editedbyI.LasieckaandB.Morton,pages1-16,Birkhauser,

自动微分范文篇2

关键词自动微分切线性模式数据相关分析统计准确率

1.引言

计算微分大致经历了从商微分，符号微分，手写代码到自动微分几个阶段。与其它几种微分方法相比，自动微分具有代码简练、计算精度高及投入人力少等优点。自动微分实现的基本出发点是：一个数据相对独立的程序对象（模式、过程、程序段、数值语句乃至数值表达式），无论多么复杂，总可以分解为一系列有限数目的基本函数（如sin、exp、log）和基本运算操作（加、减、乘、除、乘方）的有序复合；对所有这些基本函数及基本运算操作，重复使用链式求导法则，将得到的中间结果自上而下地做正向积分就可以建立起对应的切线性模式，而自下而上地做反向积分就可以建立起对应的伴随模式[1]。基于自动微分方法得到的切线性模式和伴随模式，在变分资料同化[2]、系统建模与参数辨识[3]、参数的敏感性分析[4]、非线性最优化以及数值模式的可预测性分析[5]等问题中有着十分广泛的应用。

迄今为止，已有数十所大学和研究所各自开发了能够用于求解切线性模式的自动微分系统，比较典型的有TAMC系统[6]、ADJIFOR系统[7]和ODYSSEE系统[8]。在一些特定的运用中，它们都是比较成功的，但在通用性和复杂问题的处理效率上还存在许多不足。通常，自动生成切线性模式的关键难题在于对象自身的强相关性，这给系统全局分析（如数据IO相关分析和数据依赖相关分析）和微分代码的整体优化都带来了很多困难。同时，对于程序对象不可导处的准确识别和微分处理，至今仍还没有一个统一而有效的算法。另外，最优或有效求解稀疏雅可比矩阵一直是衡量一个自动微分系统有效性的重要尺度。

统计准确率被我们视为评价一类自动微分工具及其微分模式代码可靠性与有效性的重要尺度。其基本假设是：如果对于定义域空间内随机抽样获得的至多有限个n维初始场（或网格点），微分模式输出的差分和微分逼近是成功的；那么对于定义域空间内所有可能初始场（或网格点），微分模式输出的差分和微分逼近都是成功的。微分模式统计准确率评价的具体方法是：在所有随机抽样得到的初始场（或网格点）附近，当输入扰动逐渐趋向于机器有效精度所能表示的最小正值时，模式输出的差分和微分之间应该有足够精度有效位数上的逼近。

DFT系统具有许多优点，它能够完全接受用FORTRAN77语言编写的源代码，微分代码结构清晰，其微分处理能力与问题和对象的规模及复杂性无关。它基于YACC实现，具有很强的可扩展性。DFT系统具有四个重要特色。它通过对象全局依赖相关分析，准确求解雅可比矩阵的稀疏结构，自动计算有效初始输入矩阵，从而可以用较小的代价求得整个雅可比矩阵。同时，它可以自动生成客观评价微分模式效率与可靠性的测试程序，对奇异函数做等价微分处理，并采用二元归约的方法，在语句级层次上实现微分代码优化。

2.系统概况

DFT系统主要由两部分组成：微分代码转换和微分代码评价，图2.1。微分代码转换部分接受用户输入指令并自动分析对象模式，生成切线性模式代码及其相关测试代码，后者直接构成微分代码评价系统的主体。微分代码评价是DFT系统的一个重要特色。DFT系统的开发小组认为，一个微分模式如果在可靠性、时间和存储效率上没有得到充分的验证，至少对实际应用而言，它将是毫无意义的。

原模式切线性模式

统计评价结果

图2.1DFT系统结构简图

2.1微分代码转换

DFT系统是基于YACC在UNIX环境下开发的，其结构图2.2所示。通过DFT系统产生的切线性模式代码成对出现，并在语句级程度上做了简化，可读性很强，如图2.4。

切线性模式

评价函数集

图2.2微分代码转换

微分代码转换部分从功能上分为四个部分：词法分析，语义分析，对象复杂性及数据相关分析和微分代码转换。对于一组具有复杂数据相关的程序模式对象，通常需要系统运行两遍才能得到有效而可靠的微分代码。这主要有两方面的考虑：其一，根据对象的复杂性（如最大语句长度、最大变量维数、子过程或函数数目、子过程或函数内最大变量数目等对象特征）选择合适的系统参数以求最优的运行代价；其二，模式内各子过程或函数之间以及一个子过程或函数内往往具有很强的数据相关性，需要事先保存对象的相关信息并且在考虑当前对象的属性之前必须做上下文相关分析。

图2.3PERIGEE源程序代码图2.4DFT系统生成的切线性代码

2.2微分代码评价

通常，评价一个编译系统的性能有很多方面，如处理速度、结果代码可靠性及质量、出错诊断、可扩展和可维护性等。对于一类自动微分系统来说，由于软件开发人力的局限以及对象模式的复杂多样性,通过自动转换得到的微分模式并非常常是有效而可靠的（即无论是在数学意义上还是在程序逻辑上应与期待的理想结果一致），因而在微分模式被投入实际应用前，往往需要投入一定的人力来对其做严格的分析测试。

对切线性模式做统计评价测试的主要内容可以简单叙述为：在网格化的模式定义域空间内，选择所有可能的网格点形成微分模式计算的初始场；在不同的网格点附近，随机选取至少个线性无关的初始扰动，对每个扰动输入分别进行网格点逼近，统计考察模式输出差分和微分在有效位数上的逼近程度。图2.5描述了整个测试过程，它包含网格点数据随机采样（1）和网格点数据逼近（2）两级循环。

图2.5切线性模式代码的测试过程

3.系统主要特色

DFT系统并不是一个完整的FORTRAN编译器，但它几乎可以接受和处理所有FORTRAN77编写的源模式代码，并且可以很方便地扩展并接受FORTRAN90编写的源模式代码。本节将着重介绍DFT系统（版本3.0）的以下几个重要特色。

3.1结构化的微分实现

DFT系统采用标准化的代码实现，切线性模式的扰动变量和基态值变量、微分计算语句和基态值计算语句总是成对出现，并具有清晰的程序结构。微分代码保持了原模式本身的结构和风格（如并行和向量特性、数据精度等），即语句到语句、结构到结构的微分实现。在奇异点或不可导处，DFT系统对微分扰动采取简单的清零处理，实践证明这对抑制扰动计算溢出具有重要意义，但并不影响评价测试结果。

3.2全局数据相关分析

DFT系统具有较强的数据相关分析能力，它包括全局数据IO相关分析、全局数据依赖相关分析、全局过程相关分析以及数据迭代相关分析几个不同方面。数据依赖相关与数据IO相关关系密切，但又存在根本不同。前者强调每个变量在数学关系上的依赖性；而后者描述了一个对象的输入输出特性，且具有相对性，即任何一个变量参数，无论它是独立变量还是依赖变量，在数学意义上都可等价为一个既是输入又是输出的参数来处理。

DFT系统记录所有过程参数的IO属性表，通过深度递归相关计算，准确计算每个过程参数的最终IO属性。DFT系统通过对数据相关矩阵做模二和及自乘迭代计算（An+1=An⊕An2）来完成数据的依赖相关分析，这种算法具有很好的对数收敛特性。DFT系统通过全局过程相关分析的结果，自动生成模式的局部或整体相关引用树结构（如图3.1），这对用户分析复杂数值模式和微分评价测试都具有很好的指导作用。DFT系统还具有分析局部数据迭代相关和函数迭代相关的能力，这两种形式的数据迭代相关是自动微分实现颇具挑战的难题之一。

图3.1GPSRayshooting模式的相关树结构片段

3.3自动生成测试程序

基于IO相关分析的结果，DFT系统自动生成微分测试代码，分别对切线性模式的可靠性和运行代价做统计评价测试。特别地，DFT系统还可将任何模式参数都视为输入输出参数，生成在数学意义上等价的测试代码，这样处理的不利之处在于往往需要极高的存储开销。

3.4基于语句级的代码优化

目前，DFT系统仅仅具备局地优化能力。在语句级微分实现上采用二元归约的方法对微分代码进行优化是DFT系统的一个重要特色。根据右端表达式的乘法复杂性及含变元数目的不同，DFT系统采取不同的分解策略。二元归约的方法避免了微分计算中的许多冗余计算，在一些复杂的非线性表达式的微分计算中具有最小的计算代价，同时也非常适合于微分系统的软件实现。同时，对于某些特殊的运算操作（除法、乘方）和特殊函数（如sqrt、exp），DFT系统较好地利用了基态值计算得到的中间结果，避免了微分实现中的冗余计算。

4.系统应用

运用自动微分工具得到的切线性模式，可以在无截断误差意义下求解函数的数值微分和导数、稀疏雅可比矩阵。同时这些结果在数值参数敏感性分析、非线性最优化以及其它数值理论分析中有着非常重要的应用。这里简单介绍切线性模式的几个基本应用。

4.1符号导数和微分

如果输入为数学关系式，DFT系统可以自动生成对应的微分表达式和梯度，而与数学关系式的复杂程度无关。例如我们输入关系式：

，（1）

DFT系统将自动生成其符号微分形式及其梯度形式分别为

，（2）

4.2数值导数和微分

切线性模式最基本的应用就是在一定扰动输入下求解输出变量的扰动（响应）。表4.1给出了DFT系统在对IAP9L模式、GPSRayshooting模式和GPSRaytrace模式三个数值模式做切线性化的具体应用中，一些不同计算粒度、不同引用深度和不同程序风格的核心子过程，以及它们的切线性模式在SGI2000上运行的统计评价测试结果，其中切线性模式的可靠性指标都准确到六个有效数字以上，在运行时间、存储开销和代码复杂性方面分别是原模式的两倍左右，比较接近于理想的微分代价结果（1.5倍）。除了IAP9L模式由于过于复杂仅做粗略统计外，其余模式都用非注释语句行数来表示各自的代码复杂性。

表4.1DFT系统在三个数值模式中的统计评价测试结果

性能指标

对象模式运行时间（10-3秒）存储开销（字节数）代码复杂性

原模式切线性

模式

原模式切线性

模式

原模式切线性

模式

Xyz2g2.5306.1605524110485589

IntCIRA1.5602.750133426614165

Dabel0.0350.072601202749

LSS8.30017.50669133879143

RP42.4085.10360572102238

Vgrad10.1000.21218564368282454

RefGr43.0086.0071865414373083578

LL2JK0.6261.350262252442232

RayFind462.70

×103125.4

×103985618212111179

EPSIMP1.76011.50445589101327

Hlimits0.8301.8802425774842543774

Int3sL26.9051.2082002916394584690

MAKE

NCEP1340392072292514458504584

Curvcent0.0130.038527542754

DYFRM3.800

×1037.250

×1035000*9500*161279

PHYSIC2.750

×1035.385

×10330005000*1399*

(含注释行)2826*

(含注释行)

适当设置输入扰动的初值，运用切线性模式可以简单求解输出变量对输入的偏导数。例如，对于一个含有个输入参数的实型函数

（3）

这里设，。运用DFT系统，可以得到对应的切线性模式

（4）

其中，为切线性模式的扰动输入参数。可以通过以下办法来求得偏导数：

（5）

其中。如果对于某个既是输入参数又是输出参数，可以类似以下过程引用的办法来处理。对于过程引用的情形，例如一个含有个输入参数的子过程

（6）

其中，为输入参数；，为输出参数；，既为输入参数又为输出参数。运用DFT系统，可以得到对应的切线性模式为

（7）

其中，，，分别为切线性模式的微分扰动输入、输出和输入输出参数。可以通过以下输入扰动设置并引用切线性模式（7）来求得偏导数：

a)设置；（，）；（）可以同时求得（）和（），其中。

b)设置（）；；（，）可以同时求得（）和（），其中。

4.3稀疏雅可比矩阵

运用上节讨论的方法来求解稀疏雅可比矩阵,具有极高的计算代价。例如，一个含个独立和个依赖参数的子过程，为求解整个雅可比矩阵就需要反复调用次切线性模式，当相当大时，这对许多实际的数值计算问题是不能接受的。事实上，如果雅可比矩阵的任意两列（行）相互正交，那么可以通过适当设置扰动输入值，这两列（行）的元素就可以通过一次引用切线性模式（伴随模式）完全得到。设和分别为雅可比矩阵的行宽度和列宽度，即各行和各列非零元素数目的最大值，显然有，。这里介绍几种常用的求解方法。

正向积分当时，通常采用切线性模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择右乘初始输入矩阵，可以获得接近的计算时间代价。DFT系统采用一种逐列（行）求解的方法，来有效求解右（左）乘初始输入矩阵。其基本思路是：按照某种列次序考察雅可比矩阵的各列；考察当前列中所有非零元素，并对这些非零元素所在行的行向量做类似模二和累加运算（即将非零元素视为逻辑“1”，零元素视为逻辑“0”），从而得到一个描述当前列与各行存在“某种”相关的标志向量（其元素都是“1”或“0”）；依据此标志向量，就很容易得到一个与之正交的列初始向量，其中与当前列序号对应的元素设置为“1”，而与标志向量中非零元素序号对应的元素设置为“0”，与标志向量中非零元素序号对应的元素设置为“-1”，显然，该列初始向量是唯一的，并且对应着当前右乘初始输入矩阵的最后一列；逐一考察已求解得到的列初始向量，如果某列初始向量与当前求解得到的列初始向量按下面定义的乘法（见过程4）正交，那么这两列就可以合并，即将当前列初始向量中非“-1”的元素按照对应关系分别赋值给该初始向量，并从记录中删除当前列初始向量；重复以上过程，继续按照给定列次序考察雅可比矩阵的“下一列”。不难说明，按照不同列次序求解得到的右乘初始输入矩阵可能不同。其中逐列求解右乘初始输入矩阵的过程可以简单叙述为：

1）将右乘初始输入矩阵所有元素的初值均设置为，，。。

2）如果，转6）。否则，如果雅可比矩阵的第列中的所有元素均为，，重复2）的判断。否则转3）。

3）计算标志向量。令，做如下计算：

，;

4）设为的列向量。在上定义乘法，对任意的，我们有：a）；b)如果，必有和。然后，做如下计算：

,;

，6);

2);

5）令，并做如下计算：

，;

令,。如果，转6）；否则，重复2）的判断。

6）对，，如果，则。取的前列，这样，我们就得到了一个维右乘初始输入矩阵。

这里需要说明的是，运用上面的方法求得的右乘初始输入矩阵不仅与求解雅可比矩阵的列序有关，而且与过程4）中的合并顺序也有关系。至于如何最优求解右乘初始输入矩阵，目前还很难讨论清楚。但是，大量模拟试验结果表明，运用上面自然次序求得的右乘初始输入矩阵宽度已经非常接近于其下界值。

反向积分当和时，通常采用伴随模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择左乘初始输入矩阵，可以获得接近的计算时间代价。其中左乘初始输入矩阵的求解过程完全可以按照上面的方法进行，但是在处理前必须先将雅可比矩阵转置，最后还需将得到的初始输入矩阵转置才能最终得到左乘初始输入矩阵。同时，其行宽度也已经非常接近于其下界值。

混合积分如果将切线性模式和伴随模式相结合，往往可以避免梯度向量运算中的诸多冗余计算。例如，ADJIFOR系统在求解雅可比矩阵时，在语句级微分实现中首先用伴随方法求得所有偏导数，然后做梯度向量积分；其计算时间代价与和模式的语句数目有关，而其存储代价为。具体讨论可参考文献[7]。

5．结论

切线性模式在无截断误差意义上计算函数的方向导数、梯度或雅可比矩阵，以及在模式的可预测性及参数敏感性分析、伴随模式构造等相关问题中有着广泛应用。DFT系统主要用于求解FORTRAN77语言编写的切线性模式，具有很强的全局数据相关分析能力。此外，DFT系统还具有其它几个重要特色，如结构化的微分实现、自动生成微分测试程序以及基于语句级的微分代码优化。本文简单给出了DFT系统在求解数值和符号导数和微分、稀疏雅可比矩阵中的应用。为评价一类自动微分系统，本文初步提出了统计准确率的概念。

参考文献

[1]AndreasGriewank.OnAutomaticDifferentiation.InM.IriandK.Tanabe,editors,MathematicalProgramming:

RecentDevelopmentsandApplications.KluwerAcademicPublishers,1989

[2]LeDimet,F.XandO.Talagrand,Variationalalgorithmsforanalysisandassimilationofmeteorological

observations:theoreticalaspects,Tellus,1986,38A,97-110

[3]P.Werbos,Applicationsofadvancesinnonlinearsensitivityanalysis,InsystemsModeling

andOptimization,NewYork,1982,SpringerVerlag,762-777

[4]ChristianBischof,GordonPusch,andRalfKnoesel."SensitivityAnalysisoftheMM5WeatherModelusing

AutomaticDifferentiation,"ComputersinPhysics,0:605-612,1996

[5]MuMu,etal，Thepredictabilityproblemofweatherandclimateprediction,ProgressinNatureScience,accepted.

[6]GieringR.etal.RecipesforAdjointCodeConstruction.ACMTrans.OnMath.Software.1998，24（4）：

437-474.

[7]C.Bischof,A.Carle,P.Khademi,andG.Pusch."AutomaticDifferentiation:ObtainingFastandReliable

Derivatives--Fast"inControlProblemsinIndustry,editedbyI.LasieckaandB.Morton,pages1-16,Birkhauser,

自动微分范文篇3

关键词自动微分切线性模式数据相关分析统计准确率

1.引言

计算微分大致经历了从商微分，符号微分，手写代码到自动微分几个阶段。与其它几种微分方法相比，自动微分具有代码简练、计算精度高及投入人力少等优点。自动微分实现的基本出发点是：一个数据相对独立的程序对象（模式、过程、程序段、数值语句乃至数值表达式），无论多么复杂，总可以分解为一系列有限数目的基本函数（如sin、exp、log）和基本运算操作（加、减、乘、除、乘方）的有序复合；对所有这些基本函数及基本运算操作，重复使用链式求导法则，将得到的中间结果自上而下地做正向积分就可以建立起对应的切线性模式，而自下而上地做反向积分就可以建立起对应的伴随模式[1]。基于自动微分方法得到的切线性模式和伴随模式，在变分资料同化[2]、系统建模与参数辨识[3]、参数的敏感性分析[4]、非线性最优化以及数值模式的可预测性分析[5]等问题中有着十分广泛的应用。

迄今为止，已有数十所大学和研究所各自开发了能够用于求解切线性模式的自动微分系统，比较典型的有TAMC系统[6]、ADJIFOR系统[7]和ODYSSEE系统[8]。在一些特定的运用中，它们都是比较成功的，但在通用性和复杂问题的处理效率上还存在许多不足。通常，自动生成切线性模式的关键难题在于对象自身的强相关性，这给系统全局分析（如数据IO相关分析和数据依赖相关分析）和微分代码的整体优化都带来了很多困难。同时，对于程序对象不可导处的准确识别和微分处理，至今仍还没有一个统一而有效的算法。另外，最优或有效求解稀疏雅可比矩阵一直是衡量一个自动微分系统有效性的重要尺度。

统计准确率被我们视为评价一类自动微分工具及其微分模式代码可靠性与有效性的重要尺度。其基本假设是：如果对于定义域空间内随机抽样获得的至多有限个n维初始场（或网格点），微分模式输出的差分和微分逼近是成功的；那么对于定义域空间内所有可能初始场（或网格点），微分模式输出的差分和微分逼近都是成功的。微分模式统计准确率评价的具体方法是：在所有随机抽样得到的初始场（或网格点）附近，当输入扰动逐渐趋向于机器有效精度所能表示的最小正值时，模式输出的差分和微分之间应该有足够精度有效位数上的逼近。

DFT系统具有许多优点，它能够完全接受用FORTRAN77语言编写的源代码，微分代码结构清晰，其微分处理能力与问题和对象的规模及复杂性无关。它基于YACC实现，具有很强的可扩展性。DFT系统具有四个重要特色。它通过对象全局依赖相关分析，准确求解雅可比矩阵的稀疏结构，自动计算有效初始输入矩阵，从而可以用较小的代价求得整个雅可比矩阵。同时，它可以自动生成客观评价微分模式效率与可靠性的测试程序，对奇异函数做等价微分处理，并采用二元归约的方法，在语句级层次上实现微分代码优化。

2.系统概况

DFT系统主要由两部分组成：微分代码转换和微分代码评价，图2.1。微分代码转换部分接受用户输入指令并自动分析对象模式，生成切线性模式代码及其相关测试代码，后者直接构成微分代码评价系统的主体。微分代码评价是DFT系统的一个重要特色。DFT系统的开发小组认为，一个微分模式如果在可靠性、时间和存储效率上没有得到充分的验证，至少对实际应用而言，它将是毫无意义的。

原模式切线性模式

统计评价结果

图2.1DFT系统结构简图

2.1微分代码转换

DFT系统是基于YACC在UNIX环境下开发的，其结构图2.2所示。通过DFT系统产生的切线性模式代码成对出现，并在语句级程度上做了简化，可读性很强，如图2.4。

切线性模式

评价函数集

图2.2微分代码转换

微分代码转换部分从功能上分为四个部分：词法分析，语义分析，对象复杂性及数据相关分析和微分代码转换。对于一组具有复杂数据相关的程序模式对象，通常需要系统运行两遍才能得到有效而可靠的微分代码。这主要有两方面的考虑：其一，根据对象的复杂性（如最大语句长度、最大变量维数、子过程或函数数目、子过程或函数内最大变量数目等对象特征）选择合适的系统参数以求最优的运行代价；其二，模式内各子过程或函数之间以及一个子过程或函数内往往具有很强的数据相关性，需要事先保存对象的相关信息并且在考虑当前对象的属性之前必须做上下文相关分析。

图2.3PERIGEE源程序代码图2.4DFT系统生成的切线性代码

2.2微分代码评价

通常，评价一个编译系统的性能有很多方面，如处理速度、结果代码可靠性及质量、出错诊断、可扩展和可维护性等。对于一类自动微分系统来说，由于软件开发人力的局限以及对象模式的复杂多样性,通过自动转换得到的微分模式并非常常是有效而可靠的（即无论是在数学意义上还是在程序逻辑上应与期待的理想结果一致），因而在微分模式被投入实际应用前，往往需要投入一定的人力来对其做严格的分析测试。

对切线性模式做统计评价测试的主要内容可以简单叙述为：在网格化的模式定义域空间内，选择所有可能的网格点形成微分模式计算的初始场；在不同的网格点附近，随机选取至少个线性无关的初始扰动，对每个扰动输入分别进行网格点逼近，统计考察模式输出差分和微分在有效位数上的逼近程度。图2.5描述了整个测试过程，它包含网格点数据随机采样（1）和网格点数据逼近（2）两级循环。

图2.5切线性模式代码的测试过程

3.系统主要特色

DFT系统并不是一个完整的FORTRAN编译器，但它几乎可以接受和处理所有FORTRAN77编写的源模式代码，并且可以很方便地扩展并接受FORTRAN90编写的源模式代码。本节将着重介绍DFT系统（版本3.0）的以下几个重要特色。

3.1结构化的微分实现

DFT系统采用标准化的代码实现，切线性模式的扰动变量和基态值变量、微分计算语句和基态值计算语句总是成对出现，并具有清晰的程序结构。微分代码保持了原模式本身的结构和风格（如并行和向量特性、数据精度等），即语句到语句、结构到结构的微分实现。在奇异点或不可导处，DFT系统对微分扰动采取简单的清零处理，实践证明这对抑制扰动计算溢出具有重要意义，但并不影响评价测试结果。

3.2全局数据相关分析

DFT系统具有较强的数据相关分析能力，它包括全局数据IO相关分析、全局数据依赖相关分析、全局过程相关分析以及数据迭代相关分析几个不同方面。数据依赖相关与数据IO相关关系密切，但又存在根本不同。前者强调每个变量在数学关系上的依赖性；而后者描述了一个对象的输入输出特性，且具有相对性，即任何一个变量参数，无论它是独立变量还是依赖变量，在数学意义上都可等价为一个既是输入又是输出的参数来处理。

DFT系统记录所有过程参数的IO属性表，通过深度递归相关计算，准确计算每个过程参数的最终IO属性。DFT系统通过对数据相关矩阵做模二和及自乘迭代计算（An+1=An⊕An2）来完成数据的依赖相关分析，这种算法具有很好的对数收敛特性。DFT系统通过全局过程相关分析的结果，自动生成模式的局部或整体相关引用树结构（如图3.1），这对用户分析复杂数值模式和微分评价测试都具有很好的指导作用。DFT系统还具有分析局部数据迭代相关和函数迭代相关的能力，这两种形式的数据迭代相关是自动微分实现颇具挑战的难题之一。

3.3自动生成测试程序

基于IO相关分析的结果，DFT系统自动生成微分测试代码，分别对切线性模式的可靠性和运行代价做统计评价测试。特别地，DFT系统还可将任何模式参数都视为输入输出参数，生成在数学意义上等价的测试代码，这样处理的不利之处在于往往需要极高的存储开销。

3.4基于语句级的代码优化

目前，DFT系统仅仅具备局地优化能力。在语句级微分实现上采用二元归约的方法对微分代码进行优化是DFT系统的一个重要特色。根据右端表达式的乘法复杂性及含变元数目的不同，DFT系统采取不同的分解策略。二元归约的方法避免了微分计算中的许多冗余计算，在一些复杂的非线性表达式的微分计算中具有最小的计算代价，同时也非常适合于微分系统的软件实现。同时，对于某些特殊的运算操作（除法、乘方）和特殊函数（如sqrt、exp），DFT系统较好地利用了基态值计算得到的中间结果，避免了微分实现中的冗余计算。

4.系统应用

运用自动微分工具得到的切线性模式，可以在无截断误差意义下求解函数的数值微分和导数、稀疏雅可比矩阵。同时这些结果在数值参数敏感性分析、非线性最优化以及其它数值理论分析中有着非常重要的应用。这里简单介绍切线性模式的几个基本应用。

4.1符号导数和微分

如果输入为数学关系式，DFT系统可以自动生成对应的微分表达式和梯度，而与数学关系式的复杂程度无关。例如我们输入关系式：

，（1）

DFT系统将自动生成其符号微分形式及其梯度形式分别为

，（2）

4.2数值导数和微分

切线性模式最基本的应用就是在一定扰动输入下求解输出变量的扰动（响应）。表4.1给出了DFT系统在对IAP9L模式、GPSRayshooting模式和GPSRaytrace模式三个数值模式做切线性化的具体应用中，一些不同计算粒度、不同引用深度和不同程序风格的核心子过程，以及它们的切线性模式在SGI2000上运行的统计评价测试结果，其中切线性模式的可靠性指标都准确到六个有效数字以上，在运行时间、存储开销和代码复杂性方面分别是原模式的两倍左右，比较接近于理想的微分代价结果（1.5倍）。除了IAP9L模式由于过于复杂仅做粗略统计外，其余模式都用非注释语句行数来表示各自的代码复杂性。

表4.1DFT系统在三个数值模式中的统计评价测试结果

性能指标

对象模式运行时间（10-3秒）存储开销（字节数）代码复杂性

原模式切线性

模式

原模式切线性

模式

原模式切线性

模式

Xyz2g2.5306.1605524110485589

IntCIRA1.5602.750133426614165

Dabel0.0350.072601202749

LSS8.30017.50669133879143

RP42.4085.10360572102238

Vgrad10.1000.21218564368282454

RefGr43.0086.0071865414373083578

LL2JK0.6261.350262252442232

RayFind462.70

×103125.4

×103985618212111179

EPSIMP1.76011.50445589101327

Hlimits0.8301.8802425774842543774

Int3sL26.9051.2082002916394584690

MAKE

NCEP1340392072292514458504584

Curvcent0.0130.038527542754

DYFRM3.800

×1037.250

×1035000*9500*161279

PHYSIC2.750

×1035.385

×10330005000*1399*

(含注释行)2826*

(含注释行)

适当设置输入扰动的初值，运用切线性模式可以简单求解输出变量对输入的偏导数。例如，对于一个含有个输入参数的实型函数

（3）

这里设，。运用DFT系统，可以得到对应的切线性模式

（4）

其中，为切线性模式的扰动输入参数。可以通过以下办法来求得偏导数：

（5）

其中。如果对于某个既是输入参数又是输出参数，可以类似以下过程引用的办法来处理。对于过程引用的情形，例如一个含有个输入参数的子过程

（6）

其中，为输入参数；，为输出参数；，既为输入参数又为输出参数。运用DFT系统，可以得到对应的切线性模式为

（7）

其中，，，分别为切线性模式的微分扰动输入、输出和输入输出参数。可以通过以下输入扰动设置并引用切线性模式（7）来求得偏导数：

a)设置；（，）；（）可以同时求得（）和（），其中。

b)设置（）；；（，）可以同时求得（）和（），其中。

4.3稀疏雅可比矩阵

运用上节讨论的方法来求解稀疏雅可比矩阵,具有极高的计算代价。例如，一个含个独立和个依赖参数的子过程，为求解整个雅可比矩阵就需要反复调用次切线性模式，当相当大时，这对许多实际的数值计算问题是不能接受的。事实上，如果雅可比矩阵的任意两列（行）相互正交，那么可以通过适当设置扰动输入值，这两列（行）的元素就可以通过一次引用切线性模式（伴随模式）完全得到。设和分别为雅可比矩阵的行宽度和列宽度，即各行和各列非零元素数目的最大值，显然有，。这里介绍几种常用的求解方法。

正向积分当时，通常采用切线性模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择右乘初始输入矩阵，可以获得接近的计算时间代价。DFT系统采用一种逐列（行）求解的方法，来有效求解右（左）乘初始输入矩阵。其基本思路是：按照某种列次序考察雅可比矩阵的各列；考察当前列中所有非零元素，并对这些非零元素所在行的行向量做类似模二和累加运算（即将非零元素视为逻辑“1”，零元素视为逻辑“0”），从而得到一个描述当前列与各行存在“某种”相关的标志向量（其元素都是“1”或“0”）；依据此标志向量，就很容易得到一个与之正交的列初始向量，其中与当前列序号对应的元素设置为“1”，而与标志向量中非零元素序号对应的元素设置为“0”，与标志向量中非零元素序号对应的元素设置为“-1”，显然，该列初始向量是唯一的，并且对应着当前右乘初始输入矩阵的最后一列；逐一考察已求解得到的列初始向量，如果某列初始向量与当前求解得到的列初始向量按下面定义的乘法（见过程4）正交，那么这两列就可以合并，即将当前列初始向量中非“-1”的元素按照对应关系分别赋值给该初始向量，并从记录中删除当前列初始向量；重复以上过程，继续按照给定列次序考察雅可比矩阵的“下一列”。不难说明，按照不同列次序求解得到的右乘初始输入矩阵可能不同。其中逐列求解右乘初始输入矩阵的过程可以简单叙述为：

1）将右乘初始输入矩阵所有元素的初值均设置为，，。。

2）如果，转6）。否则，如果雅可比矩阵的第列中的所有元素均为，，重复2）的判断。否则转3）。

3）计算标志向量。令，做如下计算：

，;

4）设为的列向量。在上定义乘法，对任意的，我们有：a）；b)如果，必有和。然后，做如下计算：

,;

，6);

2);

5）令，并做如下计算：

，;

令,。如果，转6）；否则，重复2）的判断。

6）对，，如果，则。取的前列，这样，我们就得到了一个维右乘初始输入矩阵。

这里需要说明的是，运用上面的方法求得的右乘初始输入矩阵不仅与求解雅可比矩阵的列序有关，而且与过程4）中的合并顺序也有关系。至于如何最优求解右乘初始输入矩阵，目前还很难讨论清楚。但是，大量模拟试验结果表明，运用上面自然次序求得的右乘初始输入矩阵宽度已经非常接近于其下界值。

反向积分当和时，通常采用伴随模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择左乘初始输入矩阵，可以获得接近的计算时间代价。其中左乘初始输入矩阵的求解过程完全可以按照上面的方法进行，但是在处理前必须先将雅可比矩阵转置，最后还需将得到的初始输入矩阵转置才能最终得到左乘初始输入矩阵。同时，其行宽度也已经非常接近于其下界值。

混合积分如果将切线性模式和伴随模式相结合，往往可以避免梯度向量运算中的诸多冗余计算。例如，ADJIFOR系统在求解雅可比矩阵时，在语句级微分实现中首先用伴随方法求得所有偏导数，然后做梯度向量积分；其计算时间代价与和模式的语句数目有关，而其存储代价为。具体讨论可参考文献[7]。

自动微分范文篇4

关键词自动微分切线性模式数据相关分析统计准确率

1.引言

计算微分大致经历了从商微分，符号微分，手写代码到自动微分几个阶段。与其它几种微分方法相比，自动微分具有代码简练、计算精度高及投入人力少等优点。自动微分实现的基本出发点是：一个数据相对独立的程序对象（模式、过程、程序段、数值语句乃至数值表达式），无论多么复杂，总可以分解为一系列有限数目的基本函数（如sin、exp、log）和基本运算操作（加、减、乘、除、乘方）的有序复合；对所有这些基本函数及基本运算操作，重复使用链式求导法则，将得到的中间结果自上而下地做正向积分就可以建立起对应的切线性模式，而自下而上地做反向积分就可以建立起对应的伴随模式[1]。基于自动微分方法得到的切线性模式和伴随模式，在变分资料同化[2]、系统建模与参数辨识[3]、参数的敏感性分析[4]、非线性最优化以及数值模式的可预测性分析[5]等问题中有着十分广泛的应用。

迄今为止，已有数十所大学和研究所各自开发了能够用于求解切线性模式的自动微分系统，比较典型的有TAMC系统[6]、ADJIFOR系统[7]和ODYSSEE系统[8]。在一些特定的运用中，它们都是比较成功的，但在通用性和复杂问题的处理效率上还存在许多不足。通常，自动生成切线性模式的关键难题在于对象自身的强相关性，这给系统全局分析（如数据IO相关分析和数据依赖相关分析）和微分代码的整体优化都带来了很多困难。同时，对于程序对象不可导处的准确识别和微分处理，至今仍还没有一个统一而有效的算法。另外，最优或有效求解稀疏雅可比矩阵一直是衡量一个自动微分系统有效性的重要尺度。

统计准确率被我们视为评价一类自动微分工具及其微分模式代码可靠性与有效性的重要尺度。其基本假设是：如果对于定义域空间内随机抽样获得的至多有限个n维初始场（或网格点），微分模式输出的差分和微分逼近是成功的；那么对于定义域空间内所有可能初始场（或网格点），微分模式输出的差分和微分逼近都是成功的。微分模式统计准确率评价的具体方法是：在所有随机抽样得到的初始场（或网格点）附近，当输入扰动逐渐趋向于机器有效精度所能表示的最小正值时，模式输出的差分和微分之间应该有足够精度有效位数上的逼近。

DFT系统具有许多优点，它能够完全接受用FORTRAN77语言编写的源代码，微分代码结构清晰，其微分处理能力与问题和对象的规模及复杂性无关。它基于YACC实现，具有很强的可扩展性。DFT系统具有四个重要特色。它通过对象全局依赖相关分析，准确求解雅可比矩阵的稀疏结构，自动计算有效初始输入矩阵，从而可以用较小的代价求得整个雅可比矩阵。同时，它可以自动生成客观评价微分模式效率与可靠性的测试程序，对奇异函数做等价微分处理，并采用二元归约的方法，在语句级层次上实现微分代码优化。

2.系统概况

DFT系统主要由两部分组成：微分代码转换和微分代码评价，图2.1。微分代码转换部分接受用户输入指令并自动分析对象模式，生成切线性模式代码及其相关测试代码，后者直接构成微分代码评价系统的主体。微分代码评价是DFT系统的一个重要特色。DFT系统的开发小组认为，一个微分模式如果在可靠性、时间和存储效率上没有得到充分的验证，至少对实际应用而言，它将是毫无意义的。

评价函数集

图2.2微分代码转换

微分代码转换部分从功能上分为四个部分：词法分析，语义分析，对象复杂性及数据相关分析和微分代码转换。对于一组具有复杂数据相关的程序模式对象，通常需要系统运行两遍才能得到有效而可靠的微分代码。这主要有两方面的考虑：其一，根据对象的复杂性（如最大语句长度、最大变量维数、子过程或函数数目、子过程或函数内最大变量数目等对象特征）选择合适的系统参数以求最优的运行代价；其二，模式内各子过程或函数之间以及一个子过程或函数内往往具有很强的数据相关性，需要事先保存对象的相关信息并且在考虑当前对象的属性之前必须做上下文相关分析。

2.2微分代码评价

通常，评价一个编译系统的性能有很多方面，如处理速度、结果代码可靠性及质量、出错诊断、可扩展和可维护性等。对于一类自动微分系统来说，由于软件开发人力的局限以及对象模式的复杂多样性,通过自动转换得到的微分模式并非常常是有效而可靠的（即无论是在数学意义上还是在程序逻辑上应与期待的理想结果一致），因而在微分模式被投入实际应用前，往往需要投入一定的人力来对其做严格的分析测试。

对切线性模式做统计评价测试的主要内容可以简单叙述为：在网格化的模式定义域空间内，选择所有可能的网格点形成微分模式计算的初始场；在不同的网格点附近，随机选取至少个线性无关的初始扰动，对每个扰动输入分别进行网格点逼近，统计考察模式输出差分和微分在有效位数上的逼近程度。图2.5描述了整个测试过程，它包含网格点数据随机采样（1）和网格点数据逼近（2）两级循环。

图2.5切线性模式代码的测试过程

3.系统主要特色

DFT系统并不是一个完整的FORTRAN编译器，但它几乎可以接受和处理所有FORTRAN77编写的源模式代码，并且可以很方便地扩展并接受FORTRAN90编写的源模式代码。本节将着重介绍DFT系统（版本3.0）的以下几个重要特色。

3.1结构化的微分实现

DFT系统采用标准化的代码实现，切线性模式的扰动变量和基态值变量、微分计算语句和基态值计算语句总是成对出现，并具有清晰的程序结构。微分代码保持了原模式本身的结构和风格（如并行和向量特性、数据精度等），即语句到语句、结构到结构的微分实现。在奇异点或不可导处，DFT系统对微分扰动采取简单的清零处理，实践证明这对抑制扰动计算溢出具有重要意义，但并不影响评价测试结果。

3.2全局数据相关分析

DFT系统具有较强的数据相关分析能力，它包括全局数据IO相关分析、全局数据依赖相关分析、全局过程相关分析以及数据迭代相关分析几个不同方面。数据依赖相关与数据IO相关关系密切，但又存在根本不同。前者强调每个变量在数学关系上的依赖性；而后者描述了一个对象的输入输出特性，且具有相对性，即任何一个变量参数，无论它是独立变量还是依赖变量，在数学意义上都可等价为一个既是输入又是输出的参数来处理。

DFT系统记录所有过程参数的IO属性表，通过深度递归相关计算，准确计算每个过程参数的最终IO属性。DFT系统通过对数据相关矩阵做模二和及自乘迭代计算（An+1=An⊕An2）来完成数据的依赖相关分析，这种算法具有很好的对数收敛特性。DFT系统通过全局过程相关分析的结果，自动生成模式的局部或整体相关引用树结构（如图3.1），这对用户分析复杂数值模式和微分评价测试都具有很好的指导作用。DFT系统还具有分析局部数据迭代相关和函数迭代相关的能力，这两种形式的数据迭代相关是自动微分

实现颇具挑战的难题之一。

3.3自动生成测试程序

基于IO相关分析的结果，DFT系统自动生成微分测试代码，分别对切线性模式的可靠性和运行代价做统计评价测试。特别地，DFT系统还可将任何模式参数都视为输入输出参数，生成在数学意义上等价的测试代码，这样处理的不利之处在于往往需要极高的存储开销。

3.4基于语句级的代码优化

目前，DFT系统仅仅具备局地优化能力。在语句级微分实现上采用二元归约的方法对微分代码进行优化是DFT系统的一个重要特色。根据右端表达式的乘法复杂性及含变元数目的不同，DFT系统采取不同的分解策略。二元归约的方法避免了微分计算中的许多冗余计算，在一些复杂的非线性表达式的微分计算中具有最小的计算代价，同时也非常适合于微分系统的软件实现。同时，对于某些特殊的运算操作（除法、乘方）和特殊函数（如sqrt、exp），DFT系统较好地利用了基态值计算得到的中间结果，避免了微分实现中的冗余计算。

4.系统应用

运用自动微分工具得到的切线性模式，可以在无截断误差意义下求解函数的数值微分和导数、稀疏雅可比矩阵。同时这些结果在数值参数敏感性分析、非线性最优化以及其它数值理论分析中有着非常重要的应用。这里简单介绍切线性模式的几个基本应用。

4.1符号导数和微分

如果输入为数学关系式，DFT系统可以自动生成对应的微分表达式和梯度，而与数学关系式的复杂程度无关。例如我们输入关系式：

，（1）

DFT系统将自动生成其符号微分形式及其梯度形式分别为

，（2）

4.2数值导数和微分

切线性模式最基本的应用就是在一定扰动输入下求解输出变量的扰动（响应）。表4.1给出了DFT系统在对IAP9L模式、GPSRayshooting模式和GPSRaytrace模式三个数值模式做切线性化的具体应用中，一些不同计算粒度、不同引用深度和不同程序风格的核心子过程，以及它们的切线性模式在SGI2000上运行的统计评价测试结果，其中切线性模式的可靠性指标都准确到六个有效数字以上，在运行时间、存储开销和代码复杂性方面分别是原模式的两倍左右，比较接近于理想的微分代价结果（1.5倍）。除了IAP9L模式由于过于复杂仅做粗略统计外，其余模式都用非注释语句行数来表示各自的代码复杂性。

4.3稀疏雅可比矩阵

运用上节讨论的方法来求解稀疏雅可比矩阵,具有极高的计算代价。例如，一个含个独立和个依赖参数的子过程，为求解整个雅可比矩阵就需要反复调用次切线性模式，当相当大时，这对许多实际的数值计算问题是不能接受的。事实上，如果雅可比矩阵的任意两列（行）相互正交，那么可以通过适当设置扰动输入值，这两列（行）的元素就可以通过一次引用切线性模式（伴随模式）完全得到。设和分别为雅可比矩阵的行宽度和列宽度，即各行和各列非零元素数目的最大值，显然有，。这里介绍几种常用的求解方法。

正向积分当时，通常采用切线性模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择右乘初始输入矩阵，可以获得接近的计算时间代价。DFT系统采用一种逐列（行）求解的方法，来有效求解右（左）乘初始输入矩阵。其基本思路是：按照某种列次序考察雅可比矩阵的各列；考察当前列中所有非零元素，并对这些非零元素所在行的行向量做类似模二和累加运算（即将非零元素视为逻辑“1”，零元素视为逻辑“0”），从而得到一个描述当前列与各行存在“某种”相关的标志向量（其元素都是“1”或“0”）；依据此标志向量，就很容易得到一个与之正交的列初始向量，其中与当前列序号对应的元素设置为“1”，而与标志向量中非零元素序号对应的元素设置为“0”，与标志向量中非零元素序号对应的元素设置为“-1”，显然，该列初始向量是唯一的，并且对应着当前右乘初始输入矩阵的最后一列；逐一考察已求解得到的列初始向量，如果某列初始向量与当前求解得到的列初始向量按下面定义的乘法（见过程4）正交，那么这两列就可以合并，即将当前列初始向量中非“-1”的元素按照对应关系分别赋值给该初始向量，并从记录中删除当前列初始向量；重复以上过程，继续按照给定列次序考察雅可比矩阵的“下一列”。不难说明，按照不同列次序求解得到的右乘初始输入矩阵可能不同。其中逐列求解右乘初始输入矩阵的过程可以简单叙述为：

1）将右乘初始输入矩阵所有元素的初值均设置为，，。。

2）如果，转6）。否则，如果雅可比矩阵的第列中的所有元素均为，，重复2）的判断。否则转3）。

3）计算标志向量。令，做如下计算：

，;

4）设为的列向量。在上定义乘法，对任意的，我们有：a）；b)如果，必有和。然后，做如下计算：

2);

5）令，并做如下计算：

，;

令,。如果，转6）；否则，重复2）的判断。

6）对，，如果，则。取的前列，这样，我们就得到了一个维右乘初始输入矩阵。

这里需要说明的是，运用上面的方法求得的右乘初始输入矩阵不仅与求解雅可比矩阵的列序有关，而且与过程4）中的合并顺序也有关系。至于如何最优求解右乘初始输入矩阵，目前还很难讨论清楚。但是，大量模拟试验结果表明，运用上面自然次序求得的右乘初始输入矩阵宽度已经非常接近于其下界值。

反向积分当和时，通常采用伴随模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择左乘初始输入矩阵，可以获得接近的计算时间代价。其中左乘初始输入矩阵的求解过程完全可以按照上面的方法进行，但是在处理前必须先将雅可比矩阵转置，最后还需将得到的初始输入矩阵转置才能最终得到左乘初始输入矩阵。同时，其行宽度也已经非常接近于其下界值。

混合积分如果将切线性模式和伴随模式相结合，往往可以避免梯度向量运算中的诸多冗余计算。

例如，ADJIFOR系统在求解雅可比矩阵时，在语句级微分实现中首先用伴随方法求得所有偏导数，然后做梯度向量积分；其计算时间代价与和模式的语句数目有关，而其存储代价为。具体讨论可参考文献[7]。

5．结论

切线性模式在无截断误差意义上计算函数的方向导数、梯度或雅可比矩阵，以及在模式的可预测性及参数敏感性分析、伴随模式构造等相关问题中有着广泛应用。DFT系统主要用于求解FORTRAN77语言编写的切线性模式，具有很强的全局数据相关分析能力。此外，DFT系统还具有其它几个重要特色，如结构化的微分实现、自动生成微分测试程序以及基于语句级的微分代码优化。本文简单给出了DFT系统在求解数值和符号导数和微分、稀疏雅可比矩阵中的应用。为评价一类自动微分系统，本文初步提出了统计准确率的概念。

参考文献

[1]AndreasGriewank.OnAutomaticDifferentiation.InM.IriandK.Tanabe,editors,MathematicalProgramming:

RecentDevelopmentsandApplications.KluwerAcademicPublishers,1989

[2]LeDimet,F.XandO.Talagrand,Variationalalgorithmsforanalysisandassimilationofmeteorological

observations:theoreticalaspects,Tellus,1986,38A,97-110

[3]P.Werbos,Applicationsofadvancesinnonlinearsensitivityanalysis,InsystemsModeling

andOptimization,NewYork,1982,SpringerVerlag,762-777

[4]ChristianBischof,GordonPusch,andRalfKnoesel."SensitivityAnalysisoftheMM5WeatherModelusing

AutomaticDifferentiation,"ComputersinPhysics,0:605-612,1996

[5]MuMu,etal，Thepredictabilityproblemofweatherandclimateprediction,ProgressinNatureScience,accepted.

[6]GieringR.etal.RecipesforAdjointCodeConstruction.ACMTrans.OnMath.Software.1998，24（4）：

437-474.

[7]C.Bischof,A.Carle,P.Khademi,andG.Pusch."AutomaticDifferentiation:ObtainingFastandReliable

Derivatives--Fast"inControlProblemsinIndustry,editedbyI.LasieckaandB.Morton,pages1-16,Birkhauser,

自动微分范文篇5

关键词自动微分切线性模式数据相关分析统计准确率

1.引言

计算微分大致经历了从商微分，符号微分，手写代码到自动微分几个阶段。与其它几种微分方法相比，自动微分具有代码简练、计算精度高及投入人力少等优点。自动微分实现的基本出发点是：一个数据相对独立的程序对象（模式、过程、程序段、数值语句乃至数值表达式），无论多么复杂，总可以分解为一系列有限数目的基本函数（如sin、exp、log）和基本运算操作（加、减、乘、除、乘方）的有序复合；对所有这些基本函数及基本运算操作，重复使用链式求导法则，将得到的中间结果自上而下地做正向积分就可以建立起对应的切线性模式，而自下而上地做反向积分就可以建立起对应的伴随模式[1]。基于自动微分方法得到的切线性模式和伴随模式，在变分资料同化[2]、系统建模与参数辨识[3]、参数的敏感性分析[4]、非线性最优化以及数值模式的可预测性分析[5]等问题中有着十分广泛的应用。

迄今为止，已有数十所大学和研究所各自开发了能够用于求解切线性模式的自动微分系统，比较典型的有TAMC系统[6]、ADJIFOR系统[7]和ODYSSEE系统[8]。在一些特定的运用中，它们都是比较成功的，但在通用性和复杂问题的处理效率上还存在许多不足。通常，自动生成切线性模式的关键难题在于对象自身的强相关性，这给系统全局分析（如数据IO相关分析和数据依赖相关分析）和微分代码的整体优化都带来了很多困难。同时，对于程序对象不可导处的准确识别和微分处理，至今仍还没有一个统一而有效的算法。另外，最优或有效求解稀疏雅可比矩阵一直是衡量一个自动微分系统有效性的重要尺度。

统计准确率被我们视为评价一类自动微分工具及其微分模式代码可靠性与有效性的重要尺度。其基本假设是：如果对于定义域空间内随机抽样获得的至多有限个n维初始场（或网格点），微分模式输出的差分和微分逼近是成功的；那么对于定义域空间内所有可能初始场（或网格点），微分模式输出的差分和微分逼近都是成功的。微分模式统计准确率评价的具体方法是：在所有随机抽样得到的初始场（或网格点）附近，当输入扰动逐渐趋向于机器有效精度所能表示的最小正值时，模式输出的差分和微分之间应该有足够精度有效位数上的逼近。

DFT系统具有许多优点，它能够完全接受用FORTRAN77语言编写的源代码，微分代码结构清晰，其微分处理能力与问题和对象的规模及复杂性无关。它基于YACC实现，具有很强的可扩展性。DFT系统具有四个重要特色。它通过对象全局依赖相关分析，准确求解雅可比矩阵的稀疏结构，自动计算有效初始输入矩阵，从而可以用较小的代价求得整个雅可比矩阵。同时，它可以自动生成客观评价微分模式效率与可靠性的测试程序，对奇异函数做等价微分处理，并采用二元归约的方法，在语句级层次上实现微分代码优化。

2.系统概况

DFT系统主要由两部分组成：微分代码转换和微分代码评价，图2.1。微分代码转换部分接受用户输入指令并自动分析对象模式，生成切线性模式代码及其相关测试代码，后者直接构成微分代码评价系统的主体。微分代码评价是DFT系统的一个重要特色。DFT系统的开发小组认为，一个微分模式如果在可靠性、时间和存储效率上没有得到充分的验证，至少对实际应用而言，它将是毫无意义的。

原模式切线性模式

统计评价结果

2.1微分代码转换

DFT系统是基于YACC在UNIX环境下开发的，其结构图2.2所示。通过DFT系统产生的切线性模式代码成对出现，并在语句级程度上做了简化，可读性很强，如图2.4。

切线性模式

评价函数集

图2.2微分代码转换

微分代码转换部分从功能上分为四个部分：词法分析，语义分析，对象复杂性及数据相关分析和微分代码转换。对于一组具有复杂数据相关的程序模式对象，通常需要系统运行两遍才能得到有效而可靠的微分代码。这主要有两方面的考虑：其一，根据对象的复杂性（如最大语句长度、最大变量维数、子过程或函数数目、子过程或函数内最大变量数目等对象特征）选择合适的系统参数以求最优的运行代价；其二，模式内各子过程或函数之间以及一个子过程或函数内往往具有很强的数据相关性，需要事先保存对象的相关信息并且在考虑当前对象的属性之前必须做上下文相关分析。

2.2微分代码评价

通常，评价一个编译系统的性能有很多方面，如处理速度、结果代码可靠性及质量、出错诊断、可扩展和可维护性等。对于一类自动微分系统来说，由于软件开发人力的局限以及对象模式的复杂多样性,通过自动转换得到的微分模式并非常常是有效而可靠的（即无论是在数学意义上还是在程序逻辑上应与期待的理想结果一致），因而在微分模式被投入实际应用前，往往需要投入一定的人力来对其做严格的分析测试。

对切线性模式做统计评价测试的主要内容可以简单叙述为：在网格化的模式定义域空间内，选择所有可能的网格点形成微分模式计算的初始场；在不同的网格点附近，随机选取至少个线性无关的初始扰动，对每个扰动输入分别进行网格点逼近，统计考察模式输出差分和微分在有效位数上的逼近程度。图2.5描述了整个测试过程，它包含网格点数据随机采样（1）和网格点数据逼近（2）两级循环。

3.系统主要特色

DFT系统并不是一个完整的FORTRAN编译器，但它几乎可以接受和处理所有FORTRAN77编写的源模式代码，并且可以很方便地扩展并接受FORTRAN90编写的源模式代码。本节将着重介绍DFT系统（版本3.0）的以下几个重要特色。

3.1结构化的微分实现

DFT系统采用标准化的代码实现，切线性模式的扰动变量和基态值变量、微分计算语句和基态值计算语句总是成对出现，并具有清晰的程序结构。微分代码保持了原模式本身的结构和风格（如并行和向量特性、数据精度等），即语句到语句、结构到结构的微分实现。在奇异点或不可导处，DFT系统对微分扰动采取简单的清零处理，实践证明这对抑制扰动计算溢出具有重要意义，但并不影响评价测试结果。

3.2全局数据相关分析

DFT系统具有较强的数据相关分析能力，它包括全局数据IO相关分析、全局数据依赖相关分析、全局过程相关分析以及数据迭代相关分析几个不同方面。数据依赖相关与数据IO相关关系密切，但又存在根本不同。前者强调每个变量在数学关系上的依赖性；而后者描述了一个对象的输入输出特性，且具有相对性，即任何一个变量参数，无论它是独立变量还是依赖变量，在数学意义上都可等价为一个既是输入又是输出的参数来处理。

DFT系统记录所有过程参数的IO属性表，通过深度递归相关计算，准确计算每个过程参数的最终IO属性。DFT系统通过对数据相关矩阵做模二和及自乘迭代计算（An+1=An⊕An2）来完成数据的依赖相关分析，这种算法具有很好的对数收敛特性。DFT系统通过全局过程相关分析的结果，自动生成模式的局部或整体相关引用树结构（如图3.1），这对用户分析复杂数值模式和微分评价测试都具有很好的指导作用。DFT系统还具有分析局部数据迭代相关和函数迭代相关的能力，这两种形式的数据迭代相关是自动微分实现颇具挑战的难题之一。

3.3自动生成测试程序

基于IO相关分析的结果，DFT系统自动生成微分测试代码，分别对切线性模式的可靠性和运行代价做统计评价测试。特别地，DFT系统还可将任何模式参数都视为输入输出参数，生成在数学意义上等价的测试代码，这样处理的不利之处在于往往需要极高的存储开销。

3.4基于语句级的代码优化

目前，DFT系统仅仅具备局地优化能力。在语句级微分实现上采用二元归约的方法对微分代码进行优化是DFT系统的一个重要特色。根据右端表达式的乘法复杂性及含变元数目的不同，DFT系统采取不同的分解策略。二元归约的方法避免了微分计算中的许多冗余计算，在一些复杂的非线性表达式的微分计算中具有最小的计算代价，同时也非常适合于微分系统的软件实现。同时，对于某些特殊的运算操作（除法、乘方）和特殊函数（如sqrt、exp），DFT系统较好地利用了基态值计算得到的中间结果，避免了微分实现中的冗余计算。

4.系统应用

运用自动微分工具得到的切线性模式，可以在无截断误差意义下求解函数的数值微分和导数、稀疏雅可比矩阵。同时这些结果在数值参数敏感性分析、非线性最优化以及其它数值理论分析中有着非常重要的应用。这里简单介绍切线性模式的几个基本应用。

4.1符号导数和微分

如果输入为数学关系式，DFT系统可以自动生成对应的微分表达式和梯度，而与数学关系式的复杂程度无关。例如我们输入关系式：

，（1）

DFT系统将自动生成其符号微分形式及其梯度形式分别为

，（2）

.2数值导数和微分

切线性模式最基本的应用就是在一定扰动输入下求解输出变量的扰动（响应）。表4.1给出了DFT系统在对IAP9L模式、GPSRayshooting模式和GPSRaytrace模式三个数值模式做切线性化的具体应用中，一些不同计算粒度、不同引用深度和不同程序风格的核心子过程，以及它们的切线性模式在SGI2000上运行的统计评价测试结果，其中切线性模式的可靠性指标都准确到六个有效数字以上，在运行时间、存储开销和代码复杂性方面分别是原模式的两倍左右，比较接近于理想的微分代价结果（1.5倍）。除了IAP9L模式由于过于复杂仅做粗略统计外，其余模式都用非注释语句行数来表示各自的代码复杂性。

表4.1DFT系统在三个数值模式中的统计评价测试结果

性能指标

对象模式运行时间（10-3秒）存储开销（字节数）代码复杂性

原模式切线性

模式

原模式切线性

模式

原模式切线性

模式

Xyz2g2.5306.1605524110485589

IntCIRA1.5602.750133426614165

Dabel0.0350.072601202749

LSS8.30017.50669133879143

RP42.4085.10360572102238

Vgrad10.1000.21218564368282454

RefGr43.0086.0071865414373083578

LL2JK0.6261.350262252442232

RayFind462.70

×103125.4

×103985618212111179

EPSIMP1.76011.50445589101327

Hlimits0.8301.8802425774842543774

Int3sL26.9051.2082002916394584690

MAKE

NCEP1340392072292514458504584

Curvcent0.0130.038527542754

DYFRM3.800

×1037.250

×1035000*9500*161279

PHYSIC2.750

×1035.385

×10330005000*1399*

(含注释行)2826*

(含注释行)

适当设置输入扰动的初值，运用切线性模式可以简单求解输出变量对输入的偏导数。例如，对于一个含有个输入参数的实型函数

（3）

这里设，。运用DFT系统，可以得到对应的切线性模式

（4）

其中，为切线性模式的扰动输入参数。可以通过以下办法来求得偏导数：

（5）

其中。如果对于某个既是输入参数又是输出参数，可以类似以下过程引用的办法来处理。对于过程引用的情形，例如一个含有个输入参数的子过程

（6）

其中，为输入参数；，为输出参数；，既为输入参数又为输出参数。运用DFT系统，可以得到对应的切线性模式为

（7）

其中，，，分别为切线性模式的微分扰动输入、输出和输入输出参数。可以通过以下输入扰动设置并引用切线性模式（7）来求得偏导数：

a)设置；（，）；（）可以同时求得（）和（），其中。

b)设置（）；；（，）可以同时求得（）和（），其中。

4.3稀疏雅可比矩阵

运用上节讨论的方法来求解稀疏雅可比矩阵,具有极高的计算代价。例如，一个含个独立和个依赖参数的子过程，为求解整个雅可比矩阵就需要反复调用次切线性模式，当相当大时，这对许多实际的数值计算问题是不能接受的。事实上，如果雅可比矩阵的任意两列（行）相互正交，那么可以通过适当设置扰动输入值，这两列（行）的元素就可以通过一次引用切线性模式（伴随模式）完全得到。设和分别为雅可比矩阵的行宽度和列宽度，即各行和各列非零元素数目的最大值，显然有，。这里介绍几种常用的求解方法。

正向积分当时，通常采用切线性模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择右乘初始输入矩阵，可以获得接近的计算时间代价。DFT系统采用一种逐列（行）求解的方法，来有效求解右（左）乘初始输入矩阵。其基本思路是：按照某种列次序考察雅可比矩阵的各列；考察当前列中所有非零元素，并对这些非零元素所在行的行向量做类似模二和累加运算（即将非零元素视为逻辑“1”，零元素视为逻辑“0”），从而得到一个描述当前列与各行存在“某种”相关的标志向量（其元素都是“1”或“0”）；依据此标志向量，就很容易得到一个与之正交的列初始向量，其中与当前列序号对应的元素设置为“1”，而与标志向量中非零元素序号对应的元素设置为“0”，与标志向量中非零元素序号对应的元素设置为“-1”，显然，该列初始向量是唯一的，并且对应着当前右乘初始输入矩阵的最后一列；逐一考察已求解得到的列初始向量，如果某列初始向量与当前求解得到的列初始向量按下面定义的乘法（见过程4）正交，那么这两列就可以合并，即将当前列初始向量中非“-1”的元素按照对应关系分别赋值给该初始向量，并从记录中删除当前列初始向量；重复以上过程，继续按照给定列次序考察雅可比矩阵的“下一列”。不难说明，按照不同列次序求解得到的右乘初始输入矩阵可能不同。其中逐列求解右乘初始输入矩阵的过程可以简单叙述为：

1）将右乘初始输入矩阵所有元素的初值均设置为，，。。

2）如果，转6）。否则，如果雅可比矩阵的第列中的所有元素均为，，重复2）的判断。否则转3）。

3）计算标志向量。令，做如下计算：

，;

4）设为的列向量。在上定义乘法，对任意的，我们有：a）；b)如果，必有和。然后，做如下计算：

,;

，6);

2);

5）令，并做如下计算：

，;

令,。如果，转6）；否则，重复2）的判断。

6）对，，如果，则。取的前列，这样，我们就得到了一个维右乘初始输入矩阵。

这里需要说明的是，运用上面的方法求得的右乘初始输入矩阵不仅与求解雅可比矩阵的列序有关，而且与过程4）中的合并顺序也有关系。至于如何最优求解右乘初始输入矩阵，目前还很难讨论清楚。但是，大量模拟试验结果表明，运用上面自然次序求得的右乘初始输入矩阵宽度已经非常接近于其下界值。

反向积分当和时，通常采用伴随模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择左乘初始输入矩阵，可以获得接近的计算时间代价。其中左乘初始输入矩阵的求解过程完全可以按照上面的方法进行，但是在处理前必须先将雅可比矩阵转置，最后还需将得到的初始输入矩阵转置才能最终得到左乘初始输入矩阵。同时，其行宽度也已经非常接近于其下界值。

混合积分如果将切线性模式和伴随模式相结合，往往可以避免梯度向量运算中的诸多冗余计算。例如，ADJIFOR系统在求解雅可比矩阵时，在语句级微分实现中首先用伴随方法求得所有偏导数，然后做梯度向量积分；其计算时间代价与和模式的语句数目有关，而其存储代价为。具体讨论可参考文献[7]。

5．结论

切线性模式在无截断误差意义上计算函数的方向导数、梯度或雅可比矩阵，以及在模式的可预测性及参数敏感性分析、伴随模式构造等相关问题中有着广泛应用。DFT系统主要用于求解FORTRAN77语言编写的切线性模式，具有很强的全局数据相关分析能力。此外，DFT系统还具有其它几个重要特色，如结构化的微分实现、自动生成微分测试程序以及基于语句级的微分代码优化。本文简单给出了DFT系统在求解数值和符号导数和微分、稀疏雅可比矩阵中的应用。为评价一类自动微分系统，本文初步提出了统计准确率的概念。

参考文献

[1]AndreasGriewank.OnAutomaticDifferentiation.InM.IriandK.Tanabe,editors,MathematicalProgramming:

RecentDevelopmentsandApplications.KluwerAcademicPublishers,1989

[2]LeDimet,F.XandO.Talagrand,Variationalalgorithmsforanalysisandassimilationofmeteorological

observations:theoreticalaspects,Tellus,1986,38A,97-110

[3]P.Werbos,Applicationsofadvancesinnonlinearsensitivityanalysis,InsystemsModeling

andOptimization,NewYork,1982,SpringerVerlag,762-777

[4]ChristianBischof,GordonPusch,andRalfKnoesel."SensitivityAnalysisoftheMM5WeatherModelusing

AutomaticDifferentiation,"ComputersinPhysics,0:605-612,1996

[5]MuMu,etal，Thepredictabilityproblemofweatherandclimateprediction,ProgressinNatureScience,accepted.

[6]GieringR.etal.RecipesforAdjointCodeConstruction.ACMTrans.OnMath.Software.1998，24（4）：

437-474.

[7]C.Bischof,A.Carle,P.Khademi,andG.Pusch."AutomaticDifferentiation:ObtainingFastandReliable

Derivatives--Fast"inControlProblemsinIndustry,editedbyI.LasieckaandB.Morton,pages1-16,Birkhauser,

自动微分范文篇6

关键词自动微分切线性模式数据相关分析统计准确率

1.引言

计算微分大致经历了从商微分，符号微分，手写代码到自动微分几个阶段。与其它几种微分方法相比，自动微分具有代码简练、计算精度高及投入人力少等优点。自动微分实现的基本出发点是：一个数据相对独立的程序对象（模式、过程、程序段、数值语句乃至数值表达式），无论多么复杂，总可以分解为一系列有限数目的基本函数（如sin、exp、log）和基本运算操作（加、减、乘、除、乘方）的有序复合；对所有这些基本函数及基本运算操作，重复使用链式求导法则，将得到的中间结果自上而下地做正向积分就可以建立起对应的切线性模式，而自下而上地做反向积分就可以建立起对应的伴随模式[1]。基于自动微分方法得到的切线性模式和伴随模式，在变分资料同化[2]、系统建模与参数辨识[3]、参数的敏感性分析[4]、非线性最优化以及数值模式的可预测性分析[5]等问题中有着十分广泛的应用。

迄今为止，已有数十所大学和研究所各自开发了能够用于求解切线性模式的自动微分系统，比较典型的有TAMC系统[6]、ADJIFOR系统[7]和ODYSSEE系统[8]。在一些特定的运用中，它们都是比较成功的，但在通用性和复杂问题的处理效率上还存在许多不足。通常，自动生成切线性模式的关键难题在于对象自身的强相关性，这给系统全局分析（如数据IO相关分析和数据依赖相关分析）和微分代码的整体优化都带来了很多困难。同时，对于程序对象不可导处的准确识别和微分处理，至今仍还没有一个统一而有效的算法。另外，最优或有效求解稀疏雅可比矩阵一直是衡量一个自动微分系统有效性的重要尺度。

统计准确率被我们视为评价一类自动微分工具及其微分模式代码可靠性与有效性的重要尺度。其基本假设是：如果对于定义域空间内随机抽样获得的至多有限个n维初始场（或网格点），微分模式输出的差分和微分逼近是成功的；那么对于定义域空间内所有可能初始场（或网格点），微分模式输出的差分和微分逼近都是成功的。微分模式统计准确率评价的具体方法是：在所有随机抽样得到的初始场（或网格点）附近，当输入扰动逐渐趋向于机器有效精度所能表示的最小正值时，模式输出的差分和微分之间应该有足够精度有效位数上的逼近。

DFT系统具有许多优点，它能够完全接受用FORTRAN77语言编写的源代码，微分代码结构清晰，其微分处理能力与问题和对象的规模及复杂性无关。它基于YACC实现，具有很强的可扩展性。DFT系统具有四个重要特色。它通过对象全局依赖相关分析，准确求解雅可比矩阵的稀疏结构，自动计算有效初始输入矩阵，从而可以用较小的代价求得整个雅可比矩阵。同时，它可以自动生成客观评价微分模式效率与可靠性的测试程序，对奇异函数做等价微分处理，并采用二元归约的方法，在语句级层次上实现微分代码优化。

2.系统概况

DFT系统主要由两部分组成：微分代码转换和微分代码评价，图2.1。微分代码转换部分接受用户输入指令并自动分析对象模式，生成切线性模式代码及其相关测试代码，后者直接构成微分代码评价系统的主体。微分代码评价是DFT系统的一个重要特色。DFT系统的开发小组认为，一个微分模式如果在可靠性、时间和存储效率上没有得到充分的验证，至少对实际应用而言，它将是毫无意义的。

微分代码转换部分从功能上分为四个部分：词法分析，语义分析，对象复杂性及数据相关分析和微分代码转换。对于一组具有复杂数据相关的程序模式对象，通常需要系统运行两遍才能得到有效而可靠的微分代码。这主要有两方面的考虑：其一，根据对象的复杂性（如最大语句长度、最大变量维数、子过程或函数数目、子过程或函数内最大变量数目等对象特征）选择合适的系统参数以求最优的运行代价；其二，模式内各子过程或函数之间以及一个子过程或函数内往往具有很强的数据相关性，需要事先保存对象的相关信息并且在考虑当前对象的属性之前必须做上下文相关分析。

2.2微分代码评价

通常，评价一个编译系统的性能有很多方面，如处理速度、结果代码可靠性及质量、出错诊断、可扩展和可维护性等。对于一类自动微分系统来说，由于软件开发人力的局限以及对象模式的复杂多样性,通过自动转换得到的微分模式并非常常是有效而可靠的（即无论是在数学意义上还是在程序逻辑上应与期待的理想结果一致），因而在微分模式被投入实际应用前，往往需要投入一定的人力来对其做严格的分析测试。

对切线性模式做统计评价测试的主要内容可以简单叙述为：在网格化的模式定义域空间内，选择所有可能的网格点形成微分模式计算的初始场；在不同的网格点附近，随机选取至少个线性无关的初始扰动，对每个扰动输入分别进行网格点逼近，统计考察模式输出差分和微分在有效位数上的逼近程度。图2.5描述了整个测试过程，它包含网格点数据随机采样（1）和网格点数据逼近（2）两级循环。

3.系统主要特色

DFT系统并不是一个完整的FORTRAN编译器，但它几乎可以接受和处理所有FORTRAN77编写的源模式代码，并且可以很方便地扩展并接受FORTRAN90编写的源模式代码。本节将着重介绍DFT系统（版本3.0）的以下几个重要特色。

3.1结构化的微分实现

DFT系统采用标准化的代码实现，切线性模式的扰动变量和基态值变量、微分计算语句和基态值计算语句总是成对出现，并具有清晰的程序结构。微分代码保持了原模式本身的结构和风格（如并行和向量特性、数据精度等），即语句到语句、结构到结构的微分实现。在奇异点或不可导处，DFT系统对微分扰动采取简单的清零处理，实践证明这对抑制扰动计算溢出具有重要意义，但并不影响评价测试结果。

3.2全局数据相关分析

DFT系统具有较强的数据相关分析能力，它包括全局数据IO相关分析、全局数据依赖相关分析、全局过程相关分析以及数据迭代相关分析几个不同方面。数据依赖相关与数据IO相关关系密切，但又存在根本不同。前者强调每个变量在数学关系上的依赖性；而后者描述了一个对象的输入输出特性，且具有相对性，即任何一个变量参数，无论它是独立变量还是依赖变量，在数学意义上都可等价为一个既是输入又是输出的参数来处理。

DFT系统记录所有过程参数的IO属性表，通过深度递归相关计算，准确计算每个过程参数的最终IO属性。DFT系统通过对数据相关矩阵做模二和及自乘迭代计算（An+1=An⊕An2）来完成数据的依赖相关分析，这种算法具有很好的对数收敛特性。DFT系统通过全局过程相关分析的结果，自动生成模式的局部或整体相关引用树结构（如图3.1），这对用户分析复杂数值模式和微分评价测试都具有很好的指导作用。DFT系统还具有分析局部数据迭代相关和函数迭代相关的能力，这两种形式的数据迭代相关是自动微分实现颇具挑战的难题之一。

3.3自动生成测试程序

基于IO相关分析的结果，DFT系统自动生成微分测试代码，分别对切线性模式的可靠性和运行代价做统计评价测试。特别地，DFT系统还可将任何模式参数都视为输入输出参数，生成在数学意义上等价的测试代码，这样处理的不利之处在于往往需要极高的存储开销。

3.4基于语句级的代码优化

目前，DFT系统仅仅具备局地优化能力。在语句级微分实现上采用二元归约的方法对微分代码进行优化是DFT系统的一个重要特色。根据右端表达式的乘法复杂性及含变元数目的不同，DFT系统采取不同的分解策略。二元归约的方法避免了微分计算中的许多冗余计算，在一些复杂的非线性表达式的微分计算中具有最小的计算代价，同时也非常适合于微分系统的软件实现。同时，对于某些特殊的运算操作（除法、乘方）和特殊函数（如sqrt、exp），DFT系统较好地利用了基态值计算得到的中间结果，避免了微分实现中的冗余计算。

4.系统应用

运用自动微分工具得到的切线性模式，可以在无截断误差意义下求解函数的数值微分和导数、稀疏雅可比矩阵。同时这些结果在数值参数敏感性分析、非线性最优化以及其它数值理论分析中有着非常重要的应用。这里简单介绍切线性模式的几个基本应用。

4.1符号导数和微分

如果输入为数学关系式，DFT系统可以自动生成对应的微分表达式和梯度，而与数学关系式的复杂程度无关。例如我们输入关系式：，（1）

DFT系统将自动生成其符号微分形式及其梯度形式分别为，（2）

4.2数值导数和微分

切线性模式最基本的应用就是在一定扰动输入下求解输出变量的扰动（响应）。表4.1给出了DFT系统在对IAP9L模式、GPSRayshooting模式和GPSRaytrace模式三个数值模式做切线性化的具体应用中，一些不同计算粒度、不同引用深度和不同程序风格的核心子过程，以及它们的切线性模式在SGI2000上运行的统计评价测试结果，其中切线性模式的可靠性指标都准确到六个有效数字以上，在运行时间、存储开销和代码复杂性方面分别是原模式的两倍左右，比较接近于理想的微分代价结果（1.5倍）。除了IAP9L模式由于过于复杂仅做粗略统计外，其余模式都用非注释语句行数来表示各自的代码复杂性。

适当设置输入扰动的初值，运用切线性模式可以简单求解输出变量对输入的偏导数。例如，对于一个含有个输入参数的实型函数（3）

这里设，。运用DFT系统，可以得到对应的切线性模式（4）

其中，为切线性模式的扰动输入参数。可以通过以下办法来求得偏导数：（5）

其中。如果对于某个既是输入参数又是输出参数，可以类似以下过程引用的办法来处理。对于过程引用的情形，例如一个含有个输入参数的子过程（6）

其中，为输入参数；，为输出参数；，既为输入参数又为输出参数。运用DFT系统，可以得到对应的切线性模式为（7）

其中，，，分别为切线性模式的微分扰动输入、输出和输入输出参数。可以通过以下输入扰动设置并引用切线性模式（7）来求得偏导数：a)设置；（，）；（）可以同时求得（）和（），其中。

b)设置（）；；（，）可以同时求得（）和（），其中。

4.3稀疏雅可比矩阵

运用上节讨论的方法来求解稀疏雅可比矩阵,具有极高的计算代价。例如，一个含个独立和个依赖参数的子过程，为求解整个雅可比矩阵就需要反复调用次切线性模式，当相当大时，这对许多实际的数值计算问题是不能接受的。事实上，如果雅可比矩阵的任意两列（行）相互正交，那么可以通过适当设置扰动输入值，这两列（行）的元素就可以通过一次引用切线性模式（伴随模式）完全得到。设和分别为雅可比矩阵的行宽度和列宽度，即各行和各列非零元素数目的最大值，显然有，。这里介绍几种常用的求解方法。

正向积分当时，通常采用切线性模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择右乘初始输入矩阵，可以获得接近的计算时间代价。DFT系统采用一种逐列（行）求解的方法，来有效求解右（左）乘初始输入矩阵。其基本思路是：按照某种列次序考察雅可比矩阵的各列；考察当前列中所有非零元素，并对这些非零元素所在行的行向量做类似模二和累加运算（即将非零元素视为逻辑“1”，零元素视为逻辑“0”），从而得到一个描述当前列与各行存在“某种”相关的标志向量（其元素都是“1”或“0”）；依据此标志向量，就很容易得到一个与之正交的列初始向量，其中与当前列序号对应的元素设置为“1”，而与标志向量中非零元素序号对应的元素设置为“0”，与标志向量中非零元素序号对应的元素设置为“-1”，显然，该列初始向量是唯一的，并且对应着当前右乘初始输入矩阵的最后一列；逐一考察已求解得到的列初始向量，如果某列初始向量与当前求解得到的列初始向量按下面定义的乘法（见过程4）正交，那么这两列就可以合并，即将当前列初始向量中非“-1”的元素按照对应关系分别赋值给该初始向量，并从记录中删除当前列初始向量；重复以上过程，继续按照给定列次序考察雅可比矩阵的“下一列”。不难说明，按照不同列次序求解得到的右乘初始输入矩阵可能不同。其中逐列求解右乘初始输入矩阵的过程可以简单叙述为：

1）将右乘初始输入矩阵所有元素的初值均设置为，，。。

2）如果，转6）。否则，如果雅可比矩阵的第列中的所有元素均为，，重复2）的判断。否则转3）。

3）计算标志向量。令，做如下计算：，;

4）设为的列向量。在上定义乘法，对任意的，我们有：a）；b)如果，必有和。然后，做如下计算：,;，6);2);

5）令，并做如下计算：，;令,。如果，转6）；否则，重复2）的判断。

6）对，，如果，则。取的前列，这样，我们就得到了一个维右乘初始输入矩阵。

这里需要说明的是，运用上面的方法求得的右乘初始输入矩阵不仅与求解雅可比矩阵的列序有关，而且与过程4）中的合并顺序也有关系。至于如何最优求解右乘初始输入矩阵，目前还很难讨论清楚。但是，大量模拟试验结果表明，运用上面自然次序求得的右乘初始输入矩阵宽度已经非常接近于其下界值。

反向积分当和时，通常采用伴随模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择左乘初始输入矩阵，可以获得接近的计算时间代价。其中左乘初始输入矩阵的求解过程完全可以按照上面的方法进行，但是在处理前必须先将雅可比矩阵转置，最后还需将得到的初始输入矩阵转置才能最终得到左乘初始输入矩阵。同时，其行宽度也已经非常接近于其下界值。

混合积分如果将切线性模式和伴随模式相结合，往往可以避免梯度向量运算中的诸多冗余计算。例如，ADJIFOR系统在求解雅可比矩阵时，在语句级微分实现中首先用伴随方法求得所有偏导数，然后做梯度向量积分；其计算时间代价与和模式的语句数目有关，而其存储代价为。具体讨论可参考文献[7]。

5．结论

切线性模式在无截断误差意义上计算函数的方向导数、梯度或雅可比矩阵，以及在模式的可预测性及参数敏感性分析、伴随模式构造等相关问题中有着广泛应用。DFT系统主要用于求解FORTRAN77语言编写的切线性模式，具有很强的全局数据相关分析能力。此外，DFT系统还具有其它几个重要特色，如结构化的微分实现、自动生成微分测试程序以及基于语句级的微分代码优化。本文简单给出了DFT系统在求解数值和符号导数和微分、稀疏雅可比矩阵中的应用。为评价一类自动微分系统，本文初步提出了统计准确率的概念。

参考文献

[1]AndreasGriewank.OnAutomaticDifferentiation.InM.IriandK.Tanabe,editors,MathematicalProgramming:

RecentDevelopmentsandApplications.KluwerAcademicPublishers,1989

[2]LeDimet,F.XandO.Talagrand,Variationalalgorithmsforanalysisandassimilationofmeteorological

observations:theoreticalaspects,Tellus,1986,38A,97-110

[3]P.Werbos,Applicationsofadvancesinnonlinearsensitivityanalysis,InsystemsModeling

andOptimization,NewYork,1982,SpringerVerlag,762-777

[4]ChristianBischof,GordonPusch,andRalfKnoesel."SensitivityAnalysisoftheMM5WeatherModelusing

AutomaticDifferentiation,"ComputersinPhysics,0:605-612,1996

[5]MuMu,etal，Thepredictabilityproblemofweatherandclimateprediction,ProgressinNatureScience,accepted.

[6]GieringR.etal.RecipesforAdjointCodeConstruction.ACMTrans.OnMath.Software.1998，24（4）：

437-474.

[7]C.Bischof,A.Carle,P.Khademi,andG.Pusch."AutomaticDifferentiation:ObtainingFastandReliable

Derivatives--Fast"inControlProblemsinIndustry,editedbyI.LasieckaandB.Morton,pages1-16,Birkhauser,

自动微分范文篇7

关键词自动微分切线性模式数据相关分析统计准确率

1.引言

计算微分大致经历了从商微分，符号微分，手写代码到自动微分几个阶段。与其它几种微分方法相比，自动微分具有代码简练、计算精度高及投入人力少等优点。自动微分实现的基本出发点是：一个数据相对独立的程序对象（模式、过程、程序段、数值语句乃至数值表达式），无论多么复杂，总可以分解为一系列有限数目的基本函数（如sin、exp、log）和基本运算操作（加、减、乘、除、乘方）的有序复合；对所有这些基本函数及基本运算操作，重复使用链式求导法则，将得到的中间结果自上而下地做正向积分就可以建立起对应的切线性模式，而自下而上地做反向积分就可以建立起对应的伴随模式[1]。基于自动微分方法得到的切线性模式和伴随模式，在变分资料同化[2]、系统建模与参数辨识[3]、参数的敏感性分析[4]、非线性最优化以及数值模式的可预测性分析[5]等问题中有着十分广泛的应用。

迄今为止，已有数十所大学和研究所各自开发了能够用于求解切线性模式的自动微分系统，比较典型的有TAMC系统[6]、ADJIFOR系统[7]和ODYSSEE系统[8]。在一些特定的运用中，它们都是比较成功的，但在通用性和复杂问题的处理效率上还存在许多不足。通常，自动生成切线性模式的关键难题在于对象自身的强相关性，这给系统全局分析（如数据IO相关分析和数据依赖相关分析）和微分代码的整体优化都带来了很多困难。同时，对于程序对象不可导处的准确识别和微分处理，至今仍还没有一个统一而有效的算法。另外，最优或有效求解稀疏雅可比矩阵一直是衡量一个自动微分系统有效性的重要尺度。

统计准确率被我们视为评价一类自动微分工具及其微分模式代码可靠性与有效性的重要尺度。其基本假设是：如果对于定义域空间内随机抽样获得的至多有限个n维初始场（或网格点），微分模式输出的差分和微分逼近是成功的；那么对于定义域空间内所有可能初始场（或网格点），微分模式输出的差分和微分逼近都是成功的。微分模式统计准确率评价的具体方法是：在所有随机抽样得到的初始场（或网格点）附近，当输入扰动逐渐趋向于机器有效精度所能表示的最小正值时，模式输出的差分和微分之间应该有足够精度有效位数上的逼近。

DFT系统具有许多优点，它能够完全接受用FORTRAN77语言编写的源代码，微分代码结构清晰，其微分处理能力与问题和对象的规模及复杂性无关。它基于YACC实现，具有很强的可扩展性。DFT系统具有四个重要特色。它通过对象全局依赖相关分析，准确求解雅可比矩阵的稀疏结构，自动计算有效初始输入矩阵，从而可以用较小的代价求得整个雅可比矩阵。同时，它可以自动生成客观评价微分模式效率与可靠性的测试程序，对奇异函数做等价微分处理，并采用二元归约的方法，在语句级层次上实现微分代码优化。

2.系统概况

DFT系统主要由两部分组成：微分代码转换和微分代码评价，图2.1。微分代码转换部分接受用户输入指令并自动分析对象模式，生成切线性模式代码及其相关测试代码，后者直接构成微分代码评价系统的主体。微分代码评价是DFT系统的一个重要特色。DFT系统的开发小组认为，一个微分模式如果在可靠性、时间和存储效率上没有得到充分的验证，至少对实际应用而言，它将是毫无意义的。

原模式切线性模式

统计评价结果

图2.1DFT系统结构简图

2.1微分代码转换

DFT系统是基于YACC在UNIX环境下开发的，其结构图2.2所示。通过DFT系统产生的切线性模式代码成对出现，并在语句级程度上做了简化，可读性很强，如图2.4。

切线性模式

评价函数集

图2.2微分代码转换

微分代码转换部分从功能上分为四个部分：词法分析，语义分析，对象复杂性及数据相关分析和微分代码转换。对于一组具有复杂数据相关的程序模式对象，通常需要系统运行两遍才能得到有效而可靠的微分代码。这主要有两方面的考虑：其一，根据对象的复杂性（如最大语句长度、最大变量维数、子过程或函数数目、子过程或函数内最大变量数目等对象特征）选择合适的系统参数以求最优的运行代价；其二，模式内各子过程或函数之间以及一个子过程或函数内往往具有很强的数据相关性，需要事先保存对象的相关信息并且在考虑当前对象的属性之前必须做上下文相关分析。

图2.3PERIGEE源程序代码图2.4DFT系统生成的切线性代码

2.2微分代码评价

通常，评价一个编译系统的性能有很多方面，如处理速度、结果代码可靠性及质量、出错诊断、可扩展和可维护性等。对于一类自动微分系统来说，由于软件开发人力的局限以及对象模式的复杂多样性,通过自动转换得到的微分模式并非常常是有效而可靠的（即无论是在数学意义上还是在程序逻辑上应与期待的理想结果一致），因而在微分模式被投入实际应用前，往往需要投入一定的人力来对其做严格的分析测试。

对切线性模式做统计评价测试的主要内容可以简单叙述为：在网格化的模式定义域空间内，选择所有可能的网格点形成微分模式计算的初始场；在不同的网格点附近，随机选取至少个线性无关的初始扰动，对每个扰动输入分别进行网格点逼近，统计考察模式输出差分和微分在有效位数上的逼近程度。图2.5描述了整个测试过程，它包含网格点数据随机采样（1）和网格点数据逼近（2）两级循环。

图2.5切线性模式代码的测试过程

3.系统主要特色

DFT系统并不是一个完整的FORTRAN编译器，但它几乎可以接受和处理所有FORTRAN77编写的源模式代码，并且可以很方便地扩展并接受FORTRAN90编写的源模式代码。本节将着重介绍DFT系统（版本3.0）的以下几个重要特色。

3.1结构化的微分实现

DFT系统采用标准化的代码实现，切线性模式的扰动变量和基态值变量、微分计算语句和基态值计算语句总是成对出现，并具有清晰的程序结构。微分代码保持了原模式本身的结构和风格（如并行和向量特性、数据精度等），即语句到语句、结构到结构的微分实现。在奇异点或不可导处，DFT系统对微分扰动采取简单的清零处理，实践证明这对抑制扰动计算溢出具有重要意义，但并不影响评价测试结果。

3.2全局数据相关分析

DFT系统具有较强的数据相关分析能力，它包括全局数据IO相关分析、全局数据依赖相关分析、全局过程相关分析以及数据迭代相关分析几个不同方面。数据依赖相关与数据IO相关关系密切，但又存在根本不同。前者强调每个变量在数学关系上的依赖性；而后者描述了一个对象的输入输出特性，且具有相对性，即任何一个变量参数，无论它是独立变量还是依赖变量，在数学意义上都可等价为一个既是输入又是输出的参数来处理。

DFT系统记录所有过程参数的IO属性表，通过深度递归相关计算，准确计算每个过程参数的最终IO属性。DFT系统通过对数据相关矩阵做模二和及自乘迭代计算（An+1=An⊕An2）来完成数据的依赖相关分析，这种算法具有很好的对数收敛特性。DFT系统通过全局过程相关分析的结果，自动生成模式的局部或整体相关引用树结构（如图3.1），这对用户分析复杂数值模式和微分评价测试都具有很好的指导作用。DFT系统还具有分析局部数据迭代相关和函数迭代相关的能力，这两种形式的数据迭代相关是自动微分实现颇具挑战的难题之一。

图3.1GPSRayshooting模式的相关树结构片段

3.3自动生成测试程序

基于IO相关分析的结果，DFT系统自动生成微分测试代码，分别对切线性模式的可靠性和运行代价做统计评价测试。特别地，DFT系统还可将任何模式参数都视为输入输出参数，生成在数学意义上等价的测试代码，这样处理的不利之处在于往往需要极高的存储开销。

3.4基于语句级的代码优化

目前，DFT系统仅仅具备局地优化能力。在语句级微分实现上采用二元归约的方法对微分代码进行优化是DFT系统的一个重要特色。根据右端表达式的乘法复杂性及含变元数目的不同，DFT系统采取不同的分解策略。二元归约的方法避免了微分计算中的许多冗余计算，在一些复杂的非线性表达式的微分计算中具有最小的计算代价，同时也非常适合于微分系统的软件实现。同时，对于某些特殊的运算操作（除法、乘方）和特殊函数（如sqrt、exp），DFT系统较好地利用了基态值计算得到的中间结果，避免了微分实现中的冗余计算。

4.系统应用

运用自动微分工具得到的切线性模式，可以在无截断误差意义下求解函数的数值微分和导数、稀疏雅可比矩阵。同时这些结果在数值参数敏感性分析、非线性最优化以及其它数值理论分析中有着非常重要的应用。这里简单介绍切线性模式的几个基本应用。

4.1符号导数和微分

如果输入为数学关系式，DFT系统可以自动生成对应的微分表达式和梯度，而与数学关系式的复杂程度无关。例如我们输入关系式：

，（1）

DFT系统将自动生成其符号微分形式及其梯度形式分别为

，（2）

4.2数值导数和微分

切线性模式最基本的应用就是在一定扰动输入下求解输出变量的扰动（响应）。表4.1给出了DFT系统在对IAP9L模式、GPSRayshooting模式和GPSRaytrace模式三个数值模式做切线性化的具体应用中，一些不同计算粒度、不同引用深度和不同程序风格的核心子过程，以及它们的切线性模式在SGI2000上运行的统计评价测试结果，其中切线性模式的可靠性指标都准确到六个有效数字以上，在运行时间、存储开销和代码复杂性方面分别是原模式的两倍左右，比较接近于理想的微分代价结果（1.5倍）。除了IAP9L模式由于过于复杂仅做粗略统计外，其余模式都用非注释语句行数来表示各自的代码复杂性。

表4.1DFT系统在三个数值模式中的统计评价测试结果

性能指标

对象模式运行时间（10-3秒）存储开销（字节数）代码复杂性

原模式切线性

模式

原模式切线性

模式

原模式切线性

模式

Xyz2g2.5306.1605524110485589

IntCIRA1.5602.750133426614165

Dabel0.0350.072601202749

LSS8.30017.50669133879143

RP42.4085.10360572102238

Vgrad10.1000.21218564368282454

RefGr43.0086.0071865414373083578

LL2JK0.6261.350262252442232

RayFind462.70

×103125.4

×103985618212111179

EPSIMP1.76011.50445589101327

Hlimits0.8301.8802425774842543774

Int3sL26.9051.2082002916394584690

MAKE

NCEP1340392072292514458504584

Curvcent0.0130.038527542754

DYFRM3.800

×1037.250

×1035000*9500*161279

PHYSIC2.750

×1035.385

×10330005000*1399*

(含注释行)2826*

(含注释行)

适当设置输入扰动的初值，运用切线性模式可以简单求解输出变量对输入的偏导数。例如，对于一个含有个输入参数的实型函数

（3）

这里设，。运用DFT系统，可以得到对应的切线性模式

（4）

其中，为切线性模式的扰动输入参数。可以通过以下办法来求得偏导数：

（5）

其中。如果对于某个既是输入参数又是输出参数，可以类似以下过程引用的办法来处理。对于过程引用的情形，例如一个含有个输入参数的子过程

（6）

其中，为输入参数；，为输出参数；，既为输入参数又为输出参数。运用DFT系统，可以得到对应的切线性模式为

（7）

其中，，，分别为切线性模式的微分扰动输入、输出和输入输出参数。可以通过以下输入扰动设置并引用切线性模式（7）来求得偏导数：

a)设置；（，）；（）可以同时求得（）和（），其中。

b)设置（）；；（，）可以同时求得（）和（），其中。

4.3稀疏雅可比矩阵

运用上节讨论的方法来求解稀疏雅可比矩阵,具有极高的计算代价。例如，一个含个独立和个依赖参数的子过程，为求解整个雅可比矩阵就需要反复调用次切线性模式，当相当大时，这对许多实际的数值计算问题是不能接受的。事实上，如果雅可比矩阵的任意两列（行）相互正交，那么可以通过适当设置扰动输入值，这两列（行）的元素就可以通过一次引用切线性模式（伴随模式）完全得到。设和分别为雅可比矩阵的行宽度和列宽度，即各行和各列非零元素数目的最大值，显然有，。这里介绍几种常用的求解方法。

正向积分当时，通常采用切线性模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择右乘初始输入矩阵，可以获得接近的计算时间代价。DFT系统采用一种逐列（行）求解的方法，来有效求解右（左）乘初始输入矩阵。其基本思路是：按照某种列次序考察雅可比矩阵的各列；考察当前列中所有非零元素，并对这些非零元素所在行的行向量做类似模二和累加运算（即将非零元素视为逻辑“1”，零元素视为逻辑“0”），从而得到一个描述当前列与各行存在“某种”相关的标志向量（其元素都是“1”或“0”）；依据此标志向量，就很容易得到一个与之正交的列初始向量，其中与当前列序号对应的元素设置为“1”，而与标志向量中非零元素序号对应的元素设置为“0”，与标志向量中非零元素序号对应的元素设置为“-1”，显然，该列初始向量是唯一的，并且对应着当前右乘初始输入矩阵的最后一列；逐一考察已求解得到的列初始向量，如果某列初始向量与当前求解得到的列初始向量按下面定义的乘法（见过程4）正交，那么这两列就可以合并，即将当前列初始向量中非“-1”的元素按照对应关系分别赋值给该初始向量，并从记录中删除当前列初始向量；重复以上过程，继续按照给定列次序考察雅可比矩阵的“下一列”。不难说明，按照不同列次序求解得到的右乘初始输入矩阵可能不同。其中逐列求解右乘初始输入矩阵的过程可以简单叙述为：

1）将右乘初始输入矩阵所有元素的初值均设置为，，。。

2）如果，转6）。否则，如果雅可比矩阵的第列中的所有元素均为，，重复2）的判断。否则转3）。

3）计算标志向量。令，做如下计算：

，;

4）设为的列向量。在上定义乘法，对任意的，我们有：a）；b)如果，必有和。然后，做如下计算：

,;

，6);

2);

5）令，并做如下计算：

，;

令,。如果，转6）；否则，重复2）的判断。

6）对，，如果，则。取的前列，这样，我们就得到了一个维右乘初始输入矩阵。

这里需要说明的是，运用上面的方法求得的右乘初始输入矩阵不仅与求解雅可比矩阵的列序有关，而且与过程4）中的合并顺序也有关系。至于如何最优求解右乘初始输入矩阵，目前还很难讨论清楚。但是，大量模拟试验结果表明，运用上面自然次序求得的右乘初始输入矩阵宽度已经非常接近于其下界值。

反向积分当和时，通常采用伴随模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择左乘初始输入矩阵，可以获得接近的计算时间代价。其中左乘初始输入矩阵的求解过程完全可以按照上面的方法进行，但是在处理前必须先将雅可比矩阵转置，最后还需将得到的初始输入矩阵转置才能最终得到左乘初始输入矩阵。同时，其行宽度也已经非常接近于其下界值。

混合积分如果将切线性模式和伴随模式相结合，往往可以避免梯度向量运算中的诸多冗余计算。例如，ADJIFOR系统在求解雅可比矩阵时，在语句级微分实现中首先用伴随方法求得所有偏导数，然后做梯度向量积分；其计算时间代价与和模式的语句数目有关，而其存储代价为。具体讨论可参考文献[7]。

5．结论

切线性模式在无截断误差意义上计算函数的方向导数、梯度或雅可比矩阵，以及在模式的可预测性及参数敏感性分析、伴随模式构造等相关问题中有着广泛应用。DFT系统主要用于求解FORTRAN77语言编写的切线性模式，具有很强的全局数据相关分析能力。此外，DFT系统还具有其它几个重要特色，如结构化的微分实现、自动生成微分测试程序以及基于语句级的微分代码优化。本文简单给出了DFT系统在求解数值和符号导数和微分、稀疏雅可比矩阵中的应用。为评价一类自动微分系统，本文初步提出了统计准确率的概念。

参考文献

[1]AndreasGriewank.OnAutomaticDifferentiation.InM.IriandK.Tanabe,editors,MathematicalProgramming:

RecentDevelopmentsandApplications.KluwerAcademicPublishers,1989

[2]LeDimet,F.XandO.Talagrand,Variationalalgorithmsforanalysisandassimilationofmeteorological

observations:theoreticalaspects,Tellus,1986,38A,97-110

[3]P.Werbos,Applicationsofadvancesinnonlinearsensitivityanalysis,InsystemsModeling

andOptimization,NewYork,1982,SpringerVerlag,762-777

[4]ChristianBischof,GordonPusch,andRalfKnoesel."SensitivityAnalysisoftheMM5WeatherModelusing

AutomaticDifferentiation,"ComputersinPhysics,0:605-612,1996

[5]MuMu,etal，Thepredictabilityproblemofweatherandclimateprediction,ProgressinNatureScience,accepted.

[6]GieringR.etal.RecipesforAdjointCodeConstruction.ACMTrans.OnMath.Software.1998，24（4）：

437-474.

[7]C.Bischof,A.Carle,P.Khademi,andG.Pusch."AutomaticDifferentiation:ObtainingFastandReliable

Derivatives--Fast"inControlProblemsinIndustry,editedbyI.LasieckaandB.Morton,pages1-16,Birkhauser,

自动微分范文篇8

关键词：发电厂；循环流化床；除渣；风室；自动控制

郭家湾电厂建设规模为2×300MW发电机组，锅炉为国产首台1065t/h亚临界双床循环流化床锅炉，每台锅炉配备6台冷渣器、3台链斗输送机、3台斗提机用于除去煤燃烧后产生的灰渣，由于双床循环流化床锅炉结构的特殊性，其除渣自动控制没有可参照的模型，投产后除渣系统一直由运行人员手动控制，根据该情况，电厂技术人员设计了除渣自动控制系统模型并进行逻辑修改和调试。

1双床循环流化床锅炉除渣系统

循环流化床锅炉是1种通过燃料燃烧将水加热后产生蒸汽的设备，双床循环流化床锅炉具有2个布风板（床），此类锅炉通常用于发电厂。双床循环流化床锅炉配备除渣系统，用于将炉膛内燃料燃烧产生的灰渣放出冷却后输送至渣仓。炉膛内燃料燃烧产生的灰渣堆积在床上，床上的灰渣通过放渣管流出到冷渣器，在冷渣器内冷却后排出，经链斗机和斗提机输送至渣仓，工艺系统如图1所示。双床循环流化床锅炉除渣系统配备6台变频式冷渣器。运行过程中，通过调节变频冷渣器的转速增加或减少放渣量来调整水冷风室压力在允许范围内，同时还要避免冷渣器出渣温度和冷渣器出水温度超限。双床循环流化床锅炉除渣系统运行时，依靠运行人员手动调节冷渣器转速来控制水冷风室压力，手动调节存在调节不准确、劳动强度大、费电、安全稳定性差等缺点。

2双床循环流化床锅炉除渣自动控制系统

根据双床循环流化床锅炉的特点，设计了一种双床循环流化床锅炉除渣自动控制系统模型，在使用过程中取得了较好的效果。以左侧冷渣系统为例，对模型进行说明（右侧冷渣系统控制模型与左侧类同），如图2所示。（1）I型选择器功能：当S1引脚为1，S2、S3引脚为0时，输出等于P1；当S2引脚为1，S1、S3引脚为0时，输出等于P2；当S3引脚为1，S1、S2引脚为0时，输出等于P3；当S1、S2、S3引脚均为0时，输出等于P0。（2）II型选择器功能：当Z引脚为1时，输出等于X1；当Z引脚为0时，输出等于X2。（3）高限块功能：当第一个引脚的输入值大于第二个引脚的设定值时，该功能块输出为1。（4）操作器：操作器可在DCS或PLC人机界面中进行操作，可以设置为自动或手动模式，当设置为自动模式时，操作器的输出为前面功能块的输出；当设置为手动模式时，操作器的输出为操作人员手动置入的数据。（5）逻辑功能：调节器的输入值为“左侧水冷风室压力目标值”与“左侧水冷风室压力实测值”的偏差，采用变比例系数、变积分时间、固定微分时间、固定微分系数的调节方式，可变比例系数和可变积分时间通过选择器获得。总操作器的输入信号为调节器的输出。当“1号冷渣器出渣温度”超过高限或“1号冷渣器出水温度”超过高限时，1号冷渣器操作器的输入为超温前总操作器输出值的1/2，否则1号冷渣器操作器的输入为总操作器的输出。当“3号冷渣器出渣温度”超过高限或“3号冷渣器出水温度”超过高限时，3号冷渣器操作器的输入为超温前总操作器输出值的1/2，否则3号冷渣器操作器的输入为总操作器的输出。当“5号冷渣器出渣温度”超过高限或“5号冷渣器出水温度”超过高限时，5号冷渣器操作器的输入为超温前总操作器输出值的1/2，否则5号冷渣器操作器的输入为总操作器的输出。根据此模型在郭家湾电厂2台机组进行组态和调试，调试过程中最终确定的调节器参数为：左侧（右侧）冷渣器有1台在自动状态时，比例放大系数Kp=0.8；积分时间Ti=40；微分时间Td=50；微分系数Kd=0.05。左侧（右侧）冷渣器有2台在自动状态时，比例放大系数Kp=0.7；积分时间Ti=55；微分时间Td=50；微分系数Kd=0.05。左侧（右侧）冷渣器有3台在自动状态时，比例放大系数Kp=0.6；积分时间Ti=70；微分时间Td=50；微分系数Kd=0.05。控制模型的特点：①采用水冷风室压力作为被调量，避免连续放渣或定期排渣对风室压力造成的波动，同时起到稳定循环流化床锅炉床压的作用。②根据锅炉两侧对应的冷渣器运行数量采用变积分和变比例控制方案，使自动控制更加精确快速。③使用总调节器分别控制3台分操作器的控制方式，使系统控制更加集中。④在控制系统中加入出渣温度和出水温度超限自动减速功能，将保护功能也融入自动控制系统。逻辑模型应用后，提高了两侧水冷风室压力的稳定性，减少了运行人员的劳动强度，各种工况下，水冷风室压力能够稳定在9.5~11.5kPa，当出水温度和出渣温度超限时，也能够立即对冷渣器自动减速保障系统安全。

3结论

自动控制模型将水冷风室压力作为被调量，采集各台冷渣器的自动状态，经过变积分、变比例调节来适应除渣系统的运行工况，同时设计了冷渣器出渣温度和出水温度越限自调整功能，实现自动控制的同时满足了安全需要。经过在郭家湾电厂的使用，证明此控制模型可以有效实现双床循环流化床锅炉除渣系统的自动调节，能够避免手动控制存在的两侧风室压力调节不准确、冷渣器出渣温度超限、冷渣器出水温度超限、劳动强度大、费电、安全稳定性差等缺点。

参考文献：

［1］佟春海.郭家湾电厂SNCR脱硝自动控制系统优化［J］.神华科技，2018，16（2）：57-58，73.

自动微分范文篇9

验证换料后的堆芯装载图；测量与核电站正常安全运行有关的物理参数，包括热态零功率下，所有棒全提时临界硼浓度和等温温度系数，控制棒的积分价值，硼微分价值；验证有关核参数的安全准则和设计准则；验证堆芯换料设计的有效性；为实施提升功率试验创造良好的条件。

2零功率物理试验前提条件

在此试验前应完成初始临界试验，包括校验RPN系统各测量量程之间的线性度和相互之间的重叠范围；通过功率量程本底噪声水平和多卜勒水平的寻找，确定不发生核加热中子通量的范围，从而建立零功率物理试验的功率水平范围；完成反应性仪的校验等。多普勒效应后RPN源量程保护定值调整已完成。反应性仪的测量精度已经验证，并能正确有效地使用。RCV和REA的调硼系统（稀释或硼化）均能正常工作。控制棒驱动机构、棒位操作和指示系统均能正常投入使用。质量安全计划上相应的操作，相关部门已签字，H点值长已签字。

3机组初始状态及试验过程中风险控制

3.1机组状态反应堆冷却剂的压力为155-2+0巴。停堆棒和功率控制棒N均全提到225步；R棒的棒位为初始临界实验结束时的棒位。一回路硼浓度为R棒170步时临界硼浓度附近。反应堆冷却剂系统的平均温度291-1+2℃，通过GCT-A进行调节。如果需要硼化或稀释，应投入2个下泄孔板及全部通断式加热器，并以恒定的速率完成，以获得反应性的线性变化。3.2反应堆功率控制反应堆在稳定临界，通量在零功率物理试验范围内（应在P6与多普勒效应点之间，L204为功率量程的1.04E-7A至7E-9A之间，以反应性仪读数为准，此时KIT内SRC/PRC显示无效，IRC显示要比反应性仪显示大10倍左右），应在反应性仪15%~90%之间。3.3反应性控制提棒和稀释时，通量的增加速率不超过10倍/min（即倍增周期大于18s），以防止反应性引入过快，造成停堆（此时IRC跳堆定值仅为2.5E-5A）。试验期间，如果出现不可控的降温或其它原因引起的反应性急剧增加，反应堆紧急停堆保护动作，执行DEC规程。3.4棒位控制允许在换料后的试验期间，慢化剂的温度系数稍许为正，R棒的插入限值可适当突破。但试验结束后，需恢复到慢化剂温度系数为负的正常运行条件，R棒恢复正常棒位要求。温度控制棒R以及停堆棒S，除非特别规定，应在“ALGN2”下移动，由RGL001QM至RGL005QM及RGL013QM和RGL014QM可读出棒位。当用R棒交替法测停堆棒组的积分价值时，在RGL003CC转到被选棒之前，必须确保被选棒的子组同步。功率控制棒，当不希望以叠步方式移动时，须使用“ALGN1”模式。在RGL004CC转到“ALGN1”位置以前，RGL002CV需转到“VALALGNT1”位置，否则将出现RGL001AA报警。在此设置下，被测棒组的棒位计数器不计数，该棒棒位指示由RGL004CC下的RGL016QM读出，这个读数仅指示棒位移动步数。在此情况下必须记录这个计数器的初始和终了值，防止RGL016QM“复零”，否则将失去棒位信息。例如：用R棒测N1棒的积分价值时，N1全提出，棒位计数器指示为225步，RGL016QM指示000，N1插入到5步时，棒位计数器指示仍为225步，RGL016QM指示780，这就是说：225-5=220步（从初始棒位插入步）。控制棒在“ALGN1”模式下移动时，需确保棒组不超过棒位的上限及下限，以防止棒位超出正常棒位。对上例，当RGL016QM为780时，N1棒必须停止插入；当RGL016QM为000时，N1棒必须停止提出。3.5蒸发器水位及GCT-A控制进行等温温度系数测量时注意控制3台SG压差以及SG水位的稳定，它经常导致试验不成功。由于新技术规范要求三台SG的GCT-A不可以同时置手动，这样就使获得4°C/min的升降温速率很困难。我们可以置一台SG的GCT-A手动，对于1号机，建议使用SG1，2号机建议使用SG3。因为1号机的温度梯度来源于RCP029MT，而2号机来源于RCP056MT，用相应环路的GCT-A可以使温度梯度更快地反映GCT-A的调节。要求4°C/min的速率是为了使燃料棒的温度更加接近于RCP冷却剂的温度，由于燃料温度不可测，只能通过冷却剂温度得到。温度变化太快会使冷却剂温度与燃料温度不一致，太慢则使试验时间延长。从运行控制角度来说，越慢越容易实现控制稳定。因此，我们可以缓慢增加温度梯度到4℃/h，每次调节，开度变化即停止，观察梯度稳定后再接着调节，调节时，瞬态时的温度偏差的变化快于温度梯度的变化（因为温度偏差来源于整定温度与三个环路平均温度最大值，而梯度仅本环路温度信号，且数值较小，变化较慢，体现不出微分的优点）。在另外两个置自动的GCT-A有开度时，调节较慢，因为降温时手动的开大时，自动的会补偿关小。待自动的全关后，手动的调节要更加缓慢。升温时，关小手动的，自动的会开启，造成调节扰动，可以事先增加开启整定值，使自动保持全关。回复时再将定值调一致。冷却时密切关注三台SG压力，避免两高一低导致安注，需要缓慢调节，且保持给水稳定，三台SG压差大约在1巴时，就同步稳定变化了。同时关注SG水位，若APA或APD调节，小阀在自动，则汽水压差在5巴左右为宜，太大则小阀会全关，失去调节水位的作用，太小则在升温时，由于SG压力增加，给水有可能克服不了汽压，造成SG水位下降而停堆，给水流量调节过大又会影响一回路温度，导致梯度太大而不满足试验要求。

4实施步骤

①全部控制棒在顶端时的临界硼浓度测量CBAROMES；②全部控制棒在顶端时的等温温度系数测量α（ISO）AROMES；③R棒组的微分和积分价值测量及硼微分价值测量；④用R棒组替代法测量其它控制棒组的积分价值；⑤确定慢化剂温度系数为负的硼浓度限值；⑥反应性仪精度校验试验。

5调整机组状态，试验结束

检查安全棒棒位，确保全提（225步），检查R棒棒位，如低于低低限（197），插入功率棒，提出R棒，调整R棒到低低限以上，最好在调节带内。检查功率棒棒位是否在零功率棒位之上，否则进行硼化提棒至零功率棒位以上。硼化结束后等待一回路均匀半小时，手工分析回路和稳压器中的硼浓度，三次取样，每次间隔15分钟，要求三次测量之差小于5ppm，稳压器与回路硼浓度之差小于20ppm。确认硼浓度已满足慢化剂温度系数为负的要求，否则通过插功率棒稀释的方法完成。一回路平均温度满足Tavg＝291.4-1+2℃。机组状态调整结束后，零功率试验主要内容已完成，机组已恢复到目标状态。

6几个公式

自动微分范文篇10

1．1PID控制原理［1，2］

常规PID控制系统原理框图如图1所示。

PID控制器是一种线性控制器，它根据给定值r（t）与实际输出构成控制偏差：

将此偏差的比例（P）、积分（I）和微分（D）通过线性组合构成控制量，对被控对象进行控制。其控制规律为：

式中，Kp为比例系数，T1为积分时间常数，TD为微分时间常数。

在PID控制中，比例项用于纠正偏差，积分项用于消除系统的稳态误差，微分项用于减小系统的超调量，增加系统稳定性。PID控制器的性能就决定于Kp、T1和TD这3个系数。如何选用这3个系数是PID控制的核心。

1．2数字PID控制算法选择

设计和调整数字PID控制器的任务就是根据被控对象和系统要求，选择合适的PID模型，将其进行离散化处理，编出计算机程序由微处理器实现，最后确定KP、T1、TD、和T，T为采样周期。微处理器控制是一种采样控制，它只能根据采样时刻的偏差值计算控制量，因此，必须对PID模型进行离散化处理。

用矩形方法数值积分代替式（3）中的积分项，对式（3）中的导数项用后向差分逼近，经推理可得到基本PID控制的位置式算法：

式中k——采样序号，k＝0，1，2，……

U（k）——第k次采样时刻输出值

E（k）——第k次采样时输入的偏差值

E（k－1）——第（k－1）次采样时刻输入的偏差值

K1——积分系数，K1＝KpT／T1

KD——微分数系，KD＝KpTD／T1

在数字控制系统中，PID控制规律是用程序来实现的，因而具有更大的灵活性。由于基本PID控制中引入了积分环节，其目的主要是为了消除静差，提高精度。但在柴油机调速过程中，突加突减负载时，会引起转速的较大波动，导致短时间内转速出现较大偏差，通过PID积分运算积累，超调量过大，系统产生振荡，严重影响发电机组输出电能的品质。为避免PID控制中积分项引起的超调，提高其调节品质，拟采用积分分离法对基本PID控制进行改进，简称变速积分PID。变速积分PID的基本思路是设法改变积分项的累加速度，使其与偏差大小相对应，偏差越大，积分越慢；反之，则越快。

式中，A、B为积分区间。

变速积分PID算法为：

式中，U1（k）为第k次采样时刻PID运算的积分部分输出值。

采用变速积分PID控制，系统具有以下特点：用比例消除大偏差，用积分消除小偏差，可完全消除积分饱和现象；各参数容易整定，易实现系统稳定，而且对A、B两参数不要求十分精确；超调量大大减小，改善了调节品质，适应性较强。

2柴油发电机组数字调速系统中PID控制参数整定［3，4］

数字PID控制参数整定的任务主要是确定数字PID的参数KP、T1、TD和T。

对于简单控制系统，可采用理论计算方法确定这些参数。但由于柴油机调速系统的工况较为复杂，其数学模型并非十分精确，在此，采用工程整定常用的扩充临界比例带法，结合经验法再对参数进行调整，得到最终的PID参数。

（1）采样周期T的选择

在数字控制系统中，采样周期T是一个比较重要的因素，采样周期的选取，应与PID参数的整定综合考虑。

首先，采样周期T的选取应满足以下要求：远小于对象扰动周期；比对象时间常数小得多；尽量缩短采样周期，以改善调节品质。

该系统中，PID调节控制过程是在定时中断状态下完成的，因此，采样周期T的大小必须保证中断服务程序的正常运行。在不影响中断程序运行的情况下，可取采样周期T＝0．1τ（τ为柴油机的纯滞后时间）。当中断程序运行时间Tz大于0．1τ时，则取T＝Tz，

（2）临界振荡周期Ts的确定

初始确定数字PID参数时，在用上述方法确定采样周期T的条件下，从调速系统的PID调节回路中，去掉数字控制器的微分控制作用和积分控制作用，只采用比例调节环节来确定系统的振荡周期Ts和临界比例系数Ks。由单片机系统自动控制比例系数KP，并逐渐增大Kp，直到系统出现持续的等幅振荡，然后由单片机系统自动记录并显示调速系统发生等幅振荡时的临界比例度δ和相应的临界振荡周期Ts。

控制度就是以模拟调节器为基础，定量衡量数字控制系统与模拟调节器对同一对象的控制效果。控制效果就是采用某一积分准则，根据系统在规定的输入下的输出响应，使用该准则取最小值时的最

如前所述，采样周期T的长短会影响系统的控制品质，同样是最佳整定，数字控制系统的品质要低于模拟系统的控制品质。即控制度总是大于1的，且控制度越大，相应的数字控制系统品质越差。

为获得与模拟控制器相当的品质，控制度选为1．05。不同控制度时，扩充临界比例带法PID参数计算公式

（4）KP、K1、KD、T的求取

根据实验所得Ks和Ts及选定的控制度，按表1计算出数字PID参数Kp、T1、TD和T。

（5）控制效果的调节

按求得的参数值在调速控制系统中运行，并观察控制效果。如控制效果达不到控制要求，可基于以下原则，根据经验法对参数做适当调整。

①增大比例系数Kp，将加快系统的响应速度，但过大会使系统产生较大超调，甚至产生振荡。

②增大积分时间T1，有利于减小超调，减少振荡，使系统更加稳定，但会增加系统过渡过程时间。

③增大微分时间常数TD有利于加快系统的响应，使超调减小，稳定性增加，但系统对扰动的抑制能力减弱，对扰动有较敏感的响应。

基于上述原则，调整PID参数时，应先比例、后积分、再微分进行调整。

参考文献：

［1］陶永华，尹怡欣，葛芦生．新型PID控制及其应用［M］．机械工业出版社，1998．

［2］王福瑞．单片微机测控系统设计大全［M］．北京航空航天大学出版社，1998．