混合整数最优控制问题松弛控制解的一种改进取整策略
An Improved Rounding Strategy for Relaxed Control Solutions of Mixed-Integer Optimal Control Problems
摘要: 混合整数最优控制问题(Mixed-Integer Optimal Control Problem, MIOCP)因其包含整数变量而更加复杂,但更能贴近实际应用需求。Sager等人(Math. Program. Vol. 133, 2012)提出了一种松弛–取整策略,将MICOP问题凸松弛为经典最优控制,再对凸松弛的最优解进行取整近似(Sum Up Rounding, SUR),得到原问题的近似解,并证明了近似误差为时间步长的同阶无穷小。然而,该近似误差估计的同阶无穷小的系数项是时间区间总长度的指数函数,当控制问题的时间区间较大时,这个误差可能会非常大。针对这一问题,本文对SUR策略进行改进,提出一个新的控制取整策略,证明了新控制策略的收敛性,并通过数值例子验证了本文的策略显著提高了MIOCP的求解精度。
Abstract: Mixed-Integer Optimal Control Problems (MIOCP) are more complex due to the inclusion of integer variables but are more aligned with practical application needs. Sager et al. (Math. Program. Vol. 133, 2012) proposed a relaxation-rounding strategy that convexly relaxes the MIOCP to a classical optimal control problem and then approximates the optimal solution of the convex relaxation (Sum Up Rounding, SUR) to obtain an approximate solution to the original problem, proving that the approximation error is of the same order as the infinitesimal of the time step. However, the coefficient of this infinitesimal error estimate is an exponential function of the total length of the time interval, and when the time interval of the control problem is large, this error can be very significant. To address this issue, this paper improves the SUR strategy and proposes a new control rounding strategy, proving the convergence of the new control strategy and verifying through numerical examples that it significantly improves the solution accuracy of MIOCP.
文章引用:谷云飞. 混合整数最优控制问题松弛控制解的一种改进取整策略[J]. 应用数学进展, 2025, 14(3): 107-116. https://doi.org/10.12677/aam.2025.143096

1. 引言

混合整数最优控制问题(MIOCP)是一类特殊的最优控制问题,其中控制函数不仅包含连续变量,还包含整数变量。具体来说,这类问题的目标是最小化(或最大化)一个性能指标,同时满足一组动态约束(通常是常微分方程)和路径约束。整数变量的引入使得问题更加复杂,但也更贴近实际应用中的需求。例如,在运输系统中,换挡操作可以被视为一个整数变量;在化工过程中,阀门的开关状态也可以用整数变量表示。

MIOCP问题的求解一直是一个挑战性的问题,受到了学者们的广泛关注。常见的求解方法可归分为两大类:启发性算法和确定性算法。前者主要是针对一些具体问题,提出一些智能优化算法进行数值求解。例如,Wang等人[1] [2]在研究一个代谢网络的结构辨识的MIOCP问题时,提出了一种竞争粒子群算法进行求解;Ahmadi等人[3]在考虑卡车行驶时间不确定性和软时间窗口的情况下,优化卡车和无人机的配送路线,以最小化总成本;Wang等人[4]提出了一种新的效率松弛模型,用于优化精馏塔的塔板数。该模型通过插值函数将整数塔板数松弛为连续变量,有效避免了非整数解或局部最优解的问题,显著提高了优化效率并降低了计算成本;Liu等人[5]提出了一种混合智能优化方法求解MIOCP问题,分别用量子退火算法求解其中的整数型控制变量和双精英量子蚁群算法求解连续型控制变量。启发性便于编程实现,但无法保证所得的解是否最优。

在确定性算法方法,Barton等人提出了一系列基于参数灵敏度函数的梯度下降算法框架[5] [6]。此类算法存在的主要问题在于整数值切换时刻,状态灵敏度轨迹会出现跳跃,因而梯度函数不具有光滑性,算法的收敛性较差。Sager等人[7] [8]提出了一种松弛–取整算法,将原问题转化为纯连续控制变量的最优控制问题,然后构造一种称为Sum Up Rounding (SUR)的取整策略,将松弛控制策略取整转化为整数值控制策略,并给出了近似误差估计,证明这一取整策略的误差是控制离散化时间步长的同阶无穷小。这种方法的优势在于构建了MIOCP问题与经典的连续最优控制问题的关系,但该项研究中仍然存在不足之处,即误差估计过于保守,误差估计中的同阶无穷小的系数随时间区间总长度呈指数增长。如果所考虑的MIOCP问题的时间区域较大,这个误差可能变得非常大。为此,本文在Sager等人所提出松弛控制策略基础上,提出一种新的取整策略,提高MIOCP的松弛–取整框架的求解精度。

