CN109491242B - 一种最优控制问题直接离散求解的网格重构方法 - Google Patents

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

Info

Publication number
CN109491242B
CN109491242B CN201811325145.3A CN201811325145A CN109491242B CN 109491242 B CN109491242 B CN 109491242B CN 201811325145 A CN201811325145 A CN 201811325145A CN 109491242 B CN109491242 B CN 109491242B
Authority
CN
China
Prior art keywords
time
grid
control
variable
value
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201811325145.3A
Other languages
English (en)
Other versions
CN109491242A (zh
Inventor
江爱朋
徐炜峰
蒋恩辉
丁强
王浩坤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hangzhou Dianzi University
Original Assignee
Hangzhou Dianzi University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hangzhou Dianzi University filed Critical Hangzhou Dianzi University
Priority to CN201811325145.3A priority Critical patent/CN109491242B/zh
Publication of CN109491242A publication Critical patent/CN109491242A/zh
Application granted granted Critical
Publication of CN109491242B publication Critical patent/CN109491242B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/0205—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric not using a model or a simulator of the controlled system
    • G05B13/024—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric not using a model or a simulator of the controlled system in which a parameter or coefficient is automatically adjusted to optimise the performance

Landscapes

  • Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Evolutionary Computation (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Feedback Control In General (AREA)

Abstract

本发明公开了一种最优控制问题直接离散求解的网格重构方法。现有网格重构方法要么给出的时间数量太大或迭代次数过多导致优化计算非常耗时，要么无法保证离散精度从而使得优化结果不够理想，并且现有方法往往难以快速准确地找到***的结构切换点。本发明方法不仅可以降低复杂最优控制问题直接离散求解变量规模，而且计算量较小，并且迭代次数少，较少的参数获得高质量的解。该方法适用大规模复杂动态优化问题的在线优化。本发明提出的最优控制问题直接离散求解的网格重构方法快速有效，不仅可以在满足精度要求的情况下最大限度地降低离散化非线性规划问题规模，降低迭代次数，而且可以快速准确定位到***结构切换点。

Description

一种最优控制问题直接离散求解的网格重构方法
技术领域
本发明属于动态优化控制技术领域，涉及一种最优控制问题直接离散求解的网格重构方法。
背景技术
化工反应过程、优化设计、动态过程***参数估计、生产过程工作点切换和过程***优化控制等过程中，存在诸多复杂最优控制问题。这类问题一般含有微分和代数方程，以及众多的轨线等式和不等式约束。对于复杂最优控制问题，传统上采用间接解法，求解的一般步骤为：(1)首先将原***扩展成Hamilton***；(2)推导出一阶最优性必要条件；(3)获得数学上的两点边值问题，进而采用相应的方法进行求解得到精确的解析解。但是对于一些大规模的并且含有约束条件的最优控制问题时，在求解过程中，则需要引入更多的乘子函数和互补条件，这也是间接法不适合求解大规模、复杂度较高的动态优化问题的原因。
迭代动态规划算法是基于动态规划计算量大，求解效率不高等缺点提出的，改进的地方主要包括两大部分：网格离散和区域缩减。网格离散的思想是：首先将整个***从时间和空间两个角度进行离散，使每个时间段的状态变量被离散成一系列网格，根据贝尔曼最优性原理的分阶段的特性，从最后一个时间段开始，根据可行控制变量进行计算，再依次从后致前进行计算。区域缩减的思想是：首先将整个时域离散成一组较为粗糙的网格从而获得精度略低的解，而后以这些低精度的解为分界点，进行区域缩减，那么缩减后的时间区间自然变得更小从而可以使用动态规划进行迭代求解。虽然迭代动态规划方法是一种具有全局收敛性的优化方法，并且不需要计算梯度等信息，但是它也有它明显的不足之处：由于离散后的时间段很细才会得到满意的最优解致使计算成本大幅度提高。
随着计算机和计算技术的发展，解决复杂最优控制问题的方法往往是采用直接法。直接法，顾名思义，与间接法相反，不需要求解最优性条件，而是直奔动态优化问题本身，直接对性能指标进行寻优。直接法的原理是将动态问题的整个时间域中的控制变量以及状态变量进行离散化，这样便能够将原动态优化问题转化成一个大规模的非线性规划问题。离散方法中采用较多的是控制变量参数化(Control Variable Parameterization，简称CVP)方法，该方法中时间网格划分得是否恰当决定了求解效率和最优控制轨迹的逼近精度。划分地疏，对最优控制轨迹的逼近程度不高；划分地密，虽然确保了离散精度满足要求，但大大增加了非线性规划问题的维数和计算成本。为此，本发明提出一种最优控制问题直接离散求解的网格重构方法以解决控制向量参数化方法逼近精度和计算时间之间的矛盾等问题。
发明内容
本发明的目的是针对以往直接离散求解方法存在的不足，提供一种最优控制问题直接离散求解的网格重构方法。
本发明基于伪Wigner-Ville分布分析快速而有效的网格重构策略，用于复杂最优控制问题直接离散求解。首先给定网格进行第一次优化迭代，快速获得控制变量的大致轨迹。然后通过伪Wigner-Ville分布分析得出不同时间网格节点瞬时频率变化对性能指标的影响，籍此对原有网格节点进行重构，包括对时间节点的消除、细化。并且结合变时间节点控制向量参数化方法的思想，将瞬时频率为极大值时对应的时间节点作为待优化参数，与控制变量一同进行求解优化，从而找到准确的最优时间切换点。
包括以下步骤：(1)该方法首先根据经验给定初始时间网格数量，初始时间网格数较少(通常为5～20个)，目的为了快速获得控制变量的大致轨迹，为下次迭代提供良好的初值条件，提高求解速度。将最优控制问题根据初始网格数量离散化为非线性规划问题。(2)求解非线性规划问题，得到当前时间网格下的最优的控制参数和目标函数值。(3)将得到的控制参数根据时间顺序形成控制参数轨迹，对该轨迹进行伪Wigner-Ville分析变换，得到控制参数轨迹的时频曲线。(4)根据规则对时间网格进行合并和细分，根据重要时间切换点的判断准则，筛选出需要优化的时间切换点。(5)判断是否满足终止条件，若满足，则输出重构后的时间网格；否则，转步骤(3)。
本发明具体包括以下步骤：
步骤(1)：采用基于控制向量参数化的直接法将式(1.1)～(1.7)所示的复杂最优控制问题离散为非线性规划问题，初始迭代次数l＝1。
Figure BDA0001858529530000031
Figure BDA0001858529530000032
Figure BDA0001858529530000033
Figure BDA0001858529530000034
x(t0)＝x0 (1.5)；
Figure BDA0001858529530000035
t0≤t≤tf (1.7)；
其中J是目标函数，由终值项φ0[x(tf),tf]和积分项
Figure BDA0001858529530000036
组成，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分别表示状态变量和控制变量的维数。
Figure BDA0001858529530000037
是x(t)的导数。φr[x(tf)]表示关于终端状态x(tf)的末端性能函数，Lr[t,x(t),u(t)]是时间t、状态变量x(t)和控制变量u(t)的复合函数，
Figure BDA0001858529530000038
表示m1个等式约束，
Figure BDA0001858529530000039
表示m2个不等式约束。x(t0)＝x0是状态变量在起始时间的初值，u和
Figure BDA00018585295300000310
分别是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维分量uj(t)就可以由各个时间子区间的值近似表示为式(1.9)：
Figure BDA0001858529530000041
其中，
Figure BDA0001858529530000042
为控制变量uj(t)在子区间[ti-1,ti)的值，T[ti-1,ti)为单位开关函数，被定义为式(1.10)：
Figure BDA0001858529530000043
在各个时间子区间内的控制变量
Figure BDA0001858529530000044
均由一系列基函数的线性组合来近似，即式(1.11)：
Figure BDA0001858529530000045
其中，
Figure BDA0001858529530000046
为Qi,j阶基函数，
Figure BDA0001858529530000047
为线性组合系数，称为控制参数。对函数
Figure BDA0001858529530000048
采用分段常量(分段零次多项式)逼近策略，则有k＝Qi,j＝1，以及
Figure BDA0001858529530000049
式(1.11)可简化为式(1.12)：
Figure BDA00018585295300000410
Figure BDA00018585295300000411
是控制参数
Figure BDA00018585295300000412
的简化形式，是式(1.1)～(1.7)需要求取的优化变量。
由此便可将一无限维的动态优化问题(1.1)～(1.7)转换为一个含有限维的控制参数
Figure BDA00018585295300000413
的非线性规划问题(1.13)～(1.19)，即:
Figure BDA00018585295300000414
Figure BDA0001858529530000051
Figure BDA0001858529530000052
Figure BDA0001858529530000053
x(t0)＝x0 (1.17)；
Figure BDA0001858529530000054
t0≤t≤tf (1.19)；
步骤(2)：使用非线性规划求解技术求解式(1.13)～(1.19)的非线性问题，得到当前时间网格下的最优的控制参数
Figure BDA0001858529530000057
和目标函数值J，令Obj1＝J，其中非线性规划求解技术为已有成熟技术。Obj1表示当前时间网格下的最小目标函数值。
步骤(3):将得到的控制参数根据时间顺序形成控制参数轨迹，对该轨迹进行伪Wigner-Ville分布分析变换，得到控制参数轨迹的时频曲线ωj(t)，以及在时间点的ti-1瞬时频率
Figure BDA0001858529530000055
其中ωj(t)表示第j维控制参数轨迹对应的频率，i＝1,2,…,N，j＝1,2,…,n。伪Wigner-Ville分布变换技术是已有成熟技术。
步骤(4):进行时间网格精细化重构策略，重新划分时间网格以确保求解精度，本部分包括三个子步骤：
子步骤1):对于相邻的时间网格节点，找出其中瞬时频率变化比较小的进行网格合并。对于相邻的时间网格[ti,ti+1]与[ti+1,ti+2]，如果控制变量uj(t)在时间点ti、ti+1和ti+2对应的瞬时频率
Figure BDA0001858529530000056
满足方程(1.20)，则将时间网格[ti,ti+1]与[ti+1,ti+2]合并成一个网格[ti,ti+2]。
Figure BDA0001858529530000061
其中，
Figure BDA0001858529530000062
为消除系数，
Figure BDA0001858529530000063
取值为0.1～0.5Hz。
Figure BDA0001858529530000064
为控制变量的变化系数，其取值满足如式(1.21)所示的规则：
Figure BDA0001858529530000065
合并后时间节点和瞬时频率都将被重新标记为式(1.22)：
Figure BDA0001858529530000066
合并后的时间网格数量记为N'。
子步骤2):根据瞬时频率的大小对时间网格进行细分，将网格均匀的划分为Δk个小区间，Δk的数目由以下式(1.23)所示经验规则决定。
Figure BDA0001858529530000067
其中
Figure BDA0001858529530000068
称为细化系数，取值由下式(1.24)决定：
Figure BDA0001858529530000069
其中
Figure BDA00018585295300000610
为控制变量的走势系数，取值由式(1.25)决定。
Figure BDA00018585295300000611
时间网格[ti-1,ti]被细化分为Δk个小区间后的时间节点将被重新标记如式(1.26)：
Figure BDA0001858529530000071
式中
Figure BDA0001858529530000072
为再次重新标记后时间节点。经过以上网格细分后，新的时间网格数量记为
Figure BDA0001858529530000073
令
Figure BDA0001858529530000074
将各个网格时间点表示为
Figure BDA0001858529530000075
子步骤3):对重要时间切换点实现精确定位，获得每个网格的最佳划分方式。根据步骤(3)经过控制参数轨迹进行伪Wigner-Ville分布变换得到的时频曲线，找到控制变量uj(t)对应的瞬时频率ωj存在极大值的时间点，若该时间点在时间网格[ti-1,ti]内，则将ti作为一个待优化变量。令新的网格数
Figure BDA0001858529530000076
将控制参数
Figure BDA0001858529530000077
和所有满足以上条件的ti作为待优化变量，使用非线性规划求解技术重新求解式(1.13-1.19)的非线性问题，得到当前时间网格下的最优的参数
Figure BDA0001858529530000078
和ti，以及新的目标函数值J，令Obj2＝J。Obj2表示当前时间网格下目标函数最小值。
步骤(5):如果l＝lmax或者(Obj1-Obj2)/Obj1≤Tol，则迭代结束，得到的ti值为满足要求的最佳时间网格节点，得到的N为满足要求额最佳网格节点数量，得到的
Figure BDA0001858529530000079
表示在此网格划分下最佳的控制参数，也就代表最佳的控制量
Figure BDA00018585295300000710
Tol表示用户允许误差，一般取值在10-4-10-8之间，lmax表示设定最大迭代次数，取值一般小于等于5。否则，如果(Obj1-Obj2)/Obj1＞Tol并且l＜lmax，则置l＝l+1，然后转到步骤(3)。
本发明不仅可以降低最优控制问题直接离散求解变量规模，而且计算量较小，并且迭代次数少，较少的参数获得高质量的解。该方法适用大规模复杂动态优化问题的在线优化：可以更有效地重构时间网格，找到准确的时间切换点，不仅计算成本低，而且计算精度更出色。
具体实施方式
一种最优控制问题直接离散求解的网格重构方法，包括如下步骤：
步骤(1)：采用基于控制向量参数化的直接法将式(1.1)～(1.7)所示的复杂最优控制问题离散为非线性规划问题，初始迭代次数l＝1。
Figure BDA0001858529530000081
Figure BDA0001858529530000082
Figure BDA0001858529530000083
Figure BDA0001858529530000084
x(t0)＝x0 (1.5)；
Figure BDA0001858529530000085
t0≤t≤tf (1.7)；
其中J是目标函数，由终值项φ0[x(tf),tf]和积分项
Figure BDA0001858529530000086
组成，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分别表示状态变量和控制变量的维数。
Figure BDA00018585295300000810
是x(t)的导数。φr[x(tf)]表示关于终端状态x(tf)的末端性能函数，Lr[t,x(t),u(t)]是时间t、状态变量x(t)和控制变量u(t)的复合函数，
Figure BDA0001858529530000087
表示m1个等式约束，
Figure BDA0001858529530000088
表示m2个不等式约束。x(t0)＝x0是状态变量在起始时间的初值，u和
Figure BDA0001858529530000089
分别是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维分量uj(t)就可以由各个时间子区间的值近似表示为式(1.9)：
Figure BDA0001858529530000091
其中，
Figure BDA0001858529530000092
为控制变量uj(t)在子区间[ti-1,ti)的值，T[ti-1,ti)为单位开关函数，被定义为式(1.10)：
Figure BDA0001858529530000093
在各个时间子区间内的控制变量
Figure BDA0001858529530000094
均由一系列基函数的线性组合来近似，即式(1.11)：
Figure BDA0001858529530000095
其中，
Figure BDA0001858529530000096
为Qi,j阶基函数，
Figure BDA0001858529530000097
为线性组合系数，称为控制参数。对函数
Figure BDA0001858529530000098
采用分段常量(分段零次多项式)逼近策略，则有k＝Qi,j＝1，以及
Figure BDA0001858529530000099
式(1.11)可简化为式(1.12)：
Figure BDA00018585295300000910
Figure BDA00018585295300000911
是控制参数
Figure BDA00018585295300000912
的简化形式，是式(1.1)～(1.7)需要求取的优化变量。
由此便可将一无限维的动态优化问题(1.1)～(1.7)转换为一个含有限维的控制参数
Figure BDA0001858529530000101
的非线性规划问题(1.13)～(1.19)，即:
Figure BDA0001858529530000102
Figure BDA0001858529530000103
Figure BDA0001858529530000104
Figure BDA0001858529530000105
x(t0)＝x0 (1.17)；
Figure BDA0001858529530000106
t0≤t≤tf (1.19)；
步骤(2)：使用非线性规划求解技术求解式(1.13)～(1.19)的非线性问题，得到当前时间网格下的最优的控制参数
Figure BDA0001858529530000107
和目标函数值J，令Obj1＝J，其中非线性规划求解技术为已有成熟技术。Obj1表示当前时间网格下的最小目标函数值。
步骤(3):将得到的控制参数根据时间顺序形成控制参数轨迹，对该轨迹进行伪Wigner-Ville分布分析变换，得到控制参数轨迹的时频曲线ωj(t)，以及在时间点的ti-1瞬时频率
Figure BDA0001858529530000108
其中ωj(t)表示第j维控制参数轨迹对应的频率，i＝1,2,…,N，j＝1,2,…,n。伪Wigner-Ville分布变换技术是已有成熟技术。
步骤(4):进行时间网格精细化重构策略，重新划分时间网格以确保求解精度，本部分包括三个子步骤：
子步骤1):对于相邻的时间网格节点，找出其中瞬时频率变化比较小的进行网格合并。对于相邻的时间网格[ti,ti+1]与[ti+1,ti+2]，如果控制变量uj(t)在时间点ti、ti+1和ti+2对应的瞬时频率
Figure BDA0001858529530000111
满足方程(1.20)，则将时间网格[ti,ti+1]与[ti+1,ti+2]合并成一个网格[ti,ti+2]。
Figure BDA0001858529530000112
其中，
Figure BDA0001858529530000113
为消除系数，
Figure BDA0001858529530000114
取值为0.1～0.5Hz。
Figure BDA0001858529530000115
为控制变量的变化系数，其取值满足如式(1.21)所示的规则：
Figure BDA0001858529530000116
合并后时间节点和瞬时频率都将被重新标记为式(1.22)：
Figure BDA0001858529530000117
合并后的时间网格数量记为N'。
子步骤2):根据瞬时频率的大小对时间网格进行细分，将网格均匀的划分为Δk个小区间，Δk的数目由以下经验规则决定。
Figure BDA0001858529530000118
其中
Figure BDA0001858529530000119
称为细化系数，取值由下式(1.24)决定：
Figure BDA0001858529530000121
其中
Figure BDA0001858529530000122
为控制变量的走势系数，取值由式(1.25)决定。
Figure BDA0001858529530000123
时间网格[ti-1,ti]被细化分为Δk个小区间后的时间节点将被重新标记如式(1.26)：
Figure BDA0001858529530000124
式中
Figure BDA0001858529530000125
为再次重新标记后时间节点。经过以上网格细分后，新的时间网格数量记为
Figure BDA0001858529530000126
令
Figure BDA0001858529530000127
将各个网格时间点表示为
Figure BDA00018585295300001211
子步骤3):对重要时间切换点实现精确定位，获得每个网格的最佳划分方式。根据步骤(3)经过控制参数轨迹进行伪Wigner-Ville分布变换得到的时频曲线，找到控制变量uj(t)对应的瞬时频率ωj存在极大值的时间点，若该时间点在时间网格[ti-1,ti]内，则将ti作为一个待优化变量。令新的网格数
Figure BDA0001858529530000128
将控制参数
Figure BDA0001858529530000129
和所有满足以上条件的ti作为待优化变量，使用非线性规划求解技术重新求解式(1.13-1.19)的非线性问题，得到当前时间网格下的最优的参数
Figure BDA00018585295300001210
和ti，以及新的目标函数值J，令Obj2＝J。Obj2表示当前时间网格下目标函数最小值。
步骤(5):如果l＝lmax或者(Obj1-Obj2)/Obj1≤Tol，则迭代结束，得到的ti值为满足要求的最佳时间网格节点，得到的N为满足要求额最佳网格节点数量，得到的
Figure BDA0001858529530000131
表示在此网格划分下最佳的控制参数，也就代表最佳的控制量
Figure BDA0001858529530000132
Tol表示用户允许误差，一般取值在10-4-10-8之间，lmax表示设定最大迭代次数，取值一般小于等于5。否则，如果(Obj1-Obj2)/Obj1＞Tol并且l＜lmax，则置l＝l+1，然后转到步骤(3)。

