CN109635330B - 一种基于直接法的复杂优化控制问题准确和快速求解方法 - Google Patents

一种基于直接法的复杂优化控制问题准确和快速求解方法 Download PDF

Info

Publication number
CN109635330B
CN109635330B CN201811325710.6A CN201811325710A CN109635330B CN 109635330 B CN109635330 B CN 109635330B CN 201811325710 A CN201811325710 A CN 201811325710A CN 109635330 B CN109635330 B CN 109635330B
Authority
CN
China
Prior art keywords
finite element
length
state variable
configuration
ith
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
CN201811325710.6A
Other languages
English (en)
Other versions
CN109635330A (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 CN201811325710.6A priority Critical patent/CN109635330B/zh
Publication of CN109635330A publication Critical patent/CN109635330A/zh
Application granted granted Critical
Publication of CN109635330B publication Critical patent/CN109635330B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Complex Calculations (AREA)

Abstract

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

Description

一种基于直接法的复杂优化控制问题准确和快速求解方法
技术领域
本发明属于动态优化控制技术领域，涉及一种基于直接法的复杂优化控制问题准确和快速求解方法。
背景技术
化工反应过程、优化设计、动态过程系统参数估计、生产过程工作点切换和过程系统优化控制等过程中，存在诸多复杂优化控制问题。这类问题一般含有微分和代数方程，以及众多的轨线等式和不等式约束。对于复杂优化控制问题，传统上采用庞特里亚金极大值原理(Pontryagin’s Maximum Principle)来得到最优变分条件，然后进行求解，但是由于存在不少不等式轨线约束，采用极大值原理往往难以处理复杂不等式轨线约束而导致求解困难。并且该方法需要引入更多的协态变量及其相应的互补条件，会不可避免地带来耗时巨大的组合问题。
动态规划算法也往往用来求解优化控制问题，但是该方法不能处理状态变量的边界约束，因此需要在目标函数增加罚函数项来消去这些约束。动态规划算法的求解效率低于基于导数的一般优化算法，且对于状态变量含有边界约束和整体含有不等式约束的复杂问题需要处理后才能求解，且效率较低，所以该方法一般仅适合求解小规模问题。
随着计算机和计算技术的发展，解决复杂优化控制问题的方法往往是采用联立法，联立法的原理是将动态问题的整个时间域中的控制变量以及状态变量进行离散化，这样便能够将原动态优化问题转化成一个大规模的非线性规划问题。离散方法中采用较多的是有限元正交配置法，一般根据经验去选取有限元的个数，并且均匀划分每个有限元的长度，为了确保离散精度满足要求，有限元的个数往往选的非常多，使得离散后非线性规划问题的变量规模很大，造成其优化求解耗时非常大。为此部分研究者给出有限元非配置点处的误差方程，并将误差方程作为约束加入到原非线性规划问题中，并由此重新划分有限元网格，但是该方法增加了原非线性规划问题的非线性程度，求解计算复杂，而且有限元非配置点处的误差方程在计算过程中与实际结果有较大偏差，也影响了优化求解结果的合理性。为此，本发明提出一种基于直接法的复杂优化控制问题准确和快速求解方法以解决上述问题。
发明内容
本发明的目的是针对以往直接离散求解方法存在的不足，提供一种基于直接法的复杂优化控制问题准确和快速求解方法。
本发明首先根据经验给定ne个均分有限元网格，然后分别采用3阶和4阶Radau配置点将原复杂优化控制问题离散化为2个非线性规划问题，接着分别对以上非线性规划问题进行求解，并根据每个有限元内非配置点处状态变量的差值、每个有限元内状态变量的非线性程度以及控制变量的跳变情况重新划分网格，直到所有有限元网格内变量均满足给定要求，这时所得到的状态变量和控制变量值为满足给定要求的最优变量值。
本发明包括以下步骤：
步骤(1)：采用基于有限元配置值的联立法将式(1.1)～(1.8)所示的复杂优化控制问题离散为非线性规划问题。
g(z(t),y(t),u(t),t,p)＝0 (1.3)；
zL≤z(t)≤zU (1.4)；
uL≤u(t)≤uU (1.5)；
yL≤y(t)≤yU (1.6)；
t0≤t≤tf (1.7)；
z(t0)＝z0 (1.8)；
这里表示标量目标函数，z(t)、y(t)和u(t)分别表示与时间t相关的微分状态变量、代数状态变量和控制变量值。t0和tf表示开始与终端时间，p表示外界环境参数。z(tf)、y(tf)和u(tf)则分别表示在终端时刻微分状态变量、代数状态变量和控制变量的值。/>表示微分状态变量z(t)对时间t的导数。f表示微分方程形式的状态或者反应函数，g表示代数方程形式的过程轨线束方程，z0表示状态变量z(t)在t0时刻的初值，zL和zU表示状态变量z(t)的下界和上界，uL和uU分别表示控制变量u(t)的下界和上界，yL和yU表示代数状态变量y(t)的下界和上界。
对于式(1.1)～(1.8)所示的复杂优化控制问题，首先将时间区间[t0,tf]均匀离散化为ne个有限元网格，ne为整数，一般取5到15之间，每个有限元网格的长度hi可表示为式(2.1)：
hi＝(tf-t0)/ne,i＝1,...,ne (2.1)；
在每个有限元内插入K个配置点，配置点的相对位置选择Radau方程的根[ρ1,ρ2,...,ρK]，则各配置点所对应的时间可表示为式(2.2)：
ti,j＝ti-1+hiρj,j＝1,...,K (2.2)；
其中ti-1表示第i个有限元的初始时刻。
在第i个有限元内，微分状态变量可表示为式(2.3)：
代数状态变量可表示为式(2.4)：
控制变量可表示为式(2.5)：
这里，zi-1,0表示z(t)在第i个有限元内的初始值，hi是第i个有限元的长度，表示在第i个有限元第q个配置点处z(t)对时间的导数值，ti-1表示第i个有限元的初始时刻。ρr表示Radau方程的第r个根，Ωq为K阶多项式，满足式(2.6)：
Ωq(0)＝0 q＝1...,K (2.6)；
式(2.7)中，Ω'q(ρr)表示多项式Ωq在ρr处对时间的导数。
yi,q和ui,q分别表示在第i个有限元第q个配置点处代数变量y(t)和控制变量u(t)的值，ψq表示在第i个有限元第q个配置点的拉格朗日函数，其形式如式(2.8)：
其中，ti,j表示在第i个有限元第j个配置点处的时间，ρq和ρj表示第q个和j个Radau方程的根，且满足式(2.9)：
考虑到微分状态变量的连续性，在下一个有限元微分状态变量的初值等于前一个有限元微分状态变量的终值，因此有式(2.10)：
zi,0表示z(t)在第i+1个有限元内的初始值，zi+1,0表示z(t)在第i+1个有限元内的最终值，Ωq(1)表示时候多项式Ωq的值。
根据以上离散策略，式(1.1)～(1.8)形式的复杂优化控制问题离散化为如式(3.1)-(3.10)所示的非线性规划问题：
0＝g(zi,j,yi,j,ui,j,p) (3.3)；
z1,0＝z0 (3.5)；
zL≤zi,k≤zU (3.6)；
zL≤zi,0≤zU (3.7)；
yL≤yi,k≤yU (3.8)；
uL≤ui,k≤uU (3.9)；
i＝1...ne,j＝1...K (3.10)；
步骤(2)：根据步骤(1)所述离散策略，给定ne个均分限元网格，对采用3个配置点和4个配置点情况下式(3.1)-(3.10)所示的优化命题进行求解，给出最优状态变量、控制变量在不同时间点的值。ne一般取值5-15，每个有限元网格的长度由式(2.1)得到；在每个有限元内插入3个配置点，配置点的相对位置选择3阶Radau方程的根[ρ1,ρ2,ρ3]；采用非线性规划求解器求解式式(3.1)-(3.10)所示的优化命题，得到3个配置点情况下所有状态变量值、控制变量值以及状态变量导数值。然后，选择 表示4阶Radau方程的根/>中的第j个根(j∈1...4)。根据式(2.3)计算每个有限元内τnc处的状态变量值，该处的值记为再次，在同样ne个有限元网格下，在每个有限元内插入4个配置点，配置点的相对位置选择4阶Radau方程的根/>然后采用非线性规划求解器求解式(3.1)-(3.10)所示的优化命题，得到ne个有限元和4个配置点下所有状态变量值、控制变量值以及状态变量导数值。
步骤(3)：根据步骤(2)所得求解信息计算在3个和4个配置点情况下每个有限元内状态变量的差值和非线性信息。对于状态变量z(t)，其在3个配置点情况下τnc处的值为在4个配置点情况下τnc处的值为/>因为该处对应4个配置点情况下/>处的值zi,j，因此在τnc处采用3个配置点和4个配置点得到的状态变量的差值可表示为式(4.1)：
这里Err(i)表示第i个有限元的τnc处状态变量的差值。
每个有限元网格内状态变量非线性度Nonl(i)可表示为式(4.2)：
Nonl(i)＝max||dzi,j|-|dzi,k||,j,k＝1,...,K (4.2)；
步骤(4)：根据以上计算信息对有限元网格进行重新划分，确定新的有限元网格数量和每个有限元网格长度。假设状态变量z(t)允许误差为σ，每个有限元内状态变量非线性度要求低于ζ，如果每个有限元内Err(i)＜σ，且Nonl(i)＜ζ，则不再重新划分有限元网格，终止计算，此时该有限元网格内变量均满足给定要求，所得到的状态变量和控制变量值为满足给定要求的最优变量值。
否则根据误差阶数与网格长度之间的关系以及变量非线性度对误差的影响，采用以下规则重新划分有限元网格：
1)在长度为hi的第i个有限元内，如果则增加两个有限元，每个有限元网格长度为hi/3。
2)在长度为hi的第i个有限元内，如果且Nonl(i)＞ζ，则增加两个有限元，每个网格长度为hi/3。如果/>但是Nonl(i)≤ζ，则仅增加一个有限元，每个有限元网格长度为hi/2。
3)在长度为hi的第i个有限元内，如果且Nonl(i)＞ζ，则增加一个有限元，每个网格长度为hi/2。如果/>且Nonl(i)≤ζ，则保持原有有限元网格长度不变。
4)在长度为hi的第i个有限元内和网格长度为hi+1的第i+1个有限元内，如果且/>则第i个有限元与第i+1个有限元合并，新的有限元网格长度为hi+hi+1。
5)在长度为hi的第i个有限元内、长度为hi+1的第i+1个有限元内和长度为hi+2的第i+2个有限元内，如果满足且/>且/>则第i个有限元与第i+1个有限元、第i+2个有限元合，新的有限元网格长度为hi+hi+1+hi+2。
6)在长度为hi的第i个有限元内，如果但是4个配置点情况下的控制变量ui,j在第j(1≤j＜4)个配置点处的值达到其上界或者下界，则增加1个有限元，前后两个有限元网格长度分别为/>和/>
步骤(5)：根据以上网格划分得到新的有限元数量和每个有限元网格的长度，重新表示为ne和hi。在每个有限元内插入3个配置点，配置点的相对位置选择3阶Radau方程的根[ρ1,ρ2,ρ3]；采用非线性规划求解器求解式(3.1)-(3.10)所示的优化命题，得到3个配置点情况下所有状态变量值、控制变量值以及状态变量导数值。然后，根据式(2.3)计算每个有限元内τnc处的状态变量值，该处的值记为这里/> 表示4阶Radau方程的根中的一个根(j∈1...4)。再次，在同样ne个有限元网格下，在每个有限元内插入4个配置点，配置点的相对位置选择4阶Radau方程的根/>然后采用非线性规划求解器求解式(3.1)-(3.10)所示的优化命题，得到ne个有限元和4个配置点下所有状态变量值、控制变量值以及状态变量导数值。然后返回步骤(3)。
本发明不仅能够保证复杂优化控制问题中的状态变量离散误差满足用户要求，给出的有限元网格数量相对较少，方法简单有效，而且很容易快速、准确找到结构切换点。
具体实施方式
一种基于直接法的复杂优化控制问题准确和快速求解方法，具体包括以下步骤：
步骤(1)：采用基于有限元配置值的联立法将式(1.1)～(1.8)所示的复杂优化控制问题离散为非线性规划问题。
g(z(t),y(t),u(t),t,p)＝0 (1.3)；
zL≤z(t)≤zU (1.4)；
uL≤u(t)≤uU (1.5)；
yL≤y(t)≤yU (1.6)；
t0≤t≤tf (1.7)；
z(t0)＝z0 (1.8)；
这里表示标量目标函数，z(t)、y(t)和u(t)分别表示与时间t相关的微分状态变量、代数状态变量和控制变量值。t0和tf表示开始与终端时间，p表示外界环境参数。z(tf)、y(tf)和u(tf)则分别表示在终端时刻微分状态变量、代数状态变量和控制变量的值。/>表示微分状态变量z(t)对时间t的导数。f表示微分方程形式的状态或者反应函数，g表示代数方程形式的过程轨线束方程，z0表示状态变量z(t)在t0时刻的初值，zL和zU表示状态变量z(t)的下界和上界，uL和uU分别表示控制变量u(t)的下界和上界，yL和yU表示代数状态变量y(t)的下界和上界。
对于式(1.1)～(1.8)所示的复杂优化控制问题，首先将时间区间[t0,tf]均匀离散化为ne个有限元网格(ne为整数，一般取5到15之间)，每个有限元网格的长度hi可表示为式(2.1)：
hi＝(tf-t0)/ne,i＝1,...,ne (2.1)；
在每个有限元内插入K个配置点，配置点的相对位置选择Radau方程的根[ρ1,ρ2,...,ρK]，则各配置点所对应的时间可表示为式(2.2)：
ti,j＝ti-1+hiρj,j＝1,...,K (2.2)；
其中ti-1表示第i个有限元的初始时刻。
在第i个有限元内，微分状态变量可表示为式(2.3)：
代数状态变量可表示为式(2.4)：
控制变量可表示为式(2.5)：
这里，zi-1,0表示z(t)在第i个有限元内的初始值，hi是第i个有限元的长度，表示在第i个有限元第q个配置点处z(t)对时间的导数值，ti-1表示第i个有限元的初始时刻。ρr表示Radau方程的第r个根，Ωq为K阶多项式，满足式(2.6)：
Ωq(0)＝0 q＝1...,K (2.6)；
式(2.7)中，Ω'q(ρr)表示多项式Ωq在ρr处对时间的导数。
yi,q和ui,q分别表示在第i个有限元第q个配置点处代数变量y(t)和控制变量u(t)的值，ψq表示在第i个有限元第q个配置点的拉格朗日函数，其形式如式(2.8)：
其中，ti,j表示在第i个有限元第j个配置点处的时间，ρq和ρj表示第q个和j个Radau方程的根，且满足式(2.9)：
考虑到微分状态变量的连续性，在下一个有限元微分状态变量的初值等于前一个有限元微分状态变量的终值，因此有式(2.10)：
zi,0表示z(t)在第i+1个有限元内的初始值，zi+1,0表示z(t)在第i+1个有限元内的最终值，Ωq(1)表示时候多项式Ωq的值。
根据以上离散策略，式(1.1)～(1.8)形式的复杂优化控制问题离散化为如式(3.1)-(3.10)所示的非线性规划问题：
0＝g(zi,j,yi,j,ui,j,p) (3.3)；
z1,0＝z0 (3.5)；
zL≤zi,k≤zU (3.6)；
zL≤zi,0≤zU (3.7)；
yL≤yi,k≤yU (3.8)；
uL≤ui,k≤uU (3.9)；
i＝1...ne,j＝1...K (3.10)；
步骤(2)：根据步骤(1)所述离散策略，给定ne个均分限元网格，对采用3个配置点和4个配置点情况下式(3.1)-(3.10)所示的优化命题进行求解，给出最优状态变量、控制变量在不同时间点的值。ne一般取值5-15，每个有限元网格的长度由式(2.1)得到；在每个有限元内插入3个配置点，配置点的相对位置选择3阶Radau方程的根[ρ1,ρ2,ρ3]；采用非线性规划求解器求解式式(3.1)-(3.10)所示的优化命题，得到3个配置点情况下所有状态变量值、控制变量值以及状态变量导数值。然后，选择 表示4阶Radau方程的根/>中的第j个根(j∈1...4)。根据式(2.3)计算每个有限元内τnc处的状态变量值，该处的值记为再次，在同样ne个有限元网格下，在每个有限元内插入4个配置点，配置点的相对位置选择4阶Radau方程的根/>然后采用非线性规划求解器求解式(3.1)-(3.10)所示的优化命题，得到ne个有限元和4个配置点下所有状态变量值、控制变量值以及状态变量导数值。
步骤(3)：根据步骤(2)所得求解信息计算在3个和4个配置点情况下每个有限元内状态变量的差值和非线性信息。对于状态变量z(t)，其在3个配置点情况下τnc处的值为在4个配置点情况下τnc处的值为/>因为该处对应4个配置点情况下/>处的值zi,j，因此在τnc处采用3个配置点和4个配置点得到的状态变量的差值可表示为式(4.1)：
这里Err(i)表示第i个有限元的τnc处状态变量的差值。
每个有限元网格内状态变量非线性度Nonl(i)可表示为式(4.2)：
Nonl(i)＝max||dzi,j|-|dzi,k||,j,k＝1,...,K (4.2)；
步骤(4)：根据以上计算信息对有限元网格进行重新划分，确定新的有限元网格数量和每个有限元网格长度。假设状态变量z(t)允许误差为σ，每个有限元内状态变量非线性度要求低于ζ，如果每个有限元内Err(i)＜σ，且Nonl(i)＜ζ，则不再重新划分有限元网格，终止计算，此时该有限元网格内变量均满足给定要求。
否则根据误差阶数与网格长度之间的关系以及变量非线性度对误差的影响，采用以下规则重新划分有限元网格：
1)在长度为hi的第i个有限元内，如果则增加两个有限元，每个有限元网格长度为hi/3。
2)在长度为hi的第i个有限元内，如果且Nonl(i)＞ζ，则增加两个有限元，每个网格长度为hi/3。如果/>但是Nonl(i)≤ζ，则仅增加一个有限元，每个有限元网格长度为hi/2。
3)在长度为hi的第i个有限元内，如果且Nonl(i)＞ζ，则增加一个有限元，每个网格长度为hi/2。如果/>且Nonl(i)≤ζ，则保持原有有限元网格长度不变。
4)在长度为hi的第i个有限元内和网格长度为hi+1的第i+1个有限元内，如果且/>则第i个有限元与第i+1个有限元合并，新的有限元网格长度为hi+hi+1。
5)在长度为hi的第i个有限元内、长度为hi+1的第i+1个有限元内和长度为hi+2的第i+2个有限元内，如果满足且/>且/>则第i个有限元与第i+1个有限元、第i+2个有限元合，新的有限元网格长度为hi+hi+1+hi+2。
6)在长度为hi的第i个有限元内，如果但是4个配置点情况下的控制变量ui,j在第j(1≤j＜4)个配置点处的值达到其上界或者下界，则增加1个有限元，前后两个有限元网格长度分别为/>和/>
步骤(5)：根据以上网格划分得到新的有限元数量和每个有限元网格的长度，重新表示为ne和hi。在每个有限元内插入3个配置点，配置点的相对位置选择3阶Radau方程的根[ρ1,ρ2,ρ3]；采用非线性规划求解器求解式(3.1)-(3.10)所示的优化命题，得到3个配置点情况下所有状态变量值、控制变量值以及状态变量导数值。然后，根据式(2.3)计算每个有限元内τnc处的状态变量值，该处的值记为这里/> 表示4阶Radau方程的根中的一个根(j∈1...4)。再次，在同样ne个有限元网格下，在每个有限元内插入4个配置点，配置点的相对位置选择4阶Radau方程的根/>然后采用非线性规划求解器求解式(3.1)-(3.10)所示的优化命题，得到ne个有限元和4个配置点下所有状态变量值、控制变量值以及状态变量导数值。然后返回步骤(3)。

