解最優控制問題結合同倫法的自適應擬譜方法

秦廷華

最優控制問題廣泛存在于科學研究和實際工程各領域,因為解析解通常難以找到,所以許多學者致力于研究處理該問題的3 類數值方法[1?2]:直接法、間接法和混合法.

擬譜方法是直接法的一種,對光滑問題具有指數收斂率[3]是其誘人的優點.因為大量的實際問題有間斷或弱間斷[4]的不光滑解,例如生產和維護最優控制問題,所以有學者關注這一優點對應的反面缺點,即不光滑點妨礙了擬譜方法的快速收斂[5].已經有各種自適應擬譜方法[5?8]可以捕捉不光滑點以改善逼近效果,它們大都依據數值解提供的后驗估計來捕捉不光滑點,然后主要用兩種手段來改善數值解精度,一是設置網格點在不光滑點可能的位置附近,二是依據估計的各區間光滑程度來調整區間內逼近多項式次數.

一些學者致力于同倫法解最優控制問題.文獻[9?10]將同倫法的基本思想簡單解釋為“構造一個與原問題有聯系但是相對容易求解的輔助問題,從求解構造的容易問題出發,通過迭代的方式逐步過渡到原來棘手的問題”.

文獻[11]研究燃料最優控制問題,先用擬譜法解光滑的最優控制問題,所得結果用于構造和估計協態,并將該協態做為間接法解題的初始猜測,然后用同倫法將光滑問題逐漸轉變為不光滑的原問題,在此過程中用間接法解題.文獻[11]的思路還被文獻[12]用于研究時間最優控制問題.文獻[13]也采用類似的思路,研究的問題模型更為特殊,但是所得的協態正好為零,這為同倫法和間接法解這類特殊的題帶來便利.

文獻[10]用同倫法將光滑的最優控制問題逐漸轉變為不光滑的原問題,產生的若干最優控制問題都用直接法中的自適應控制參數化方法(Adaptive control parametrization method)求解.文獻[14]研究月面上升最優控制問題,通過調整問題本身的參數得到易于求解的問題，然后用擬譜法和同倫法求解.

本文的思路來自以下三點:

1)文獻[8]在約束界限內尋找數列收斂到約束界限,在數值解等于數列各元素時求根,利用數值解的收斂性,在根中尋找數列收斂到弱間斷點.本文采用該思路,而且用同倫法思想放寬約束界限,這又帶來兩個好處,一是可以在原問題約束界限之外尋找數列收斂到約束界限,二是放寬約束界限可能增加問題的光滑性,避免直接處理不光滑的原問題.

2)文獻[1]提到:Grigoriev 用同倫法和間接法“先放松對推力幅值的限制進行求解,再慢慢減少最大推力幅值進行求解,用上一步的解作為下一步求解的初值,直到得出符合推力幅值約束的結果”.本文采用同樣做法,但使用的是同倫法和直接法.

3)前述其他同倫法文獻與文獻[14]不同,均調整人為加入的參數,這些參數會改變原問題形式,例如文獻[10]需要構造一個易解的最優控制問題與人為加入的參數一起合并到原問題中,其他文獻往往也要改變目標函數的形式,文獻[14]僅調整問題本身的參數顯然更為簡單.本文將原問題自身的約束界限做為同倫參數予以調整,避免了改變原問題形式.

本文方法的主要思想是:對約束上下界先放松再慢慢收緊,用擬譜法解由此產生的多個最優控制問題,并將上一個最優控制問題的解做為下一個求解的初值,直到得出符合約束上下界的結果;與此同時,用數列收斂到約束上下界,在約束方程等于數列各元素時求根,在這些根中尋找數列收斂到不光滑點,據此實現自適應調整網格點分布和逼近多項式次數.

1 問題描述和離散

1.1 問題描述

本文主要研究Bolza 型最優控制問題(問題B).首先引入若干記號,令Mx,Mu,Me和Mh是正整數,a和b是實數且a

其中,控制函數u(t):是弱間斷[4]函數或Bang-Bang 控制,x(t):是狀態函數,表示x關于自變量t的導數;0e=(0,···,0)T∈已知的函數F:和h:關于變元x和u連續可微.

于是,式(1)～(4)構成的問題B 被修改為式(1)～(3)和式(5)構成的問題B(ch),其中ch即本文的同倫參數,用于調整約束界限.

1.2 問題離散

分割區間(a,b)為K個子區間Ik=(tk?1,tk),a=t0

問題B(ch)可轉化為K段連續最優控制問題,每段用Chebyshev 擬譜方法離散,得到一個非線性規劃問題,即如下形式的離散問題BM(δM,ch)

其中,對于?i ∈{1,···,Mx}和?j ∈{1,···,Mu}有

2 收斂性分析

