动态优化又称为最优控制理论，自庞特里亚金(Pontryagin)极大值原理[1]和动态规划理论[2]提出后开始快速发展，并成功应用到航空航天、生物工程等诸多领域．由于各领域中受控系统日益复杂、非线性程度不断提高及实际约束的限制，基于理论分析的动态优化问题解析法难以甚至无法获得理想的系统控制律．因此，随着计算机性能的提升，间接法和直接法等数值求解方法[3]备受关注．其中直接法应用最为广泛，该方法主要通过直接转录将动态优化问题转化为非线性规划问题进行求解．当动态优化问题中的状态方程为微分代数方程表示的显式表达式时，可基于直接法和自动微分技术对动态优化问题进行快速求解．但在一些工程应用中，通过工业软件构建的动态系统仿真模型无法提供状态方程的显式表达式，如基于Modelica语言开发的多领域物理系统建模和仿真平台MWorks[4]中，其功能模块接口可集成由用于执行工程系统模拟的高级建模环境(AMESim)、机械系统动力学自动分析系统(ADAMS)和可视化仿真工具(SIMULINK)等平台构建的仿真模型，这些仿真模型的功能模型单元是二进制或C代码实现，无法提取显式表达的状态方程，这种情况下的动态优化问题可视为状态方程隐式表达的动态优化问题．虽然结合有限差分技术可对状态方程隐式表达的动态优化问题进行优化求解，但须要对仿真模型进行大量估值，大大增加了优化求解时间．为高效求解状态方程隐式表达的动态优化问题，现有方法主要利用响应面技术构建状态方程右端式的代理模型以显式化表达状态方程[5-6]．但在高维复杂的动态优化问题中，在整个状态空间上用单个模型近似状态方程右端式显得尤为困难，不可避免会遇到模型训练耗时长、收敛速度慢等问题．因此，可根据合适的空间划分策略将整个状态空间划分成多个局部子空间，进而在各子空间中分别构建近似模型．现有的空间划分策略主要体现在基于树结构的机器学习方法中，分类与回归树算法[7]通过穷举搜索法在每个维度上寻找最佳划分点．平滑和非平滑分段多项式回归树[8]和广义无偏交互检测与估计[9]算法通过计算残差统计值选择须要进行划分的维度，而基于模型的递归划分算法[10]则基于参数的不确定性来确定划分的维度．与上述基于局部优化算法进行维度划分的策略所不同的是，决策树算法[11]构建全局最优分类与回归树，分段回归方法和数学规划树[12]方法则采用数学规划方法选取划分点的位置．虽然这些机器学习算法主要用来对训练数据集进行划分，但也可用于对状态空间的划分．Qiao等[13]受这些划分策略启发，提出一种基于灵敏度指标和局部区域的区间长度信息的状态空间划分方法．然而，这些划分方法只基于维度信息对空间进行划分，子空间之间的非光滑边界使得代理模型导数具有较大的不连续性，影响了动态优化问题的求解精度．针对上述问题，为降低划分空间的非光滑边界对代理模型连续性及动态优化问题求解精度造成的影响，提出了一种新的空间划分策略，根据状态方程右端式的非线性程度将状态空间划分为多个局部区域．鉴于神经网络模型结构简洁，显式表达式便于提取，利用神经网络模型构建各个局部区域的代理模型，进而基于神经网络模型采用直接法对状态方程隐式表达的动态优化问题进行求解．最后，通过两个工程算例对所提方法的性能进行了验证．1 基础理论和方法1.1　动态优化问题及其数值求解动态优化问题[5]的目标是获得使动态系统性能J达到最优的控制变量u*(t)．状态方程显式表达的动态优化问题的一般形式如下    minu(t),tf  J=f(x(t0),t0,x(tf),tf)+∫t0tfL(x(t),u(t),t)dt,s.t.  x˙=f(x(t),u(t),t);       Ψ(x(t0),t0,x(tf),tf)=0;        g(x(t),u(t),t)≤0;uL≤u(t)≤uU,式中：性能指标J由迈耶尔(Mayer)项ϕ和拉格朗日(Lagrange)项组成；L为Lagrange项中的被积函数；t∈[t0,tf]为时间区间，t0为初始时刻，tf为终端时刻；x(t0)和x(tf)为系统的初始状态和终端状态；x(t)和u(t)为状态变量和控制变量在时刻t处的值；x˙=f(∙)为状态方程约束，其显式特性主要体现在微分变量导数x˙在状态方程左端，右端式f(∙)的数学表达式明确给出．除状态方程约束外，动态优化问题还受限于边界约束Ψ和路径约束g．这三种约束中，状态方程约束和路径约束是连续约束，在[t0,tf]整个时间段内都必须满足；而边界约束是离散约束，只须在时刻t0和tf上满足即可．最后，控制变量u(t)限制于区间约束[uL,uU]中．1.2　高斯混合模型聚类算法高斯混合模型聚类算法采用概率模型来表达聚类原型，是机器学习领域中被广泛应用的一种聚类方式[14]．高斯混合模型聚类算法的假设前提是所有数据样本X服从k个多元高斯分布组合而成的混合分布，该混合分布表示为Px | μ,δ2=∑i=1kαi∙ϕ(x | μi,δi2)，式中ϕ(x | μi,δi2)为第i个高斯分布的概率密度函数，ϕ(x | μi,δi2)=12πδiexp-(x-μi)22δi2，其中，αi为混合系数，αi≥0且∑αi=1；μi和δi2为第i个高斯分布的期望值和方差．因此，高斯混合模型聚类算法的实质是：设定聚类簇数k后，基于训练样本X，通过极大似然估计法推导出混合分布中每个高斯成分的期望值、方差及混合系数，共3k个未知参数．在训练高斯混合模型过程中，当聚类簇数k未指定时，可通过赤池信息准则(Akaike information criterion，AIC)或者贝叶斯信息准则(Bayesian information criteria，BIC)来确定合适的聚类簇数．当基于训练完成的高斯混合模型训练对新样本点xnew进行分类时，将该样本带入到k个高斯分布中求出属于每个类别的概率，概率最大者即为新样本点所归属的类．1.3　基于最大曲率最小点距的序列采样方法基于最大曲率最小点距的序列采样方法(sequence sampling method based on maximum curvature minimum point distance)主要根据曲率，也即局部区域的弯曲程度进行采样；同时，最大化任意两样本点之间的最小距离以避免在极端弯曲的区域附近集中采样．因此，基于最大曲率最小点距的序列采样方法通过最大化代理模型在新样本点xnew处的曲率K(xnew)和新样本点xnew与现有样本点之间最小距离dmin(xnew)的乘积以获得新样本点xnew．根据微分几何原理，曲率K可通过海森矩阵的特征值进行估算，当海森矩阵有q个特征值k1,k2,…,kq，点x处的曲率可通过下式K(x)=∑i=1qki(x)2进行估算．基于最大曲率最小点距的序列采样方法的具体应用及其可行性在文献[15]进行了具体说明和验证．2 基于高斯混合模型聚类的分区神经网络构建方法当基于空间划分策略在整个状态空间上构建状态方程右端式的近似模型时，须要将状态空间划分成多个局部区域．当划分状态空间的非光滑边界经过状态方程非线性程度较高、梯度较大的区域时，如图1(a)所示，状态空间上有波峰A，B和波谷C，D，两条边界将状态空间划分为四个局部区域Ω1,Ω2,Ω3,Ω4，且分别将每个波峰波谷分割成两部分，并归属到相邻的两个局部区域中，将导致建模复杂度提高．同时，四个局部区域的代理模型构建完成后，波峰A所在区域梯度较大，以点A1和A2为代表的边界两边附近的点将通过两个不同的局部模型进行估值，会造成两点之间的近似结果存在较严重的不连续性．当基于该模型求解动态优化问题时，这种不连续性会导致代理模型提供的雅可比信息矩阵与真实值误差较大，影响动态优化问题的求解精度．10.13245/j.hust.238733.F001图1不同空间划分策略的示意图本研究提出一种基于高斯混合模型聚类的空间划分策略(如图1(b)所示)，该策略使得边界经过状态方程非线性程度较低、梯度较小的区域，避免了对非线性程度较高区域的分割，降低了建模复杂度；同时，边界附近的点虽然通过不同局部模型估值，但由于边界两边梯度较小，对估值结果的连续性影响较小，进而保证了基于代理模型的隐式状态方程动态优化问题的求解精度．2.1　基于高斯混合模型聚类的空间划分策略基于高斯混合模型的空间划分及多局部神经网络模型构建的具体算法流程如下．步骤1 初始化相关参数，如高斯混合模型聚类的最大簇数，局部神经网络模型参数、训练方法、训练次数、目标精度等．步骤2 基于最大曲率最小点距的序列采样方法获取初始样本点集S0，S0主要特征为聚集在状态方程右端式非线性程度高的区域，且非线性程度越高样本点越密集．步骤3 构建基于S0的高斯混合模型M，M={Ci | Ci=Cluster(S0),i=1,2,…,n}，式中：Cluster表示采用聚类算法对S0进行聚类；聚类簇数n可根据模型反馈的赤池信息准则或贝叶斯信息准则自适应确定．步骤4 采用并行仿真串行化采样方法[16]在MWorks平台中，通过仿真模型获取大批量便宜估值的样本点集S0'．步骤5 基于已构建的高斯混合模型对S0'进行分类，S={S1,S2,…,Sn}，同时将状态空间划分为n个不规则的局部区域Ω={Ω1,Ω2,…,Ωn}，S={Si | Si=Classify(S0'),i=1,2,…,n},Ω={Ωi | Ωi=Classify(S0'),i=1,2,…,n}.步骤6 将每类样本点作为训练集分别训练局部模型，获得n个局部模型LN={n1,n2,…,nn}，LN={ni | ni=N(Si),i=1,2,…,n}，式中N(Si)表示采用神经网络模型对样本集进行训练．通过组合局部神经网络模型生成隐式状态方程右端式的高精度全局神经网络模型．当基于全局神经网络模型求解动态优化问题的过程时，通过高斯混合模型对输入向量进行归类，确定了输入向量的类别也即确定了该向量所属的局部区域，进而调用该局部神经网络ni近似计算该输入向量所对应的雅可比矩阵信息．2.2　示例说明为了更好阐述基于高斯混合模型聚类的空间划分策略，以下式为例，对上述流程进行具体说明，f(x1,x2)=12exp(-x12-(x2-2.5)2)+12exp(-x12-(x2+2.5)2)-12exp(-(x1+2.5)2-x22)-12exp(-(x1-2.5)2-x22).首先，通过基于最大曲率最小点距的序列采样方法获取初始样本点．如图2所示，该函数在设计域上有四个非线性程度较高的区域，初始样本点主要集中在这些区域上．然后，通过高斯混合模型聚类算法对初始样本点进行聚类．如图3所示，初始样本点主要被划分为四类，且每类均涵盖了一个非线性程度较高的区域．进而，为清晰表示高斯混合模型对状态空间的划分，本示例通过拉丁超立方采样方法获取批量样本点．基于高斯混合模型对样本点进行归类，归类结果如图4所示，整个状态空间被划分为四个不规则的局部区域，分别为Ω1,Ω2，10.13245/j.hust.238733.F002图2初始样本点分布图10.13245/j.hust.238733.F003图3初始样本点分类示意图10.13245/j.hust.238733.F004图4基于高斯混合模型的空间划分结果Ω3,Ω4．最后，在局部区域中分别构建神经网络模型，神经网络的节点结构均为[2,5,5,1]．为进行比较，使用相同数量的样本点构建基于维度划分的空间划分模型(策略1)，两种不同划分策略下各模型精度见表1．10.13245/j.hust.238733.T001表1不同划分策略下各局部模型精度策略Ω1Ω2Ω3Ω4策略119.507 042.506 010.130 07.590 3策略25.500 70.781 72.031 56.368 110-4由图2和图4可知：本研究提出的基于高斯混合模型聚类的空间划分策略当对示例函数的状态空间进行划分时，边界经过状态方程非线性程度较低、梯度较小的区域，降低了局部模型的建模复杂度．如表1所示，在样本点数量相同的情况下，基于高斯混合模型聚类的空间划分策略(策略2)所构建的局部模型精度更高．3 工程算例3.1　2自由度机械臂最短时间控制问题2自由度机械臂最短时间控制问题[17]是一个经典的动态优化问题，该问题的目标是寻找最优控制变量u*(t)使得机械臂以最短时间从指定起点到达指定终点．2自由度机械臂系统的仿真模型可在MWorks平台上通过可视化拖放方式进行构建，仿真模型中隐含了该系统的状态方程．同时，该系统的非线性动力学模型为minu1,u2J=tf，s.t.  v˙1=sinα94(cosα)v12+2v22+43(u1-u2)-32(cosα)u2/ 3136+94(sin2α);v˙2=-sinα72v12+94(cosα)v22-73u2+32(cosα)(u1-u2)/ 3136+94(sin2α);α˙=v2-v1;β˙=v1;v1(0)=0;v2(0)=0;α(0)=0.5;β(0)=0;v1(tf)=0;v2(tf)=0;α(tf)=0.5;β(tf)=0.522;-1≤u1≤1;  -1≤u2≤1,式中：v1和v2为肩部和肘部连杆的速度；α和β为相应的角度；u1和u2为控制变量．基于动力学模型对本算例中的动态优化问题进行求解，得到目标函数的精确解为2.982 3 s．根据系统中速度与角度之间的关系，角加速度α˙和β˙的右端式较为简单，无须近似，因此只近似加速度v˙1和v˙2的右端式．通过基于高斯混合模型的空间划分策略将状态空间划分为多个局部区域，并对局部模型进行训练，分别得到两个均包含2个局部神经网络模型的代理模型，其中局部神经网络的节点结构均为[5,15,15，1]．得到代理模型后，调用SNOPT求解器对该问题进行求解，求解过程中进行了6次离散时间网格自适应更新，最后一次更新后的网格点数目为71，求得该动态优化问题目标函数的近似解为2.982 7 s．因此，近似解与精确解中最优目标函数值的误差为0.013 4%．同时，近似解和精确解中部分状态轨迹和控制曲线如图5所示(ve，αe和ue分别为状态分量和控制分量近似解与精确解之间的误差)，可见近似解虽是次优，但与精确解的误差在允许范围内．10.13245/j.hust.238733.F005图5机械臂系统优化问题的精确解与近似解同时，本研究构建了2自由度机械臂系统状态方程右端式的单个全局模型和基于维度划分的多局部模型．其中基于维度划分的多局部模型主要根据某一维度对状态空间进行划分，为便于比较，也将状态空间划分为两个局部区域．基于不同划分策略的代理模型方法求解得到的近似结果见表2(策略3表示不采用空间划分策略，只构建一个全局模型；模型精度为多个局部模型的平均精度)，可见本研究提出的基于高斯混合模型聚类的空间划分策略所训练的模型精度最高，且基于该模型求解得到的近似解与精确解的误差最小．10.13245/j.hust.238733.T002表2机械臂系统优化问题在不同划分策略下的模型精度及求解结果参数策略1策略2策略3模型精度/10-43.460 42.743 53.060 0求解结果/s2.983 12.982 72.982 83.2　悬挂式滑翔机终端水平距离最大化问题悬挂式滑翔机终端水平距离最大化问题[18]的优化目标是通过控制气动升力系数CL，使得滑翔机在终端时间tF时刻的水平位移xF最大．该问题的仿真模型可在MWorks平台上通过可视化拖放方式进行构建，状态方程隐含于仿真模型中．同时，滑翔机平面运动的状态方程为x˙=vx;y˙=vy;v˙x=(-Lsin η-Dcos η)/m;v˙y=(Lcos η-Dsin η)/m-g,式中：x为水平距离；y为海拔高度；vx为水平速度；vy为垂直速度；L=CLρSvr2/2；D=CDρSvr2/2； CD=C0+kCL2； X=(x/R-2.5)2； ua(x)=uM(1-X)exp(-X)； Vy=vy-ua(x)； vr=vx2+Vy2；sin η=Vy/vr，cos η=vx/vr．模型中的常量为：uM=2.5,  m=100 kg,  R=100,  S=14 m2,C0=0.034,  ρ=1.13 kg/m3,k=0.069 662,  g=9.806 65 m/s2.同时，该动态优化问题的约束主要包括状态变量的边界约束和控制变量的区间约束为：0≤CL≤1.4,  x(0)=0 m,  x(tF)取任意值,y(0)=1 000 m,  y(tF)=900 m,vx(0)= 13.227 567 5 m/s, vx(tF)= 13.227 567 5 m/s,vy(0)=-1.287 500 52 m/s, vy(tF)=-1.287 500 52 m/s.  根据速度与距离之间的关系，状态方程中x˙和y˙的右端式不须要近似，只须近似v˙x和v˙y的右端式．通过本研究提出的方法，得到两个分别包含5个和4个局部神经网络模型的代理模型，其中局部神经网络的节点结构均为[5,20,20,1]．基于近似模型对该问题进行求解，求解过程中进行了8次离散时间网格自适应更新，最后一次更新后的网格点数目为805，并得到该滑翔机最大水平位移xF为1 248.068 5 m，终端时间x(tF)为98.421 8 s．而通过滑翔机平面运动的状态方程求解得到的最大水平距离为1 248.031 m，终端时间为98.437 s．因此，近似解与精确解中最优目标函数值的误差为0.003 0%．滑翔机部分状态轨迹和控制曲线的近似解和精确解如图6所示，可见近似解和精确解之间的误差在允许范围内．10.13245/j.hust.238733.F006图6滑翔机系统优化问题的精确解与近似解同时，本研究构建了滑翔机系统状态方程右端式的单个全局模型和基于维度划分策略的多局部模型．各划分策略所得到的近似模型精度及基于近似模型求解得到的结果见表3，可见本算例中局部模型的精度要高于全局模型，且基于高斯混合模型聚类的空间划分策略训练得到的模型精度最高，求解得到的近似解也与精确解的误差最小．10.13245/j.hust.238733.T003表3滑翔机系统优化问题在不同划分策略下的模型精度及求解结果参数策略1策略2策略3模型精度/10-58.319 14.147 821.485 0求解结果/m1 248.137 31 248.068 51 248.311 24 结语为实现基于代理模型的方法高效求解状态方程隐式表达的动态优化问题，本研究提出了一种基于高斯混合模型聚类的空间划分策略，根据状态方程非线性程度将状态空间划分为多个不规则的局部区域，并采用神经网络构建每个局部区域的代理模型．通过组合所有局部区域的代理模型获得状态方程右端式的全局神经网络模型，进而基于全局神经网络模型，采用直接法对动态优化问题进行求解．两个工程算例中，基于所提方法得到最优目标函数的近似值与精确值的误差分别为0.013 4%和0.003 0%，这表明基于高斯混合模型聚类的空间划分策略降低了局部区域的建模复杂度，提高了局部模型的精度；同时，降低了代理模型在边界两边不同局部区域之间估值的不连续性，保证了基于代理模型方法求解动态优化问题的精度．