Claims (1)

1.一种基于直接法的复杂优化控制问题准确和快速求解方法；其特征在于：首先给定ne个均分有限元网格，然后分别采用3阶和4阶Radau配置点将原复杂优化控制问题离散化为2个非线性规划问题，接着分别对以上非线性规划问题进行求解，并根据每个有限元内非配置点处状态变量的差值、每个有限元内状态变量的非线性程度以及控制变量的跳变情况重新划分网格，直到所有有限元网格内变量均满足给定要求，这时所得到的状态变量和控制变量值为满足给定要求的最优变量值；
具体包括以下步骤：
步骤(1)：采用基于有限元配置值的联立法将式(1.1)～(1.8)所示的复杂优化控制问题离散为非线性规划问题；
g(z(t),y(t),u(t),t,p)＝0 (1.3)；
zL≤z(t)≤zU (1.4)；
uL≤u(t)≤uU (1.5)；
yL≤y(t)≤yU (1.6)；
t0≤t≤tf (1.7)；
z(t0)＝z0 (1.8)；
这里表示标量目标函数，z(t)、y(t)和u(t)分别表示与时间t相关的微分状态变量、代数状态变量和控制变量值；t0和tf表示开始与终端时间，p表示外界环境参数；z(tf)、y(tf)和u(tf)则分别表示在终端时刻微分状态变量、代数状态变量和控制变量的值；/>表示微分状态变量z(t)对时间t的导数；f表示微分方程形式的状态或者反应函数，g表示代数方程形式的过程轨线束方程，z0表示状态变量z(t)在t0时刻的初值，zL和zU表示状态变量z(t)的下界和上界，uL和uU分别表示控制变量u(t)的下界和上界，yL和yU表示代数状态变量y(t)的下界和上界；
对于式(1.1)～(1.8)所示的复杂优化控制问题，首先将时间区间[t0,tf]均匀离散化为ne个有限元网格，ne为整数，取5到15之间，每个有限元网格的长度hi表示为式(2.1)：
hi＝(tf-t0)/ne,i＝1,...,ne (2.1)；
在每个有限元内插入K个配置点，配置点的相对位置选择Radau方程的根[ρ1,ρ2,...,ρK]，则各配置点所对应的时间表示为式(2.2)：
ti,j＝ti-1+hiρj,j＝1,...,K (2.2)；
其中ti-1表示第i个有限元的初始时刻；
在第i个有限元内，微分状态变量表示为式(2.3)：
代数状态变量表示为式(2.4)：
控制变量表示为式(2.5)：
这里，zi-1,0表示z(t)在第i个有限元内的初始值，hi是第i个有限元的长度，表示在第i个有限元第q个配置点处z(t)对时间的导数值，ti-1表示第i个有限元的初始时刻；ρr表示Radau方程的第r个根，Ωq为K阶多项式，满足式(2.6)：
Ωq(0)＝0 q＝1,...,K (2.6)；
式(2.7)中，Ω'q(ρr)表示多项式Ωq在ρr处对时间的导数；
yi,q和ui,q分别表示在第i个有限元第q个配置点处代数变量y(t)和控制变量u(t)的值，ψq表示在第i个有限元第q个配置点的拉格朗日函数，其形式如式(2.8)：
其中，ti,j表示在第i个有限元第j个配置点处的时间，ρq和ρj表示第q个和j个Radau方程的根，且满足式(2.9)：
考虑到微分状态变量的连续性，在下一个有限元微分状态变量的初值等于前一个有限元微分状态变量的终值，因此有式(2.10)：
zi,0表示z(t)在第i+1个有限元内的初始值，zi+1,0表示z(t)在第i+1个有限元内的最终值，Ωq(1)表示时候多项式Ωq的值；
根据以上离散策略，式(1.1)～(1.8)形式的复杂优化控制问题离散化为如式(3.1)-(3.10)所示的非线性规划问题：
0＝g(zi,j,yi,j,ui,j,p) (3.3)；
z1,0＝z0 (3.5)；
zL≤zi,k≤zU (3.6)；
zL≤zi,0≤zU (3.7)；
yL≤yi,k≤yU (3.8)；
uL≤ui,k≤uU (3.9)；
i＝1,...,ne,j＝1,...,K (3.10)；
步骤(2)：根据步骤(1)所述离散策略，给定ne个均分限元网格，对采用3个配置点和4个配置点情况下式(3.1)-(3.10)所示的优化命题进行求解，给出最优状态变量、控制变量在不同时间点的值；ne取值为5-15，每个有限元网格的长度由式(2.1)得到；在每个有限元内插入3个配置点，配置点的相对位置选择3阶Radau方程的根[ρ1,ρ2,ρ3]；采用非线性规划求解器求解式式(3.1)-(3.10)所示的优化命题，得到3个配置点情况下所有状态变量值、控制变量值以及状态变量导数值；然后，选择 表示4阶Radau方程的根/>中的第j个根，j∈1,...,4；根据式(2.3)计算每个有限元内τnc处的状态变量值，该处的值记为/>再次，在同样ne个有限元网格下，在每个有限元内插入4个配置点，配置点的相对位置选择4阶Radau方程的根/>然后采用非线性规划求解器求解式(3.1)-(3.10)所示的优化命题，得到ne个有限元和4个配置点下所有状态变量值、控制变量值以及状态变量导数值；
步骤(3)：根据步骤(2)所得求解信息计算在3个和4个配置点情况下每个有限元内状态变量的差值和非线性信息；对于状态变量z(t)，其在3个配置点情况下τnc处的值为在4个配置点情况下τnc处的值为/>因为该处对应4个配置点情况下/>处的值zi,j，因此在τnc处采用3个配置点和4个配置点得到的状态变量的差值表示为式(4.1)：
这里Err(i)表示第i个有限元的τnc处状态变量的差值；
每个有限元网格内状态变量非线性度Nonl(i)表示为式(4.2)：
Nonl(i)＝max||dzi,j|-|dzi,k||,j,k＝1,...,K (4.2)；
步骤(4)：根据以上计算信息对有限元网格进行重新划分，确定新的有限元网格数量和每个有限元网格长度；假设状态变量z(t)允许误差为σ，每个有限元内状态变量非线性度要求低于ζ，如果每个有限元内Err(i)＜σ，且Nonl(i)＜ζ，则不再重新划分有限元网格，终止计算，此时该有限元网格内变量均满足给定要求，所得到的状态变量和控制变量值为满足给定要求的最优变量值；
否则根据误差阶数与网格长度之间的关系以及变量非线性度对误差的影响，采用以下规则重新划分有限元网格：
1)在长度为hi的第i个有限元内，如果则增加两个有限元，每个有限元网格长度为hi/3；
2)在长度为hi的第i个有限元内，如果且Nonl(i)＞ζ，则增加两个有限元，每个网格长度为hi/3；如果/>但是Nonl(i)≤ζ，则仅增加一个有限元，每个有限元网格长度为hi/2；
3)在长度为hi的第i个有限元内，如果且Nonl(i)＞ζ，则增加一个有限元，每个网格长度为hi/2；如果/>且Nonl(i)≤ζ，则保持原有有限元网格长度不变；
4)在长度为hi的第i个有限元内和网格长度为hi+1的第i+1个有限元内，如果且/>则第i个有限元与第i+1个有限元合并，新的有限元网格长度为hi+hi+1；
5)在长度为hi的第i个有限元内、长度为hi+1的第i+1个有限元内和长度为hi+2的第i+2个有限元内，如果满足且/>且/>则第i个有限元与第i+1个有限元、第i+2个有限元合并，新的有限元网格长度为hi+hi+1+hi+2；
6)在长度为hi的第i个有限元内，如果但是4个配置点情况下的控制变量ui,j在第j个配置点处的值达到其上界或者下界，1≤j＜4，则增加1个有限元，前后两个有限元网格长度分别为/>和/>
步骤(5)：根据以上网格划分得到新的有限元数量和每个有限元网格的长度，重新表示为ne和hi；在每个有限元内插入3个配置点，配置点的相对位置选择3阶Radau方程的根[ρ1,ρ2,ρ3]；采用非线性规划求解器求解式(3.1)-(3.10)所示的优化命题，得到3个配置点情况下所有状态变量值、控制变量值以及状态变量导数值；然后，根据式(2.3)计算每个有限元内τnc处的状态变量值，该处的值记为这里/> 表示4阶Radau方程的根/>中的一个根，j∈1，...，4；再次，在同样ne个有限元网格下，在每个有限元内插入4个配置点，配置点的相对位置选择4阶Radau方程的根/>然后采用非线性规划求解器求解式(3.1)-(3.10)所示的优化命题，得到ne个有限元和4个配置点下所有状态变量值、控制变量值以及状态变量导数值；然后返回步骤(3)。
CN201811325710.6A 2018-11-08 2018-11-08 一种基于直接法的复杂优化控制问题准确和快速求解方法 Active CN109635330B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811325710.6A CN109635330B (zh) 2018-11-08 2018-11-08 一种基于直接法的复杂优化控制问题准确和快速求解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811325710.6A CN109635330B (zh) 2018-11-08 2018-11-08 一种基于直接法的复杂优化控制问题准确和快速求解方法

