一种最优控制问题直接离散求解的网格重构方法与流程

本发明属于动态优化控制技术领域，涉及一种最优控制问题直接离散求解的网格重构方法。

背景技术：

化工反应过程、优化设计、动态过程系统参数估计、生产过程工作点切换和过程系统优化控制等过程中，存在诸多复杂最优控制问题。这类问题一般含有微分和代数方程，以及众多的轨线等式和不等式约束。对于复杂最优控制问题，传统上采用间接解法，求解的一般步骤为：(1)首先将原系统扩展成hamilton系统；(2)推导出一阶最优性必要条件；(3)获得数学上的两点边值问题，进而采用相应的方法进行求解得到精确的解析解。但是对于一些大规模的并且含有约束条件的最优控制问题时，在求解过程中，则需要引入更多的乘子函数和互补条件，这也是间接法不适合求解大规模、复杂度较高的动态优化问题的原因。

迭代动态规划算法是基于动态规划计算量大，求解效率不高等缺点提出的，改进的地方主要包括两大部分：网格离散和区域缩减。网格离散的思想是：首先将整个系统从时间和空间两个角度进行离散，使每个时间段的状态变量被离散成一系列网格，根据贝尔曼最优性原理的分阶段的特性，从最后一个时间段开始，根据可行控制变量进行计算，再依次从后致前进行计算。区域缩减的思想是：首先将整个时域离散成一组较为粗糙的网格从而获得精度略低的解，而后以这些低精度的解为分界点，进行区域缩减，那么缩减后的时间区间自然变得更小从而可以使用动态规划进行迭代求解。虽然迭代动态规划方法是一种具有全局收敛性的优化方法，并且不需要计算梯度等信息，但是它也有它明显的不足之处：由于离散后的时间段很细才会得到满意的最优解致使计算成本大幅度提高。

随着计算机和计算技术的发展，解决复杂最优控制问题的方法往往是采用直接法。直接法，顾名思义，与间接法相反，不需要求解最优性条件，而是直奔动态优化问题本身，直接对性能指标进行寻优。直接法的原理是将动态问题的整个时间域中的控制变量以及状态变量进行离散化，这样便能够将原动态优化问题转化成一个大规模的非线性规划问题。离散方法中采用较多的是控制变量参数化(controlvariableparameterization，简称cvp)方法，该方法中时间网格划分得是否恰当决定了求解效率和最优控制轨迹的逼近精度。划分地疏，对最优控制轨迹的逼近程度不高；划分地密，虽然确保了离散精度满足要求，但大大增加了非线性规划问题的维数和计算成本。为此，本发明提出一种最优控制问题直接离散求解的网格重构方法以解决控制向量参数化方法逼近精度和计算时间之间的矛盾等问题。

技术实现要素：

本发明的目的是针对以往直接离散求解方法存在的不足，提供一种最优控制问题直接离散求解的网格重构方法。

本发明基于伪wigner-ville分布分析快速而有效的网格重构策略，用于复杂最优控制问题直接离散求解。首先给定网格进行第一次优化迭代，快速获得控制变量的大致轨迹。然后通过伪wigner-ville分布分析得出不同时间网格节点瞬时频率变化对性能指标的影响，籍此对原有网格节点进行重构，包括对时间节点的消除、细化。并且结合变时间节点控制向量参数化方法的思想，将瞬时频率为极大值时对应的时间节点作为待优化参数，与控制变量一同进行求解优化，从而找到准确的最优时间切换点。

包括以下步骤：(1)该方法首先根据经验给定初始时间网格数量，初始时间网格数较少(通常为5～20个)，目的为了快速获得控制变量的大致轨迹，为下次迭代提供良好的初值条件，提高求解速度。将最优控制问题根据初始网格数量离散化为非线性规划问题。(2)求解非线性规划问题，得到当前时间网格下的最优的控制参数和目标函数值。(3)将得到的控制参数根据时间顺序形成控制参数轨迹，对该轨迹进行伪wigner-ville分析变换，得到控制参数轨迹的时频曲线。(4)根据规则对时间网格进行合并和细分，根据重要时间切换点的判断准则，筛选出需要优化的时间切换点。(5)判断是否满足终止条件，若满足，则输出重构后的时间网格；否则，转步骤(3)。

