摘 要 目前国际上对动态优化问题中的状态变量路径约束已有一些研究, 但专门处理控制变量路径约束的方法却鲜见报道.本文首先介绍两种分别基于三角函数变换、约束算子截断来处理控制变量路径约束的方法, 然后提出一种基于光滑化的二次罚函数方法. 光滑化罚函数方法不仅能够处理控制变量路径约束, 而且还能同时处理关于状态变量的路径约束. 最后使用目前流行的控制变量参数化(Control variable parameterization, CVP) 策略对最终获得的、不再含控制变量路径约束的动态优化问题求解. 实例测试一展现了三种方法各自的特点; 实例测试二表明了光滑罚函数方法的有效性和优越性.
关键词 动态优化, 路径约束, 控制变量, 控制变量参数化引用格式
Methods to Deal with Control Variable Path Constraints inDynamic Optimization Problems
Abstract Many researchers have reported their methods to handle the state variable path constraints in dynamicoptimization problems. However, very few of them mentioned the techniques to deal with the control variable pathconstraints. This paper ¯rstly introduces two ways to eliminate the control variable path constraints, then presents asmoothed quadratic penalty function method which can optimize the dynamic problems with control and state variablepath constraints at the same time. The ¯nal transformed problem is solved by control variable parameterization (CVP)strategy in the study. Test on Case 1 demonstrates the characteristics of the three techniques; test on Case 2 indicatesthe e®ectiveness and the superiority of the proposed smoothed penalty function method.
Key words Dynamic optimization, path constraint, control variable, control variable parameterization (CVP)Citation
　　在工业过程的优化领域里, 动态优化算法的研究成为当前的热点和趋势[1]. 自从动态优化的概念提出以来, 由于其巨大的应用价值, 已经引起了广泛关注[2¡5], 并在化学工程、系统生物工程、清洁能源、企业级管理优化等方面逐渐得到应用[6¡9].
　　工业过程在严格意义上是状态变量随时间的演进、空间的转移而发生改变的动态过程. 动态过程由微分方程或差分方程描述, 称为动态模型. 动态优化就是对动态模型中的操作变量实施控制, 使得过程的性能指标达到最优. 实际的工业过程常常需要对某些状态变量进行约束(例如操作过程中压力不能超过指定值).如果约束只在过程的某一时刻要求满足, 则称为点约束; 如果在一段时间内都要求满足, 则称为路径约束.
　　当动态优化问题含有关于控制变量和状态变量的路径约束时, 其求解算法是一个难点[10]. 目前国际上对动态优化问题中的状态变量路径约束已有一些研究[11¡12], 但专门处理控制变量路径约束的方法却鲜见报道. 本文首先介绍两种分别基于三角函数变换、约束算子截断来处理控制变量路径约束的方法, 然后新提出一种基于光滑化的二次罚函数方法. 最后使用控制变量参数化(Control variableparameterization, CVP)策略对最终获得的、不再含控制变量路径约束的动态优化问题求解. 在实例测试中展示了这三种方法各自的特点, 以及所提出的光滑化罚函数方法的优点.
　　1 带控制变量路径约束的动态优化问题描述
　　动态优化就是对动态过程模型中的操作变量实4期 胡云卿等: 处理动态优化问题中控制变量路径约束的方法 441施控制, 使性能指标达到最优. 在动态优化问题中, 动态模型往往是由一系列常微分方程(Ordinarydi®erential equations, ODE)或微分代数方程(Dif-ferential and algebraic equations, DAE) 共同组成的微分代数方程组, 并作为优化目标的约束条件而存在. Pantelides[13]指出高阶DAE方程组都可以通过降阶变为一阶DAE方程组, 而一阶DAE方程组与ODE方程组可以等价转换. 最终得到带路径约束(包括关于控制变量和状态变量)的一阶动态优化问题的数学模型(问题1)如下:min J= ©[x(tf); tf] +Ztft0L[x(t); u u(t); t]dts:t: _ x(t) =f[x(t); u u(t); t]x(t0) =x0a(t)·u(t)·b(t)'[x(t); t]·0t2[t0; tf]其中, 目标函数J 由终端时刻tf 的性能泛函©[x(tf); tf] 和固定时间段t 2[t0; tf] 内的性能泛函Rtft0L[x(t); u u(t); t]dt 共同组成, x(t) 2Rn是n维状态变量, x0 是状态变量x(t) 在t0 时刻的初始值, u(t)2Rr是r 维控制变量, f 表示连续可微的n维向量函数, b(t)和a(t)分别是控制变量路径约束的上下界, 且满足ai(t)·bi(t) (i = 1;2;¢ ¢ ¢; r),'2Rq(q·n)表示状态变量不等式路径约束.
　　2 算法描述
　　2.1 三角函数变换法
　　三角函数变换法由张兵等[14]提出. 该方法通过引入中间变量z(t) 2Rr和三角函数代换, 将原来的控制变量u(t)转化为关于中间变量z(t)的连续、有界三角函数形式. 对于上下界分别为b(t)和a(t)的控制变量u(t) , 按照下式进行转换:ui(t) = 0:5[bi(t)¡ai(t)]fsin[zi(t)] + 1g+ai(t); i = 1;2;¢ ¢ ¢; r(1)注意z(t)与u(t)的不同在于z(t)的取值范围无限制, 优化时将z(t)作为决策变量, 优化完成后再用最优决策变量z¤(t)通过关系式(1) 获得u¤(t). 这样就避开了关于控制变量u(t)的路径约束. 当z(t)超出[¡¼; ¼]时,容易发生振荡,不利于收敛. 因此需要对相邻两点之间的步长加以限制, 使得在t 2[t0; tf]内均满足:°z(k+1)i(t)¡z(k)i(t)°·¼M; i = 1;2;¢ ¢ ¢; r (2)其中, k为迭代次数, M为正整数, 且实测表明最佳取值范围为M2[5;8]
　　.三角函数变换法是一种有效的处理控制变量路径约束的方法, 但无法处理状态变量的路径约束, 需要和专门处理状态变量路径约束的方法结合使用.
　　2.2 约束算子法
　　如果控制变量受到的路径约束为a(t)·u(t)·b(t), 其中b(t) 和a(t) 是关于t 的函数并且满足ai(t)·bi(t) (i = 1;2;¢ ¢ ¢; r),则在CVP策略中可以被离散化下列形式:ai(tj)·ui(tj)·bi(tj); j = 1;2;¢ ¢ ¢; N (3)定义约束算子c 2Rr, 迭代过程中对每次优化的结果都进行如下约束:ciui(tj) =8><>:ai(tj); ui(tj)·a(tj)ui(tj); a(tj)< ui(tj)< b(tj)bi(tj); ui(tj)¸b(tj)(4)该方法被称为约束算子法[15]. 约束算子法便于程序实现, 计算量很小, 便于处理复杂的控制容许集. 但该方法也只能用来处理控制变量的路径约束,无法同时处理状态变量路径约束.
　　2.3 光滑化罚函数方法
　　基于惩罚函数的思想, 问题1可以表述成无约束形式(问题2):^J = J+½Ztft0q Xi=1maxf0; 'i[x(t); t]g2dt+½Ztft0r Xi=1maxf0; ai(t)¡ui(t)g2dt+½Ztft0r Xi=1maxf0; ui(t)¡bi(t)g2dt(5)其中½是惩罚因子, ½Rtft0Pqk=1maxf0; 'i[x(t); u u(t);t]g2dt 是对状态变量路径约束的惩罚项, ½Rtft0Pri=1maxf0; ai(t) ¡ ui(t)g2dt 和½Rtft0Pri=1maxf0; ui(t)¡bi(t)g2dt 是对控制变量路径约束的惩罚项. 罚函数方法的原理是: 当½! 1时, 无约束问题2的目标函数值^J将收敛到原来带约束问题1的最优解J¤.
　　用惩罚函数法求解带约束的动态优化问题, 最大的难点在于非光滑的max运算符的处理[14]. 本文引入一种max运算符号的光滑化处理方法[16]: 设gi[x(t); u u(t); t] 为路径约束(gi 既可以是关于控制变量的路径约束也可以是关于状态变量的路径约束),442 自 动 化 学 报 39卷引进光滑化函数Á(gi; ®) =12(gi +pg2i+®) (6)它具有lim®!0+Á(gi; ®) ¼maxf0; gi[x(t); u u(t); t]g的良好性质. 所以, 当®!0+时得到^J的光滑化表达式为J=J+½Ztft0q+2r Xi=1Á(gi; ®)2dt¼^J (7)
　　CVP是求解动态优化问题的一种重要策略, 在该策略中, 控制变量u(t) 被离散为N段(称为控制参数, 记作u2RrN), 而状态变量仍保持连续形式. CVP 策略中包括两个主要步骤: 当状态变量x(k)(t) 已知时, 可以通过数值优化方法求得新的控制参数u(k+1); 当控制参数更新为u(k+1)时, 又可以通过解ODE方程组或DAE方程组获得新的状态变量x(k+1)(t). 交替执行这两个主要步骤最终可以获得最优控制变量u¤(t)的近似离散值u¤. 关于CVP策略的原理和实现步骤请参考Vassiliadis[17].利用CVP策略对问题1中的控制变量进行离散化, 得到问题3:min J= ©[x(tf); tf] +Ztft0L[x(t); u uu; t]dts:t: _ x(t) =f[x(t); u uu; t]x(t0) =x0a(tj)·uj ·b(tj)'[x(t); t]·0t2[t0; tf]
　　根据式(5) 中的二次惩罚函数法将关于控制变量和状态变量的路径约束增广到目标函数, 则问题3可以化为问题4:min^J=J+½Ztft0q Xi=1maxf0; 'i[x(t); t]g2dt+½Ztft0r Xi=1maxf0; ai(t)¡uig2dt+½Ztft0r Xi=1maxf0; ui ¡bi(t)g2dts:t: _ x(t) =f[x(t); u uu; t]x(t0) =x0t2[t0; tf]再根据式(6) 对问题4的目标函数进行光滑化处理, 进一步得到问题5:min¹J=J+½Ztft0q+2r Xi=1Ái2dt¼^Js:t: _ x(t) =f[x(t); u uu; t]x(t0) =x0t2[t0; tf]其中, Ái =Áfgi[x(t); u uu; t]; ®g. 值得注意的是: 1)问题5中不再含有路径约束; 2) 所有约束实质上构成关于x(t)的ODE方程组, 而目标函数则是关于控制参数u、惩罚因子½以及光滑化因子®的数学规划问题, 可以直接利用CVP策略求解.将问题4的解记作¹ u¤, 由于问题4的目标函数是对问题3目标函数的近似, 所以¹ u¤也是对u¤的近似, 其近似程度与惩罚因子½以及光滑化因子®有关. 下面证明¹ u¤与u¤的误差关系并给出½和®的一种有效更新方式.
　　定义1.对于问题3,如果控制参数u®满足gi[u®; x(u®; t); t]·®; i= 1;2;¢ ¢ ¢ ; q+2r (8)则称u®是问题3的®- 可行解.定理1. 设® > 0, 对于光滑化函数Áfgi[x(t); u uu; t]; ®g则有0< Áfgi[x(t); u uu; t]; ®g2¡maxfgi[x(t); u uu; t];0g2<2®2(9)证明.Áfgi[x(t); u uu; t]; ®g2¡maxfgi[x(t); u uu; t];0g2=Áfgi[x(t); u uu; t]; ®g2¡(0; gi[x(t); u uu; t]<0gi[x(t); u uu; t]2; gi[x(t); u uu; t]¸0=8><>:12gi[x(t); u uu; t]2+12gi[x(t); u uu; t]£qgi[x(t); u uu; t]2+ 4®2+®2; gi[x(t); u uu; t]<012gi[x(t); u uu; t]2+12gi[x(t); u uu; t]£qgi[x(t); u uu; t]2+ 4®2+®2; gi[x(t); u uu; t]¸0当gi[x(t); u uu; t]<0且® >0时, 可得:12gi[x(t); u uu; t]2+12gi[x(t); u uu; t]qgi[x(t); u uu; t]2+4®2+®2>12gi[x(t); u uu; t]2+12gi[x(t); u uu; t]qgi[x(t); u uu; t]2= 0当gi[x(t); u uu; t]¸0且® >0时, 可得:¡12gi[x(t); u uu; t]2+12gi[x(t); u uu; t]qgi[x(t); u uu; t]2+4®2+®2>¡12gi[x(t); u uu; t]2+12gi[x(t); u uu; t]qgi[x(t); u uu; t]2= 04期 胡云卿等: 处理动态优化问题中控制变量路径约束的方法 443因此:Áfgi[x(t); u uu; t]; ®g2¡maxfgi[x(t); u uu; t];0g2>0 ¤另一方面, 当gi[x(t); u uu; t] <0且® >0时,12gi[x(t); u uu; t]2+12gi[x(t); u uu; t]qgi[x(t); u uu; t]2+ 4®2+®2<2®2可等价表示为gi[x(t); u uu; t]2<2®2¡gi[x(t); u uu; t]qgi[x(t); u uu; t]2+ 4®2由于2®2¡gi[x(t); u uu; t]qgi[x(t); u uu; t]2+ 4®2>0,不等号两边同时平方, 得到:0< ®2+gi[x(t); u uu; t]2¡gi[x(t); u uu; t]£qgi[x(t); u uu; t]2+ 4®2容易验证上式恒成立; 当gi[x(t); u uu; t] ¸0且® >0时, 可等价表示为gi[x(t); u uu; t]qgi[x(t); u uu; t]2+ 4®2<2®2+gi[x(t); u uu; t]2由于2®2+gi[x(t); u uu; t]2>0,容易验证上式也恒成立; 因此:Áfgi[x(t); u uu; t]; ®g2¡maxfgi[x(t); u uu; t];0g2<2®2¤定理2. 设½ >0,且® >0对于光滑化函数Áfgi[x(t); u uu; t]; ®g则有0<¹J(uu; ½; ®)¡^J(uu; ½)<2½(q+2r)®2(tf ¡t0)(10)证明. 根据定理1,如果® >0,则有:0< Áfgi[x(t); u uu; t]; ®g2¡maxfgi[x(t); u uu; t];0g2<2®2对下标i 求和, 得到:0<q+2r Pi=1Áfgi[x(t); u uu; t]; ®g2¡q+2r Pi=1maxfgi[x(t); u uu; t]; t];0g2<2(q+2r)®2由于½ >0,可以进一步得到:0 < ½Ztft0q+2r Xi=1Áfgi[x(t); u uu; t]; ®g2dt¡½Ztft0q+2r Xi=1maxfgi[x(t); u uu; t]; t];0g2dt <2½(q+2r)®2(tf ¡t0)上式可以等价写为0<¹J(uu; ½; ®)¡^J(uu; ½)<2½(q+2r)®2(tf ¡t0)¤定理3. 设u¤是问题3的解, ¹ u¤是问题5的解, 则0<¹J(¹ u¤; ½; ®)¡^J(u¤; ½)<2½(q+2r)®2(tf ¡t0)(11)证明. 根据定理2,可知:0<¹J(¹ u¤; ½; ®)¡^J(¹ u¤; ½)<2½(q+2r)®2(tf ¡t0)以及0<¹J(u¤; ½; ®)¡^J(u¤; ½)<2½(q+2r)®2(tf ¡t0)由于u¤是问题3的解, 可得:¹J(¹ u¤; ½; ®)¡^J(¹ u¤; ½)·¹J(¹ u¤; ½; ®)¡^J(u¤; ½)¹ u¤是问题5的解, 可得:¹J(¹ u¤; ½; ®)¡^J(u¤; ½)·¹J(u¤; ½; ®)¡^J(u¤; ½)联立以上两式可得:0 <¹J(¹ u¤; ½; ®)¡^J(¹ u¤; ½)·¹J(¹ u¤; ½; ®)¡^J1(u¤; ½)·¹J(u¤; ½; ®)¡^J(u¤; ½)<2½(q+2r)®2(tf ¡t0)因此:0<¹J(¹ u¤; ½; ®)¡^J(u¤; ½)<2½(q+2r)®2(tf ¡t0)¤定理4. 设u¤是问题3的解, ¹ u¤是问题5的解, 进一步, 设u¤是问题3的®- 可行解且¹ u¤是问题3的®- 可行解, 则有:¡3 +p52½(q+2r)®2(tf ¡t0)< J(¹ u¤)¡J(u¤)<2½(q+2r)®2(tf ¡t0) (12)证明. 由于u¤是3的可行解, 可知:q+2r Xi=1maxfgi[x(u¤;t); u uu; t];0g2= 0又由于¹ u¤是问题3的®- 可行解, 从式(6) 和式(8)可知:0< Áfgi[x(¹ u¤; t);¹ u¤; t]; ®g2·3 +p52®2444 自 动 化 学 报 39卷因此有:0< ½Ztft0q+2r Xi=1Áfgi[x(¹ u¤; t);¹ u¤; t]; ®g2dt·3+p52½(q+2r)®2(tf ¡t0)从定理3可知:0< J(¹ u¤) +½Ztft0q+2r Xi=1Áfgi[x(¹ u¤; t);¹ u¤; t]; ®g2dt¡J(u¤)¡½Ztft0q+2r Xi=1maxfgi[x(u¤; t);u u¤; t];0g2dt <2½(q+2r)®2(tf ¡t0)所以有:¡3+p52½(q+2r)®2(tf ¡t0)< J(¹ u¤)¡J(u¤)<2½(q+2r)®2(tf ¡t0)¤定理3说明: 当½(q+2r)®2(tf ¡t0)足够小时,将问题5的解¹ u¤(通过CVP策略可以求得)代入问题4,将使得问题4和问题5的目标函数值近似相等. 定理4说明: 当½(q+2r)®2(tf ¡t0) 足够小且¹ u¤是问题3的®- 可行解时, 问题3和5的目标函数值近似相等.Fiacco等已经证明[15], 可以通过逐渐增加½值的方式求解一系列问题4, 使获得的¹ u¤逐渐成为问题3的可行解. 根据定理3进一步可知, 当½(q+2r)®2(tf ¡t0) 足够小时, 通过逐渐增加½值的方式求解一系列问题5,获得的¹ u¤也将逐渐逼近wenti 3的可行解.由以上分析, 可以导出光滑化二次罚函数法的实现结构, 如图1所示, 其中k是循环次数.图1 光滑罚函数方法的数值实现结构Fig. 1 Implementation structure of the smoothedpenalty method实际计算时, 设J(¹ u¤)和J(u¤)的容许误差为", 且循环过程中½的增大系数½和®的缩小系数d满足cd= 1的关系, 根据式(12) 可知: 当光滑化因子按照式(13)更新时, 最终将满足jJ(¹ u¤)¡J(u¤)j · ".®(k)=s2"(3 +p5)½(k)(q+ 2r)(tf ¡t0)(13)光滑化罚函数方法可以同时处理多个控制变量和状态变量的路径约束, 适合大规模复杂动态优化问题求解. 但这种方法也有缺点: 需要更新惩罚因子和光滑化因子, 并求解一系列的问题5,所以耗时相对较多.
　　3 实例测试
　　本文采用的CVP策略中, 控制变量u(t) 的离散化段数为100 (即N= 100),使用带预测校正的4步4阶Adams法求解状态微分方程组, 使用共轭梯度法求解关于控制参数的优化问题.
　　3.1 一个有解析解的实例问题6的数学模型
　　如下:min J=Z10[x(t)¡12u(t)]dts:t: _ x(t) =¡x(t) +u(t)x(0) = 1:0ju(t)j · 1:0t2[0;1]许多最优控制类书籍都给出了其理论解, 但无数值解. 问题6的理论最优目标函数值保留六位有效数字后为J¤= 0:457388,对应的最优控制策略为u¤(t) =(¡1; 0·t·lne21; lne2·t·1(其中lne2¼0:3)用离散段数为N= 100、误差容限为"= 10¡6的CVP策略求解, 采用三角函数变换法对控制变量的路径约束ju(t)j · 1:0进行处理, 得到的最优目标函数值J¤= 0:457474, 精确到了千分位. 相应的最优控制轨迹和最优状态轨迹分别如图2 (a)和图2 (b)所示. 从图2 (a)可以看出, 控制变量在t 2[0;1]内很好地满足了路径约束, 但在t = lne2时刻, 控制变量的跳变过程并非瞬时完成, 而是需要一小段过渡时间. 此外, 如果跳变过程以t = lne2为对称中心, 其两端的控制曲线也不完全对称.用离散段数为N= 100、误差容限为"= 10¡6的CVP策略求解, 其中对控制变量的路径约束ju(t)j · 1:0采用约束算子法进行处理, 得到的最优目标函数值J¤= 0:457389, 精度提高到了十万分位. 相应的最优控制轨迹和最优状态轨迹分别如图3 (a)和图3 (b)所示. 从图3 (a)可以看出, 控制变量在内很好地满足了路径约束ju(t)j · 1:0. 在t = lne2时刻, 控制变量的跳变过程比较匀称, 需要的过渡时间也比图2(a)更短.4期 胡云卿等: 处理动态优化问题中控制变量路径约束的方法 445图2 问题6,第2.1节方法Fig. 2 Problem 6, Method in Section 2.1采用光滑罚函数方法进行处理, 得到的目标函数表达式为J = J+½Ztft0Áf[u(t)¡1]; ®g2dt+½Ztft0Áf¡[u(t)¡1]; ®g2dt(14)用离散段数为N= 100、误差容限为"= 10¡6的CVP策略求解, 初始惩罚因子设置为½(0)= 0:1,增大系数为c = 5. 当惩罚因子½= 312:5,光滑化因子®= 10¡8时算法停止, 得到的最优目标函数值J¤= 0:458111. 图4和图5显示了当惩罚因子½= 312:5时, 最优控制轨迹和最优状态轨迹随光滑化因子®变化的情况. 从图3 (a)可以看出, 当算法停止时, 控制变量在t 2[0;1]内满足路径约束ju(t)j · 1:0.在t= lne2时刻, 控制变量的跳变过程非常匀称, 但需要较长的过渡时间.表1列出了在AMD Turion64/1G DDR2计算平台上, 分别用三种控制变量路径约束的处理方法处理问题5得到的结果. 其中过渡时间指u¤(t)的值从¡1变化到1所需的时间(在N= 100时精确到0.025 s).图3 问题6,第2.2节方法Fig. 3 Problem 6, Method in Section 2.2表1 问题6的求解结果对比Table 1 The comparison results of Problem 6控制约束处理方法 过渡时间(s) 目标函数值时间消耗(s)三角函数变换法 0.1 (0.025£4) 0.457474 0.06约束算子法 0.05 (0.025£2) 0.457389 0.03光滑化罚函数法 0.075 (0.025£3) 0.457615 0.563.2 间歇反应饱和度控制问题管式间歇反应器中发生并行反应: A→B和A→C, 反应速率分别为k1、k2, B是目标产物, C是副产物, 目标是控制饱和度使得产物B在最终时刻的产量最大. 用x1(t)、x2(t)分别表示A、B的浓度分率, u(t)表示饱和度, 即控制变量, u(t) =k1L=v(L为反应器长度, v为阻塞流率),该过程为问题7,可描述为446 自 动 化 学 报 39卷图4 控制曲线的惩罚过程(问题6,第2.3节方法)Fig. 4 The penalty process of control pro¯le (Problem 6, Method in Section 2.3)图5 状态曲线的惩罚过程(问题6,第2.3节方法)Fig. 5 The penalty process of state pro¯le (Problem 6, Method in Section 2.3)4期 胡云卿等: 处理动态优化问题中控制变量路径约束的方法 447图6 控制曲线的惩罚过程(问题7,第2.3节方法)Fig. 6 The penalty process of control pro¯le (Problem 7, Method in Section 2.3)max x2(tf)s:t: _ x1(t) =¡[u(t) +u2(t)2]x1(t)_ x2(t) =u(t)x1(t)x1(0) = 1:0x2(0) = 0:00:0·u(t)·5:0x1(t)¸0:2t2[0;1]这是一个既包含控制变量路径约束又包含状态变量路径约束的动态优化问题, 只能用本文所述三种方法中的光滑化罚函数法求解, 得到的目标函数表达式为J=J+½Ztft0Áf[u(t)¡5]; ®g2dt+½Ztft0Áf¡u(t); ®g2dt+½Ztft0Áf[¡x1(t) + 0:2]; ®g2dt (15)仍然采用离散段数为N= 100、误差容限为" = 10¡6的CVP策略求解, 初始惩罚因子设置为½(0)= 1,增大系数为c = 10. 当惩罚因子½ = 100 000,光滑化因子®= 3:5682£10¡8时算法停止, 得到的最优目标函数值J¤= 0:536243.图6和图7显示了控制轨迹和状态轨迹随惩罚因子和光滑化因子变化的情况. 当算法停止时, 从图6可以看出, 控制变量在t 2[0;1]内满足路径约束0:0·u(t)·5:0;从图7可以看出, 控制变量在t2[0;1]内满足路径约束x1(t)¸0:2.448 自 动 化 学 报 39卷图7 状态曲线的惩罚过程(问题7,第2.3节方法)Fig. 7 The penalty process of state pro¯le (Problem 7, Method in Section 2.3)由于计算机平台、控制变量参数化策略中的子算法和离散段数不同, 优化过程的消耗时间并不能作为比较算法优劣的依据, 只列出作为参考, 但仍然能体现优化过程是高效的. 表2列出了Bloss &Biegler[10]以及本文对问题6的求解结果. 可以看出, 0.5400 是目前文献中的最优值, 而本文用光滑罚函数方法获得的最优目标函数值与之也很接近.表2 问题7的求解结果对比Table 2 The comparison results of Problem 7研究者 约束处理方法 目标函数值时间消耗(s)Bloss & Biegler 外点罚函数法 0.5362 5.3Bloss & Biegler Kreisselmeier-Steinhauser 0.5400* 17.5本本本 文文文 (N= 100) 光光光 滑滑滑化化化 罚罚罚 函函函 数数数 法法法 0.536243 2.32
　　4 结论
　　本文介绍了三角函数变换、约束算子、光滑化罚函数这三种处理动态优化问题中控制变量路径约束的方法. 其中三角函数变换法、约束算子法较为简单, 但仅能处理关于控制变量的路径约束; 光滑化罚函数法能处理同时包含关于控制变量和状态变量路径约束的动态优化问题, 但需要一定时间得到高精度结果. 对于问题6,在相同的CVP策略下由光滑化罚函数法获得的目标函数值0.457615最接近理论值; 对于问题7 | 间歇式管式反应器的饱和度控制问题, 利用光滑罚函数法同时处理关于控制变量和状态变量的路径约束, 获得的最优目标值0.536243同样与目前文献中的最优值0.5400非常接近, 说明了该方法的有效性, 及其在使用范围上的优越性.
　　References1 Biegler L T, Grossmann I E. Retrospective on optimiza-tion. Computers and Chemical Engineering, 2004, 28(8):1169¡11922 Asgari S A, Pishvaie M R. Dynamic optimization in chemicalprocesses using region reduction strategy and control vec-tor parameterization with an ant colony optimization algo-rithm.Chemical Engineering and Technology, 2008, 31(4):507¡5123 Balku S, Yuceer M, Berber R. Control vector parameteriza-tion approach in optimization of alternating aerobic-anoxicsystems. Optimal Control Applications and Methods, 2009,30(6): 573¡5844 Gugat M, Herty M. The smoothed-penalty algorithm forstate constrained optimal control problems for partial di®er-ential equations.Optimization Methods and Software, 2010,25(4): 573¡5994期 胡云卿等: 处理动态优化问题中控制变量路径约束的方法 4495 Vega M P, Mancini M C, Calcada L A. Multi-objective dy-namic optimization of ¯xed bed dryers: simulation and ex-periments. In: Proceedigns of the 19th European Sympo-sium on Computer Aided Process Engineering. Amsterdam:Elsevier Science, 2009. 147¡1526 Hirmajer T, Balsa-Canto E, Banga J R. DOTcvpSB, a soft-ware toolbox for dynamic optimization in systems biology.Bmc Bioinformatics, 2009, 107 Fikar M, Kovacs Z, Czermak P. Dynamic optimization ofbatch dia¯ltration processes.Journal of Membrane Science,2010,355(1¡2): 168¡1748 Ni B C, Sourkounis C. Stochastic dynamic optimization forwind energy converters.International Journal for Compu-tation and Mathematics in Electrical and Electronic Engi-neering, 2011, 30(1): 265¡2799 Biegler L T. Large-scale nonlinear programming: an inte-grating framework for enterprise-wide dynamic optimiza-tion. In: Proceedings of the 17th European Symposium onComputer Aided Process Engineering. Amsterdam: ElsevierScience, 2007. 575¡58210 Bloss K F, Biegler L T, Schiesser W E. Dynamic process op-timization through adjoint formulations and constraint ag-gregation.Industrial and Engineering Chemistry Research,1999,38(2): 421¡43211 Bell M L, Sargent R W H. Optimal control of inequalityconstrained DAE systems.Computers and Chemical Engi-neering, 2000, 24(11): 2385¡240412 Luus R. Handling inequality constraints in optimal con-trol by problem reformulation. Industrial and EngineeringChemistry Research, 2009, 48(21): 9622¡963013 Pantelides C C. The mathematical modelling of transientsystem using di®erential-algebraic equations. Computersand Chemical Engineering, 1988, 1214 Zhang Bing, Chen De-Zhao, Wu Xiao-Hua. Graded op-timization strategy and its application to chemical dy-namic optimization with ¯xed boundary.Journal of Chemi-cal Industry and engineering Society of China, 2005, 56(7):1276¡1280(张兵, 陈德钊, 吴晓华. 分级优化用于边值固定的化工动态优化问题. 化工学报, 2005, 56(7): 1276¡1280)15 Fiacco A V, McCormick G P.Nonlinear Programming: Se-quential Unconstrained Minimization Techniques. Philadel-phia: SIAM, 199016 Chen X, Nashed Z, Qi L. Smoothing methods andsemismooth methods for nondi®erentiable operator equa-tions. SIAM Journal on Numerical Analysis, 2001, 38(4):1200¡121617 Vassiliadis V S. Computational Solution of Dynamic Opti-mization Problems with General Di®erential-algebraic Con-straints [Ph. D. dissertation], University of London, UK,1993