Publications (2)

Publication Number Publication Date
CN109635330A CN109635330A (zh) 2019-04-16
CN109635330B true CN109635330B (zh) 2024-01-19

Family

ID=66067542

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811325710.6A Active CN109635330B (zh) 2018-11-08 2018-11-08 一种基于直接法的复杂优化控制问题准确和快速求解方法

Country Status (1)

Country Link
CN (1) CN109635330B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110109430B (zh) * 2019-04-30 2020-09-22 杭州电子科技大学 一种间歇式啤酒发酵装置优化控制系统
CN110173589B (zh) * 2019-04-30 2020-08-04 杭州电子科技大学 一种基于开关式压电阀的智能阀门定位系统

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107391885A (zh) * 2017-08-29 2017-11-24 西北工业大学 基于有限体积法的剪切滑移动网格方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CA2484649A1 (en) * 2002-05-07 2003-11-20 Markov Processes International, Llc A method and system to solve dynamic multi-factor models in finance

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107391885A (zh) * 2017-08-29 2017-11-24 西北工业大学 基于有限体积法的剪切滑移动网格方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
An improved finite element meshing strategy for optimal control of chemical process;Aipeng Jiang等;Proceedings of the 36th Chinese Control Conference July 26-28, 2017, Dalian, China;第4351-4356页 *
一种结构动态分析中有限元网格细分方法;程耀东等;浙江大学学报(自然科学版);第33-38页 *
一类非线性系统的变负荷预测控制;陈杨;邵之江;游江洪;万娇娜;钱积新;;化工学报(第08期);第2253-2257页 *
关于二次型目标函数算法优化仿真研究;邢丽娟等;计算机仿真;第300-303+312页 *