本发明具体包括以下步骤：

步骤(1)：采用基于控制向量参数化的直接法将式(1.1)～(1.7)所示的复杂最优控制问题离散为非线性规划问题，初始迭代次数l＝1。

x(t0)＝x0(1.5)；

t0≤t≤tf(1.7)；

其中j是目标函数，由终值项φ0[x(tf),tf]和积分项组成，f[t,x(t),u(t)]是微分方程，t为时间，t0和tf分别为起始时间和终端时间。u(t)是(n×1)维的控制变量，n大于1则表示u(t)为矢量形式的控制变量；x(t)是(m×1)维的状态变量，m大于1则表示x(t)为矢量形式的状态变量；m，n分别表示状态变量和控制变量的维数。是x(t)的导数。φr[x(tf)]表示关于终端状态x(tf)的末端性能函数，lr[t,x(t),u(t)]是时间t、状态变量x(t)和控制变量u(t)的复合函数，表示m1个等式约束，表示m2个不等式约束。x(t0)＝x0是状态变量在起始时间的初值，u和分别是u(t)的下界和上界。

对于式(1.1)～(1.7)所示的复杂最优控制问题，首先将整个控制时域[t0，tf]均匀分成式(1.8)所示的n个时间子区间[ti-1,ti](i＝1,2,…,n)，分段数n通常取值为5～20，目的是为了快速获得控制变量的大体轨迹，同时也为第二次迭代提供良好的初值。

t0＜t1＜…＜tn-1＜tn＝tf(1.8)；

其中离散化后的时间节点ti(i＝1,2,…,n)都是固定值。在整个控制时域上，n维的控制变量u(t)的第j维分量u^j(t)就可以由各个时间子区间的值近似表示为式(1.9)：

其中，为控制变量u^j(t)在子区间[ti-1,ti)的值，t[ti-1,ti)为单位开关函数，被定义为式(1.10)：

在各个时间子区间内的控制变量均由一系列基函数的线性组合来近似，即式(1.11)：

其中，为qi,j阶基函数，为线性组合系数，称为控制参数。对函数采用分段常量(分段零次多项式)逼近策略，则有k＝qi,j＝1，以及式(1.11)可简化为式(1.12)：

是控制参数的简化形式，是式(1.1)～(1.7)需要求取的优化变量。

由此便可将一无限维的动态优化问题(1.1)～(1.7)转换为一个含有限维的控制参数的非线性规划问题(1.13)～(1.19)，即:

x(t0)＝x0(1.17)；

t0≤t≤tf(1.19)；

步骤(2)：使用非线性规划求解技术求解式(1.13)～(1.19)的非线性问题，得到当前时间网格下的最优的控制参数和目标函数值j，令obj1＝j，其中非线性规划求解技术为已有成熟技术。obj1表示当前时间网格下的最小目标函数值。

步骤(3):将得到的控制参数根据时间顺序形成控制参数轨迹，对该轨迹进行伪wigner-ville分布分析变换，得到控制参数轨迹的时频曲线ω^j(t)，以及在时间点的ti-1瞬时频率其中ω^j(t)表示第j维控制参数轨迹对应的频率，i＝1,2,…,n，j＝1,2,…,n。伪wigner-ville分布变换技术是已有成熟技术。

步骤(4):进行时间网格精细化重构策略，重新划分时间网格以确保求解精度，本部分包括三个子步骤：

子步骤1):对于相邻的时间网格节点，找出其中瞬时频率变化比较小的进行网格合并。对于相邻的时间网格[ti,ti+1]与[ti+1,ti+2]，如果控制变量u^j(t)在时间点ti、ti+1和ti+2对应的瞬时频率满足方程(1.20)，则将时间网格[ti,ti+1]与[ti+1,ti+2]合并成一个网格[ti,ti+2]。

