求解非光滑最優控制問題的自適應網格優化

陳 琦，王中原，常思江，舒敬榮

（1．南京理工大學能源與動力工程學院，江蘇南京210094；2．陸軍軍官學院二系，安徽合肥230031）

求解非光滑最優控制問題的自適應網格優化

陳 琦1，王中原1，常思江1，舒敬榮2

（1．南京理工大學能源與動力工程學院，江蘇南京210094；2．陸軍軍官學院二系，安徽合肥230031）

針對傳統直接配點法在求解非光滑最優控制問題時存在離散誤差大、精度低的問題，提出了一種自適應直接配點法。利用局部分段插值多項式逼近最優解，將最優控制問題離散為非線性規劃問題，并給出了離散誤差估計方法，根據離散誤差的大小確定區間內節點的加密量，提出了自適應網格優化算法，利用該算法將大部分節點配置在非光滑區域以降低離散誤差。最后通過仿真算例將所提算法與傳統直接配點法和文獻中的擬譜自適應算法分別進行比較，驗證了所提算法的高精度和有效性。

最優控制問題；非光滑；直接配點法；網格優化；自適應算法

0 引 言

近年來，最優控制問題的求解引起了人們的廣泛關注，在實際的應用中，許多最優控制問題往往是不光滑的，例如狀態受約束的非對稱航天器最優控制問題［1］、月球著陸問題［2］、橋式吊車最優控制問題［3］等。在文獻［1］研究的非對稱航天器最優控制問題中，原本光滑的最優解在對其中一個角速度施加不等式約束后，相應的控制力矩會出現不光滑現象。文獻［2］研究了在月球著陸過程中的燃料最優控制問題，其數值仿真表明，在存在推力約束的條件下，控制量會存在切換現象。文獻［3］研究了橋式吊車的最優控制問題，由于對控制量施加了上下界約束，導致最終的狀態量不光滑且控制量不連續。

針對非光滑最優控制問題，國內外學者進行了大量的研究，并提出了較多的方法。文獻［4－5］提出了hp自適應偽譜法，設計了一種多重區間分解策略，通過判斷離散誤差的分布情況，來確定區間分割點位置及子區間內的配點數。文獻［6］采用hp自適應偽譜法研究了無人作戰飛機的軌跡規劃問題。文獻［7］同樣以hp自適應偽譜法為基礎對N脈沖軌道進行了優化設計。為了進一步提高hp自適應偽譜法對非光滑問題的求解效率，文獻［8］提出了分段點最佳化思想，提高了分段點收斂至不光滑點的速度。文獻［9］以Chebyshev偽譜法為基礎，引入方波脈沖函數，提出了一種復合Chebyshev有限差分法，該算法具有使用方便、精度高等優勢。文獻［10］提出了一種分段擬譜法，通過設置分點，將時間區域在非光滑區域進行劃分，獲得了較高的求解精度。文獻［11］以Chebyshev偽譜法為基礎提出了一種自適應算法，通過計算相鄰2個配置點上數值解導數之差，來確定新分點的位置，從而提高了算法的精度。上述方法針對求解非光滑問題所做的改進均以偽譜法為基礎，偽譜法作為一種全局離散方法，存在諸如微分矩陣計算量大、離散后的非線性規劃問題（nonlinear programming problem，NLP）的雅克比矩陣較為稠密等缺點［12－13］，且隨著節點數目的增加，這些缺點會更加突出。相比于偽譜法，直接配點法作為一種局部離散方法，其本身無需計算微分矩陣，主要采用分段插值多項式近似狀態變量和控制變量，因此離散后的NLP規模很小，雅克比矩陣也更為稀疏，求解效率很高［12－13］，近年來得到了廣泛的應用［14－17］。但該方法的計算精度低于偽譜法，且精度降低在求解非光滑最優控制問題時更為突出。若能對該方法進行改進，在保證計算效率高的同時，盡可能地降低離散誤差，使其具有和偽譜法相當的計算精度，那么對于擴大該算法的使用范圍將非常有意義。