Also Published As

Publication number Publication date
CN109635330A (zh) 2019-04-16

Similar Documents

Publication Publication Date Title
CN110110419B (zh) 一种基于多目标学习的tbm掘进参数预测方法
CN104698842B (zh) 一种基于内点法的lpv模型非线性预测控制方法
CN109635330B (zh) 一种基于直接法的复杂优化控制问题准确和快速求解方法
CN107272409B (zh) 一种基于迭代学习的直线伺服系统振动抑制方法
JP5227254B2 (ja) プロセスモデルの状態量のリアルタイム計算方法およびシミュレータ
CN103472723A (zh) 基于多模型广义预测控制器的预测控制方法及系统
CN105892296B (zh) 一种工业加热炉系统的分数阶动态矩阵控制方法
Braga et al. Discretization and event triggered digital output feedback control of LPV systems
CN109491242B (zh) 一种最优控制问题直接离散求解的网格重构方法
Neçaibia et al. Self-tuning fractional order PI λ D µ controller based on extremum seeking approach
JPWO2020227383A5 (zh)
CN111553118B (zh) 基于强化学习的多维连续型优化变量全局优化方法
CN110032706A (zh) 一种低阶时滞系统的两阶段参数估计方法及系统
CN102880046A (zh) 一种化工多变量过程解耦预测函数控制方法
CN109581864A (zh) 参数自整定的mimo异因子偏格式无模型控制方法
CN109143853B (zh) 一种石油炼化过程中分馏塔液位自适应控制方法
CN103399488B (zh) 基于自学习的多模型控制方法
CN114509949A (zh) 一种机器人预定性能的控制方法
CN108009362B (zh) 一种基于稳定性约束rbf-arx模型的系统建模方法
CN107256290A (zh) 一种基于残差摄动法的边界条件变化对流场影响差量的计算方法
CN109782586A (zh) 参数自整定的miso异因子紧格式无模型控制方法
CN107908110B (zh) 基于控制网格精细化的管式反应器动态优化系统
Srinivasan et al. Dynamic optimization under uncertainty via NCO tracking: A solution model approach
CN110888323A (zh) 一种用于切换系统智能优化的控制方法
Costello et al. Real-time optimization when plant and model have different sets of inputs

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