混沌優化算法求解Kepler方程

楊玲玲， 馬 良

（上海理工大學 管理學院，上海 200093）

1 問題的提出

1619年Kepler首次推導出了Kepler方程，該方程在物理、數學領域，特別是在經典天體力學研究中都起著重要作用.

在天體力學和軌道計算中，求解Kepler方程

是經常遇到的問題，即在給定平近點角M（0≤M≤π）和偏心率e后去求偏近點角x［1］.在橢圓軌道情況下，0≤e≤1.由于方程中含有超越函數sinx，因此，Kepler方程是個超越方程，這意味著x無法用代數方法求出，只能用數值分析或級數求解，這就要求Kepler方程在實際求解中應盡可能地接近真值.

超越方程常用迭代方法求解，目前應用較多的數值迭代方法是不動點迭代法或牛頓迭代法，但前者收斂速度慢，后者在選取初值上要求較高，而且每次都要計算導數值［2］.

在實際應用中，由于對上千萬年的天體軌道作積分需要反復求解Kepler方程，計算量大且耗時，因此，研究快速求解Kepler方程的方法具有重要的意義［3］.本文試圖使用一種智能優化算法——混沌優化算法求解Kepler方程，探索從新的角度求解此類問題.

2 現有解法

盡管Kepler方程的提出已近4個世紀，但因其重要性而一直是研究的重要課題.每10～20年就會提出一個新的求解方法.Ng［4］在1979年提出一立方收斂的方法.

1983年Danby和Burkardt提出了具有立方收斂性的Halley方法.

同時，他們基于該方法提出了四次方和五次方的收斂性方法［5］.Colwell，Battin和 Vallado分別提出的方法也是以迭代或級數展開為基礎的經典方法［6］.

3 算法設計

混沌運動具有遍歷性、隨機性及規律性等特點，它能在一定范圍內，按其自身規律不重復地遍歷所有狀態［7］.李兵和蔣慰孫［8］選用Logistic映射做算法

式中，μ為控制參量，取μ＝4.

設0≤a0≤1，n＝0，1，2，…，當μ＝4時，系統完全處于混沌狀態，運用載波的方式將混沌狀態引入至優化變量中，同時將遍歷范圍放大到優化變量的取值范圍，然后利用混沌變量進行搜索［8］.文獻［8］中介紹了算法的基本步驟.

針對Kepler方程設計下面的算法.

Step 1 等價轉換式（1）.

將求解式（1）的問題轉化為求x，使得目標函數F（x）＝（x－esinx－M）2的值最小.

Step 2 確定x的取值范圍.

已知0≤x≤2π，進一步來說，在文獻［9］中得出結論

因為f′（x）＝1－ecosx≥0恒成立，0≤e≤1，函數f（x）單調遞增，同時

故可得出M≤x≤M＋e.f（x）的圖像如圖1所示.

圖1 f（x）的圖形Fig.1 Graph of f（x）

Step 3 初始化.

設置k＝1，a0，最大迭代次數n.

Step 4 根據式（4），產生混沌變量ai.通過式（6），利用載波方式將混沌變量放大到相應優化變量在Step 2中所得到的取值范圍.

Step 5 如果F（xi）＜F＊，那么F＊＝F（xi），x＊＝xi；否則，放棄xi.

Step 6 如經若干次搜索后保持不變，則按式（7）進行第二次載波.

式中，x＊為當前最優解；α為較小的調節常數，本文中取小于1的常數.

Step 7 如果F［x（k）］＜F＊，那么F＊＝F［x（k）］，x＊＝x（k）；否則，放棄x（k）.

Step 8 當達到最大迭代次數時，終止搜索，輸出最優解x＊，F＊.

不難估算，該算法的時間復雜度為O（n），n為最大迭代次數.

4 數值測試

為檢驗求解效果，進行了一系列數值測試計算，并與已有方法進行比較.實驗所用的硬件為Pentium（R）4CPU 2.80GHz，256MB內存的微機，軟件為 Windows XP操作系統和Matlab7.0運行環境.

在現有文獻中，大多數計算實例是在M＝0和e＝1附近進行的，因為，在這一點有f′＝f″＝0，在運用傳統方法時容易導致解法不收斂，本文首先求解了這類情況，現列出部分結果.表1和表2分別列出了迭代次數為1 000和10 000時，在臨界情況（零平近點角和等于1的偏心率，即M＝0，e＝1）附近的測試結果，混沌初值a0＝0.201 7.

表1 最大迭代次數為1 000時的結果Tab.1 Results of n＝1000

從表1和表2中可以看出，迭代次數越大，求解精度越好.由于算法的時間復雜度較低，算法的運行時間基本可以忽略不計.

文獻［10］中指出，用改進的牛頓法求解時可能收斂至一個或多個錯誤的值，例如，當M＝1.347，e＝0.66，迭 代 序 列 在 0.448 902，0.354 721，0.469 946和0.362 751這4個值上循環，而此時的正確解是1.958 11［10］.所以，本文也對該情況進行了求解，求解結果為

表2 最大迭代次數為10 000時的結果Tab.2 Results of n＝10 000

本文還測試了文獻［3］和文獻［10］中的一些算例，結果如表3所示.

表3 部分算例測試結果Tab.3 Part of test results

測試結果表明，混沌優化算法對于Kepler方程求解的有效性，該算法較低的時間復雜度提升了運算的效率，且免除了多次求導的復雜運算，避免了在臨界狀態容易不收斂的缺點.

5 結 論

常見的求解Kepler方程的方法是簡單迭代法、牛頓法等經典方法，本文運用智能優化方法——混沌優化算法求解該問題.算法的時間復雜度較低，具有較高的運算效率，同時避免了傳統方法中多次求導等復雜計算，且取得了滿意的結果.對于Kepler方程，本文所給出的混沌優化算法是一種簡單而又快速地獲得有效解的途徑.

［1］周慶林，黃天衣.Kepler方程的解法［J］.天文學報，1988，29（1）：106－112.

［2］束雄英，李紅.Kepler方程的一種迭代加速［J］.航空計算技術，2005，35（1）：41－44.

［3］王玉詔，鐘雙英，孫威，等.Kepler方程的六階迭代解法［J］.江西科學，2009，27（6）：790－792.

［4］Colwell P.Solving Kepler’s equation over three centuries［M］.Richmond：Willmann－Bell，1993.

［5］Danby J M A，Burkardt T M.The solution of Kepler’s equation［J］.Celestial Mechanics，1983，31（2）：317－328.

［6］Davis J J.Sequential solution to Kepler’equation［J］.Celest Mech Dyn Astr，2010，108（7）：59－72.

［7］苗東升.系統科學大學講稿［M］.北京：中國人民大學出版社，2007.

［8］李兵，蔣慰孫.混沌優化方法及其應用［J］.控制理論與應用，1997，14（4）：613－615.

［9］Smith G R.A simple，efficient starting value for the iterative solution of Kepler’s equation［J］.Celestial Mechanics，1979，19（2）：163－166.

［10］岳錦海，黃天衣.橢圓Kepler方程求解中的非線性現象［J］.南京大學學報（自然科學版），1998，34（1）：21－28.