有關直接配點法的改進，文獻［18－19］提出了密度函數法，以密度函數來量化狀態量和控制量的劇烈變化程度，并以此為基準重新分配離散節點。在密度函數選取合適的條件下，該方法可較快地將非光滑區域進行節點加密。然而針對不同的問題，選取合適的密度函數并不容易。文獻［20－21］均以小波多尺度分析的方法構建了相應的節點調整策略。該方法主要以小波系數的幅值來確定非光滑區域，基于小波系數和二分點的對應關系，確定下一次迭代步所使用的節點并進行序列優化。與以上方法類似，本文也是通過調整離散節點的位置來提高求解精度。但在非光滑區域位置的確定方式上，不同于密度函數或是小波系數的方法，本文給出了一種離散誤差的估計公式，然后根據該公式得到的離散誤差的大小來確定非光滑區域的位置，構造了一種自適應網格更新策略，使節點能夠根據離散誤差的大小自動調整，從而避免了密度函數的挑選以及較為復雜的小波推導過程，在保證計算效率的基礎上，提高了直接配點法對非光滑最優控制問題的求解精度，以期為求解非光滑最優控制問題提供一種新的思路。最后通過數值仿真并與文獻中其他算法的對比驗證了所提算法的準確性和有效性。

1 直接配點法

考慮一般Bolza型最優控制問題：

式中，狀態變量x（t）∈Rn；控制變量u（t）∈Rm；f為動力學函數；e為終端約束函數；h為路徑約束函數；t0和tf分別為初始和終端時間；Φ為終端性能指標函數；L為性能指標被積函數。

在區間［t0，tf］內，將連續的時間等分為N段，三階Hermite－Simpson方法通過構造3次插值多項式近似子區間內的狀態變量，記第k個子區間為［tk，tk＋1］（k＝0，1，…，N－1），該區間內的3次插值多項式可表示為

邊界條件為

記hk＝tk＋1－tk，fk＝f［x（tk），u（tk），tk］，fk＋1＝f［x（tk＋1），u（tk＋1），tk＋1］，Hermite－Simpson方法還需要利用區間中點處的信息，為此令利用邊界條件及區間中點處的導數信息，可得

從式（4）解出a（k）0，a（k）1，a（k）2和a（k）3并帶入到式（2）

將tk＋1處的邊值條件帶入式（5）

記xk＝x（tk），xk＋1＝x（tk＋1），由式（6）可得

式（7）即為Hermite－Simpson defect向量，該向量構成了離散形式的系統微分方程約束。令，其中，式（1）中性能指標函數中的積分項利用復合Simpson公式求解，這樣連續的最優控制問題式（1）便可在節點處轉換為離散格式為

式中，k＝0，…，N－1。從式（8）可看出，連續的無限維最優控制問題被離散為一般的NLP問題，對于該問題，可采用非線性最優化算法進行求解。此外，在轉換過程中加入了一個松弛因子δ，該因子的引入在一定程度上保證了離散后的NLP問題可行解的集合是非空的。

由于直接配點法將時間區域進行了N等分，因此便產生了均勻分布的離散節點，這些節點不會根據離散信息而改變，當最優控制問題在某些區域內變化劇烈甚至產生不光滑現象時，這種均勻節點的配置方式將會產生較大的離散誤差，如果能在這些區域內適當增加節點的數目，那么可有效地減小離散誤差。

2 網格更新算法

2．1 離散誤差的估計

節點更新算法需要利用離散誤差的信息，然而在一般情況下，最優控制問題的解事先是未知的，因此離散誤差的大小只能采用特殊的方法進行估計，本小節主要給出離散誤差的估計公式。求解式（8）得到每個節點上的狀態量xi及其導數值＝fi。為了充分利用節點上的數值及其導數信息，構造Hermite插值函數（t）近似狀態量x（t），有

式中，Ns為節點個數；αi，βi（i＝1，2，…，Ns）為Hermite插值基函數。每個插值基函數為2 Ns－1次多項式，并滿足如下條件［18］