針對最優控制連續的情況,本節首先證明離散問題的收斂性(定理1),這證明對弱間斷的最優控制顯然也適用,據此設計了第3 節算法1 的步驟2;然后,第2.1 節證明在同倫參數收斂的情況下可以得到原問題最優解(定理2).對于觸及約束界限的弱間斷點,第2.2 節證明可以捕捉這類不光滑點(定理3).在第3 節的算法1 中,步驟3 和步驟4 正是以定理2 和定理3 為依據.

為了論述方便,除了沿用前邊的記號之外,在后文中令N 表示自然數集.

2.1 離散問題的收斂性

假設1.對于?j ∈N,存在Mj=(M1,j,···,MK,j),其中Mk,j ∈N,是遞增序列,k=1,···,K,向量Mj的元素都足夠大;離散問題序列為最優解序列,一致收斂的極限為其中都屬于C[a,b].

定理1.如果假設1 成立,問題B(ch)的最優解是

證明.注意到ch是非負常數組成的向量,問題B(ch)的式(5)可等價寫為h(x(t),u(t))?ch≤0h.此時,式(1)～(3)和(5)構成的問題B(ch)形式上變為式(1)～(4)構成的問題B,而問題B 在文獻[15]定理5.4.1 中已有證明.

不難看出,假設1 和定理1 是文獻[8]中假設1和定理1 的推廣,當ch為零向量時,它們與文獻[8]一致.

假設2.非負常數組成的向量序列滿足問題序列由定理1得到的最優解中存在序列且一致收斂于其中q(t)和u∞(t)都屬于C[a,b].

定理2.如果假設2 成立,問題B 的最優解是

證明.首先證明u∞(t)和x∞(t)是可行解.

注意到問題B 中已經假設函數h關于變元x和u連續可微.于是,根據假設2 和式(5)易知

即u∞(t)和x∞(t)滿足問題B 的式(4).類似地,可證u∞(t)和x∞(t)也滿足問題B 的式(2)和式(3),所以,得證u∞(t)和x∞(t)是問題B 的可行解.

令u?(t)和x?(t)是問題B 的最優解,由上述證明可知

注意到問題B 與問題B(ch)的區別僅在于式(4)與式(5)不同,再由假設2 可以得知,對于?j ∈N,問題B 的最優解u?(t)和x?(t)必然是問題B的可行解,所以有

注意到問題B 中已經假設函數E和F關于變元x和u連續可微.根據假設2,易知

由式(6)～(8)知,J[x?,u?]=J[x∞,u∞].

2.2 弱間斷點的捕捉

對于問題B 和B(ch)的向量函數h(x(t),u(t)),令hi(x(t),u(t))表示其第i個分量,i=1,···,Mh,又令不難看出,C[a,b],相應地有

假設3.問題B 的中,至少存在一個和相應兩個集合滿足

其中,i ∈{1,···,Mh}.

定理3.在假設2 和假設3 成立的情況下,對于存在數列使得以下三個等式成立.

其中,Sj,i是在[a,b]內的根組成的集合,的閉包.

證明.令

3 誤差指示量和算法

采用文獻[8]的誤差指示量,給出定義1,與文獻[8]不同,本文定義1 的網格由式(13)的同倫參數決定.

定義1.令L,K ∈N,K≥L≥2,對于網格

后文算法1 的步驟4 將用誤差指示量EIj度量兩套網格點彼此的某種距離,根據柯西收斂原理,該距離可判斷是否有網格點足夠接近不光滑點.

的CGL 點被算法1 的步驟4 用來近似tj,i,如此得到新網格點替換原有網格點,實現網格自適應調整.

算法1.結合同倫法的自適應擬譜方法

步驟1.設置j=K1=1.Nmin=4,NInitial=8,=33,δ=0,θ=0.1,TolEI==Tol,Mj=NInitial,根據算例情況輸入參數

步驟2.解得數值解和目標值J(j).

否則,改Mj各元素為原值的2 倍,以此為新的Mj,以剛才算出的數值解為初值,重新執行步驟2.

步驟3.如果為零向量,則上次解所得為最終解,程序停止,否則繼續.

步驟4.如果Kj >1,則設置

如果Kj=1,則記錄在各CGL 點的最大值,i=1,···,Mh,以此組成向量設置

對于i=1,···,Mh,適當增加CGL 點,使之與當前使用的CGL 點一起組成集合PCGL,找出PCGL中滿足式(13)的點為新增網格點.用端點a和b與新增網格點一起組成新網格,并得到新網格數Kj+1.

如果Kj=1,則設置

設置EIj=+∞.如果min{Kj+1,Kj}≥2,則用步驟2 的網格和步驟4 產生的網格計算當前指標j對應的誤差指示量EIj,計算方法見定義1.

如果EIj≤TolEI和成立,則設置為零向量.

設置j=j+1.設置Mj為Kj個Nmin組成的向量,轉到步驟2.

