自動微分技術(shù)及其在軌跡優(yōu)化中的應(yīng)用

王晨曦，李新國

(西北工業(yè)大學航天學院，陜西西安 710072)

引言

飛行器軌跡優(yōu)化問題多采用依賴于導(dǎo)數(shù)信息數(shù)值優(yōu)化算法來求解，特別是對于得到廣泛應(yīng)用的直接法，其基本思想是將求解無限維的函數(shù)優(yōu)化問題轉(zhuǎn)化為有限維參數(shù)非線性規(guī)劃問題。在采用非線性規(guī)劃算法求解時，需要計算目標函數(shù)及約束對參數(shù)的微分信息，微分信息計算的精度對非線性規(guī)劃的精度及收斂性有較大影響。因此如何快速、準確、高效地計算微分信息，對于軌跡優(yōu)化設(shè)計具有重要意義。

現(xiàn)有的導(dǎo)數(shù)計算方法種類繁多，但計算的效率、精度卻存在差別，常見的方法有:有限差分法、符號微分法、復(fù)變量法和自動微分方法［1］，其中自動微分法具有很高的計算精度，理論上可以達到機器精度，同時又方便編程實現(xiàn)，近年來在多學科設(shè)計優(yōu)化［1］、飛行器氣動外形優(yōu)化［2］等數(shù)值優(yōu)化領(lǐng)域得到廣泛應(yīng)用。本文將自動微分技術(shù)與內(nèi)點非線性規(guī)劃方法相結(jié)合，用于求解直接軌跡優(yōu)化問題。

1 自動微分技術(shù)

自動微分(Automatic Differentiation，AD)又稱算法微分，其理論和實現(xiàn)方法在近幾年得到迅速發(fā)展。當前，各領(lǐng)域的工程設(shè)計問題越來越依賴于計算機，通常都是將問題的數(shù)學模型轉(zhuǎn)化為計算機程序，然后基于計算程序來實現(xiàn)基于仿真的設(shè)計及優(yōu)化。以計算機語言描述的數(shù)學模型不同于簡單的數(shù)學表達式，而是成千上萬行的程序代碼，不可能依賴于手工求導(dǎo)，必須尋求適用于程序微分信息計算的新途徑，自動微分技術(shù)恰好是這樣一種技術(shù)。

自動微分的實質(zhì)是基于鏈式法則(Chain Rule)自動生成計算程序微分信息或微分計算程序的技術(shù)。其基本思想是無論描述函數(shù)的計算機程序多復(fù)雜，它本質(zhì)上都是執(zhí)行一系列的初等計算、初等函數(shù)及邏輯運算的組合。對這些初等運算迭代運用鏈式求導(dǎo)法則，可以自動并精確地得到目標函數(shù)的任意階微分信息，與其他微分計算方法相比具有明顯優(yōu)勢［3］。

1.1 自動微分模式

當通過鏈式法則獲取了復(fù)雜函數(shù)(程序)的微分計算關(guān)系之后，需要將初等函數(shù)微分進行連續(xù)累積得到復(fù)雜函數(shù)的微分信息，根據(jù)累積方向的不同可以分為前向模式(Forward Mode)和逆向模式(Reverse Mode)兩種。

顧名思義，前向模式是按照從自變量(Independent Variable)到因變量(Dependent Variable)的方向應(yīng)用鏈式法則逐次計算每個計算單元(初等函數(shù))的局部偏導(dǎo)數(shù)，得到所需導(dǎo)數(shù)信息;逆向模式與前向模式相對應(yīng)，是按照從因變量到自變量反方向逐次計算微分信息。首先前向執(zhí)行函數(shù)代碼列表，并記錄中間變量^ti，然后通過^ti傳遞因變量與自變量之間的導(dǎo)數(shù)關(guān)系，反向累積獲得微分信息(見表1)。

表1 自動微分的兩種模式

兩種計算模式各有特點，逆向模式計算量與設(shè)計變量的數(shù)目無關(guān)，適用于設(shè)計變量數(shù)目較多的情況;而前向模式適用于函數(shù)數(shù)目多于設(shè)計變量數(shù)目的情況?？筛鶕?jù)對應(yīng)問題的特點選擇合適的計算模式。以下面的算例說明這兩種計算模式［4］。

1.2 自動微分方法和工具