式中，Ci（t）為3次B－樣條函數的基。

考慮區間［tk，tk＋hk］，根據系統微分方程，可知

式（13）的計算需要準確的狀態量x（t）和控制量u（t），但這些均未知。為此，利用上文的和u）來近似求解x（tk＋hk），有

對式（16）兩邊取絕對值

定義區間［tk，tk＋hk］上的絕對局部離散誤差

式中

δi，k表示第i個狀態量在第k區間內的離散誤差。進而可定義最大相對局部離散誤差

2．2 節點增加量的確定

設εtol為預先設定的誤差閾值，且，那么需對第k個區間內的節點進行加密處理，設增加的節點數為M′k，增加節點后的離散誤差為。為盡可能提高算法效率，令，因此得

式中，p為離散算法的計算精度。對于本文的Hermite－Simpson方法，對應p＝4．0［22］。求解式（21）可得區間k內需增加的節點數為

考慮到M′k為整數，且過大的M′k會使離散后的NLP產生稠密的雅克比矩陣，加劇求解難度，為此對M′k進行如下處理：

式中，M1為預先設定的常數，用于限制M′k的增長幅值。

2．3 自適應網格更新算法步驟

自適應網格更新算法步驟如下。

步驟1 初始化。給定正整數M0，M1，Nmaxiter以及Niter＝1；設置εtol的值，并在區間［t0，tf］內均勻布置M0個節點。

步驟2 求解NLP問題。采用直接配點法在節點上離散最優控制問題式，得到NLP問題，利用內點法（interior point optimization，IPOPT）［23］求解該問題，得到xN和uN。

步驟3 估計離散誤差。根據式（9）和式（13），構造插值函數?x（t）和?u（t）近似xN和uN，結合式（18）～式（20）求解每個區間內的離散誤差。

（2）節點更新次數未達到最大值，即Niter＜。

則根據式（23）計算M′k個節點加入到區間k中，否則不增加節點，循環本步操作以遍歷k值。

3 數值驗證

3．1 能量最優控制問題

為了驗證本文所提算法的準確性，本小節采用該算法求解文獻［24］中具有解析解的非光滑最優控制問題：

自適應直接配點法的計算參數取為M0＝5，M1＝5，＝10，εtol＝1×10－5，δ＝1×10－10。該最優控制問題的解析解為

由解析解可看出，控制量在t＝0．3和t＝0．7處出現不光滑現象。圖1分別將自適應算法和采用均勻節點的傳統直接配點法的數值計算結果與解析結果進行了對比，從中可以明顯地發現，由自適應算法求解出的控制量在邊界點及t＝0．3和t＝0．7點均能很好地逼近解析解，而傳統直接配點法則產生了較大的誤差。圖2為狀態變量的對比結果，可以發現在t＝0．3和t＝0．7附近，自適應算法配置了較多的節點，且節點處的狀態變量的值與解析解十分吻合，而傳統直接配點法則不能較好地逼近解析解。

圖1 控制量的數值解與解析解對比

圖2 狀態量的數值解與解析解的對比

圖3為每步迭代過程中的節點配置情況，從中可以看出，自適應算法在迭代的過程中能夠自動地增加兩端及t＝0．3和t＝0．7附近的節點數目，對比圖1可知，自適應算法所增加節點的位置正是傳統算法不能較好逼近真值的區域，這說明自適應算法能夠自動識別出離散誤差較大的局域，并對其進行加密節點處理，這一措施可以有效地提高求解精度。圖4對比了2種算法的收斂效果，從中可以看出，當節點數目為30時，自適應算法的最大離散誤差為10－4，而傳統算法的離散誤差仍大于3×10－4，若使其誤差接近當10－4，所需節點數為50個，但此時自適應算法的離散誤差已經小于10－5，可見自適應算法的收斂速度要明顯快于傳統算法，并且隨著節點數的增加，這種差異愈發明顯。從以上分析可得出，在同樣節點數目的情況下，通過合理地安排節點的配置可以有效地降低離散誤差。以上的仿真對比結果驗證了本文所提算法在處理非光滑最優控制問題時的準確性，這為今后算法的有效使用提供了依據。

