曹小群
(國防科學(xué)技術(shù)大學(xué)計算機學(xué)院,長沙 410073)
(2012年11月11日收到;2012年12月27日收到修改稿)
近年來混沌控制和同步已經(jīng)成為非線性科學(xué)研究中的重要方向之一,在已知精確參數(shù)的前提下相應(yīng)的理論和方法已經(jīng)得到廣泛發(fā)展[1-3].在實際中,由于混沌系統(tǒng)的復(fù)雜性,某些參數(shù)經(jīng)常無法測量或確定;或者出于保密通信等特殊原因,系統(tǒng)的某些參數(shù)不可知.此時,要實現(xiàn)對混沌系統(tǒng)的控制或同步,上述方法有相當?shù)木窒扌?因此參數(shù)估計是混沌控制與同步中首先必須解決的課題,具有重要的現(xiàn)實意義.文獻[4]引入基于反饋的同步和自適應(yīng)控制方法,對若干混沌系統(tǒng)進行了參數(shù)估計;文獻[5]通過最小化平均同步誤差,對一個給定的混沌動態(tài)系統(tǒng)進行了參數(shù)估計;文獻[6]通過參數(shù)自適應(yīng)方法對目標系統(tǒng)的參數(shù)進行估計,達到了廣義同步的目的;文獻[7]提出了未知參數(shù)辨識觀測器的概念,并對Lorenz系統(tǒng)的參數(shù)進行了辨識;文獻[8]提出了一種基于變分原理的混沌系統(tǒng)參數(shù)估計方法,并對Lorenz系統(tǒng)和超混沌Chen系統(tǒng)的參數(shù)進行了估計.近年來一系列智能優(yōu)化算法已經(jīng)成功應(yīng)用于混沌系統(tǒng)參數(shù)估計中,如改進粒子群優(yōu)化算法[9]、混沌螞蟻群算法[10]、混合差分進化算法[11]和混合交叉進化算法[12]等.
非線性映射是混沌控制和同步領(lǐng)域中一類重要的混沌系統(tǒng),其未知參數(shù)估計也是一個必須解決的關(guān)鍵問題.非線性映射參數(shù)估計是離散動力系統(tǒng)研究中典型的反問題,和正問題不同,它一般通過觀測數(shù)據(jù)來反求系統(tǒng)的參數(shù)或初/邊值.因為觀測中不可避免地包含有噪聲以及系統(tǒng)固有的非線性,所以反問題求解過程中經(jīng)常出現(xiàn)不適定現(xiàn)象:解不一定存在、即使解存在也不惟一、在解存在惟一條件下也不穩(wěn)定(即解不連續(xù)依賴于觀測數(shù)據(jù))[13,14].因此非線性映射系統(tǒng)參數(shù)估計問題的求解通常采用特殊方法[8-13].本文基于二階離散變分原理[14-18]來估計非線性映射中的未知參數(shù).首先引入拉格朗日乘子向量,將以非線性離散系統(tǒng)為約束的最優(yōu)化問題轉(zhuǎn)換為無約束最優(yōu)化問題,并利用離散變分方法導(dǎo)出了伴隨方程和目標泛函梯度表達式;其次采用二階離散變分方法導(dǎo)出了二階伴隨方程和Hessian矩陣-向量乘積的通用公式;然后設(shè)計了二階離散變分方法估計非線性映射中未知參數(shù)的算法流程,并以此估計了Hyperhen′on映射[19]和二維拋物映射[20]中的未知參數(shù).仿真結(jié)果表明:該方法具有很高的估計精度,同時具有較好的抗噪聲性能.
考慮如下n維非線性映射系統(tǒng)的未知參數(shù)估計問題:
其中x0和xk∈Rn分別表示n維初值和狀態(tài)向量,M(xk,θ):Rn×Rm→Rn表示非線性映射的模式函數(shù)向量,θ =(θ1,θ2,···,θm)T∈ Rm表示需要估計的未知物理參數(shù)向量.目標是利用非線性映射的觀測數(shù)據(jù)準確辨識出不可測的未知參數(shù)向量θ,在觀測噪聲滿足高斯分布的情況下,該問題可通過對下面具有最小二乘意義的目標泛函求極小值而解決:
其中yko和Rk(k=1,2,3,···,N)分別表示映射系統(tǒng)的觀測向量和觀測誤差協(xié)方差矩陣,〈·,·〉表示歐氏空間的內(nèi)積.顯然,非線性映射的未知參數(shù)估計問題變成在離散狀態(tài)方程(1)式約束下的目標泛函最優(yōu)化問題.但是,由于非線性映射物理系統(tǒng)的強非線性和不穩(wěn)定性將導(dǎo)致目標泛函J(θ)中可能存在多個局部極值點,從而使得僅僅依靠一階導(dǎo)數(shù)(即目標泛函梯度 ?θJ=(?θ1J,?θ2J,···,?θmJ)T)信息很難準確估計未知參數(shù)值,進而使得如何精確計算目標泛函的高階信息(例如,與Hessian矩陣有關(guān)的信息)就成為一個重要問題.同時考慮到未知參數(shù)向量θ是通過非線性映射系統(tǒng)與目標泛函間接相關(guān)的,因此要準確地得到J關(guān)于θ的二階導(dǎo)數(shù)信息是非常困難的.本文利用二階離散變分方法給出了切線性方程、二階伴隨方程和準確計算Hessian矩陣(?θ2J)-向量乘積的公式,既可以避免Hessian矩陣直接求解的巨大計算代價,又可以加快最優(yōu)化算法的收斂速度和提高未知參數(shù)的估計精度.
首先引入n維拉格朗日乘子向量λk∈Rn(k=0,1,2,···,N),將以非線性映射系統(tǒng)為約束條件的最優(yōu)化問題(2)式轉(zhuǎn)化為下面(3)式表示的無約束最優(yōu)化問題:
Jk(θ)的一階變分可由下面的表達式給出:
其中
另外(3)式的一階變分為
上式中δxk和δλk分別表示非線性映射狀態(tài)向量和拉格朗日乘子向量的一階變分.由(1)式顯然有〈δλk+1,xk+1-M(xk,θ)〉=0.而目標泛函 J(θ)和向量函數(shù)M(xk,θ)的一階變分按定義可以進一步寫為
其中DMx k∈Rn×n和DθM∈Rm×n分別表示模式向量函數(shù)M(xk,θ)關(guān)于xk和θ的Jacobian矩陣.將(8)式中第二式代入(7)式右邊第三項,同時利用恒等式和伴隨性質(zhì)〈λ,A x〉=〈ATλ,x〉可得
將(8)式中第一式與(6)式和(9)式代入(7)式并合并同類項.考慮到一階變分δxk,δλk和δθ可以取任意值,且δx0=0,δxN/=0,進而可導(dǎo)出下面的方程組:
(10)式中的第一式稱為伴隨方程,λk表示k時次的伴隨狀態(tài)向量;從表達式易知,伴隨方程積分之前必須求得非線性映射的狀態(tài)向量序列xk,而且需按照順序 k=N,N-1,···,1 進行逆向積分.(10)式中第二式給出了伴隨方程的初值條件λN+1=0,因為是逆向積分,所以初值條件在末端時刻k=N+1給出.(10)式中的第三式是目標泛函關(guān)于未知參數(shù)向量梯度的表達式,計算梯度之前需要確定狀態(tài)向量xk和伴隨向量λk.在本小節(jié)中利用一階變分方法獲得了伴隨方程和目標泛函梯度的公式,通過計算目標泛函梯度向量的值,進而利用基于一階導(dǎo)數(shù)信息的最優(yōu)化算法(共軛梯度方法和擬牛頓方法)可以估計未知參數(shù)[8,13].但是,在混沌系統(tǒng)參數(shù)估計中,完全依賴一階導(dǎo)數(shù)信息的最優(yōu)化算法通常收斂速度較慢,對于具有較高維數(shù)的強非線性混沌系統(tǒng)尤其如此.因此為了加快收斂速度,有必要研究能同時利用導(dǎo)數(shù)和精確Hessian矩陣信息的非線性映射未知參數(shù)估計方法.
為了導(dǎo)出Hessian矩陣-向量乘積及相關(guān)的二階伴隨方程的通用表達式,首先對(1)式求一階和二階變分,則分別得到相應(yīng)的切線性系統(tǒng)(TLS):
和二階切線性系統(tǒng):
(11)式和(12)式中的δxk和δθ表示一階變分或切線性向量,而δ2xk和δ2θ表示二階變分或二階切線性向量.
其次對(10)式中的伴隨方程直接求一階變分,并將伴隨向量λk的一階變分δλk用ζk表示(即ζk=δλk),則可得
ζk稱為二階伴隨向量,考慮到λN+1=0,進而可確定二階伴隨量的末端值為ζN+1=0.下面推導(dǎo)Jacobian矩陣Dhxk,DMxk和DθM的一維變分表達式.
因為
所以有
同理可得非線性模式函數(shù)向量Jacobian矩陣的一維變分表達式:
將(15)和(16)式代入(12)和(13)式,可得:
(17)式被稱為非線性映射的二階切線性模式,刻畫的是二階切線性向量從時刻k=0到k=N的演化;相應(yīng)地,(18)式稱為二階伴隨模式,刻畫的是二階伴隨向量從時刻k=N到k=0的演化.對(17)和(18)式左右兩端分別用-λk+1和δxk求內(nèi)積,并將兩式所得到的結(jié)果相加,則可得:
為了簡化(19)式,利用伴隨性質(zhì)后有以下等式成立:
另一方面,由一階切線性模式和一階伴隨模式可導(dǎo)出以下等式:
將(20)式和(21)式都代入(19)式,進行化簡后可得:
將鏈式法則應(yīng)用到(4)式上,可以獲得Jk(θ)的二階變分表達式:
另外對 δJ(θ)= 〈?θJ,δθ〉兩邊進行一階變分運算,可以導(dǎo)出下式:
將(10)式中目標泛函梯度表達式代入(24)式,同時對(23)式左右兩邊從k=1到N求和,再考慮到進而得到:
對(22)式左右兩邊從k=0到N-1求和,并與(25)式相加,化簡后得到
考慮到一階變分δθ的任意性,從而有
(27)式是目標泛函關(guān)于未知物理參數(shù)的Hessian矩陣-向量乘積的表達式,其顯式給出了可以改善最優(yōu)化算法性能的二階導(dǎo)數(shù)信息,與Hessian矩陣相乘的向量δθ表示在每次最優(yōu)化迭代過程中未知參數(shù)向量逼近物理參數(shù)真值時的增量.在本小節(jié)中利用二階變分方法給出了切線性方程、二階伴隨方程以及計算目標泛函關(guān)于未知物理參數(shù)的Hessian矩陣-向量乘積的公式.從(27)式可知,在計算?θ2Jδθ的值之前需要確定狀態(tài)向量xk,一階切線性向量δxk,伴隨向量λk和二階伴隨向量ζk的值,因此需要依次求解非線性映射正模式(1),切線性模式(11),一階伴隨模式(10)和二階伴隨模式(18).其中,非線性映射系統(tǒng)和切線性模式需要進行正向積分,一/二階伴隨模式需要進行逆向積分.
在導(dǎo)出目標泛函梯度?θJ和Hessian矩陣-向量乘積?θ2Jδθ的公式之后,可以選擇合適的牛頓型最優(yōu)化算法對未知物理參數(shù)向量進行迭代估計.本文中選用的是牛頓-共軛梯度方法(Newton-conjugate gradient method,NCG)[21,22].在NCG方法中使用共軛梯度方法近似求解牛頓方程?2J(θk)sk=-?J(θk),當牛頓方程的殘余量足夠小(即滿足(28)式)或者遇到負曲率方向時,對共軛梯度方法能進行足夠精確地截斷.
在(28)式中ηk∈(0,1),確定方向矢量sk后,使用簡單Armijo線搜索方法[21]能計算出步長大小αk.因為完整Hessian矩陣的直接計算和存儲代價太高,所以不適合快速算法.而NCG方法的優(yōu)點是只需要使用Hessian矩陣-向量乘積的值,因此本文3.2部分中利用二階離散變分方法給出的(27)式正好能滿足此要求.二階離散變分方法估計非線性映射系統(tǒng)的未知參數(shù)的具體計算流程如圖1所示,具體步驟說明如下.
第1步 設(shè)定未知物理參數(shù)向量的初始猜測值;
第2步 利用參數(shù)向量的新估計值,對非線性映射模式(1)進行正向積分,獲得系統(tǒng)狀態(tài)量的演化軌跡xk并加以存貯;
圖1 二階變分方法估計非線性映射參數(shù)的流程圖
第3步 利用xk和yok根據(jù)(2)式計算目標泛函 J(θ)的值;
第4步 利用xk和yok和伴隨初值條件,對伴隨模式從k=N到k=1進行逆向積分,求得伴隨向量的演化軌跡λk并加以存貯;
第5步 利用xk和λk根據(jù)(10)式中第三式求目標泛函關(guān)于各個未知物理參數(shù)的梯度向量值;
第6步 如果滿足程序終止條件(如達到所要求的收斂精度或預(yù)先指定的最大迭代次數(shù)),則終止程序,同時獲得非線性系統(tǒng)的參數(shù)估計值;若不滿足,則繼續(xù)進行下面的步驟;
第7步 利用演化軌跡xk和切線性量初值條件,對切線性模式(11)進行正向積分,獲得切線性量或擾動量的演化軌跡δxk并加以存貯;
第8步 利用xk,δxk,λk和二階伴隨初值條件ζN=0,對二階伴隨模式(18)從k=N到k=1進行反向積分,求得二階伴隨向量的演化軌跡ζk并加以存貯;
第9步 利用xk,δxk,λk和ζk按照(27)式計算Hessian矩陣與未知參數(shù)向量增量之間的乘積;
第10步 將目標泛函梯度值和Hessian矩陣-向量乘積值輸入到基于NCG方法的最優(yōu)化程序中,計算出方向矢量和最優(yōu)步長,進而利用θk+1=θk+αksk對未知物理參數(shù)向量進行迭代更新.利用新估計值從第二步開始進行新一輪的迭代循環(huán).
首先以典型的Hyperheno′n映射為例,研究二階離散變分方法估計非線性映射系統(tǒng)中未知物理參數(shù)的有效性.Hyperheno′n映射的控制方程可表示如下:
當a=1.76和b=0.1時,該系統(tǒng)處于超混沌狀態(tài)[19].狀態(tài)向量表示為x=(x,y,z)T,θ=(a,b)T表示未知物理參數(shù)向量.目標是在獲得x的觀測量yko=(xok,yok,zok)T的情況下,辨識未知參數(shù)向量θ.通過比較(29)式與(1)式,易知非線性模式函數(shù)向量可表示為它的一階變分為
(30)式中 δxk=(δxk,δyk,δzk)T和 δθ =(δa,δb)T分別表示xk和θ的一階變分.根據(jù)(11)式可得切線性方程組:
首先引入與xk對應(yīng)的伴隨向量序列λk=(λk,pk,qk)T,并將 xk,yko,λk,θ,(DMxk)T和 (DθM)T等的表達式代入(10)式中的第一和第三式,同時考慮到伴隨初始條件,進而得到Hyperheno′n超混沌映射系統(tǒng)的伴隨方程和目標泛函關(guān)于未知參數(shù)的梯度公式,分別如(32)和(33)式所示:
其次將λk,
和Dxk代入(18)式中,則得到Hyperheno′n映射系M統(tǒng)的二階伴隨方程:
在仿真試驗中,觀測數(shù)據(jù)的采集過程如下:首先設(shè)定物理參數(shù)的真實值,讓Hyperheno′n映射自由演化,在經(jīng)歷暫態(tài)過程之后任取混沌狀態(tài)的一點作為初值,并以此時刻作為初始時刻;其次任其迭代演化至第N步,可獲得Hyperheno′n映射的離散混沌狀態(tài)序列xk(k=0,1,2,···,N).在本試驗中,初始條件設(shè)定為混沌狀態(tài)中一點:x0=0.8350073,y0=-0.8732880和z0=1.6236088,抽取前面4個狀態(tài)量作為觀測數(shù)據(jù)yko,(k=1,2,···,N);然后可在所選狀態(tài)量上疊加高斯隨機數(shù)作為觀測噪聲,噪聲的均值取為0,標準偏差取為σo.因為在這種情況下觀測誤差協(xié)方差矩陣Rk為3×3大小的單位矩陣,所以離散形式的目標泛函表示如下:
根據(jù)伴隨方程、目標泛函及其梯度表達式,利用共軛梯度算法能對非線性映射的未知物理參數(shù)進行最優(yōu)估計,這里稱之為一階變分方法;而根據(jù)一/二階伴隨方程、切線性方程、目標泛函、梯度及Hessian矩陣-向量乘積的表達式,利用NCG算法也能對非線性映射的未知物理參數(shù)進行最優(yōu)估計,稱之為二階變分方法.在兩個方法中,未知參數(shù)的迭代初值都設(shè)定為:θ0=(0,0,0)T.
圖2和圖3中分別給出了在采用一階和二階變分方法時目標泛函對數(shù)值和梯度模對數(shù)值隨迭代次數(shù)的下?降圖.?理論上當目標泛函關(guān)于未知參數(shù)的梯度模時,未知參數(shù)向量取最優(yōu)估計值.但本文中為了減小計算代價,采用通常的迭代收斂標準,即前后兩次?迭代的?目?標泛函梯?度模變化小于預(yù)先指定的閾值此時的計算誤差小于機器誤差,可以忽略.從圖2和圖3中可以看出,就目標泛函及其梯度值的變化而言,二階變分方法的下降速度明顯快于一階變分方法.由表1進一步可知,前者達到收斂只需要經(jīng)歷8次迭代,而后者需要72次迭代;而且二階變分方法達到收斂后的目標泛函值和梯度值比一階變分方法都要小數(shù)個量級.圖4和圖5表示的是物理參數(shù)a和b的估計值在最小化過程中隨迭代次數(shù)的變化以及逼近真實值的程度,從圖中易知,采用二階變分方法時估計值收斂到真值的速度明顯高于一階變分方法.分析和總結(jié)圖2—5中所示的試驗結(jié)果,可得如下結(jié)論:由于二階離散變分方法能夠同時利用目標泛函梯度信息和準確的Hessian矩陣信息,因此目標泛函能夠快速下降,進而達到收斂.表1中給出了在無觀測噪聲情況下一/二階離散變分方法對Hyperhen′on映射系統(tǒng)中未知物理參數(shù)進行估計的比較結(jié)果,容易得出結(jié)論:二階離散變分方法用非常少的迭代次數(shù)獲得了更高的估計精度.隨著觀測誤差增加,雖然未知參數(shù)的估計精度下降,但在相同誤差水平情況下二階變分方法無論是在目標泛函迭代收斂速度還是在參數(shù)估計精度上都優(yōu)于一階變分方法,因篇幅限制具體情況不再贅述.
在本部分中以二維拋物映射系統(tǒng)為例,進一步說明利用二階離散變分方法估計非線性映射中未知參數(shù)的有效性.二維拋物映射可用離散方程表示如下:
圖2 Hyperhen′on映射目標泛函值變化圖
圖3 Hyperhen′on映射泛函梯度值變化圖
圖4 Hyperhen′on映射參數(shù)a辨識曲線圖
圖5 Hyperhen′on映射參數(shù)b辨識曲線圖
表1 無觀測噪聲時Hyperhen′on映射的參數(shù)估計結(jié)果
(37)式表示的是一個二維雙參數(shù)離散迭代映射,兩個控制參數(shù)a和b必須為正實數(shù).a具有全局動態(tài)幅度調(diào)節(jié)功能,稱之為幅度調(diào)節(jié)參數(shù);而b具有動力學(xué)特性調(diào)節(jié)功能,稱之為特性調(diào)節(jié)參數(shù)[20].當a=1.00和b=1.72時,二維拋物映射系統(tǒng)呈現(xiàn)混沌狀態(tài).設(shè)定狀態(tài)向量x=(x,y)T和未知物理參數(shù)向量θ=(a,b)T,目標是在獲得x的觀測序列yko=(xok,yok)T的情況下,估計未知參數(shù)向量θ.
首先通過比較(37)式與 (1)式,易知M(xk,θ)=[ay2k,byk(1-xk)]T,其一階變分為
根據(jù)(11)式可得切線性方程組:
其次引入與xk對應(yīng)的伴隨向量序列λk=(λk,pk)T,并將 xk,yko,λk,θ,(DMxk)T和 (DθM)T等的表達式代入(10)式中的第一和第三式,同時考慮到伴隨初始條件,進而得到二維拋物映射系統(tǒng)的伴隨方程和目標泛函梯度公式,分別如(40)和(41)式所示:
中,則得到二維拋物映射的二階伴隨方程:
在仿真試驗中,觀測數(shù)據(jù)的生成采用與4.1部分中相同的方法.模式初始狀態(tài)量(x0,y0)取為(0.0147653,-0.1697946),其是二維拋物映射混沌系統(tǒng)狀態(tài)中一點.考慮到觀測誤差協(xié)方差矩陣是2×2大小的單位矩陣,目標泛函表示如下:
根據(jù)伴隨方程、目標泛函及其梯度表達式,利用共軛梯度算法能對二維拋物映射的未知物理參數(shù)進行最優(yōu)估計,稱為一階變分方法;而根據(jù)一/二階伴隨方程、切線性方程、目標泛函、梯度及Hessian矩陣-向量乘積的表達式,利用NCG算法能夠?qū)ΧS拋物映射的未知物理參數(shù)進行最優(yōu)估計,稱為二階變分方法.兩種方法中未知參數(shù)向量的迭代初值都設(shè)定為θ0=(0,0)T.最小化迭代收斂采用與4.1部分中相同的準則.
圖6 二維拋物映射目標泛函值變化圖
圖7 二維拋物映射泛函梯度值變化圖
在二維拋物映射系統(tǒng)的未知物理參數(shù)估計中,圖6和7中分別給出了在采用一階和二階變分方法時目標泛函對數(shù)值和梯度模對數(shù)值隨迭代次數(shù)的下降圖.從圖6和7中易知:就目標泛函及其梯度模對數(shù)值的變化而言,二階變分方法的下降速度明顯快于一階變分方法.由表2進一步可得:前者達到收斂只經(jīng)歷17次迭代,而后者需要85次迭代.同時,二階變分方法達到收斂后的目標泛函值和梯度模值比一階變分方法都要小數(shù)個量級.圖8和圖9展示的是二維拋物映射系統(tǒng)中參數(shù)a和b的估計值隨迭代次數(shù)的變化及其逼近真實值的程度,從圖中易知:在利用二階變分方法估計兩個未知物理參數(shù)時,收斂速度都快于一階變分方法,而對參數(shù)a估計的收斂速度差異更明顯.表2中分別給出了在無觀測噪聲情況下一/二階離散變分方法對二維拋物映射系統(tǒng)中未知參數(shù)進行估計的比較結(jié)果:一階變分方法對參數(shù)a和b的估計精度分別是小數(shù)點后3位和4位,而二階變分方法對兩個參數(shù)的估計精度都能達到小數(shù)點后7位.根據(jù)上面的分析歸納出如下結(jié)論:因為二階離散變分方法能夠同時利用準確的目標泛函梯度和Hessian矩陣信息,所以能夠用更少的迭代次數(shù)獲得更高的估計精度,明顯優(yōu)于一階變分方法.上述結(jié)論在不同觀測誤差水平情況下仍然成立,由于篇幅所限不再贅述.
圖8 二維拋物映射參數(shù)a的辨識曲線圖
圖9 二維拋物映射參數(shù)b的辨識曲線圖
表2 無觀測噪聲時二維拋物映射的參數(shù)估計結(jié)果
本文提出一種估計非線性映射未知物理參數(shù)向量的二階離散變分方法.首先針對xk+1=M(xk,θ)形式的非線性離散混沌系統(tǒng),通過拉格朗日乘子將系統(tǒng)的狀態(tài)控制方程引入到目標泛函中,將映射方程約束最優(yōu)化問題轉(zhuǎn)換為無約束最優(yōu)化問題;其次利用變分方法給出了一/二階伴隨方程、目標泛函梯度和Hessian矩陣-向量乘的顯式表達式;然后設(shè)計了估計非線性映射未知參數(shù)的新算法,并以此對Hyperhen′on映射和二維拋物映射中的未知參數(shù)進行了估計.數(shù)值仿真結(jié)果表明了二階離散變分方法的有效性和優(yōu)點:因為能夠同時利用準確的目標泛函梯度信息和與Hessian矩陣相關(guān)的二階信息,所以二階變分方法在無論有無觀測噪聲的情況下都能夠用較少的迭代次數(shù)獲得更高的估計精度,明顯優(yōu)于一階變分方法.該方法可推廣應(yīng)用于其他類型非線性映射和混沌系統(tǒng)的未知物理參數(shù)估計中.
[1]Liu FC,Liang X M 2005 Acta Phys.Sin.54 4584(in Chinese)[劉福才,梁曉明2005物理學(xué)報54 4584]
[2]Wang X Y,Wu X J2006 Acta Phys.Sin.55 605(in Chinese)[王興元,武相軍2006物理學(xué)報55 605]
[3]Wu Z Q,Tan F X,Wang SX 2006 Acta Phys.Sin.55 1651(in Chinese)[吳忠強,譚拂曉,王紹仙2006物理學(xué)報55 1651]
[4]Maybhate A,Amritkar RE 1999 Phys.Rev.E 59 284
[5]Parlitz U 1996 Phys.Rev.Lett.76 1232
[6]He M F,Mu Y M,Zhao L Z 2000 Acta Phys.Sin.49 830(in Chinese)[賀明峰,穆云明,趙立中2000物理學(xué)報49 830]
[7]Guan X P,Peng H P,Li L X,Wang Y Q 2001 Acta Phys.Sin.50 26(in Chinese)[關(guān)新平,彭海朋,李麗香,王益群2001物理學(xué)報50 26]
[8]Cao X Q,Song JQ,Zhang W M,Zhao J,Zhang L L 2011 Acta Phys.Sin.60 070511(in Chinese)[曹小群,宋君強,張衛(wèi)民,趙軍,張理論2011物理學(xué)報60 070511]
[9]Gao F,Tong H Q 2006 Acta Phys.Sin.55 577(in Chinese)[高飛,童恒慶2006物理學(xué)報55 577]
[10]Li L X,Peng H P,Yang Y X,Wang X D 2007 Acta Phys.Sin.56 51(in Chinese)[李麗香,彭海朋,楊義先,王向東2007物理學(xué)報56 51]
[11]Wang JY,Huang D X 2008 Acta Phys.Sin.57 2755(in Chinese)[王鈞炎,黃德先2008物理學(xué)報57 2755]
[12]Long W,Jiao JJ 2012 Acta Phys.Sin.61 110507(in Chinese)[龍文,焦建軍2012物理學(xué)報61 110507]
[13]Cao X Q,Zhang W M,Song JQ,Zhu X Q,Zhao J 2012 Acta Phys.Sin.61 020507(in Chinese)[曹小群,張衛(wèi)民,宋君強,朱小謙,趙軍2012物理學(xué)報61 020507]
[14]Huang SX,Wu R S 2001 Mathematical Physics Problems in Atmosphere Science(Beijing:Meteorology Press)(in Chinese)[黃思訓(xùn),伍榮生2001大氣科學(xué)中的數(shù)學(xué)物理問題(北京:氣象出版社)]
[15]He JH 2008 Int.J.Modern.Phys.B.22 3487
[16]He JH 2000 Appl.Math.Mech.21 797
[17]He JH 2001 Int.J.Nonlin.Sci.Numer.2 309
[18]He JH,Lee EWM 2009 Phys.Lett.A 373 1644
[19]Chen ED 2005 Ph.D.Dissertation(Dalian:Dalian University of Science and Technology)(in Chinese)[陳爾東2005博士學(xué)位論文(大連:大連理工大學(xué))]
[20]Meng JD,Bao B C,Xu Q 2011 Acta Phys.Sin.60 010504(in Chinese)[孟繼德,包伯成,徐強2011物理學(xué)報60 010504]
[21]Kelley C T 1999 Iterative Methods for Optimization(Philadelphia:SIAM)p95
[22]Nocedal J,Wright SJ2006 Numerical Optimization(2nd Ed.)(Verlag,Berlin,Heidelberg,New York:Springer)p198