其中，为消除系数，取值为0.1～0.5hz。为控制变量的变化系数，其取值满足如式(1.21)所示的规则：

合并后时间节点和瞬时频率都将被重新标记为式(1.22)：

合并后的时间网格数量记为n'。

子步骤2):根据瞬时频率的大小对时间网格进行细分，将网格均匀的划分为δk个小区间，δk的数目由以下式(1.23)所示经验规则决定。

其中称为细化系数，取值由下式(1.24)决定：

其中为控制变量的走势系数，取值由式(1.25)决定。

时间网格[ti-1,ti]被细化分为δk个小区间后的时间节点将被重新标记如式(1.26)：

式中为再次重新标记后时间节点。经过以上网格细分后，新的时间网格数量记为令将各个网格时间点表示为

子步骤3):对重要时间切换点实现精确定位，获得每个网格的最佳划分方式。根据步骤(3)经过控制参数轨迹进行伪wigner-ville分布变换得到的时频曲线，找到控制变量u^j(t)对应的瞬时频率ω^j存在极大值的时间点，若该时间点在时间网格[ti-1,ti]内，则将ti作为一个待优化变量。令新的网格数将控制参数和所有满足以上条件的ti作为待优化变量，使用非线性规划求解技术重新求解式(1.13-1.19)的非线性问题，得到当前时间网格下的最优的参数和ti，以及新的目标函数值j，令obj2＝j。obj2表示当前时间网格下目标函数最小值。

步骤(5):如果l＝lmax或者(obj1-obj2)/obj1≤tol，则迭代结束，得到的ti值为满足要求的最佳时间网格节点，得到的n为满足要求额最佳网格节点数量，得到的表示在此网格划分下最佳的控制参数，也就代表最佳的控制量tol表示用户允许误差，一般取值在10^-4-10^-8之间，lmax表示设定最大迭代次数，取值一般小于等于5。否则，如果(obj1-obj2)/obj1＞tol并且l＜lmax，则置l＝l+1，然后转到步骤(3)。

本发明不仅可以降低最优控制问题直接离散求解变量规模，而且计算量较小，并且迭代次数少，较少的参数获得高质量的解。该方法适用大规模复杂动态优化问题的在线优化：可以更有效地重构时间网格，找到准确的时间切换点，不仅计算成本低，而且计算精度更出色。

具体实施方式

一种最优控制问题直接离散求解的网格重构方法，包括如下步骤：

步骤(1)：采用基于控制向量参数化的直接法将式(1.1)～(1.7)所示的复杂最优控制问题离散为非线性规划问题，初始迭代次数l＝1。

x(t0)＝x0(1.5)；

t0≤t≤tf(1.7)；

其中j是目标函数，由终值项φ0[x(tf),tf]和积分项组成，f[t,x(t),u(t)]是微分方程，t为时间，t0和tf分别为起始时间和终端时间。u(t)是(n×1)维的控制变量，n大于1则表示u(t)为矢量形式的控制变量；x(t)是(m×1)维的状态变量，m大于1则表示x(t)为矢量形式的状态变量；m，n分别表示状态变量和控制变量的维数。是x(t)的导数。φr[x(tf)]表示关于终端状态x(tf)的末端性能函数，lr[t,x(t),u(t)]是时间t、状态变量x(t)和控制变量u(t)的复合函数，表示m1个等式约束，表示m2个不等式约束。x(t0)＝x0是状态变量在起始时间的初值，u和分别是u(t)的下界和上界。

对于式(1.1)～(1.7)所示的复杂最优控制问题，首先将整个控制时域[t0，tf]均匀分成式(1.8)所示的n个时间子区间[ti-1,ti](i＝1,2,…,n)，分段数n通常取值为5～20，目的是为了快速获得控制变量的大体轨迹，同时也为第二次迭代提供良好的初值。

t0＜t1＜…＜tn-1＜tn＝tf(1.8)；

其中离散化后的时间节点ti(i＝1,2,…,n)都是固定值。在整个控制时域上，n维的控制变量u(t)的第j维分量u^j(t)就可以由各个时间子区间的值近似表示为式(1.9)：