圖3 每步迭代過程中的節點配置

圖4 算法收斂過程

3．2 航天器最優控制問題

選取文獻［1］中的航天器最優控制問題，采用本文的算法進行求解有

式中，ω1（t），ω2（t）和ω3（t）為航天器角速度；I1，I2和I3為航天器轉動慣量；u1，u2和u3為控制變量。計算參數εtol＝1×10－6，其他參數的選取與前文一致。

計算結果如圖5～圖9及表1所示。圖5為狀態量變化曲線，可見求解得到的結果滿足約束條件，同時可看出，對ω1（t）施加約束后，自適應算法將大部分節點配置在靠近約束的區域。圖6～圖8為最優控制量變化曲線，為了驗證本文算法的正確性，圖中將自適應算法和傳統算法求解的控制量分別與文獻［1］中的結果分別進行了對比。從圖6可看出，控制量u1在t＝40s附近產生了不光滑現象，自適應算法與傳統算法所得的結果均能吻合文獻［1］中的結果。圖7和圖8分別為u2和u3變化曲線，從中可以看出，自適應算法所得的結果與文獻［1］中的結果相符，但是傳統算法得到的結果則不能吻合文獻［1］中的結果，產生了較大的偏差。由此驗證了自適應算法的求解精度要高于傳統算法。

圖5 狀態變量變化曲線

圖6 控制變量u1變化曲線

圖7 控制變量u2變化曲線

圖8 控制變量u3變化曲線

圖9 自適應算法和傳統算法計算誤差對比（N＝28）

表1 本文算法與文獻［11］算法對比

圖9為在同樣節點數目的情況下（28個），2種算法計算誤差的對比曲線。可以看出，傳統方法使用了均勻分布的離散節點，且在端點及t＝40s附近產生了較大的離散誤差。而自適應算法采用非均勻的節點配置方式，將大部分節點設置在t＝40s附近，極大地降低了離散誤差，從而驗證了本文算法的高精度。

針對該航天器最優控制問題，文獻［11］以自適應擬譜法進行了求解。為了進一步考察本文所提的自適應直接配點法在每步迭代過程中的性能，表1將本文方法在求解航天器最優控制問題時的每步迭代結果與文獻［11］中的迭代結果進行了對比。從中可看出，本文算法需迭代7次，大于文獻［11］中的5次迭代，但所需的節點數目卻由176個減小為28個，降低了84．1%，節點數目的降低有助于減小離散后的NLP的規模，從而提高優化求解效率。此外，從表1也可看出，在相同迭代次數時，本文方法得到的目標值要比文獻［11］中得到的目標值大，這是由于雖然迭代過程中節點位置在不斷地調整，數量也在不斷地增加，但整體來講其數量還是偏少的（迭代7步后節點數才達到17），位置也不是最優的，而且直接配點法本身的離散精度也要低于擬譜法。因此，在節點數目偏少且位置并非最優的情況下，其得到的目標值便大于文獻［11］中得到的目標值。但是這種情況在后續的迭代中得到了改善，從表1中可看出，從第6步開始，隨著節點數目的增加及節點位置的優化，目標值在不斷減小，最終得到了比文獻［11］更小的目標值，這表明了本文算法在求解精度上具有一定的優勢。

4 結 論

本文對非光滑最優控制問題的數值解法進行了研究，提出了一種自適應直接配點法。通過構造分段多項式近似狀態變量和控制變量，將最優控制問題離散為非線性規劃問題，并給出了離散誤差的估計方法及節點加密量的計算公式，根據離散誤差的大小自適應地配置節點。該算法能夠自動地確定出最優控制問題中的非光滑區域，并對其進行加密處理，能較好地解決傳統直接配點法計算精度有限、收斂速度過慢等問題。此外，通過將數值算例的結果與已有文獻中的擬譜自適應算法進行比較，驗證了所提算法僅需較少的節點便可獲得相當的計算精度。本文算法具有一定的通用性，可為此類問題的求解提供參考。

