2.845

2023影响因子

(CJCR)

  • 中文核心
  • EI
  • 中国科技核心
  • Scopus
  • CSCD
  • 英国科学文摘

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

具有不等式路径约束的微分代数方程系统的动态优化

孙燕 张弛 路兴龙 王靖戈 付俊

孙燕, 张弛, 路兴龙, 王靖戈, 付俊. 具有不等式路径约束的微分代数方程系统的动态优化. 自动化学报, 2019, 45(5): 897-905. doi: 10.16383/j.aas.c180302
引用本文: 孙燕, 张弛, 路兴龙, 王靖戈, 付俊. 具有不等式路径约束的微分代数方程系统的动态优化. 自动化学报, 2019, 45(5): 897-905. doi: 10.16383/j.aas.c180302
SUN Yan, ZHANG Chi, LU Xing-Long, WANG Jing-Ge, FU Jun. Dynamic Optimization of Differential-algebraic Equations With Inequality Path Constraints. ACTA AUTOMATICA SINICA, 2019, 45(5): 897-905. doi: 10.16383/j.aas.c180302
Citation: SUN Yan, ZHANG Chi, LU Xing-Long, WANG Jing-Ge, FU Jun. Dynamic Optimization of Differential-algebraic Equations With Inequality Path Constraints. ACTA AUTOMATICA SINICA, 2019, 45(5): 897-905. doi: 10.16383/j.aas.c180302

具有不等式路径约束的微分代数方程系统的动态优化

doi: 10.16383/j.aas.c180302
基金项目: 

国家自然科学基金 61590922

中央高校基本科研专项资金 N150802001

国家自然科学基金 61473063

国家自然科学基金 61503041

详细信息
    作者简介:

    孙燕   东北大学流程工业综合自动化国家重点实验室硕士研究生.主要研究方向为动态优化, 不等式路径约束.E-mail:suny618@foxmail.com

    张弛   东北大学流程工业综合自动化国家重点实验室博士研究生.主要研究方向为切换系统的动态优化, 全局优化.E-mail:deville136@hotmail.com

    路兴龙   东北大学流程工业综合自动化国家重点实验室博士研究生.主要研究方向为自适应动态规划, 非线性控制理论, 故障检测与诊断.E-mail:xinglonglu@stumail.neu.edu.cn

    王靖戈   东北大学硕士研究生.主要研究方向为多智能体系统一致性问题, 切变网络, 事件驱动控制策略, 非凸受限.E-mail:wangjingge@hotmail.com

    通讯作者:

    付俊   东北大学流程工业综合自动化国家重点实验室教授.主要研究方向为动态优化, 切换系统, 非线性控制.本文通信作者.E-mail:junfu@mail.neu.edu.cn

Dynamic Optimization of Differential-algebraic Equations With Inequality Path Constraints

Funds: 

Natural Science Foundations of China 61590922

the Fundamental Research Funds for the Central Universities N150802001

Natural Science Foundations of China 61473063

Natural Science Foundations of China 61503041