自動微分程序的實現(xiàn)方法分為兩種:源代碼轉(zhuǎn)換法(Source Transformation Method，STM)和算子重載法(Operator Overloading Method，OOM)。

STM通過對程序代碼進行分析，利用編譯技術(shù)，將原計算程序轉(zhuǎn)換為可進行函數(shù)導(dǎo)數(shù)計算的新函數(shù)(見圖1)，它可用于任何程序語言［5］。STM需要針對源代碼程序語言的特性，設(shè)計專門的語法、詞法分析器，由于具有針對性，因而可獲得優(yōu)化的計算程序，具有高的計算效率，但是對于某些語法靈活的程序語言(如CC++)，設(shè)計自動微分代碼轉(zhuǎn)換的AD工具是一項非常有挑戰(zhàn)性的工作。目前，針對C語言，通過代碼轉(zhuǎn)換策略實現(xiàn)的AD工具僅有美國Argonne國家實驗室開發(fā)的ADIC?；诖a轉(zhuǎn)化實現(xiàn)的軟件包有:ADIC，ADIFOR，OpenAD和TAPENADE等。

圖1 自動微分的代碼轉(zhuǎn)化實現(xiàn)方式

OOM是利用計算機程序語言的特性，采用算子重載技術(shù)實現(xiàn)微分的同步計算［5］。目前許多高級計算機程序語言均支持這一技術(shù)，如C++，F(xiàn)ortran-90，Ada和Matlab等。采用算子重載法無需改變原有程序的程序代碼，而是通過聲明特定類型的變量實現(xiàn)自動微分，如圖2所示。例如本文采用的ADOL-C工具包［6］，只需將期望獲得微分信息的自變量和因變量聲明為“active”變量，即可在程序運行時生成包含微分信息的“tape”文件，與STM技術(shù)實現(xiàn)的自動微分相比更加便于使用，不生成微分計算程序，但是這種方法無法通過編譯實現(xiàn)微分計算過程的優(yōu)化?；谒阕又剌d實現(xiàn)的軟件包有:ADC，ADF，ADOLC，F(xiàn)ADBADU++，CppAD和MAD等。本文采用ADOL-C作為自動微分工具。

圖2 自動微分的算子重載實現(xiàn)方式

2 軌跡優(yōu)化方法

軌跡優(yōu)化本質(zhì)上是一種無限維的函數(shù)最優(yōu)化問題，直接軌跡優(yōu)化方法的基本思想是將函數(shù)優(yōu)化問題轉(zhuǎn)化為有限維參數(shù)優(yōu)化問題。本文首先利用直接配點法將其轉(zhuǎn)化為非線性規(guī)劃問題，然后結(jié)合自動微分和內(nèi)點非線性規(guī)劃算法進行求解。

2.1 最優(yōu)控制問題

最優(yōu)控制問題的求解是找到滿足一定約束的最優(yōu)控制變量u*(t)及其對應(yīng)的最優(yōu)狀態(tài)(軌跡)x*(t)，并使目標函數(shù)最小。其中x(t)∈Rnx和u(t)∈Rnu為狀態(tài)函數(shù)和控制函數(shù)向量，目標函數(shù)為:

2.2 直接配點法

直接配點法同時將狀態(tài)變量和控制變量離散化后實現(xiàn)最優(yōu)控制問題向參數(shù)優(yōu)化問題的轉(zhuǎn)化，將離散點處的狀態(tài)及控制量同時作為優(yōu)化設(shè)計參數(shù)。直接配點法的核心內(nèi)容是通過一定的轉(zhuǎn)化方式，將動力學約束轉(zhuǎn)化為非線性規(guī)劃的等式約束。根據(jù)不同的積分近似準則，直接配點法分為Trapezoidal法、Hermite-Simpson法［7］等局部配點法和 Legendre偽譜法［8］、Chebyshev 偽譜法［9］、Guass 偽譜法［10］等全局配點法。

通過直接配點法可將最優(yōu)控制問題轉(zhuǎn)化為下面的非線性規(guī)劃問題(有關(guān)直接配點法的詳細描述參見文獻［7-10］):

本文采用C++編制了Trapezoidal法、Hermite-Simpson法、Legendre偽譜法的計算程序。

2.3 內(nèi)點非線性規(guī)劃