Claims (1)

1.一种最优控制问题直接离散求解的网格重构方法；其特征在于：基于伪Wigner-Ville分布分析快速而有效的网格重构策略，用于复杂最优控制问题直接离散求解；首先给定网格进行第一次优化迭代，快速获得控制变量的大致轨迹；然后通过伪Wigner-Ville分布分析得出不同时间网格节点瞬时频率变化对性能指标的影响，籍此对原有网格节点进行重构，包括对时间节点的消除、细化；并且结合变时间节点控制向量参数化方法的思想，将瞬时频率为极大值时对应的时间节点作为待优化参数，与控制变量一同进行求解优化，从而找到准确的最优时间切换点；
包括以下步骤：(1)该方法首先根据经验给定初始时间网格数量，初始时间网格数为5～20个，目的为了快速获得控制变量的大致轨迹，为下次迭代提供良好的初值条件，提高求解速度；将最优控制问题根据初始网格数量离散化为非线性规划问题；(2)求解非线性规划问题，得到当前时间网格下的最优的控制参数和目标函数值；(3)将得到的控制参数根据时间顺序形成控制参数轨迹，对该轨迹进行伪Wigner-Ville分析变换，得到控制参数轨迹的时频曲线；(4)根据规则对时间网格进行合并和细分，根据重要时间切换点的判断准则，筛选出需要优化的时间切换点；(5)判断是否满足终止条件，若满足，则输出重构后的时间网格；否则，转步骤(3)；
具体包括以下步骤：
步骤(1)：采用基于控制向量参数化的直接法将式(1.1)～(1.7)所示的复杂最优控制问题离散为非线性规划问题，初始迭代次数l＝1；
Figure FDA0003202025360000011
Figure FDA0003202025360000012
Figure FDA0003202025360000013
Figure FDA0003202025360000014
x(t0)＝x0 (1.5)；
Figure FDA0003202025360000021
t0≤t≤tf (1.7)；
其中J是目标函数，由终值项φ0[x(tf),tf]和积分项
Figure FDA0003202025360000022
组成，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分别表示状态变量和控制变量的维数；
Figure FDA0003202025360000027
是x(t)的导数；φr[x(tf)]表示关于终端状态x(tf)的末端性能函数，Lr[t,x(t),u(t)]是时间t、状态变量x(t)和控制变量u(t)的复合函数，
Figure FDA0003202025360000023
表示m1个等式约束；
Figure FDA0003202025360000024
表示m2个不等式约束；x(t0)＝x0是状态变量在起始时间的初值，u和
Figure FDA0003202025360000025
分别是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维分量uj(t)就可以由各个时间子区间的值近似表示为式(1.9)：
Figure FDA0003202025360000026
其中，
Figure FDA0003202025360000031
为控制变量uj(t)在子区间[ti-1,ti)的值，T[ti-1,ti)为单位开关函数，被定义为式(1.10)：
Figure FDA0003202025360000032
在各个时间子区间内的控制变量
Figure FDA0003202025360000033
均由一系列基函数的线性组合来近似，即式(1.11)：
Figure FDA0003202025360000034
其中，
Figure FDA0003202025360000035
为Qi,j阶基函数，
Figure FDA0003202025360000036
为线性组合系数，称为控制参数；对函数
Figure FDA0003202025360000037
采用分段常量逼近策略，即分段零次多项式逼近策略，则有k＝Qi,j＝1，以及
Figure FDA0003202025360000038
式(1.11)可简化为式(1.12)：
Figure FDA0003202025360000039
Figure FDA00032020253600000310
是控制参数
Figure FDA00032020253600000311
的简化形式，是式(1.1)～(1.7)需要求取的优化变量；
由此便可将一无限维的动态优化问题(1.1)～(1.7)转换为一个含有限维的控制参数
Figure FDA00032020253600000316
的非线性规划问题(1.13)～(1.19)，即:
Figure FDA00032020253600000312
Figure FDA00032020253600000313
Figure FDA00032020253600000314
Figure FDA00032020253600000315
x(t0)＝x0 (1.17)；
Figure FDA0003202025360000041
t0≤t≤tf (1.19)；
步骤(2)：使用非线性规划求解技术求解式(1.13)～(1.19)的非线性问题，得到当前时间网格下的最优的控制参数
Figure FDA0003202025360000042
和目标函数值J，令Obj1＝J，其中非线性规划求解技术为已有成熟技术；Obj1表示当前时间网格下的最小目标函数值；
步骤(3):将得到的控制参数根据时间顺序形成控制参数轨迹，对该轨迹进行伪Wigner-Ville分布分析变换，得到控制参数轨迹的时频曲线ωj(t)，以及在时间点的ti-1瞬时频率
Figure FDA0003202025360000048
其中ωj(t)表示第j维控制参数轨迹对应的频率，i＝1,2,…,N，j＝1,2,…,n；伪Wigner-Ville分布变换技术是已有成熟技术；
步骤(4):进行时间网格精细化重构策略，重新划分时间网格以确保求解精度，本部分包括三个子步骤：
子步骤1):对于相邻的时间网格节点，找出其中瞬时频率变化比较小的进行网格合并；对于相邻的时间网格[ti,ti+1]与[ti+1,ti+2]，如果控制变量uj(t)在时间点ti、ti+1和ti+2对应的瞬时频率
Figure FDA0003202025360000043
满足方程(1.20)，则将时间网格[ti,ti+1]与[ti+1,ti+2]合并成一个网格[ti,ti+2]；
Figure FDA0003202025360000044
其中，
Figure FDA0003202025360000045
为消除系数，
Figure FDA0003202025360000046
取值为0.1～0.5Hz；
Figure FDA0003202025360000047
为控制变量的变化系数，其取值满足如式(1.21)所示的规则：
Figure FDA0003202025360000051
合并后时间节点和瞬时频率都将被重新标记为式(1.22)：
Figure FDA0003202025360000052
合并后的时间网格数量记为N'；
子步骤2):根据瞬时频率的大小对时间网格进行细分，将网格均匀的划分为Δk个小区间，Δk的数目由以下式(1.23)所示经验规则决定；
Figure FDA0003202025360000053
其中
Figure FDA0003202025360000054
称为细化系数，取值由下式(1.24)决定：
Figure FDA0003202025360000055
其中
Figure FDA0003202025360000056
为控制变量的走势系数，取值由式(1.25)决定；
Figure FDA0003202025360000057
时间网格[ti-1,ti]被细化分为Δk个小区间后的时间节点将被重新标记如式(1.26)：
Figure FDA0003202025360000061
式中
Figure FDA0003202025360000062
为再次重新标记后时间节点；经过以上网格细分后，新的时间网格数量记为
Figure FDA0003202025360000063
令
Figure FDA0003202025360000064
将各个网格时间点表示为
Figure FDA0003202025360000065
子步骤3):对重要时间切换点实现精确定位，获得每个网格的最佳划分方式；根据步骤(3)经过控制参数轨迹进行伪Wigner-Ville分布变换得到的时频曲线，找到控制变量uj(t)对应的瞬时频率ωj存在极大值的时间点，若该时间点在时间网格[ti-1,ti]内，则将ti作为一个待优化变量；令新的网格数
Figure FDA0003202025360000066
将控制参数
Figure FDA0003202025360000067
和所有满足以上条件的ti作为待优化变量，使用非线性规划求解技术重新求解式(1.13-1.19)的非线性问题，得到当前时间网格下的最优的参数
Figure FDA00032020253600000610
和ti，以及新的目标函数值J，令Obj2＝J；Obj2表示当前时间网格下目标函数最小值；
步骤(5):如果l＝lmax或者(Obj1-Obj2)/Obj1≤Tol，则迭代结束，得到的ti值为满足要求的最佳时间网格节点，得到的N为满足要求额最佳网格节点数量，得到的
Figure FDA0003202025360000068
表示在此网格划分下最佳的控制参数，也就代表最佳的控制量
Figure FDA0003202025360000069
Tol表示用户允许误差，取值在10-4-10-8之间，lmax表示设定最大迭代次数，取值小于等于5；否则，如果(Obj1-Obj2)/Obj1＞Tol并且l＜lmax，则置l＝l+1，然后转到步骤(3)。
CN201811325145.3A 2018-11-08 2018-11-08 一种最优控制问题直接离散求解的网格重构方法 Active CN109491242B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811325145.3A CN109491242B (zh) 2018-11-08 2018-11-08 一种最优控制问题直接离散求解的网格重构方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811325145.3A CN109491242B (zh) 2018-11-08 2018-11-08 一种最优控制问题直接离散求解的网格重构方法