2. 混合整数最优控制问题的松弛–取整相关结果

混合整数最优控制问题

考虑如下最优控制问题:

(P1) min u( t ),v( t ) J=Φ( x( T ) )

s.t.{ x ˙ =f( x( t ),u( t ),v( t ) ),t[ 0, t f ] x( 0 )= x 0 (1)

c( x( t ),u( t ) )0,t[ 0, t f ]. (2)

其中, x( t ) R n x 是状态变量, u( t )U R n u 是连续型控制变量, v( t )Ω:={ v 1 , v 2 ,, v n ω } R n v 是离散型控制变量, J 是目标性能, c: R n x × R n u R n c Φ: R n x R 均为连续可微函数。满足约束(1)和(2)的三元组 ( x( ),u( ),v( ) ) 称为问题(P1)的一个可行轨迹,此时,称 x( ) 为可行状态轨迹, ( u( ),v( ) ) 为可行控制轨迹或可行控制对。

由于问题(P1)含有连续型和离散型两类决策变量,故该问题被称为混合整数最优控制问题。对于离散型控制变量 v( t ) 的每个取值 v i ,通过引入二元控制函数 ω i ( ):[ 0, t f ]{ 0,1 } ,则问题(P1)可等价表示如下:

(P2) min u( t ),ω( t ) J=Φ( x( T ) )

s.t.{ x ˙ = i=1 n ω f( x( t ),u( t ), v i ) ω i ( t ) ,t[ 0, t f ], x( 0 )= x 0 (3)

i=1 n ω ω i ( t ) =1 (4)

ω i ( t ){ 0,1 },i=1,2,, n ω , (5)

c( x( t ),u( t ) )0,t[ 0, t f ].  (6)

进一步地,将二元控制函数 ω i ( t ){ 0,1 } 松弛为 α i ( t )[ 0,1 ] ,得到如下凸松弛最优控制问题:

(P3) min u( t ),α( t ) J=Φ( x( T ) )

s.t.{ x ˙ = i=1 n ω f( x( t ),u( t ), v i ) α i (t) ,t[ 0, t f ], x( 0 )= x 0 (7)

i=1 n ω α i ( t ) =1 (8)

α i ( t )[ 0,1 ],i=1,2,, n ω , (9)

c( x( t ),u( t ) )0,t[ 0, t f ].  (10)

问题(P3)是仅含有连续型控制变量的经典最优控制问题,可运用经典最优控制的数值算法求解,设其

最优解为 ( u * ( t ), α * ( t ) ) ,其中 α * ( t )=( α 1 * ( t ), α 2 * ( t ),, α n ω * ( t ) ) 。当然,问题(P3)与问题(P1)不再等价。为

此,需要回到如下问题:

(I) 如何由问题(P3)的最优解构造原问题(P1)的(近似)最优解;

(II) 所构造的(近似)最优解是否收敛于原问题(P1)的最优解。

在文献[7]中,Sager等人给出了一种称为Sum Up Rounding (SUR)的取整策略,由问题(P3)的最优解构造原问题(P1)的近似最优解,并证明了这种构造方法的收敛性。下面对这一方面进行叙述。

考虑给定的可测函数 α j :[ 0, t f ][ 0,1 ] ,其中 j=1,, n ω ,以及一个时间分割 0= t 0 < t 1 << t m = t f 。定义 Δ t i := t i+1 t i ,并用 Δt 表示最大时间步长,即:

Δt:= max i=0,,m1 Δ t i .

然后,定义函数 ω( ):[ 0, t f ] { 0,1 } n ω 如下:

ω j ( t )= p j,i ,t[ t i , t i+1 )i=0,1,,m1,j=1,2,, n ω . (11)

其中,

p j,i ={ 1,  p ^ j,i p ^ k,i kjjk: p ^ j,i > p ^ k,i 0,  (12)

p ^ j,i = 0 t i+1 α j ( τ )dτ k=0 i1 p j,k Δ t k (13)

显然,所构造的函数 ω( t )=( ω 1 ( t ), ω 2 ( t ),, ω nω ( t ) ) 满足问题(P2)的约束(4)和约束(5)。

关于 ω( ) α( ) 满足如下关系:

引理1 [7]。任给可测函数 α( ):[ 0, t f ] [ 0,1 ] n ω ,通过公式(11)~(13)所构造的函数 ω( ) ,则有如下不等式成立:

0 t α( τ )ω( τ )dτ 0.5Δt.

引理1说明公式(11)~(13)的取整策略所得的整数值控制函数 ω( ) Δt0 时可以任意逼近实值控制函数 α( ) 。这两个控制函数所对应的约束动力系统(3)的解轨迹之间还满足如下收敛性质。

引理2 [7]。设 ( x,α, u * ) 是松弛问题(P3)的一个可行轨迹, ω( ) α( ) 通过公式(11)~(13)在给定时间分割上所构造, y( ) 是在控制对 ( ω, u * ) 下求解初值问题(3)得到的状态轨迹。假设对于给定的可测控制 u * 和所有 v i Ω ,存在常数 L,M R + ,使得函数 f( x( t ), u * ( t ), v i ) 满足:

f( y( t ), u * ( t ), v i )f( x( t ), u * ( t ), v i ) L y( t )x( t )

f( x( ), u * ( ), v i ) M.

其中, 分别表示欧式范数和 L 范数。那么, ( y,ω, u * )( ) 是问题(P2)的一个可行轨迹,并且对于所有 t[ 0, t f ] ,有

y( t )x( t ) ( M+Ct )c( n ω ) e Lt Δt (14)

其中, c( n ω )= n ω 1

注:引理2对所有可行轨迹 ( x,α,u ) 都成立,因此也对全局和局部最优轨迹成立。

3. 一种新的松弛–取整算法与其收敛性

从引理2可以看出,在相同的连续控制 u( ) 下,任给 α( ) ,通过SUR取整策略构造的 ω( ) ,两者对应的系统(3)的状态轨迹 x( t ) y( t ) 之间误差随着时间 t 呈指数增长。因此,为了保证误差足够小, Δt 要充分小。换而言之,SUR策略中的时间分割节点数随着考虑问题的时间规模呈指数增长。另一方面,从引理1可以看到,SUR策略能够实现 α( ) 与其取整策略 ω( ) 的累计误差代数和不超过 0.5Δt 。因此,这种取整策略偏向于限制控制输入误差,而非系统状态轨迹的误差。然而,从实际运用的角度,一般更关注系统状态之间的误差。为此,本文从系统状态误差最小化的角度出发,提出一种新的控制输入取整策略,具体如下:

给定控制对 ( α, u * ) ,设约束系统(3)在某一微分方程数值解框架下进行求解,得到的数值解为 { x( t i ) } i=1 m ,对应的步长为 Δt 。记 e j 为第 j 个分量取1的 n ω 维单位向量。定义 ω ˜ ( ):[ 0, t f ] { 0,1 } n ω 如下:

ω ˜ j ( t )= p ˜ j,i ,t[ t i , t i+1 ), i=0,1,,m1,j=1,2,, n ω . (15)

其中 p ˜ i =( p ˜ i 1 , p ˜ i 2 ,, p ˜ j n ω ) 由如下公式给出:

p ˜ i =argmin{ x( t i )z( t i ;e ) |e{ e 1 , e 2 ,, e n ω } } (16)

这里, z( t i ,e ) 为以 z( t i1 ) 为前一个节点值,在给定微分方程数值解框架下得到的第 i 个节点的计算值,且 z( t 0 )= x 0 。我们将公式(15)和公式(16)所定义的取整策略称为轨迹极小化跟踪取整策略,简记为MTR策略。

可以证明, ( u * ( ), ω ˜ ( ) ) 对应的系统(3)的状态轨迹 z( t ) Δt0 时也是收敛于松弛控制 ( u * ( ),ω( ) ) 对应的状态轨迹 x( t ) 。下面证明这一收敛性结果。

定理1. 设 ( ω( ), u * ( ) ) 是问题(P2)的可行控制对, ω ˜ ( ) 由公式(15)和(16)所定义, z( t ) ( ω ˜ ( ), u * ( ) ) 对应的系统(3)的状态轨迹。假设对于给定的可测控制 u * 和所有 v i Ω ,存在常数 L,M R + ,使得函数 f( x( t ), u * ( t ), v i ) 满足:

f( y( t ), u * ( t ), v i )f( x( t ), u * ( t ), v i ) L y( t )x( t )

f( x( t ), u * ( t ), v i ) M.

那么,对于所有 t[ 0, t f ] ,有

z( t )x( t ) ( ( M+C t i )c( n ω ) e L t i +2M )Δt (17)

证明:首先,由公式(16)可知,在 t= t i 处有

z( t i )x( t i ) y( t i )x( t i ) ( M+C t i )c( n ω ) e L t i Δt.

又对于 t[ t i1 , t i ) ,有

z( t )=z( t i1 )+ t i1 t f( z( t ), u * ( t ), ω ˜ ( t ) )dt

x( t )=x( t i1 )+ t i1 t f( x( t ), u * ( t ),α( t ) )dt

故有

z( t )x( t ) = z( t i1 )x( t i1 )+ t i1 t f( z( t ), u * ( t ), ω ˜ ( t ) )f( x( t ), u * ( t ),α( t ) )dt z( t i1 )x( t i1 ) + t i1 t f( z( t ), u * ( t ), ω ˜ ( t ) )f( x( t ), u * ( t ),α( t ) )dt z( t i1 )x( t i1 ) + t i1 t f( z( t ), u * ( t ), ω ˜ ( t ) )dt + t i1 t f( x( t ), u * ( t ),α( t ) )dt z( t i1 )x( t i1 ) +2M( t t i1 ) ( M+C t i 1 )c( n ω ) e L t i1 Δt+2MΔt=( ( M+C t i )c( n ω ) e L t i +2M )Δt

证毕。

注:对比引理2和定理1的误差界估计,似乎本文构造的取整策略所对应的系统(3)的状态轨迹与松弛控制的状态轨迹之间的误差要比文献[7]的误差界更大。但误差界估计仅仅给出的是一种上界估计,而本文的估计结果是基于文献[7]进行推导的,因此所估计误差界有可能是偏于保守的。这一点将从下一节的数值实验中可以得到验证。

4. 数值结果

为了验证本文所提出的改进策略的有效性,我们通过MATLAB进行了数值实验。由于本文关注的是取整策略的有效性,因此我们假设所给的策略都是问题(P3)的可行策略。正如前面所指出的,如果我们的取整策略对于任意可行的松弛控制输入都具有收敛性,那么对最优的松弛控制输入自然也是收敛的。因此,下文的数值算例中,仅需比较不同取整策略所得的整数值控制输入对应的系统(3)状态轨迹近似解与连续型的松弛控制对应的状态轨迹之间的误差。

例1考虑如下初值问题[7]

(IVP1) { x ˙ 1 = x 1 α 1 +( x 1 + x 2 ) α 2 +( x 1 x 2 ) α 3 , x ˙ 2 =( x 1 +2 x 3 ) α 1 +( x 1 2 x 3 ) α 2 +( x 1 + x 2 ) α 3 ,t[ 0,1 ] x ˙ 3 = x 1 2 + x 2 2 , x 1 ( 0 )=0.5, x 2 ( 0 )=0.5, x 3 ( 0 )=0.

这里,取连续控制输入 α 1 ( t )=sint, α 2 ( t )=0.40.4sint, α 2 ( t )=0.60.6sint 。显然, α i ( t )[ 0,1 ] i=1 3 α i ( t ) =1 。利用文献[7]的SUR策略对 α( ) 进行取整所得控制策略记为 ω( ) ,利用本文的取整策略(公式(14)和公式(15))所得控制策略记为 ω ˜ ( ) 。将这三个控制策略分别代入初值问题(IVP1),分别通过欧拉法求解对应的初值问题,所得系统状态轨迹的如图1~3所示。这里,取步长 Δt= 10 3 。可以看出,利用SUR策略和本文的MTR策略,所得整数值控制策略对应的系统状态轨迹都非常接近连续控制输入 α( ) 对应系统状态轨迹。

例2考虑如下初值问题:

(IVP2) { x ˙ 1 =sin( x 1 ) α 1 +cos( x 1 + x 2 ) α 2 + x 1 α 3 , x ˙ 2 =sin( x 1 +2 x 2 ) α 1 +( x 1 2 x 2 ) α 2 +( x 1 + x 2 ) α 3 ,t[ 0,10 ] x ˙ 3 =sin( x 1 + x 2 x 1 2 + x 2 2 ), x 1 ( 0 )=0.5, x 2 ( 0 )=0.5, x 3 ( 0 )=0.

这里,取连续控制输入 α 1 ( t )=sin t 2 , α 2 ( t )=0.40.4sin t 2 , α 2 ( t )=0.60.6sin t 2 。同样地,将松弛控制 α( ) 以及由SUR策略和MTR策略分别所得的 ω( ) ω ˜ ( ) 代入初值问题(IVP2),通过欧拉法求解对应的初值问题,所得系统状态轨迹的如图4~6所示。这里,依然取步长 Δt= 10 3 。可以看出,与SUR策略所得整数值控制策略相比较,本文的MTR策略所得整数值控制策略对应的系统状态轨迹更接近连续控制输入 α( ) 对应系统状态轨迹。

Figure 1. Comparison of the calculated values of state component x1 for initial value problem (IVP1) corresponding to two different strategies

1. 两种不同策略对应初值问题(IVP1)的状态分量x1的计算值比较

Figure 2. Comparison of the calculated values of state component x2 for initial value problem (IVP1) corresponding to two different strategies

2. 两种不同策略对应初值问题(IVP1)的状态分量x2的计算值比较

Figure 3. Comparison of the calculated values of state component x3 for initial value problem (IVP1) corresponding to two different strategies

3. 两种不同策略对应初值问题(IVP1)的状态分量x3的计算值比较

Figure 4. Comparison of the calculated values of state component x1 for initial value problem (IVP2) corresponding to two different strategies

4. 两种不同策略对应初值问题(IVP2)的状态分量x1的计算值比较

Figure 5. Comparison of the calculated values of state component x2 for initial value problem (IVP2) corresponding to two different strategies

5. 两种不同策略对应初值问题(IVP2)的状态分量x2的计算值比较

Figure 6. Comparison of the calculated values of state component x3 for initial value problem (IVP2) corresponding to two different strategies

6. 两种不同策略对应初值问题(IVP2)的状态分量x3的计算值比较

5. 结论

本文针对混合整数最优控制问题在松弛–取整数值求解框架下需要解决的松弛控制取整问题,提出了一种新的误差估计策略,并通过具体的数值例子验证新策略在误差上的精度上优于传统的SUR策略。然而,本文在进行误差限的理论估计推导过程中,仍然沿用了文献[7]的方法,因此所得误差估计与该文献的结果是同阶的。从数值结果来看,这个误差估计结果可能是比较保守的。在下一步的工作中,我们将尝试用新的方法来对本文所提出的方法的误差进行重新估计。

参考文献

[1] Wang, J., Xu, H., Teo, K.L., Sun, J. and Ye, J. (2020) Mixed-Integer Minimax Dynamic Optimization for Structure Identification of Glycerol Metabolic Network. Applied Mathematical Modelling, 82, 503-520.
https://doi.org/10.1016/j.apm.2020.01.042
[2] Liu, Z. and Li, S. (2021) A Quantum Computing Based Numerical Method for Solving Mixed-Integer Optimal Control Problems. Journal of Systems Science and Complexity, 34, 2428-2469.
https://doi.org/10.1007/s11424-020-9278-6
[3] Ahmadi, M. and Tavakkoli-Moghaddam, R. (2024) The Multi-Visit Drone-Assisted Routing Problem with Soft Time Windows and Stochastic Truck Travel Times. Transportation Research Part B: Methodological, 183, 102894.
[4] Wang, Y. and Jia, S. (2022) A New Efficiency Relaxation Model for the Optimization of Distillation Column Stages. Chemical Engineering Science, 264, 117073.
[5] Galán, S. and Barton, P.I. (2002) Parametric Sensitivity Functions for Hybrid Discrete/Continuous Systems. Mathematical Programming, 93, 347-388.
[6] Barton, P.I. and Lee, C.K. (2002) Modeling, Simulation, Sensitivity Analysis, and Optimization of Hybrid Systems. ACM Transactions on Modeling and Computer Simulation, 12, 256-289.
https://doi.org/10.1145/643120.643122
[7] Sager, S., Bock, H.G. and Diehl, M. (2010) The Integer Approximation Error in Mixed-Integer Optimal Control. Mathematical Programming, 133, 1-23.
https://doi.org/10.1007/s10107-010-0405-3
[8] Sager, S. and Zeile, C. (2020) On Mixed-Integer Optimal Control with Constrained Total Variation of the Integer Control. Computational Optimization and Applications, 78, 575-623.
https://doi.org/10.1007/s10589-020-00244-5