秦廷华

最优控制问题广泛存在于科学研究和实际工程各领域,因为解析解通常难以找到,所以许多学者致力于研究处理该问题的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拟谱法在各种研究中可更换为其他直接法,如此修改后的算法可以发挥相应直接法特有的长处.本文算法捕捉触碰到约束界限的两类不光滑点,后续研究将考虑如何捕捉其他位置的弱间断点和间断点.