Publications (2)

Publication Number Publication Date
CN109491242A CN109491242A (zh) 2019-03-19
CN109491242B true CN109491242B (zh) 2021-10-08

Family

ID=65695381

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811325145.3A Active CN109491242B (zh) 2018-11-08 2018-11-08 一种最优控制问题直接离散求解的网格重构方法

Country Status (1)

Country Link
CN (1) CN109491242B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110119085B (zh) * 2019-04-03 2022-03-15 杭州电子科技大学 一种Manutec R3型工业机器人动态优化***
CN110173589B (zh) * 2019-04-30 2020-08-04 杭州电子科技大学 一种基于开关式压电阀的智能阀门定位***
CN110109430B (zh) * 2019-04-30 2020-09-22 杭州电子科技大学 一种间歇式啤酒发酵装置优化控制***
CN111324035A (zh) * 2019-11-21 2020-06-23 浙江大学 一种高超声速飞行器轨迹优化自适应最优控制器

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101887260A (zh) * 2010-06-30 2010-11-17 浙江大学 一种自适应同步策略的工业过程最优控制***及方法
CN101915911A (zh) * 2010-08-26 2010-12-15 哈尔滨工业大学 基于相消积累空时谱的空间任意构型分布式sar动目标参数估计方法
CN102798891A (zh) * 2012-08-22 2012-11-28 电子科技大学 基于短时分数阶傅里叶变换的地震信号时频分解方法
CN104618279A (zh) * 2015-01-22 2015-05-13 沈阳大学 基于快速独立分量分析的Wigner-Ville分布交叉项消除的信号处理方法
WO2015175090A1 (en) * 2014-05-13 2015-11-19 Wisconsin Alumni Research Foundation Method and apparatus for rapid acquisition of elasticity data in three dimensions
CN106157339A (zh) * 2016-07-05 2016-11-23 华南理工大学 基于低秩顶点轨迹子空间提取的动画网格序列压缩算法
CN107703899A (zh) * 2017-11-13 2018-02-16 浙江大学 一种基于经验模态分解动态优化的催化剂混合反应控制装置

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10204643B2 (en) * 2016-03-31 2019-02-12 OmniSpeech LLC Pitch detection algorithm based on PWVT of teager energy operator
US10573065B2 (en) * 2016-07-29 2020-02-25 Activision Publishing, Inc. Systems and methods for automating the personalization of blendshape rigs based on performance capture data

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101887260A (zh) * 2010-06-30 2010-11-17 浙江大学 一种自适应同步策略的工业过程最优控制***及方法
CN101915911A (zh) * 2010-08-26 2010-12-15 哈尔滨工业大学 基于相消积累空时谱的空间任意构型分布式sar动目标参数估计方法
CN102798891A (zh) * 2012-08-22 2012-11-28 电子科技大学 基于短时分数阶傅里叶变换的地震信号时频分解方法
WO2015175090A1 (en) * 2014-05-13 2015-11-19 Wisconsin Alumni Research Foundation Method and apparatus for rapid acquisition of elasticity data in three dimensions
CN104618279A (zh) * 2015-01-22 2015-05-13 沈阳大学 基于快速独立分量分析的Wigner-Ville分布交叉项消除的信号处理方法
CN106157339A (zh) * 2016-07-05 2016-11-23 华南理工大学 基于低秩顶点轨迹子空间提取的动画网格序列压缩算法
CN107703899A (zh) * 2017-11-13 2018-02-16 浙江大学 一种基于经验模态分解动态优化的催化剂混合反应控制装置

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于高斯混合模型与改进网格搜索法的轴承故障诊断;陈远帆;《重庆理工大学学报》;20160331;第30卷(第3期);第34-40页 *