內(nèi)點法［11］的基本思想是將有約束問題轉(zhuǎn)換為一系列障礙問題 (Barrier Problem)，在可行域內(nèi)迭代求解障礙問題極值點的方法。

早期的內(nèi)點法均采用精確罰函數(shù)法保證算法的收斂性，即將約束函數(shù)作為懲罰項，新的評價函數(shù)為目標函數(shù)與懲罰函數(shù)兩項之和。近期由Fletcher等提出了過濾方法代替評價函數(shù)來保證求解NLP問題的收斂性，其核心思想是減小性能指標函數(shù)值或減小約束函數(shù)違背值的搜索點均將被算法接受，而非僅接受減小評價函數(shù)值的點。這使得內(nèi)點法逐漸成為高效、穩(wěn)定的求解大規(guī)模NLP問題的一種方法。內(nèi)點法數(shù)學描述如下:

對于式(9)中描述的不等式約束 g(X)和h(X)，為描述方便，令

式中，X為nx維實數(shù)向量;c為nc維實數(shù)向量。引入nx維松弛向量s將不等式約束轉(zhuǎn)換為等式約束，即c(X)=ˉ(X)－s=0，s≤0，則優(yōu)化問題可表示為僅包含等式約束的優(yōu)化問題:

內(nèi)點法將上述問題轉(zhuǎn)換為求解一系列障礙問題的極值:

在內(nèi)點法的求解過程中，障礙參數(shù)μ＞0逐漸減小并趨于零，障礙項用于防止設(shè)計變量X(i)達到邊界，從而保證其始終在可行域內(nèi)。應(yīng)用Homology方法可得到問題(12)的主-對偶極值點條件:

式中，λ∈Rnc和z∈Rnx分別為式(11)中等式約束及邊界約束對應(yīng)的拉格朗日乘子向量;:=diag(X)(表示取 X對角線元素構(gòu)成對角矩陣)，同理:=diag(z);e為單位向量。當式(13)中的Homology參數(shù)μ=0時，式(13)和X，z≥0就可以構(gòu)成原問題式(12)的KKT條件。內(nèi)點法首先求解μ為固定值下障礙問題的近似解，然后減小μ繼續(xù)求解，直到得到μ趨于0(趨近原問題的解)。在求解固定μ的障礙問題中，該算法應(yīng)用了衰減牛頓法對原-對偶方程(13)迭代求解，通過式(13)簡化的線性化形式構(gòu)造搜索方向。采用線性搜索過濾法構(gòu)造搜索步長，并提供二階修正方法對步長進行校正。

本文采用的內(nèi)點法求解器是由美國卡耐基梅隆大學AndreasWaechter等人［11］采用C++開發(fā)的軟件包IPOPT(Interior Point OPTimizer)。

3 飛行器軌跡優(yōu)化設(shè)計

3.1 飛行器數(shù)學模型

建立飛行器數(shù)學模型時做以下假設(shè):(1)假設(shè)地球為均質(zhì)不旋轉(zhuǎn)圓球;(2)忽略附加哥氏力的影響;(3)采用指數(shù)大氣模型和平方反比重力模型。基于上述假設(shè)建立如下動力學模型:

式中，m為飛行器質(zhì)量;L，D分別為飛行器所受升力和阻力;H為飛行高度;γ為航跡角，ψ為航向角，北向為0°，東向為90°;θ為經(jīng)度，φ為緯度;α為攻角，σ為傾側(cè)角;g為重力加速度;Re為地球半徑。氣動力模型詳細參見文獻［7］。

目標函數(shù)為:

軌跡優(yōu)化的目標函數(shù)為最大橫程，計算過程中是對其等效的目標函數(shù)進行計算。

3.2 計算與分析

本文采用文獻［4］中航天飛機的算例進行計算，并與文獻中的計算結(jié)果進行比較，以驗證算法的正確性。在進行自動微分和有限差分比較時，采用了ADOL-C中前向自動微分和中心有限差分。研究中采用的編譯環(huán)境及計算機配置分別為:Visual 2008(SP1)和Intel(R)Core2 Duo CPU-E8400(3.0 GHz)處理器，內(nèi)存2.00 G。