［1］Jaddu H．Direct solution of nonlinear optimal control problems using quasi linearization and Chebyshev polynomials［J］．Journal of the Franklin Institute，2002，339（4／5）：479－498．

［2］Shamsi M．A modified pseudospectral scheme for accurate solution of bang－bang optimal control problems［J］．Optimal Control Application and Methods，2011，32（6）：668－680．

［3］Ross I M，Fahroo F．A direct method for solving nonsmooth optimal control problems［C］∥Proc．of the International Federation of Automatic Control World Conference，2002：3860－3864．

［4］Darby C L，Hager W W，Rao A V．An hp－adaptive pseudospectral method for solving optimal control problems［J］．Optimal Control Application and Methods，2010，32（4）：476－502．

［5］Darby C L，Hager W W，Rao A V．Direct trajectory optimization using a variable low－order adaptive pseudospectral method［J］．Journal of Spacecraft and Rockets，2011，48（3）：433－445．

［6］Liu H M，Ding D L，Huang C Q，et al．UCAV low observable attacking trajectory planning based on adaptive pseudospectral method［J］．Systems Engineering and Electronics，2013，35（1）：78－84．（劉鶴鳴，丁達理，黃長強，等．基于自適應偽譜法的UCAV低可探測攻擊軌跡規劃研究［J］．系統工程與電子技術，2013，35（1）：78－84．）

［7］Han W H，Yang X，Jiang M Z．Optimal design of N－impulse trajectories based on hp－adaptive pseudospectral method［J］．Systems Engineering and Electronics，2013，35（12）：2566－2571．（韓威華，楊新，姜萌哲．基于hp自適應偽譜法的N脈沖軌道優化設計［J］．系統工程與電子技術，2013，35（12）：2566－2571．）

［8］Liu Y B，Zhu H W，Huang X N，et al．Optimal mesh segmentation algorithm for pseudospectral methods for non－smooth optimal control problems［J］．Systems Engineering and Electronics，2013，35（11）：2396－2399．（劉淵博，朱恒偉，黃小念，等．偽譜法求解非光滑最優控制問題的網格優化［J］．系統工程與電子技術，2013，35（11）：2396－2399．）

［9］Marzban H R，Hoseini S M．A composite Chebyshev finite difference method for nonlinear optimal control problems［J］．Communication in Nonlinear Science and Numerical Simula－

tion，2013，18（6）：1347－1361．

［10］Zhang W．A numerical method for solving nonsmooth optimal control problems［J］．Applied Mathematics Journal of Chinese Universities，2009，24（2）：207－220．（張穩．非光滑最優控制問題的一種數值解法［J］．高校應用數學學報，2009，24（2）：207－220．）

［11］Qin T H，Ma H P．Adaptive algorithm for weakly discontinuous solutions of optimal control problems［J］．Control and Decision，2013，28（2）：313－316．（秦廷華，馬和平．最優控制問題弱間斷解的一個自適應算法［J］．控制與決策，2013，28（2）：313－316．）

［12］Huntington G T，Rao A V．Comparison of global and local collocation methods for optimal control［J］．Journal of Guidance，Control，and Dynamics，2008，31（2）：432－436．

［13］Huntington G T，Rao A V．A Comparison between global and local orthogonal collocation methods for solving optimal control problems［C］∥Proc．of the AIAA American Control Conference，2007：1950－1957．

［14］Herman A L，Conway B A．Direct optimization using collocation based on high－order Gauss－Lobatto quadrature rules［J］．Journal of Guidance，Control，and Dynamics，1996，19（3）：592－599．

［15］Hull D G．Conversion of optimal control problems into parameter optimization problems［J］．Journal of Guidance，Control，and Dynamics，1997，20（1）：57－60．