Also Published As

Publication number Publication date
CN109491242A (zh) 2019-03-19

Similar Documents

Publication Publication Date Title
CN109491242B (zh) 一种最优控制问题直接离散求解的网格重构方法
Peng et al. A new unbiased stochastic derivative estimator for discontinuous sample performances with structural parameters
CN107563550B (zh) 一种基于pmu的配电网实时分布式状态估计及pmu的优化配置方法
CN109839825B (zh) 一种稀土萃取过程组分含量的预测控制方法及***
CN111182564B (zh) 一种基于lstm神经网络的无线链路质量预测方法
CN106156434B (zh) 基于局部时滞重构的滑动窗时间差-高斯过程回归建模方法
Pei et al. The Improved GM (1, N) Models with Optimal Background Values: a Case Study of Chinese High-tech Industry.
CN112018758A (zh) 基于数字孪生的含高比例新能源交直流混联***建模方法
CN112381279B (zh) 一种基于vmd和bls组合模型的风电功率预测方法
CN112113146B (zh) 供水管网管道粗糙系数和节点需水量同步自适应校核方法
Wang et al. H∞ fuzzy PID control for discrete time-delayed TS fuzzy systems
CN109684723B (zh) 一种二维结构内部声学性能分析方法
CN107909202B (zh) 一种基于时间序列的油井产液量集成预测方法
Sun et al. Short-term power load prediction based on VMD-SG-LSTM
CN109635330B (zh) 一种基于直接法的复杂优化控制问题准确和快速求解方法
CN116449779A (zh) 基于Actor-Critic结构的汽车车身喷涂用环境数据分析方法
CN106647250A (zh) 基于离线优化/在线查表方式的双层结构预测控制方法
CN114626207A (zh) 一种构建面向工业负荷的谐波发射水平的通用模型的方法
Wu et al. U-model-based adaptive control for a class of stochastic non-linear dynamic plants with unknown parameters
CN110658722B (zh) 一种基于gap的自均衡多模型分解方法及***
CN114329765A (zh) 一种高超声速飞行器气动外形自动优化设计***
CN111950123A (zh) 一种陀螺仪误差系数曲线拟合预测方法及***
Li et al. Closed-loop identification of systems using hybrid Box–Jenkins structure and its application to PID tuning
Rauh et al. Quantification of Time-Domain Truncation Errors for the Reinitialization of Fractional Integrators
CN112748660A (zh) 存储器、加热炉出口温度控制方法、装置和设备

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
EE01 Entry into force of recordation of patent licensing contract

Application publication date: 20190319

Assignee: Hangzhou LIANTENG Network Technology Co.,Ltd.

Assignor: HANGZHOU DIANZI University

Contract record no.: X2022330000004

Denomination of invention: A grid reconstruction method for direct discrete solution of optimal control problems

Granted publication date: 20211008

License type: Common License

Record date: 20220106

EE01 Entry into force of recordation of patent licensing contract