4 數值算例

算例用MATLAB 編程,筆記本電腦AMD A4-5000 處理器,1.5 GHz 主頻,8 GB 內存.為了與算法1 進行對比,使用Chebyshev 擬譜方法計算K=1時的問題BM(δ,000h).所有問題都用SNOPT[16]求解,并用文獻[17]的快速變換加快計算速度,使用的參數如表1 所示.

表1 算法1 解全部算例使用的參數Table 1 The parameters of Algorithm 1 in all examples

例1.有解析解的弱間斷最優控制問題[8].

由文獻[8]知,這個最優控制問題的弱間斷點是t=2?ln(4.5)和t=2?ln(2.5),最優目標值為J?=59+14 ln 2?(81/2)ln 3+(25/4)ln 5?14e2≈?69.177535595535176.使用算法1,參數的設置相當于將問題中的0≤u(t)≤2 改為?1≤u(t)≤3.

表2 和表3 是算法1 和Chebyshev 擬譜法解例1 的結果,表4 是三種方法解例1 的結果.從表4可以看出,算法1 的精度優于另兩種方法;算法1 的時間分別是Chebyshev 擬譜方法和文獻[8]方法的26%(11.4/44.59)和59%(11.4/19.39).表4 最后一行數據來自文獻[8],由于使用電腦不同,該數據僅供參考.

例2.考慮一個有解析解的Bang-Bang 最優控制問題[18].

由文獻[18]知,u(t)的間斷點為t=1?ln 2,并且可以推知最優目標值為J?=0.5?2e?1+ln 2≈0.457388298217061.使用算法1,參數的設置相當于問題的|u(t)|≤1 被改為|u(t)|≤2.

表5 和表6 是算法1 和Chebyshev 擬譜法解例2 的結果,表7 是三種方法解例2 的結果.從表7可以看出,算法1 的誤差和耗費的時間均小于另外兩種方法.表7 最后一行數據來自文獻[18],由于使用電腦不同,該數據僅供參考.

例3.關于生產和維護的Bang-Bang 最優控制問題[19].

表2 算法1 解例1 的結果Table 2 The results of Example 1 by Algorithm 1

表3 Chebyshev 擬譜方法解例1 的結果Table 3 The results of Example 1 by the Chebyshev pseudospectral method

表4 三種方法解例1 的結果Table 4 The results of Example 1 by three methods

表5 算法1 解例2 的結果Table 5 The results of Example 2 by Algorithm 1

表6 Chebyshev 擬譜方法解例2 的結果Table 6 The results of Example 2 by the Chebyshev pseudospectral method

表7 三種方法解例2 的結果Table 7 The results of Example 2 by three methods

用5 000 個等距網格點時,文獻[19]方法算出Bang-Bang 控制m(t)間斷點為t=0.65691,最優目標值為J=?26.705,相應的控制函數數值解可見文獻[19]的圖5.1(b).不難看出,該圖與本文算法1 得到的圖1 基本一致.

圖1 表8 中Tol=5×10?2 對應的數值解Fig.1 Numerical solutions for Tol=5×10?2 in Table 8

表8 和表9 是算法1 和Chebyshev 擬譜方法解例3 的結果,表10 是三種方法解例3 的結果.從表10 可以看出,在目標值精度大致相同的情況下,算法1 的計算時間和配置點數明顯少于常用的Chebyshev 擬譜方法,尤其與文獻[19]對比時,算法1 節省了大量點數,使得算法1 處理的問題規模更小,有利于節省計算時間.

表8 算法1 解例3 的結果Table 8 The results of Example 3 by Algorithm 1

表9 Chebyshev 擬譜方法解例3 的結果Table 9 The results of Example 3 by the Chebyshev pseudospectral method

表10 三種方法解例3 的結果Table 10 The results of Example 3 by three methods

5 結論

本文針對弱間斷最優控制問題和Bang-Bang最優控制問題,提出一種結合同倫法的自適應擬譜方法,證明了數值解的收斂性和該方法捕捉弱間斷點的能力,據此設計的算法對含有弱間斷點和Bang-Bang 控制間斷點的三個算例都有良好表現.本文方法主要有以下三個優點:

1)本文結合同倫法與直接法,避免了同倫法與間接法結合時以及單獨使用某一方法時的缺點.

2)本文中同倫法只調整約束界限的值,不需為使用同倫法明顯調整問題形式,而問題形式的明顯調整有時是很難想到的.

3)數值實驗表明,與其他幾種數值方法相比,本文方法在時間和精度上有明顯優勢.

各種直接法有各自的優點,本文的Chebyshev擬譜法在各種研究中可更換為其他直接法,如此修改后的算法可以發揮相應直接法特有的長處.本文算法捕捉觸碰到約束界限的兩類不光滑點,后續研究將考慮如何捕捉其他位置的弱間斷點和間斷點.