［16］Tu L H，Yuan J P，Luo J J，et al．Lunar soft landing rapid trajectory optimization using direct collocation method［J］．Chinese Space Science and Technology，2008，（4）：19－24．（涂梁輝，袁建平，羅建軍，等．基于直接配點法的月球軟著陸軌道快速優化［J］．中國空間科學技術，2008，（4）：19－24．）

［17］Geiger B R，Horn J F，Delullo A M，et al．Optimal path planning of UAVs using direct collocation with nonlinear programming［C］∥Proc．of the AIAA Guidance，Navigation，and Control Conference and Exhibit，2006：1－13．

［18］Zhao Y M，Tsiotras P．A density－function based mesh refinement algorithm for solving optimal control problems［C］∥Proc．of the AIAA Infotech and Aerospace Conference，2009．

［19］Zhao Y M，Tsiotras P．Density functions for mesh refinement in numerical optimal control［J］．Journal of Guidance，Control，and Dynamics，2011，34（1）：271－277．

［20］Jain S，Tsiotras P．Trajectory optimization using multiresolution techniques［J］．Journal of Guidance，Control，and Dynamics，2008，31（5）：1424－1436．

［21］Feng Z W，Zhang Q B，Tang Q G，et al．Node adaptive refinement for trajectory optimization based on second－generation wavelets［J］．Journal of Aerospace Power，2013，28（7）：1659－1665．（豐志偉，張青斌，唐乾剛，等．基于二代小波的軌跡優化節點自適應加密［J］．航空動力學報，2013，28（7）：1659－1665．）

［22］Kress R．Numerical analysis［M］．New York：Wiley，1998．

［23］Wachter A，Biegler L T．On the implementation of a primaldual interior point filter line search algorithm for large－scale nonlinear programming［J］．Mathematical Programming，2006，106（1）：25－57．

［24］Benson D A．A Gauss pseudospectral transcription for optimal control［D］．Cambridge：Massachusetts Institute of Technology，2005．

E－mail：qiychan＠126．com

王中原（1958－），男，研究員，博士，主要研究方向為飛行器飛行控制理論與技術。

E－mail：zywang＠njust．edu．cn

常思江（1983－），男，講師，博士，主要研究方向為彈箭飛行與控制技術。

E－mail：ballistics＠126．com

舒敬榮（1974－），男，講師，博士，主要研究方向為外彈道學。

E－mail：shujr1974＠126．com

Adaptive mesh refinement for solving non－smooth optimal control problems

CHEN Qi1，WANG Zhong－yuan1，CHANG Si－jiang1，SHU Jing－rong2
（1．School of Energy and Power Engineering，Nanjing University of Science and Technology，Nanjing 210094，China；2．Department 2，Army Officer Academy of PLA，Hefei 230031，China）

Due to the large discrete errors and low accuracy of the conventional direct collocation method for solving non－smooth optimal control problems，an adaptive direct collocation method is presented．The optimal control problem is transcribed into a nonlinear programming problem by using local piecewise interpolation polynomials to approximate the optimal solution．The estimation method of discrete errors is also presented，and an adaptive mesh refinement algorithm is used to refine the grid by adding nodes to the segments in which the optimal solution is non－smooth，the algorithm is repeated until a user－specified error tolerance is met．Finally，the simulation results demonstrate the utility and efficiency of the proposed method by comparing it with the conventional direct collocation method and the adaptive pseudospectral algorithm respectively．

optimal control problem；non－smooth；direct collocation method；mesh refinement；adaptive algorithm

TP 273．1

A

10．3969／j．issn．1001－506X．2015．06．23

陳 琦（1989－），男，博士研究生，主要研究方向為飛行器軌跡優化、制導與控制。

1001－506X（2015）06－1377－07

2014－07－14；

2014－10－23；網絡優先出版日期：2014－11－20。

網絡優先出版地址：http：／／www．cnki．net／kcms／detail／11．2422．TN．20141120．2110．009．html

國家自然科學基金（11272356）；中國博士后科學基金（2013M541676）資助課題