其中，为控制变量u^j(t)在子区间[ti-1,ti)的值，t[ti-1,ti)为单位开关函数，被定义为式(1.10)：

在各个时间子区间内的控制变量均由一系列基函数的线性组合来近似，即式(1.11)：

其中，为qi,j阶基函数，为线性组合系数，称为控制参数。对函数采用分段常量(分段零次多项式)逼近策略，则有k＝qi,j＝1，以及式(1.11)可简化为式(1.12)：

是控制参数的简化形式，是式(1.1)～(1.7)需要求取的优化变量。

由此便可将一无限维的动态优化问题(1.1)～(1.7)转换为一个含有限维的控制参数的非线性规划问题(1.13)～(1.19)，即:

x(t0)＝x0(1.17)；

t0≤t≤tf(1.19)；

步骤(2)：使用非线性规划求解技术求解式(1.13)～(1.19)的非线性问题，得到当前时间网格下的最优的控制参数和目标函数值j，令obj1＝j，其中非线性规划求解技术为已有成熟技术。obj1表示当前时间网格下的最小目标函数值。

步骤(3):将得到的控制参数根据时间顺序形成控制参数轨迹，对该轨迹进行伪wigner-ville分布分析变换，得到控制参数轨迹的时频曲线ω^j(t)，以及在时间点的ti-1瞬时频率其中ω^j(t)表示第j维控制参数轨迹对应的频率，i＝1,2,…,n，j＝1,2,…,n。伪wigner-ville分布变换技术是已有成熟技术。

步骤(4):进行时间网格精细化重构策略，重新划分时间网格以确保求解精度，本部分包括三个子步骤：

子步骤1):对于相邻的时间网格节点，找出其中瞬时频率变化比较小的进行网格合并。对于相邻的时间网格[ti,ti+1]与[ti+1,ti+2]，如果控制变量u^j(t)在时间点ti、ti+1和ti+2对应的瞬时频率满足方程(1.20)，则将时间网格[ti,ti+1]与[ti+1,ti+2]合并成一个网格[ti,ti+2]。

其中，为消除系数，取值为0.1～0.5hz。为控制变量的变化系数，其取值满足如式(1.21)所示的规则：

合并后时间节点和瞬时频率都将被重新标记为式(1.22)：

合并后的时间网格数量记为n'。

子步骤2):根据瞬时频率的大小对时间网格进行细分，将网格均匀的划分为δk个小区间，δk的数目由以下经验规则决定。

其中称为细化系数，取值由下式(1.24)决定：

其中为控制变量的走势系数，取值由式(1.25)决定。

时间网格[ti-1,ti]被细化分为δk个小区间后的时间节点将被重新标记如式(1.26)：

式中为再次重新标记后时间节点。经过以上网格细分后，新的时间网格数量记为令将各个网格时间点表示为

子步骤3):对重要时间切换点实现精确定位，获得每个网格的最佳划分方式。根据步骤(3)经过控制参数轨迹进行伪wigner-ville分布变换得到的时频曲线，找到控制变量u^j(t)对应的瞬时频率ω^j存在极大值的时间点，若该时间点在时间网格[ti-1,ti]内，则将ti作为一个待优化变量。令新的网格数将控制参数和所有满足以上条件的ti作为待优化变量，使用非线性规划求解技术重新求解式(1.13-1.19)的非线性问题，得到当前时间网格下的最优的参数和ti，以及新的目标函数值j，令obj2＝j。obj2表示当前时间网格下目标函数最小值。

步骤(5):如果l＝lmax或者(obj1-obj2)/obj1≤tol，则迭代结束，得到的ti值为满足要求的最佳时间网格节点，得到的n为满足要求额最佳网格节点数量，得到的表示在此网格划分下最佳的控制参数，也就代表最佳的控制量tol表示用户允许误差，一般取值在10^-4-10^-8之间，lmax表示设定最大迭代次数，取值一般小于等于5。否则，如果(obj1-obj2)/obj1＞tol并且l＜lmax，则置l＝l+1，然后转到步骤(3)。