秦廷華

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

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

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

文獻[11]研究燃料最優(yōu)控制問題,先用擬譜法解光滑的最優(yōu)控制問題,所得結(jié)果用于構(gòu)造和估計協(xié)態(tài),并將該協(xié)態(tài)做為間接法解題的初始猜測,然后用同倫法將光滑問題逐漸轉(zhuǎn)變?yōu)椴还饣脑瓎栴},在此過程中用間接法解題.文獻[11]的思路還被文獻[12]用于研究時間最優(yōu)控制問題.文獻[13]也采用類似的思路,研究的問題模型更為特殊,但是所得的協(xié)態(tài)正好為零,這為同倫法和間接法解這類特殊的題帶來便利.

文獻[10]用同倫法將光滑的最優(yōu)控制問題逐漸轉(zhuǎn)變?yōu)椴还饣脑瓎栴},產(chǎn)生的若干最優(yōu)控制問題都用直接法中的自適應控制參數(shù)化方法(Adaptive control parametrization method)求解.文獻[14]研究月面上升最優(yōu)控制問題,通過調(diào)整問題本身的參數(shù)得到易于求解的問題，然后用擬譜法和同倫法求解.

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

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

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

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

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

1 問題描述和離散

1.1 問題描述

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

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

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

1.2 問題離散

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

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

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

2 收斂性分析

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

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

2.1 離散問題的收斂性

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

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

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

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

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

定理2.如果假設2 成立,問題B 的最優(yōu)解是

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

注意到問題B 中已經(jīng)假設函數(shù)h關于變元x和u連續(xù)可微.于是,根據(jù)假設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 的最優(yōu)解,由上述證明可知

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

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

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

2.2 弱間斷點的捕捉

對于問題B 和B(ch)的向量函數(shù)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 成立的情況下,對于存在數(shù)列使得以下三個等式成立.

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

證明.令

3 誤差指示量和算法

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

定義1.令L,K ∈N,K≥L≥2,對于網(wǎng)格

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

的CGL 點被算法1 的步驟4 用來近似tj,i,如此得到新網(wǎng)格點替換原有網(wǎng)格點,實現(xiàn)網(wǎng)格自適應調(diào)整.

算法1.結(jié)合同倫法的自適應擬譜方法

步驟1.設置j=K1=1.Nmin=4,NInitial=8,=33,δ=0,θ=0.1,TolEI==Tol,Mj=NInitial,根據(jù)算例情況輸入?yún)?shù)

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

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

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

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

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

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

如果Kj=1,則設置

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

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

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

4 數(shù)值算例

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

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

例1.有解析解的弱間斷最優(yōu)控制問題[8].

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

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

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

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

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

例3.關于生產(chǎn)和維護的Bang-Bang 最優(yōu)控制問題[19].

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

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