首先對研究中編制的直接配點法程序進行了驗證，采用3.1中的模型及約束條件，對三種配點法方法進行了驗證，離散點個數(shù)采用了逐步細化策略，微分方法采用自動微分。三種方法(Trapezoidal法(TRA)、Hermite-Simpson法(H-S)、Legendre偽譜法(LPS))計算結(jié)果與文獻［7］計算結(jié)果的比較如表2所示，結(jié)果表明三種算法正確、可行，計算結(jié)果與文獻中非常接近。圖3～圖6給出了其中一組優(yōu)化計算的結(jié)果。

圖3 高度、速度隨時間變化曲線

表2 多種算法優(yōu)化結(jié)果比較

圖4 航跡角、航向角隨時間變化曲線

圖5 飛行軌跡地面投影(緯度-經(jīng)度)

圖6 攻角、傾側(cè)角隨時間變化曲線

為了比較自動微分和有限差分在軌跡優(yōu)化中的應(yīng)用，采用了不同離散方法(三種)進行計算，非線性規(guī)劃的最大允許迭代次數(shù)為1 000，誤差精度為5e－6;不進行網(wǎng)格細分的計算中離散節(jié)點個數(shù)為60。表3中給出了不進行網(wǎng)格細分的計算結(jié)果，結(jié)果表明采用ADM后，在收斂時間、穩(wěn)定性方面得到了改善，同時可看到目標函數(shù)、約束函數(shù)的評價次數(shù)極大地減少，這對于目標及約束函數(shù)計算量大的問題極為有利(由于本文在氣動力及大氣模型計算時采用了擬合模型，降低了目標及約束的計算代價，所以在計算時間方面差別不太明顯)。其中LPS中采用FDM計算方案達到最大迭代次數(shù)時沒有獲得滿足精度要求的解，而采用ADM的方案獲得了收斂解，這表明高精度的自動微分法有利于大規(guī)模非線性規(guī)劃問題的收斂。

表3 ADM和FDM用于DCNLP軌跡優(yōu)化結(jié)果

4 結(jié)束語

在基于直接配點的軌跡優(yōu)化設(shè)計中，優(yōu)化計算的效率受微分信息的影響很大，快速、精確地計算所需的微分信息對問題的求解至關(guān)重要。本文采用ADM進行微分計算，相對于傳統(tǒng)FDM對目標函數(shù)、約束函數(shù)的計算次數(shù)大大減少，同時由于ADM可獲得僅受機器精度限制、而與計算步長無關(guān)的高精度微分信息，使得非線性規(guī)劃算法的收斂速度有所提高，這對飛行器的軌跡快速優(yōu)化、快速生成具有重要意義，有著良好的應(yīng)用前景。

［1］ 顏力，陳小前，王振國.飛行器MDO中靈敏度計算的自動微分方法［J］.國防科技大學學報，2006，28(2):13-16.

［2］ 潘雷，谷良賢，龔春林.改進自動微分方法及其在飛行器氣動外形優(yōu)化中的應(yīng)用［J］.西北工業(yè)大學學報，2007，25(3):398-401.

［3］ 蔣占四，吳義忠，蔣慧.敏度分析的數(shù)值方法比較研究［J］.計算機與數(shù)字工程，2009，37(5):1-5.

［4］ Kedem G.Automatic differentiation of computer programs［J］.ACM Transactions on Mathematical Software，1980，6(2):150-165.

［5］ 張春暉，程強，曹建文.針對C語言的自動微分系統(tǒng)及其應(yīng)用［J］.計算機應(yīng)用研究，2009，26(1):155-159.

［6］ Griewank A，Juedes D.ADOL-C:a package for the automatic differentiation of algorithms written in C/C++［J］.ACM Transactions on Mathematical Software，1996，22(2):131-167.

［7］ Betts J T.Practical methods for optimal control and estimation using nonlinear programming(2rd)［M］.Philadelphia:SIAM Press，2010.

［8］ Fahroo F，Ross IM.Advances in pseudospectralmethods for optimal control［R］.AIAA 2008-7309，2008.

［9］ Fahroo F，Ross I M.Direct trajectory optimization by a Chebyshev pseudospectral method［J］.Journal of Guidance，Control，and Dynamics，2002，25(1):160-166.

［10］ Garg D，Patterson M A，HagerW W，et al.An overview of three pseudospectralmethods for the numerical solution of optimal control problems［R］.AAS-2009-332，2009.

［11］ Walther Andrea，Griewank Andreas.On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming［J］.Mathematical Programming，2006，106(1):25-57.