More Information
    Author Bio:

    Master student at the State Key Laboratory of Synthetical Automation for Process Industries, Northeastern University. Her research interest covers dynamic optimization and inequality path constraint

    Ph. D. candidate at the State Key Laboratory of Synthetical Automation for Process Industries, Northeastern University. His research interest covers dynamic optimization of switched systems and global optimization

    Ph. D. candidate at the State Key Laboratory of Synthetical Automation for Process Industries, Northeastern University. His research interest covers adaptive dynamic programming, nonlinear control theory, and fault detection and diagnosis

    Master student at Northeastern University. His research interest covers multi-agent systems, event-triggered mechanism, and nonconvex constraints

    Corresponding author: FU Jun Professor at the State Key Laboratory of Synthetical Automation for Process Industries, Northeastern University. His research interest covers dynamic optimization, switching system, and nonlinear control. Corresponding author of this paper
  • 摘要: 针对具有不等式路径约束的微分代数方程(Differential-algebraic equations,DAE)系统的动态优化问题,通常将DAE中的等式路径约束进行微分处理,或者将其转化为点约束或不等式约束进行求解.前者需要考虑初值条件的相容性或增加约束,在变量间耦合度较高的情况下这种转化求解方法是不可行的;后者将等式约束转化为其他类型的约束会增加约束条件,增加了求解难度.为了克服该缺点,本文提出了结合后向差分法对DAE直接处理来求解上述动态优化问题的方法.首先利用控制向量参数化方法将无限维的最优控制问题转化为有限维的最优控制问题,再利用分点离散法用有限个内点约束去代替原不等式路径约束,最后用序列二次规划(Sequential quadratic programming,SQP)法使得在有限步数的迭代下,得到满足用户指定的路径约束违反容忍度下的KKT(Karush Kuhn Tucker)最优点.理论上证明了该算法在有限步内收敛.最后将所提出的方法应用在具有不等式路径约束的微分代数方程系统中进行仿真,结果验证了该方法的有效性.
  • 动态优化是约束中含有微分或差分方程的一类数学规划问题[1-3], 其广泛存在于电力系统[4-5]、石油化工[6-7]、生物工程[8-9]、清洁能源[10-11]等.目前在化工过程中基于常微分方程模型的动态优化问题不仅涉及微分方程, 还包括代数方程约束, 这样的系统统称为微分代数系统.此外, 实际的流程工业过程中, 需要某些状态变量或控制变量的函数在其全部运行过程中或部分运行时间内不能超过约束的设定值, 来保证工业过程的质量和安全, 这就需要考虑具有不等式路径约束的动态优化问题.

    动态优化问题(Dynamic optimization problem, DOP)的数值求解算法通常分为直接法[12-13]和间接法[14-15].间接法通过求解原问题的最优性条件(即必要条件), 间接地获得原问题的最优解.但是对于相对复杂的问题, 比如包含不等式约束的最优控制问题, 由于状态变量或者控制变量在约束中存在, 求解极为困难[16].此外, 复杂的非线性系统最优性条件很难被定义, 且两点边值问题的求解过程中收敛域可能很小[17].直接法是通过对控制变量离散化, 即控制向量参数化[18] (Control vector parameterization, CVP), 或控制与状态变量同时离散化, 即正交配置(Orthogonal collocation, OC), 将无限维的最优控制问题转化为有限维的非线性最优化问题.由于OC方法的计算复杂度较高, 本文采用CVP方法.

    针对具有不等式路径约束的微分代数方程(Differential-algebraic equations, DAE)系统的动态优化问题, 采用直接或间接法对其进行转化后, 仍需要对等式路径约束和不等式路径约束进行处理.胡云卿等将等式路径约束进行微分[19], 将原问题转化为具有不等式的常微分方程(Ordinary differential equations, ODE)动态优化问题, 再使用基于ODE的动态优化方法进行处理, 但是DAE到ODE的转化过程中需要考虑初值条件的相容性, 并且无法在理论上证明有关解的存在性和唯一性[20], 此外, 当约束中变量的幂较高或者变量间耦合程度较高的情况下是不可行的. Pontryagin等将等式路径约束转化为两点边值问题[21]进行求解, 该方法需要增加约束条件, 增加了求解难度. Fabien将等式约束转化为终点约束[22], 这样会导致问题的维数变得很大, 需要对转化后的问题进行降维处理.所以间接地处理DAE中的等式路径约束会引起一系列的问题, 因此本文采用后向差分法(Backward differentiation formula, BDF)[23]直接求解DAE.

    对于不等式路径约束, 文献[24]采用多重打靶(Multiple shooting, MS)法, 将不等式路径约束转化为分段的点约束, 每次打靶均需要给出对应打靶分段内该点约束涉及的控制变量和各状态变量合理的范围, 各变量段间的连续还需要一个合适的系数保证, 不然很难得到理想的结果.文献[25]对不等式路径约束的积极与否进行判断, 若不等式路径约束为积极约束, 则结合微分方程在DAE求解器中进行求解, 否则不考虑该不积极约束; 不等式路径约束为积极约束时, DAE阶次可能很高, 所以该方法对DAE求解器要求较高. Jacobson等提出的松弛变量法[26], 通过增加松弛变量将不等式路径约束转化为相应的等式路径约束, 但是这种方法难以处理不等式路径约束数比控制变量数多的问题; Vassiliadis等[27]通过离散化将约束项在问题时间约束段内选择有限个内点进行离散化, 通过惩罚函数法将约束违反的积分设为0来保证约束的无违反, 但是此方法可能会在NLP问题中引进较多的约束条件; Floudas提出的凸函数近似法[28]主要应用于静态系统, Chachuat等[29]将其改进应用在动态系统的最优控制上, 此方法求得的数值解虽然可以严格满足不等式路径约束, 但算法复杂度较高; Rehbock等提出的精确罚函数法[30]能够成功求解具有不等式路径约束的复杂最优控制问题[31], 该方法虽然在算法复杂度上有所改善但是不能严格满足不等式约束, 且约束违反程度不能确定; Liu等[32-33]提出了一种新颖的光滑化精确惩罚函数法, 将不等式路径约束转化为一个光滑化惩罚项加到目标函数中, 将问题转变成一个无约束优化问题; 由于目标函数表达式会因不等式约束的违反程度随时间而改变, 所以如何获取问题转化后目标函数的梯度信息是一个难点. 2005年Chen在之前的CVP-OC混合方法基础上进一步提出了有限收敛法[34], 问题初始求解时不考虑不等式路径约束的最优控制问题, 将初始所得解做为初始解以此优化原问题, 这样可以使得优化过程得以优化.然而该方法受限于带有不等式路径约束的ODE动态优化问题, 对于带有不等式路径约束的非0阶DAE动态系统的适用性有待探究.

    因此, 本文针对具有不等式路径约束的DAE系统的动态优化问题, 采用只对控制变量进行参数化的CVP方法, 将无限维的DOP转化为有限维的DOP, 仍然保留系统的动态特性.通过分点离散化方法将不等式路径约束转化为有限的内点约束进行处理, 设计了在指定的路径约束违反容忍度下通过有限步迭代获得KKT (Karush Kuhn Tucker)最优点的算法.然后利用BDF方法求解DAE方程组, 再采用序列二次规划法(Sequential quadratic programming, SQP)获得满足不等式路径约束的最优解.最后将该算法应用在工业化工问题中, 仿真结果表明了该算法可以获得最优控制轨迹并有着良好的路径约束效果及收敛性.

    本文的主要贡献如下:

    1) 采用后向差分法直接求解DAE, 避免将DAE化为ODE的过程中会产生的解的相容性、维数增加等前文所述问题. 2)所采用的处理路径约束的方法不会改变目标函数的结构形式.并且可以在指定路径约束违反的容忍度的条件下, 经过有限步迭代得到最优解.

    本文结构如下:第1节对具有不等式路径约束的DAE系统的最优控制问题进行描述.第2节阐述本文所采用的求解动态优化问题和处理不等式路径约束的算法.第3节将算法应用于工业化工问题.最后对本文所做工作进行了总结并提出了对未来可以改进的方面进行了展望.

    本文考虑的具有不等式路径约束的DAE系统的最优控制问题, 根据Biegler[35]给出的动态优化问题描述形式添加不等式路径约束条件, 记为P1, 给出数学描述如下:

    $ \begin{align} &\min\limits_{\boldsymbol{x}(t), \boldsymbol{v}(t)}J =\phi[\boldsymbol{x}(t_f), \boldsymbol{ y}(t_f)]+ \nonumber\\&\qquad\int_{t_0}^{t_f}\ell[\boldsymbol{x}(t), \boldsymbol{y}(t), \boldsymbol{v}(t)]{\rm d}t\end{align} $

    (1)

    s. t.

    $ \dot{\boldsymbol{x}}(t) = f[t, \boldsymbol{x}(t), \boldsymbol{y}(t), \boldsymbol{v}(t)] $

    (2)

    $ g[t, \boldsymbol{x}(t), \boldsymbol{y}(t), \boldsymbol{v}(t)]=0 $

    (3)

    $ h[t, \boldsymbol{x}(t), \boldsymbol{y}(t), \boldsymbol{v}(t)] \leq 0 $

    (4)

    $ \boldsymbol{x}(t_0) = \boldsymbol{x}_0 $

    (5)

    $ \boldsymbol{u} _l \leq \boldsymbol{v}(t) \leq \boldsymbol{u}_u $

    (6)

    $ t_0 \leq t \leq t_f $

    (7)

    其中, $\boldsymbol{x}(t)\in{\bf R}^{n_x}$是状态变量, $\boldsymbol{y}(t)\in{\bf R}^{n_y}$是代数变量, $\boldsymbol{v}(t)\in{\bf R}^{n_u}$是控制变量且上下边界约束分别为$\boldsymbol{u} _l$和$\boldsymbol{u}_u$, 这些变量都是时间标量$t\in[t_0, t_f]$的函数; 动态过程开始时刻$t_0$和结束时刻$t_f$都是固定的; 动态系统特性由常微分方程(2)组成; 其初始可行状态为$\boldsymbol{x}_0$; 动态优化问题的最优解需要满足等式路径约束(3)、不等式路径约束(4);目标函数式(1)由终值项$\phi[\boldsymbol{x}(t_f), \boldsymbol{y}(t_f)]$和被积函数为$\ell[t, \boldsymbol{x}(t), \boldsymbol{y}(t), \boldsymbol{v}(t)]$的积分项组成.其中, 等式路径约束(3)和常微分方程(2)组成DAE.为不失一般性, 假设这里的DAE对各变量均有连续二阶偏导数, 且index为1.

    这些DAE系统具有更大的灵活性:可以是线性的也可以是非线性的; 不等式路径约束(4)可以有多个, 若有多个不等式路径约束, 其处理方法同(4).

    该最优控制问题可以简单描述为:在可行初始状态(5)条件下, 寻找最优的控制轨迹$\boldsymbol{v}^{*}(t)$, 使得式(2)和式(3)描述的系统满足约束(4), 且使得目标函数(1)的值最小.

    本文针对具有不等式路径约束的DAE系统的最优控制问题, 首先利用CVP技术将无限维的DOP转化为有限维的DOP, 将原问题转化为仍保留系统动态特性的具有不等式路径约束的DOP.对于系统的动态方程, 采用BDF直接求解, 避免了将DAE转化为ODE或两点边值问题.利用逐点离散法处理不等式路径约束, 将不等式约束在违反处转化为有限的点约束, 在有限步的迭代下, 得到满足用户指定的路径约束违反容忍度的KKT最优点.算法的主要结构如图 1所示.

    图 1  算法主要结构图
    Fig. 1  Structure diagram of the algorithm

    本节首先简要介绍了控制向量参数化方法及其处理后得到的有限维的DOP.然后在逐点离散法中, 路径约束被转化为点约束的形式, 同时给出了不等式路径约束的DAE系统最优控制问题的具体求解算法.并给出了分点离散法收敛性的证明.

    CVP是求解动态优化问题的一种重要策略. CVP方法的核心思想[36]是:利用基函数将控制变量进行参数化, 将原无限维最优控制问题转化为一个含有动态过程的有限维优化问题后, 再用传统的DOP求解方法进行求解.即将控制向量$v(t)$被离散为$n$段而状态变量保持连续, 每个时间分段内用带参数的多项式(如常数、线性或二次多项式).关于控制变量参数化方法的收敛性, 文献[37-38]都进行了相应的证明.

    首先不考虑不等式路径约束, 将动态优化问题的时间域$[t_0, t_f]$划分为$N$个控制阶段,

    $ \begin{equation} \begin{aligned} t_0\le t_1\le \cdots \le t_N \equiv t_f \end{aligned} \end{equation} $

    (8)

    在某一时间分段$[t_{k-1}, t_k], k=1, 2, \cdots, N$内, 根据函数逼近理论, 控制变量可以表示为低阶多项式参数依赖的表达形式,

    $ \begin{equation} \begin{aligned} v_{jk}(t)\approx \!\bar{v}_{jk}(\!t)=\sum\limits_{i=0}^{M}u_{ijk}\zeta_{ijk}^{\!(\!M\!)}(\!t), t_{k-1}\!\leq t \!\leq t_k \end{aligned} \end{equation} $

    (9)

    其中, $j$表示原问题有$1, 2, \cdots, j, \cdots, n_u$个控制变量, $\zeta_{ijk}^{(M)}(t)$为多项式基函数, $M$为基函数阶数.实际计算和应用中控制向量的参数化可以有多种参数依赖形式, 如图 2所示常用的有分段常量逼近策略(如图 2实线)、分段线性逼近策略(如图 2虚线)、分段二次多项式逼近策略(如图 2点虚线)等.由于计算机采用数字信号, 非常符合分段常量策略形式, 同时这种策略简便易行、生成参数数目较少.因此本文采用分段常量逼近策略进行控制变量参数化, 即采用$\zeta_{ij}^k(t)=1$为基函数.

    图 2  控制变量分段示例
    Fig. 2  Illustration of control variable profiles

    从而整个时间段内的控制变量可以表示为:

    $ \begin{equation} \begin{aligned} \boldsymbol{v}(t)\approx \sum\limits_{k=1}^{N}\boldsymbol{u}_{k}X_{[t_{k-1}, t_k)}(t), \ \ \ \ t\in[t_0, t_f] \end{aligned} \end{equation} $

    (10)

    其中, $X_{[t_{k-1}, t_k)}(t)$为开关函数, $k=1, 2, \cdots, N$, 具体形式如下:

    $ \begin{equation} \begin{aligned} X_{[t_{k-1}, t_k)}(t)=\left\{ \begin{aligned} 1, \ \ \ \text{若} \ \ t\in[t_{k-1}, t_k) \\ 0, \ \ \ \text{若} \ \ t\not\in[t_{k-1}, t_k) \\ \end{aligned} \right. \end{aligned} \end{equation} $

    (11)

    进而每一分段内, 控制变量可以表示为:

    $ \begin{equation} \begin{aligned} v_j(t)\approx u_{jk}\ \ \ \ t\in[t_{k-1}, t_k] \end{aligned} \end{equation} $

    (12)

    利用上述中CVP思想对原问题(1)~(7)进行处理, 得到如下有限维动态优化问题, 记为P2:

    $ \min\limits_{\boldsymbol{x}(t), \boldsymbol{U}}J =\phi[\boldsymbol{x}(t_f), \boldsymbol{y}(t_f)]+\int_{t_0}^{t_f}\ell[\boldsymbol{x}(t), \boldsymbol{y}(t), \boldsymbol{u}]\rm dt $

    (13)

    s.t.

    $ \dot{\boldsymbol{x}}(t) = f[t, \boldsymbol{x}(t), \boldsymbol{y}(t), \boldsymbol{u}] $

    (14)

    $ g[t, \boldsymbol{x}(t), \boldsymbol{y}(t), \boldsymbol{u}]=0 $

    (15)

    $ h[t, \boldsymbol{x}(t), \boldsymbol{y}(t), \boldsymbol{u}] \leq 0 $

    (16)

    $ \boldsymbol{x}(t_0) = \boldsymbol{x}_0 $

    (17)

    $ \boldsymbol{u}_l \leq \boldsymbol{u} \leq \boldsymbol{u}_u $

    (18)

    $ t_0 \leq t \leq t_f $

    (19)

    其中, $\boldsymbol{u}\in \boldsymbol{U}$是$\boldsymbol{ v}(t)$离散化后的有限维控制向量, 并且$\boldsymbol{U}\in {\bf R}^N$是个非空$N$维向量.其他变量的定义同问题P1.

    对于参数化后的DOP问题可使用成熟的基于梯度的优化算法进行求解.本文首先采用BDF直接对微分代数方程组进行求解, 进而得到带有有限维控制向量的目标函数, 再用BFGS迭代法计算拉格朗日函数的Hessian矩阵的正定拟牛顿近似值, 利用SQP算法便可得到最优的控制参数组合.

    2.2.1   算法步骤

    本文利用文献[34]中的分点离散法来处理不等式路径约束, 记分段点集合$P_N=\{ t_0, t_1, \cdots, $ $t_{N-1}, t_N\equiv t_f\}, t_{i-1}<t_i, i=1, 2, \cdots, N$.对于DOP问题中的式(13), 若某时间段上存在不等式路径约束的违反, 则在该段内将路径约束离散为有限多个点约束:

    $ \begin{equation} \begin{aligned} h[t_i, \boldsymbol{x}(t_i), \boldsymbol{y}(t_i), \boldsymbol{u}_i]\leq 0 \end{aligned} \end{equation} $

    (20)

    其中, $t_j, j=1, 2, \cdots, M$是时间域内$M$个违反时间段上的离散形式.给出求解具有不等式路径约束的DAE系统最优控制问题求解具体算法如下:

    初始化. 设置初始控制参数$\boldsymbol{ u}^{(0)}$、初始化不等式路径约束违反点集合$T_{*}^{(0)}\subset T$; 有限且任意小的容忍度参数$\eta >0$; 分段参数$N\in N^{*}$; 迭代次数$p=0$; 迭代精度$\delta>0$;

    步骤1. 先解决不带有不等式路径约束的DOP问题P2, i. e.解式(13)$\, \sim\, $(15)、式(17)$\, \sim\, $(19)构成的动态优化问题;

    步骤2. 利用BDF法求解式(14), (15)和式(17)构成的DAE初值问题, 得到$\boldsymbol{x}(t)$的值;

    步骤3. 在容忍度$\eta$下检查$h(t)\leq \eta$, $t\in[t_0, t_f]$的违反情况.将每个CVP分段中违反量最大的点对应的点约束加到DOP的约束中, 并且检查已有违反点的违反情况, 若不存在违反则去除该点, 更新$T_{*}^{(p)}$;

    步骤4. 使用BFGS迭代法计算拉格朗日函数的Hessian矩阵的正定拟牛顿近似值获取关于参数向量$\boldsymbol{u}$的梯度信息;

    步骤5. 通过SQP方法, 求解QP子问题, 生成新的迭代步$\hat{\boldsymbol{u}}_{k+1}^{(p)}$;

    步骤6. 若在容忍度$\eta$下$h(t_i)\leq \eta$, 其中$t_i\in T_{*}^{(p)}$, 且$\| J(\hat{\boldsymbol{ u}}_{k+1}^{(p)})-J(\hat{\boldsymbol{u}}_{k}^{(p)})\|<\delta$, 则算法结束; 否则重新求解带现有点约束的NLP问题, $p=p+1$, 继续步骤2.

    若问题P2中的不等式路径约束存在违反, 则会循环迭代求解带不同点约束的DOP问题, 通过SQP方法调整控制变量$\boldsymbol{u}$, 使得点约束的违反程度有效降低, 直至收敛至容忍度$\eta$内.不等式路径约束到点约束的转换是收敛的(详细内容参考第2.2.2节), 因此通过有限步迭代, 能找到满足指定容忍度的KKT最优点.

    2.2.2   收敛性证明

    在原DOP问题经过CVP处理后的$N$个控制段上, 对于分点离散法的收敛性分析, 本文给出以下证明:

    对于DAE问题P1的不等式路径约束均有以下假设:

    假设1. $\ h(t)$至少是一阶可导的.

    假设2. $\ h(t)$对于$t$是Lipschitz连续的, 即存在一个常数$K$, 使得对所有$(t)$和$(\hat{t})$, 在域$D=\{t_0\leq t\leq t_f, \left| \boldsymbol{x} \right| \le \infty, \left| \boldsymbol{u} \right| \le \infty \}$上都有下式成立:

    $ \begin{eqnarray} \left| h(t)-h(\hat{t}) \right|\leq K \left| t-\hat{t} \right| \end{eqnarray} $

    (21)

    假设3. 若$h(t)$对于$t$不是Lipschitz连续的, 记不连续点集合为$P_{non}$且$P_{non}\in P_M$, 则$h(t)$满足假设2, 其中$t\notin P_{non}$.

    定理1. 若问题P1的任意一个不等式路径约束$h(t)\leq 0$都满足上述假设, 则对于该$h(t)$存在有限的点集$P_M=\{\hat{t}_0, \hat{t}_1, \cdots, \hat{t}_{N(\eta)}\equiv t_f\}$使得$h(t)\leq \eta: \eta > 0, t\in[t_0, t_f]$.

    证明.   原DOP问题经CVP处理后, 时间域$T\in[t_0$, $t_f]$, $i=0, 1, \cdots, N-1$会被划分为$N$段相等或者不等的子时间段$T_i=[t_i, t_{i+1}], i=0, 1, \cdots, N-1$.在每个子时间段$T_i$上求解DAE初值问题时, 该子时间段将被离散为$M$个小时间分段$\hat{T_k}=[\hat{t}_k, \hat{t}_{k+1}], k=0, 1, \cdots, M-1$, 其中$\hat{t}_0=t_i$, $\hat{t}_M=t_{i+1}$.对于任意一个小时间分段$[\hat{t}_k, \hat{t}_{k+1})$, 取$\hat{t}<\hat{t}_{k+1}, \hat{t}\in \mathring{U}(\hat{t}_{k+1})$.由$h(t)$在$[\hat{t}_{k}, \hat{t}]$上连续, 在$(\hat{t}_{k}, \hat{t})$内可导, 那么在$(\hat{t}_{k}, \hat{t})$内至少存在一点$\xi(\hat{t}_{k}< \xi<\hat{t})$, 使等式$h(\hat{t})-h(\hat{t}_{k})=h' (\xi)(\hat{t}-\hat{t}_{k})$成立, 当$\left| \hat{t}_{k+1}-\hat{t}_{k} \right|<\varepsilon $时, 对等式进一步处理有:

    $ \begin{align} &\left| h(\hat{t})-h(\hat{t}_{k}) \right|=\left| h'(\xi)\right|\cdot\left| \hat{t}-\hat{t}_{k} \right| \le& \nonumber\\ &\qquad K\left| \hat{t}-\hat{t}_{k} \right| \le K\left| \hat{t}_{k+1}-\hat{t}_{k} \right|\le \eta, \quad \hat{t}\in \hat{T_i} \end{align} $

    (22)

    去绝对值进行移项, 即

    $ \begin{eqnarray} \begin{aligned} h(\hat{t}_{k})-\eta\le h(\hat{t})\le h(\hat{t}_{k})+\eta, \ &\hat{t}\in \hat{T_i} \end{aligned} \end{eqnarray} $

    (23)

    又$h(\hat{t}_{k})\le 0$, 则由上式可得到:

    $ \begin{equation} \begin{aligned} h(\hat{t})\le h(\hat{t}_{k})+\eta\le \eta, \ \ &\hat{t}\in \hat{T_i} \end{aligned} \end{equation} $

    (24)

    综上, 在整个时间域$[t_0, t_f]$上的每个子时间段$\hat{T_i}, i=0, 1, \cdots, N(\eta)-1$上都有上式成立, 因此有

    $ \begin{equation} \begin{aligned} h(t)\le \eta, \ \ &t\in[t_0, t_f] \end{aligned} \end{equation} $

    (25)

    其中, $N(\eta)$是$\eta$的函数, $\eta$是式(22)满足的条件.

    对于有限的$\eta$, $N(\eta)$也是有限的.

    对于问题P2的每一个不等式路径约束均可此类证明.

    本节采用两个DAE系统的最优控制问题对提出的算法进行仿真研究.仿真环境为: Intel i7-4790, 3.6 GHz CPU和DDMM1/1 600 MHZ 8 GB内存.在仿真中求解DAE方程组求解精度设为$10^{-6}$; 使用SQP方法解NLP问题, 求解精度设为$10^{-4}$; 内点约束的收敛容忍度为$\eta=10^{-6}$.

    由Guun等提出的催化剂混合问题[39], 一直以来被研究者研究用于动态优化问题的研究[40-41].问题可描述为:在长度为$l_f$的管式反应器中, 两种物质A、B在管内发生反应.反应式为$A\leftrightarrow B\rightarrow C$.其中B是副产物, C是目标产物, 反应速率随催化剂沿着管长的混合率而变化.为了得到更多的目标产物C, 该动态优化控制问题的求解目标是获得催化剂最佳配比方案.其数学模型如下式所示:

    $ \begin{align} &\max\limits_{u(l)} \ J = x_3(l_f)\nonumber\\ \rm s.\, t.\ & \notag\\ &\dot{x}_1(l) = u(l)[10x_2(l)-x_1(l)] \nonumber\\ \ & \dot{x}_2(l) = u(l)[x_1(l)-10x_2(l)]-[1-u(l)]x_2(l) \nonumber\\ \ & 1-x_1(l)- x_2(l)-x_3(l)=0\nonumber\\ \ & x_2(l)\leq 0.075 \end{align} $

    (26)

    其中, 状态变量$x_1(l)$、$x_2(l)$、$x_3(l)$分别表示为A、B、C的摩尔分率, 初始状态为$x_1(0)=1$, $x_2(0)=0$, $x_3(0)=0$; 控制变量$u(l)$表示催化剂沿管长的混合率, 满足的约束为$0\leq u(l) \leq 1$; $l$在这里表示管长.这里取$l_0=0, l_f=12$, 仿真中取$u_0=0.17$.

    根据本文给出的具有不等式路径约束的DAE系统最优控制问题的求解算法, 若不等式在$[l_0, l_f]$上存在违反, 则将其转化为内点约束$x_2(l_i)-0.075\leq\eta$.图 3图 4是控制变量取为$N=24$时得到的控制轨迹和状态轨迹. 图 5为$x_2(l)$的状态轨迹.可以看出在整个长度域$l\in[l_0, l_f]$内都满足不等式路径约束, 说明本文提出的方法在路径约束处理方面有着很好的效果.

    图 3  控制曲线
    Fig. 3  Control profile
    图 4  路径约束下的状态曲线
    Fig. 4  The status curve with path constraint
    图 5  路径约束下$x_2$的状态曲线
    Fig. 5  The status curve of $x_2$ with path constraint

    仿真中分别取$N=24$、29、34, 得到问题最优目标函数和求解时间如表 1.由表 1中可以看出, 本文方法求解的时间与文献[18]相比, 省时至少20 %, 精确度方面会因为参数有所浮动, 但是目标函数最差在$N=29$时, 本文所得目标函数也已经达到文献[18]的99.97 %.

    表 1  催化剂混合问题测试结果
    Table 1  Results of catalyst mixing problem
    测试对象分段数$N$目标函数$J$计算时间(s)
    240.4754609111.2188
    本文方法290.4759645190.1719
    340.4763620348.0000
    240.4755665138.9219
    文献[18]方法290.4761231290.1563
    340.4765138460.9375
    下载: 导出CSV 
    | 显示表格

    在文献[42]中描述的一个青霉素分批补料发酵问题(Fed-batch penicillin fermentation), 该问题具有很强的非线性和奇异性.该动态优化问题描述如下:

    $ \begin{eqnarray} \begin{aligned} \min\limits_{S(t)} \ J = &\int_{t_0}^{t_f}(-\theta X(t)V(t)+0.0103P(t)V(t)+ \\ &\nonumber 0.0744\mu(t) X(t)+0.00102X(t)V(t) +\\&\nonumber6\, 913.58){\rm d}t\\ \rm s.\, t.\ & \dot{X}(t)= \mu(t) X(t)-\frac {FX(t)}{V(t)}\\ \ & \dot{P}(t)= \theta X(t)-0.01P(t)-\frac {FP(t)}{V(t)}\\ \ & \dot{V}(t)= F \\ \ & \mu(t)=\frac {0.11S}{S+0.006X(t)} \\ \ & \theta=\frac{0.004}{1+\frac {0.0001}{S}+\frac {S}{0.1}}\\ \ & X(t)\leq 41 \end{aligned} \end{eqnarray} $

    其中, 目标函数$J$是总收益. $X$是生物量浓度$(\rm g/L)$, $P$是现有青霉素产物数量(每升活性), 具体增长率$\mu$, 具体产品形成速率$\theta$和反应堆的体积$V(\rm L)$都是系统的状态变量.底物浓度$S(\rm g/L)$是控制变量.进料流量$F$为$1\, 666.67\, \rm L/h$.该模型初始状态为$[X(0), P(0), V(0), \mu(0)]=[1, 0, 250\, 000, 0]$.起始时刻$[t_0, t_f]=[0, 150]$.控制变量$S\in[0.001, 0.5]$.

    采用文献[18]中处理不等式路径约束的方法与本文方法对比, 仿真中取$N=12$、14、16分别得到最优目标函数和求解时间如表 2所示.由表 2中可以看出, 本文方法求解的时间约为文献[18]的十分之一, 精确度方面会因为参数有所浮动, 但是目标函数最差在$N=12$时, 本文所得目标函数也已经达到文献[18]的99.3 %.

    表 2  青霉素分批补料发酵问题测试结果
    Table 2  Results of fed-batch penicillin fermentation
    测试对象分段数$N$目标函数$J$计算时间(s)
    本文方法121 108 269187.5313
    141 121 255301.5781
    161 115 655315.1250
    文献[18]方法121 113 7001 496.0
    141 102 6002 578.9
    161 102 4004 109.8
    下载: 导出CSV 
    | 显示表格

    仿真时设置$u_0=0.5$.取控制变量$N=12$时, 图 6为得到的控制轨迹, 图 7是状态轨迹, 其中状态变量$X(t)\leq41$.从图中可以看出本文所提方法的有效性.

    图 6  控制曲线
    Fig. 6  Control profile
    图 7  路径约束下的状态曲线
    Fig. 7  The status curve with path constraint

    综上两个仿真实验, 本文方法允许路径约束违反在用户指定范围之内, 在问题允许存在误差或者实际上确实存在误差的情况下, 以一定的准确度换取求解速度, 满足应用需求.

    本文针对带有不等式路径约束的DAE系统的动态优化问题, 提出了一种直接求解该类型问题的框架, 首先, 通过CVP将无限维的动态优化问题转化为有限维的动态优化问题; 然后, 利用分点离散方法, 设计了一个可以在一定容忍度下满足路径约束的算法, 并在理论上论证了该算法可以在有限步数的迭代下收敛; 最后, 利用催化剂混合问题和压力限定批量反应堆问题仿真研究验证了本文所提方法的有效性.

    本文对控制变量进行的是等分离散化, 但对局部变化较大的控制信号, 等分离散化不能准确地反应真实的控制信号.因此本文下一步的工作将研究根据控制信号曲线自适应地进行CVP离散化.


  • 本文责任编委 何海波
  • 图  1  算法主要结构图

    Fig.  1  Structure diagram of the algorithm

    图  2  控制变量分段示例

    Fig.  2  Illustration of control variable profiles

    图  3  控制曲线

    Fig.  3  Control profile

    图  4  路径约束下的状态曲线

    Fig.  4  The status curve with path constraint

    图  5  路径约束下$x_2$的状态曲线

    Fig.  5  The status curve of $x_2$ with path constraint

    图  6  控制曲线

    Fig.  6  Control profile

    图  7  路径约束下的状态曲线

    Fig.  7  The status curve with path constraint

    表  1  催化剂混合问题测试结果

    Table  1  Results of catalyst mixing problem

    测试对象分段数$N$目标函数$J$计算时间(s)
    240.4754609111.2188
    本文方法290.4759645190.1719
    340.4763620348.0000
    240.4755665138.9219
    文献[18]方法290.4761231290.1563
    340.4765138460.9375
    下载: 导出CSV

    表  2  青霉素分批补料发酵问题测试结果

    Table  2  Results of fed-batch penicillin fermentation

    测试对象分段数$N$目标函数$J$计算时间(s)
    本文方法121 108 269187.5313
    141 121 255301.5781
    161 115 655315.1250
    文献[18]方法121 113 7001 496.0
    141 102 6002 578.9
    161 102 4004 109.8
    下载: 导出CSV
  • [1] 陈美蓉, 郭一楠, 巩敦卫, 杨振.一类新型动态多目标鲁棒进化优化方法.自动化学报, 2017, 43 (11): 2014-2032 http://www.aas.net.cn/CN/abstract/abstract19176.shtml

    Chen Mei-Rong, Guo Yi-Nan, Gong Dun-Wei, Yang Zhen. A novel dynamic multi-objective robust evolutionary optimization method. Acta Automatica Sinica, 2017, 43 (11): 2014-2032 http://www.aas.net.cn/CN/abstract/abstract19176.shtml
    [2] 陈龙, 刘全利, 王霖青, 赵珺, 王伟.基于数据的流程工业生产过程指标预测方法综述.自动化学报, 2017, 43(6): 944-954 http://www.aas.net.cn/CN/abstract/abstract19072.shtml

    Chen Long, Liu Quan-Li, Wang Lin-Qing, Zhao Jun, Wang Wei. Data-driven prediction on performance indicators in process industry: a survey. Acta Automatica Sinica, 2017, 43 (6): 944-954 http://www.aas.net.cn/CN/abstract/abstract19072.shtml
    [3] 丁进良, 杨翠娥, 陈立鹏, 柴天佑.基于参考点预测的动态多目标优化算法. 自动化学报, 2017, 43 (2): 313-320 http://www.aas.net.cn/CN/abstract/abstract19009.shtml

    Ding Jin-Liang, Yang Cui-E, Chen Li-Peng, Chai Tian-You. Dynamic multi-objective optimization algorithm based on reference point prediction. Acta Automatica Sinica, 2017, 43 (2): 313-320 http://www.aas.net.cn/CN/abstract/abstract19009.shtml
    [4] Bienstock D. Optimal control of cascading power grid failures. In: Proceedings of the 2010 IEEE Conference on Decision and Control and European Control Conference, New York, USA: IEEE, 2010. 2166-2173 https://www.researchgate.net/publication/48170873_Optimal_adaptive_control_of_cascading_power_grid_failures
    [5] Chomat M, Schreier L, Bendl J. Optimal control of input rectifier in voltage-Source inverter supplied from unbalanced power grid. In: Proceedings of the 2006 IEEE International Symposium on Industrial Electronics, New York, USA: IEEE, 2006. 1042-1045 https://www.researchgate.net/publication/224687371_Optimal_Control_of_Input_Rectifier_in_Voltage-Source_Inverter_Supplied_from_Unbalanced_Power_Grid
    [6] Liu P, Li G, Liu X. Fast engineering optimization: a novel highly effective control parameterization approach for industrial dynamic processes. ISA Transactions, 2015, 58 : 248-254 doi: 10.1016/j.isatra.2015.06.006
    [7] Jie X, Huang Y, Lou H H. A probability distribution estimation based method for dynamic optimization. AIChE Journal, 2010, 53 (7): 1805-1816 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=a19510f734e0e72ed71850edef587b7c
    [8] Hirmajer T, Balsa-Canto E, Banga J R. DOTcvpSB, a software toolbox for dynamic optimization in systems biology. BMC Bioinformatics, 2009, 10 (1): 199 doi: 10.1186/1471-2105-10-199
    [9] Bredies K, Lorenz D A, Maass P. An optimal control problem in medical image processing. In: Proceedings of the IFIP-TC7 Conference, Turin, Italy: DBLP, 2005: 249-259 doi: 10.1007%2F0-387-33882-9_23
    [10] Barve H A, Banavar R N. Energy-optimal control of a particle in a dielectrophoretic system. In: Proceedings of the 2010 IEEE Conference on Decision & Control. New York, USA: IEEE, 2010: 3353-3358 https://www.researchgate.net/publication/224108918_Energy-optimal_control_of_a_particle_in_a_dielectrophoretic_system?ev=auth_pub
    [11] Barve H A, Banavar R N. Energy-optimal control of a particle in a dielectrophoretic system. In: Proceedings of the 2010 IEEE Conference on Decision & Control. New York, USA: IEEE, 2010: 3353-3358 https://www.researchgate.net/publication/224108918_Energy-optimal_control_of_a_particle_in_a_dielectrophoretic_system?ev=auth_pub
    [12] 孙勇, 张卯瑞, 梁晓玲.求解含复杂约束非线性最优控制问题的改进Gauss伪谱法. 自动化学报, 2013, 39 (5): 672-678 http://www.aas.net.cn/CN/abstract/abstract17921.shtml

    Sun Yong, Zhang Mao-Rui, Liang Xiao-Ling. Improved gauss pseudospectral method for solving nonlinear optimal control problem with complex constraints. Acta Automatica Sinica, 2013, 39 (5): 672-678 http://www.aas.net.cn/CN/abstract/abstract17921.shtml
    [13] Biegler L T. An overview of simultaneous strategies for dynamic optimization. Chemical Engineering & Processing Process Intensification, 2007, 46 (11): 1043-1053 https://www.sciencedirect.com/science/article/abs/pii/S0255270107001122
    [14] Biegler L T. Nonlinear programming: concepts, algorithms, and applications to chemical processes. Society for Industrial and Applied Mathematics, 2010. doi: 10.1038-emboj.2011.267/
    [15] 彭海军, 高强, 吴志刚, 钟万勰.求解最优控制问题的混合变量变分方法及其航天控制应用. 自动化学报, 2011, 37 (10): 1248-1255 http://www.aas.net.cn/CN/abstract/abstract17614.shtml

    Peng Hai-Jun, Gao Qiang, Wu Zhi-Gang, Zhong Wan-Xie. A mixed variable variational method for optimal control problems with applications in aerospace control. Acta Automatica Sinica, 2011, 37 (10): 1248-1255 http://www.aas.net.cn/CN/abstract/abstract17614.shtml
    [16] Chachuat B, Mitsos A, Barton P I. Optimal design and steady-state operation of micro power generation employing fuel cells. Chemical Engineering Science, 2005, 60 (16): 4535-4556 doi: 10.1016/j.ces.2005.02.053
    [17] Betts J T, Huffman W P. Application of sparse nonlinear programming to trajectory optimization. Journal of Guidance Control Dynamics, 1992, 15 (1): 198-206 doi: 10.2514/3.20819
    [18] Fu J, Faust J M M, Chachuat B, Mitsosc1 A. Local optimization of dynamic programs with guaranteed satisfaction of path constraints. Automatica, 2015, 62 (C): 184-192 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=f45099105990174d5d267c2b523fe5d9
    [19] 胡云卿, 刘兴高, 薛安克.带不等式路径约束最优控制问题的惩罚函数法. 自动化学报, 2013, 39 (12): 1996-2001 http://www.aas.net.cn/CN/abstract/abstract18238.shtml

    Hu Yun-Qing, Liu Xing-Gao, Xue An-Ke. A penalty method for solving inequality path constrained optimal control problems. Acta Automatica Sinica, 2013, 29 (12): 1996-2001 http://www.aas.net.cn/CN/abstract/abstract18238.shtml
    [20] 钱积新, 宋春跃, 王可心, 陈扬.非线性预测控制, 北京:科学出版社, 2015.

    Qian Ji-Xin, Song Chun-Yue, Wang Ke-Xin, Cheng Yang. Nonlinear Predictive Control, Beijing: Science Press, 2015.
    [21] Pontryagin L S, Boltyanskii V G, Gamkrelidze R V, Mishchenko E F. The mathematical theory of optimal processes. Interscience, 1962.
    [22] Fabien B C. A technique for the direct optimization of dynamic systems described by differential-algebraic equations. Optimal Control Applications & Methods, 2007, 29 (6): 445-466 http://cn.bing.com/academic/profile?id=6f2365ff0f36a2cdd559444ad66b3356&encoded=0&v=paper_preview&mkt=zh-cn
    [23] Brenan K E, Campbell S L, Petzold L R. Numerical solution of initial-value problems in differential-algebraic equations. North-Holland, 1989.
    [24] Styczen K, Drag P. A modified multipoint shooting feasible-SQP method for optimal control of DAE systems. In: Proceedings of Computer Science and Information Systems, New York, USA: IEEE, 2011. 477-484 https://www.researchgate.net/publication/220726423_A_modified_multipoint_shooting_feasible-SQP_method_for_optimal_control_of_DAE_systems
    [25] Feehery W F, Barton P I. Dynamic simulation and optimization with inequality path constraints. Computers & Chemical Engineering, 1996, 20 (3): 169-176 http://d.old.wanfangdata.com.cn/Periodical/kzllyyy-e201202012
    [26] Jacobson D, Lele M. A transformation technique for optimal control problems with a state variable inequality constraint. IEEE Transactions on Automatic Control, 1969, 14 (5): 457-464 doi: 10.1109/TAC.1969.1099283
    [27] Vassiliadis V S, Sargent R W H, Pantelides C C. Solution of a class of multistage dynamic optimization problems. 2. problems with path constraints. Industrial & Engineering Chemistry Research, 1994, 10 (9): 2122-2133 doi: 10.1021-ie00033a015/
    [28] Floudas C A, Stein O. The adaptive convexification algorithm: a feasible point method for semi-infinite programming. SIAM Journal on Optimization, 2008, 18 (4): 1187-1208 doi: 10.1137/060657741
    [29] Chachuat B. Nonlinear and dynamic optimization: from theory to practice. Automatic Control Laboratory EPFL, 2007, 107 (Spring): 192-193 http://d.old.wanfangdata.com.cn/Periodical/sxdsjyrs201701024
    [30] Rehbock V, Teo K L, Jennings L S, Lee C S. An exact penalty function approach to all-time-step constrained discrete-time optimal control problems. Applied Mathematics & Computation, 1992, 49 (2-3): 215-230 http://cn.bing.com/academic/profile?id=17831c1e7ccab4a3cd9152bfb3160291&encoded=0&v=paper_preview&mkt=zh-cn
    [31] Loxton R C, Teo K L, Rehbock V, Yiu K F C. Optimal control problems with a continuous inequality constraint on the state and the control. Automatica, 2009, 45 (10): 2250-2257 doi: 10.1016/j.automatica.2009.05.029
    [32] Liu X, Hu Y, Feng J, Liu K. A novel penalty approach for nonlinear dynamic optimization problems with inequality path constraints. IEEE Transactions on Automatic Control, 2014, 59 (10): 2863-2867 doi: 10.1109/TAC.2014.2317293
    [33] Liu P, Li X, Liu X, Hu Y. An improved smoothing technique-based control vector parameterization method for optimal control problems with inequality path constraints. Optimal Control Applications & Methods, 2016, 38 (4): 586-600 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=e6355d35d514a3baab87112e5634ac3e
    [34] Chen T W C, Vassiliadis V S. Inequality path constraints in optimal control: a finite iteration ε-convergent scheme based on pointwise discretization. Journal of Process Control, 2005, 15 (3): 353-362 doi: 10.1016/j.jprocont.2004.04.002
    [35] Biegler L T. An overview of simultaneous strategies for dynamic optimization. Chemical Engineering & Processing Process Intensification, 2007, 46 (11): 1043-1053 doi: 10.1016-j.cep.2006.06.021/
    [36] Vassiliadis V. Computational solution of dynamic optimization problems with general differential-algebraic constraints by. Journal of Guidance Control & Dynamics, 1993, 15 (2): 457-460
    [37] Teo K L, Goh C J, Wong K H. A unified computational approach to optimal control problems. Longman Scientific and Technical, 1991: 2763-2774 http://cn.bing.com/academic/profile?id=0f2e06643f437ebed6f98391ba8ba19a&encoded=0&v=paper_preview&mkt=zh-cn
    [38] Martin R, Teo K L. Optimal control of drug administration in cancer chemotherapy. World Scientific, 1994. http://cn.bing.com/academic/profile?id=ee8ed3510986556119907e09e7f82d9e&encoded=0&v=paper_preview&mkt=zh-cn
    [39] Gunn D J, Thomas W J. Mass transport and chemical reaction in multifunctional catalyst systems. Chemical Engineering Science, 1965, 20 (2): 89-100 http://cn.bing.com/academic/profile?id=51cbe371f25dc6d294db3f75816b0e9f&encoded=0&v=paper_preview&mkt=zh-cn
    [40] Hirmajer T, Fikar M, Balsa-Canto E, Banga J R. DOTcvp: dynamic optimization toolbox with control vector parameterization approach. 2008.
    [41] Irizarry R. A generalized framework for solving dynamic optimization problems using the artificial chemical process paradigm: applications to particulate processes and discrete dynamic systems. Chemical Engineering Science, 2005, 60 (21): 5663-5681 doi: 10.1016/j.ces.2005.05.028
    [42] Huang Y J, Reklaitis G V, Venkatasubramanian V. Model decomposition based method for solving general dynamic optimization problems. Computers & Chemical Engineering, 2002, 26 (6): 863-873 http://cn.bing.com/academic/profile?id=42713eb59ef18938f24aad61d71e904c&encoded=0&v=paper_preview&mkt=zh-cn
  • 期刊类型引用(3)

    1. 付俊,彭燕,刘彦辉. 具有未知参数的非线性系统动态优化. 控制与决策. 2023(08): 2223-2230 . 百度学术
    2. 张磊,郑丹,卢天林,陈庆,叶婧,熊致知. 考虑工业园区CHP机组动态约束的低碳经济调度. 电力建设. 2022(12): 122-130 . 百度学术
    3. 陈淼,于继明,张静,吴洪兵. 矿仓煤渣装运无人操作抓斗摆动抑制研究. 煤矿机械. 2021(09): 155-158 . 百度学术

    其他类型引用(5)

  • 加载中
  • 图(7) / 表(2)
    计量
    • 文章访问数:  2580
    • HTML全文浏览量:  381
    • PDF下载量:  525
    • 被引次数: 8
    出版历程
    • 收稿日期:  2018-05-14
    • 录用日期:  2018-10-06
    • 刊出日期:  2019-05-20

    目录

    /

    返回文章
    返回