1 前言

带时间延迟最优控制问题(time-delay optimal control problem，TDOCP)的求解是实现时滞工业过程对象最优控制的一个核心^[1-2]。在当前最优控制问题数值求解算法中，控制变量参数化(control variable parameterization，CVP)方法因具有求解精度高、离散化后非线性规划(nonlinear programming，NLP)问题规模较小等优点，在工业过程控制领域得到了众多学者的青睐^[3-5]。近年来，国内外学者提出了诸多改进方法用于TDOCP问题的求解。比如，Gui等^[6]基于CVP方法推导了带时间延迟系统的协态方法，并在锌冶炼过程TDOCP问题中进行仿真试验，结果显示优化计算后锌粉原料可以有效减少；Yu等^[7]研究了混合时间尺度变换(hybrid time-scaling，HTS)HTS-CVP方法进行时间延迟尺度转化，从而实现了TDOCP问题的有效求解；Lin等^[8]提出了2种用于计算TDOCP问题目标函数梯度信息的方法，实现了时滞非线性系统的参数预估。此外，Jajarmi等^[9]基于庞特里亚金极大值原理将TDOCP问题转换为耦合两点边值问题，并推导了有限差分求解方法；文献[10]和文献[11]分别提出了单时间转化和终值时间转换算法用于终值时间不固定TDOCP问题的求解。

对于传统CVP方法，动态系统时间延迟的存在让状态变量和梯度信息的求解变得困难，进而影响到TDOCP问题的高效求解，因此状态变量转化与梯度信息计算一直是研究的热点。然而，以传统CVP方法为基础的目标函数和约束条件梯度求取方法^[6-8]在算法复杂性和效率方面存在一定的不足，比如，协态方法通过构造Hamilton函数得到协态系统进行梯度求取，但迟滞系统和协态系统因为初始值不同并不能同时求解，需要引入插值方法进行近似；时间尺度转化方法引入时间尺度变换可以实现时间节点的精确控制，但也增加了系统梯度信息求解的复杂度。因此，本文提出了一种求解TDOCP问题的扩展控制向量参数化方法。首先，采用控制变量参数化技术和光滑化罚函数法将TDOCP问题近似转化为不含路径约束的有限维NLP问题；进一步，引入扩展变分系统通过辅助状态变量分析时间延迟状态变量扩展变分系统求解方法，以此实现迟滞系统与变分系统在统一初始时刻下求解；同时，针对扩展变分系统直接推导出目标函数和约束函数对于控制参数的梯度信息求解公式，为基于梯度信息的NLP求解器提供高效梯度信息求取方法；最后，在两个典型工业过程TDOCP问题上进行测试以验证所提方法的正确性和有效性。

2 问题描述

考虑控制变量和状态变量同时存在时间延迟的工业过程最优控制问题，其优化过程为对延迟系统中的控制变量实施控制，使目标性能函数达到最优，该类问题数学描述如下(问题P1)：

$\begin{array}{l} \min \;{g_0} = {\mathit{\Phi}_0}(\mathit{\boldsymbol{x}}(T|\mathit{\boldsymbol{u}})) + \int_0^T {{L_0}(\mathit{\boldsymbol{x}}(t - h|\mathit{\boldsymbol{u}}),\mathit{\boldsymbol{x}}(t|\mathit{\boldsymbol{u}}),\mathit{\boldsymbol{u}}(t))} {\rm{d}}t\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\mathit{\boldsymbol{\dot x}}(t) = \mathit{\boldsymbol{f}}[\mathit{\boldsymbol{x}}(t - h),\mathit{\boldsymbol{x}}(t),\mathit{\boldsymbol{u}}(t - h),\mathit{\boldsymbol{u}}(t),t],\;\;t \in [0,T]\\ \;\;\;\;\;\;\mathit{\boldsymbol{x}}(t) = \phi(t),\;\;t \in [ - h,0],\;\;\mathit{\boldsymbol{x}}(0){\rm{ = }}{\mathit{\boldsymbol{x}}_0};\;\;\mathit{\boldsymbol{u}}(t) = \varphi (t),\;\;t \in [ - h,0];\;\;{\mathit{\boldsymbol{u}}_{\min }} \le \mathit{\boldsymbol{u}}(t) \le {\mathit{\boldsymbol{u}}_{\max }}\\ \;\;\;\;\;\;{g_i} = {\mathit{\Phi }_i}(\mathit{\boldsymbol{x}}(T|u)) + \int_0^T {{L_i}(\mathit{\boldsymbol{x}}(t - h|u),\mathit{\boldsymbol{x}}(t|u),\mathit{\boldsymbol{u}}(t))} {\rm{d}}t = 0,\;\;i \in E\\ \;\;\;\;\;\;{g_j} = {\mathit{\Phi }_j}(\mathit{\boldsymbol{x}}(T|u)) + \int_0^T {{L_j}(\mathit{\boldsymbol{x}}(t - h|u),\mathit{\boldsymbol{x}}(t|u),\mathit{\boldsymbol{u}}(t))} {\rm{d}}t \ge 0,\;\;j \in I \end{array} $

(1)

式中：$\mathit{\boldsymbol{x}}(t) \in {\mathbb{R}^r}$为状态变量，$\mathit{\boldsymbol{u}}(t) \in {\mathbb{R}^n}$为控制变量(其中r和n代表维度)，h为给定时间延迟且0<h<T；目标函数${g_0}$由终值项$ \mathit{\Phi}_{0}{}_{0}(·)$和积分项$ {\displaystyle {\int }_{0}^{T}{L}_{0}(·)}{\rm{d}}t$构成；$\mathit{\boldsymbol{f}}(·)$和$\phi (t)$表示时间延迟动态系统，由给定连续微分方程组表示，${\mathit{\boldsymbol{x}}_0}$是$\mathit{\boldsymbol{x}}(t)$在零时刻初始可行状态，$\varphi (t)$是延迟控制变量给定分段常量方程；${g_i}(i \in E), {\rm{ }}{g_j}(j \in I)$表示动态系统约束条件，$E$和$I$分别为规范等式约束和规范不等式约束序号集；${\mathit{\boldsymbol{u}}_{\max }}$和${\mathit{\boldsymbol{u}}_{\min }}$是控制变量上下边界约束。

3 算法描述 3.1 控制变量参数化

针对问题(P1)，将控制时域$[{t_0}, {t_f}]$划分成$N$个子区间，记各区间时间节点为：${t_0} < {t_1} < \cdots < {t_N} = {t_f}$，在时间段$ [{t}_{k-1}, {t}_{k}]{，}_{}(k=1, 2, \cdots , N)$内将控制变量${u_j}(t)$, $(j = 1, 2, ..., {n_u}{\rm{)}}$近似表示为一组参数的函数：

$u_j^k(t) \approx \tilde u_j^k(t) = \sum\limits_{r = 1}^{{M_{jk}} + 1} {{\sigma _{jkr}}} \phi _{jkr}^{({M_{jk}})}(t){\rm{ }}\;\;t \in [{t_{k - 1}}, {t_k})$

(2)

式中：${\sigma _{jkr}}$为线性组合系数，为$u_j^k(t)$的控制参数，$\phi _{jkr}^{({M_{jk}})}(t)$为基函数，${M_{jk}}$为基函数的阶。本文选取分段常量函数(piecewise constant function，PCF)进行基函数$\phi _{jkr}^{({M_{jk}})}(t)$逼近，得到$u_j^k(t)$为

${u_j}(t) \approx {\tilde u_j}(t) = \sum\limits_{k = 1}^N {{\sigma _{j, k}}{\chi _{\left[ {{t_{k - 1}}, {t_k}} \right)}}(t)} , {\rm{ }}t \in [{t_{k - 1}}, {t_k})$

(3)

${\chi _{_{\left[ {{t_{k - 1}}, {t_k}} \right)}}}(t) = \left\{ \begin{gathered} 1, {\rm{ }}t \in \left[ {{t_{k - 1}}, {t_k}} \right), {\rm{ }} \\ 0, {\rm{ }}t \notin \left[ {{t_{k - 1}}, {t_k}} \right). \\ \end{gathered} \right.$

(4)

式中：${\sigma _{j, k}}$为控制参数，${\chi _{\left[ {{t_{k - 1}}, {t_k}} \right)}}(t)$为分段常量函数。从而，记第$j$个控制分量的参数依次为${\sigma _{j, 1}}, {\sigma _{j, 2}}, \ldots , {\sigma _{j, N}}$，定义向量${\mathit{\boldsymbol{\sigma}} _k} = {[{\sigma _{1, k}}, {\sigma _{2, k}}, \ldots , {\sigma _{{n_u}, k}}]^{\rm{T}}}$，则有$\mathit{\boldsymbol{\sigma}} = {[\mathit{\boldsymbol{\sigma}} _1^{\rm{T}}, \mathit{\boldsymbol{\sigma}} _2^{\rm{T}}, \ldots , \mathit{\boldsymbol{\sigma}} _n^{\rm{T}}]^{\rm{T}}}$；进一步定义${t_\rm{delay}} = t - h$，则问题(P1)经过控制变量参数化后可以转化为如下优化控制参数为$\mathit{\boldsymbol{\sigma}} $的最优控制问题(问题P2)：

$\begin{gathered} \min \;{g_0}(\mathit{\boldsymbol{\sigma}} ) = {\mathit{\Phi} _0}(\mathit{\boldsymbol{x}}(T|\mathit{\boldsymbol{\sigma}} )) + \sum\limits_{k = 1}^N {\int_{{t_{k - 1}}}^{{t_k}} {{L_0}(\mathit{\boldsymbol{x}}({t_\rm{delay}}|\mathit{\boldsymbol{\sigma}} ), x(t|\mathit{\boldsymbol{\sigma}} ), \mathit{\boldsymbol{\sigma}} ){\rm{d}}t} } \\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\mathit{\boldsymbol{\dot x}}(t) = \left\{ \begin{gathered} \mathit{\boldsymbol{f}}[\mathit{\boldsymbol{x}}({t_\rm{delay}}), \mathit{\boldsymbol{x}}(t), {\mathit{\boldsymbol{\sigma}} _j}, {\mathit{\boldsymbol{\sigma}} _k}, t], \;\;t \in [{t_{k - 1}}, {t_k}], \;\;{t_\rm{delay}} \in [{t_{j - 1}}, {t_j}] \\ \mathit{\boldsymbol{f}}[\phi (t), \mathit{\boldsymbol{x}}(t), \varphi (t), {\mathit{\boldsymbol{\sigma}} _k}, t], \;\;t < h, \;\;t \in [{t_{k - 1}}, {t_k}] \\ \end{gathered} \right. \\ \;\;\;\;\;\mathit{\boldsymbol{x}}(t) = \phi (t), \;\;t \in [ - h, 0];\;\;\mathit{\boldsymbol{x}}(0) = {\mathit{\boldsymbol{x}}_0}, \;\;\mathit{\boldsymbol{u}}(t) = \varphi (t);\;\;t \in [ - h, 0], \;\;{\mathit{\boldsymbol{u}}_{\min }} \leqslant \mathit{\boldsymbol{\sigma}} \leqslant {\mathit{\boldsymbol{u}}_{\max }} \\ \;\;\;\;\;{g_i}(\mathit{\boldsymbol{\sigma}} ) = {\mathit{\Phi} _i}(\mathit{\boldsymbol{x}}(T|\mathit{\boldsymbol{\sigma}} )) + \sum\limits_{k = 1}^N {\int_{{t_{k - 1}}}^{{t_k}} {{L_i}(\mathit{\boldsymbol{x}}({t_\rm{delay}}|\mathit{\boldsymbol{\sigma}} ), x(t|\mathit{\boldsymbol{\sigma}} ), \mathit{\boldsymbol{\sigma}} ){\rm{d}}t} } = 0, \;\;i \in E \\ \;\;\;\;\;{g_j}(\mathit{\boldsymbol{\sigma}} ) = {\mathit{\Phi} _j}(\mathit{\boldsymbol{x}}(T|\mathit{\boldsymbol{\sigma}} )) + \sum\limits_{k = 1}^N {\int_{{t_{k - 1}}}^{{t_k}} {{L_j}(\mathit{\boldsymbol{x}}({t_\rm{delay}}|\mathit{\boldsymbol{\sigma}} ), \mathit{\boldsymbol{x}}(t|\mathit{\boldsymbol{\sigma}} ), \mathit{\boldsymbol{\sigma}} ){\rm{d}}t} } \leqslant 0, \;\;j \in I \\ \end{gathered} $

(5)

有关控制变量参数化后问题(P2)与最优控制问题(P1)的收敛性可由定理1保证，鉴于定理1已经在文献[12]中进行了详细证明，本文省略其证明过程，给出具体定理描述如下：

定理1.设${\mathit{\boldsymbol{u}}^ * }(t)$是问题(P1)的最优解，在一定控制时域分段数N $(N \geqslant 1)$下，假设${\mathit{\boldsymbol{\sigma}} ^ * }$是问题(P2)的最优解，${\mathit{\boldsymbol{\tilde u}}^ * }(t)$是由${\mathit{\boldsymbol{\sigma}}^ * }$在分段常量函数逼近下得到的控制变量，随着分段数$N \to \infty $，如果最大时间间隔$[{t_{k - 1}}, {t_k}]$$(k = 1, 2, \cdots , N)$趋近于0，则有：

$\mathop {\lim }\limits_{N \to \infty } {g_0}({\mathit{\boldsymbol{\tilde u}}^{\rm{*}}}(t)) = {g_0}({\mathit{\boldsymbol{\tilde u}}^ * }(t))$

(6)

经过控制变量参数化后，问题(P2)的目标函数值将收敛到问题(P1)的目标函数值，但经过控制变量参数化处理后，问题(P2)依旧包含不等式路径约束，这使得问题(P2)求解变得困难。因此，结合文献[13]提出的光滑化罚函数法，本文引入max函数，将不等式路径约束${g_j}(\mathit{\boldsymbol{\sigma}} ) \leqslant 0, \;\;j \in I$转化为如式(7)：

$\max \{ {g_j}(\mathit{\boldsymbol{\sigma}} ), 0\} = 0{\rm{\;\;\;}}j \in I$

(7)

通过引入文献[14]提出的光滑化函数，式(7)可以通过定理2近似得到，其证明过程可参阅文献[12]。

定理2.对于不等式路径约束$\max \{ {g_j}(\mathit{\boldsymbol{\sigma}} ), 0\} = 0$，$j \in I$，当光滑化因子$\varepsilon \to {0^ + }$时，

$\max \{ {g_j}(\mathit{\boldsymbol{\sigma}} ), 0\} \approx p\{ {g_j}(\mathit{\boldsymbol{\sigma}} ), \varepsilon \} = {g_j}(\mathit{\boldsymbol{\sigma}} ) + \sqrt {{g_j}{{(\mathit{\boldsymbol{\sigma}} )}^2} + 4{\varepsilon ^2}} )/2$

(8)

进一步，引入路径约束惩罚因子$\rho $将式(8)增广到目标函数中，则问题(P2)中不等式路径约束将转化到扩展后新目标函数中：

$\min \;{\bar g_0}(\mathit{\boldsymbol{\sigma}} ) = {\mathit{\Phi}_0}(x(T|\mathit{\boldsymbol{\sigma}} )) + \sum\limits_{k = 1}^N {\int_{{t_{k - 1}}}^{{t_k}} {{L_0}(x({t_\rm{delay}}|\mathit{\boldsymbol{\sigma}} ), x(t|\mathit{\boldsymbol{\sigma}} ), \mathit{\boldsymbol{\sigma}} ){\rm{d}}t} } + \rho \sum\limits_{j \in I} {\sum\limits_{k = 1}^N {\int_{{t_{k - 1}}}^{{t_k}} {p\{ {g_j}(\mathit{\boldsymbol{\sigma}} ), \varepsilon \} {\rm{d}}t} } } $

(9)

设$\varsigma $为误差容限值，且$\varsigma = \left| {{{\bar g}_0}(\mathit{\boldsymbol{\mathit{\boldsymbol{\sigma}}}} ) - {g_0}(\mathit{\boldsymbol{\sigma}} )} \right|$，则光滑化函数$p\{ {g_j}(\sigma ), \varepsilon \} $中光滑化因子$\varepsilon $的选取公式为：$\varepsilon=\zeta /\left[6(1+\sqrt{5})\left(t_{f}-t_{0}\right)\right]$，惩罚因子通过$\rho \geqslant {\rho ^ * }$选取，其中${\rho ^ * } = {{{10}/ \varsigma }}$，以此保证不等式路径约束近似精度高于问题(P2)求解精度。限于篇幅，光滑化近似函数与$\max \{ {g_j}(\mathit{\boldsymbol{\sigma}} ), 0\} $函数的误差分析请参阅文献[15]。

3.2 时间延迟状态变量扩展变分系统

分析可知，由于时间延迟微分方程组的存在，问题(P2)中目标函数和约束函数对于控制参数$\mathit{\boldsymbol{\sigma}} $的梯度信息并不能直接得到，进而不利于基于梯度信息非线性优化算法的精确求解。基于此，本文将常规变分系统扩展到时间延迟状态变量变分系统上进行梯度信息求取。首先，定义辅助状态变量：

$\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\left( t \right) = \partial \mathit{\boldsymbol{x}}\left( t \right)/\partial \mathit{\boldsymbol{\sigma , \boldsymbol{\varGamma} }}\left( {{t_{{\rm{delay}}}}} \right) = \partial \mathit{\boldsymbol{x}}\left( {{t_{{\rm{delay}}}}} \right)/\partial \mathit{\boldsymbol{\sigma }} $

(10)

则带时间延迟状态变量变分系统可由定理3和定理4推导得到。

定理3.假设$\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} (t)$、$\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} ({t_{delay}})$为时间延迟系统$f[x({t_{delay}}), x(t), {\mathit{\boldsymbol{\sigma}} _j}, {\mathit{\boldsymbol{\sigma}} _k}, t]$的辅助状态变量，则$\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} (t)$可由以下扩展变分系统求得：

$\begin{array}{l} \mathit{\boldsymbol{ \boldsymbol{\dot \varGamma} }}\left( t \right) = \partial \mathit{\boldsymbol{f}}\left[ {\mathit{\boldsymbol{x}}\left( {{t_{{\rm{delay}}}}} \right),\mathit{\boldsymbol{x}}\left( t \right),{\mathit{\boldsymbol{\sigma }}_j},{\mathit{\boldsymbol{\sigma }}_k}} \right]/\partial \mathit{\boldsymbol{x}}\left( t \right) \times \mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\left( t \right) \times \partial \mathit{\boldsymbol{f}}\left[ {\mathit{\boldsymbol{x}}\left( {{t_{{\rm{delay}}}}} \right),\mathit{\boldsymbol{x}}\left( t \right),{\mathit{\boldsymbol{\sigma }}_j},{\mathit{\boldsymbol{\sigma }}_k},t} \right]/\partial \mathit{\boldsymbol{x}}\left( {{t_{{\rm{delay}}}}} \right) \times \mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\left( {{t_{{\rm{delay}}}}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\partial \mathit{\boldsymbol{f}}\left[ {\mathit{\boldsymbol{x}}\left( {{t_{{\rm{delay}}}}} \right),\mathit{\boldsymbol{x}}\left( t \right),{\mathit{\boldsymbol{\sigma }}_j},{\mathit{\boldsymbol{\sigma }}_k},t} \right]/\partial \mathit{\boldsymbol{\sigma }}\\ \mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\left( 0 \right) = {\bf{0}} \end{array} $

(11)

证明.由问题(P2)可知，$x(t) = {x_0}{\rm{ + }}\int_{{t_0}}^{{t_f}} {f[x({t_\rm{delay}}), x(t), {\mathit{\boldsymbol{\sigma}} _j}, {\mathit{\boldsymbol{\sigma}} _k}, t] {\rm{d}}t} , \;\;t > 0$，因此有：

$\frac{\partial \boldsymbol{x}(t)}{\partial \mathit{\boldsymbol{ \boldsymbol{\sigma} }}}=\frac{\boldsymbol{x}_{0}}{\partial \mathit{\boldsymbol{ \boldsymbol{\sigma} }}}+\int_{0}^{t_{f}}\left\{\partial \boldsymbol{f}\left[\boldsymbol{x}\left(t_{\text {delay }}\right), \boldsymbol{x}(t), \mathit{\boldsymbol{ \boldsymbol{\sigma} }}_{j}, \mathit{\boldsymbol{ \boldsymbol{\sigma} }}_{k}, t\right] / \partial \boldsymbol{x}(t) \times \mathit{\boldsymbol{ \boldsymbol{\varGamma} }}(t)+\right.\\ \left.\partial \boldsymbol{f}\left[\boldsymbol{x}\left(t_{ {delay}}\right), \boldsymbol{x}(t), \mathit{\boldsymbol{ \boldsymbol{\sigma} }}_{j}, \mathit{\boldsymbol{ \boldsymbol{\sigma} }}_{k}, t\right] / \partial \boldsymbol{x}\left(t_{\text {delay }}\right) \times \mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\left(t_{\text {delay }}\right)+\partial \boldsymbol{f}\left[\boldsymbol{x}\left(t_{\text {delay }}\right), \boldsymbol{x}(t), \mathit{\boldsymbol{ \boldsymbol{\sigma} }}_{j}, \mathit{\boldsymbol{ \boldsymbol{\sigma} }}_{k}, t\right] / \partial \mathit{\boldsymbol{ \boldsymbol{\sigma} }}\right\} \mathrm{d} t $

(12)

式(12)对时间t求导得到：

$\frac{{\rm{d}}}{{{\rm{d}}t}}\left\{ {\frac{{\partial x(t)}}{{\partial \mathit{\boldsymbol{\sigma}} }}} \right\} = \frac{{\partial f[x({t_\rm{delay}}), x(t), {\mathit{\boldsymbol{\sigma}} _j}, {\mathit{\boldsymbol{\sigma}} _k}, t]}}{{\partial x(t)}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} (t) + \frac{{\partial f[x({t_\rm{delay}}), x(t), {\mathit{\boldsymbol{\sigma}} _j}, {\mathit{\boldsymbol{\sigma}} _k}, t]}}{{\partial x({t_\rm{delay}})}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} ({t_\rm{delay}}) + \frac{{\partial f[x({t_\rm{delay}}), x(t), {\mathit{\boldsymbol{\sigma}} _j}, {\mathit{\boldsymbol{\sigma}} _k}, t]}}{{\partial \mathit{\boldsymbol{\sigma}} }}$

(13)

可知，$\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} (t)$为式(13)的解，又${\mathit{\boldsymbol{x}}_0}$为已知，则$\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}\left( 0 \right) = \partial {\mathit{\boldsymbol{x}}_0}/\partial \mathit{\boldsymbol{\sigma = }}{\bf{0}}$，得证。

定理4.定义${\tilde g_l}(\mathit{\boldsymbol{\sigma}} ), \;l \in \{ 0\} \cup E$为目标函数、等式约束的Canonical函数，其中${\mathit{\tilde \Phi }_l}(x(T|\mathit{\boldsymbol{\sigma }}))$表示终端约束项，${\tilde L_l}(x({t_\rm{delay}}|\mathit{\boldsymbol{\sigma}} ), x(t|\mathit{\boldsymbol{\sigma}} ), \mathit{\boldsymbol{\sigma}} )$表示积分项，则${\tilde g_l}(\mathit{\boldsymbol{\sigma}} )$关于控制参数$\mathit{\boldsymbol{\sigma}} $的梯度信息是：

$\begin{gathered} \frac{{\partial {{\tilde g}_l}(\mathit{\boldsymbol{\sigma}} )}}{{\partial \mathit{\boldsymbol{\sigma}} }} = \frac{{\partial {\mathit{\tilde \Phi }_l}(x(T|\mathit{\boldsymbol{\sigma}} ))}}{{\partial x}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} (T) + \frac{{\partial {\mathit{\tilde \Phi }_l}(x(T|\mathit{\boldsymbol{\sigma}} ))}}{{\partial \mathit{\boldsymbol{\sigma}} }} + \sum\limits_{k = 1}^N {\int_{{t_{k - 1}}}^{{t_k}} {\{ \frac{{\partial {{\tilde L}_l}(x({t_\rm{delay}}|\mathit{\boldsymbol{\sigma}} ), x(t|\mathit{\boldsymbol{\sigma}} ), \mathit{\boldsymbol{\sigma}} )}}{{\partial x}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} (t)} } + \\ \frac{{\partial {{\tilde L}_l}(x({t_\rm{delay}}|\mathit{\boldsymbol{\sigma}} ), x(t|\mathit{\boldsymbol{\sigma}} ), \mathit{\boldsymbol{\sigma}} )}}{{\partial x({t_\rm{delay}}|\mathit{\boldsymbol{\sigma}} )}}\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} ({t_\rm{delay}}) + \frac{{\partial {{\tilde L}_l}(x({t_\rm{delay}}|\mathit{\boldsymbol{\sigma}} ), x(t|\mathit{\boldsymbol{\sigma}} ), \mathit{\boldsymbol{\sigma}} )}}{{\partial \mathit{\boldsymbol{\sigma}} }}\} {\rm{d}}t \\ \end{gathered} $

(14)

证明.采用链式法则对函数${\tilde g_l}(\mathit{\boldsymbol{\sigma}} )$关于$\mathit{\boldsymbol{\sigma}} $求导可得：

$\begin{gathered} \frac{{\partial {{\tilde g}_l}(\mathit{\boldsymbol{\sigma}} )}}{{\partial \mathit{\boldsymbol{\sigma}} }} = \frac{{\partial {\mathit{\tilde \Phi }_l}(x(T|\mathit{\boldsymbol{\sigma}} ))}}{{\partial x}}\frac{{\partial x}}{{\partial \mathit{\boldsymbol{\sigma}} }} + \frac{{\partial {\mathit{\tilde \Phi }_n}(x(T|\mathit{\boldsymbol{\sigma}} ))}}{{\partial \mathit{\boldsymbol{\sigma}} }} + \sum\limits_{k = 1}^N {\int_{{t_{k - 1}}}^{{t_k}} {\{ \frac{{\partial {{\tilde L}_l}(x({t_\rm{delay}}|\mathit{\boldsymbol{\sigma}} ), x(t|\mathit{\boldsymbol{\sigma}} ), \mathit{\boldsymbol{\sigma}} )}}{{\partial x}}\frac{{\partial x}}{{\partial \mathit{\boldsymbol{\sigma}} }}} } + \\ \frac{{\partial {{\tilde L}_l}(x({t_\rm{delay}}|\mathit{\boldsymbol{\sigma}} ), x(t|\mathit{\boldsymbol{\sigma}} ), \mathit{\boldsymbol{\sigma}} )}}{{\partial x({t_\rm{delay}}|\mathit{\boldsymbol{\sigma}} )}}\frac{{\partial x({t_\rm{delay}}|\mathit{\boldsymbol{\sigma}} )}}{{\partial \mathit{\boldsymbol{\sigma}} }} + \frac{{\partial {{\tilde L}_l}(x({t_\rm{delay}}|\mathit{\boldsymbol{\sigma}} ), x(t|\mathit{\boldsymbol{\sigma}} ), \mathit{\boldsymbol{\sigma}} )}}{{\partial \mathit{\boldsymbol{\sigma}} }}\} {\rm{d}}t \\ \end{gathered} $

(15)

将式(10)入式(15)，即可得到定理4式(14)，证毕。

3.3 算法步骤

结合3.1节和3.2节算法推导，给出时间延迟最优控制问题扩展控制变量参数化算法(extended control variable parameterization, Extended-CVP)实现过程，其算法流程如图 1所示，主要步骤为：

步骤1.输入带时间延迟最优控制问题(P1)，设置控制时域分段数$N$，给定初始控制参数${\mathit{\boldsymbol{\sigma}} _0}$；

步骤2.在分段数$N$下，采用PCF函数对控制变量进行离散化，得到参数化后非线性优化问题(P2)；

步骤3.求解时间延迟状态变量扩展变分系统，得到辅助状态变量$\mathit{\boldsymbol{ \boldsymbol{\varGamma} }} (t)$的值；

步骤4.采用式(14)得到${\tilde g_l}(\mathit{\boldsymbol{\sigma}} ), \;l \in \{ 0\} \cup E$函数关于控制参数$\mathit{\boldsymbol{\sigma}} $的扩展梯度信息；

步骤5.选择NLP求解器求解问题(P2)，得到优化后最优控制参数${\mathit{\boldsymbol{\sigma}} ^ * }$和目标函数${g_0}({\mathit{\boldsymbol{\sigma}} ^ * })$，输出优化结果。

4 仿真测试 4.1 闭环连续搅拌釜反应器最优控制

该连续搅拌釜反应器生产模型由Soliman和Ray^[16]提出，其生产过程如图 2所示，图中${x_1}(t)$和${x_2}(t)$表示产物浓度和催化剂浓度，${x_3}(t - r)$表示反应物出口温度，反应过程主要控制反应器热交换量${u_1}(t)$和催化剂补料${u_2}(t)$，其中催化剂补料被分成了$\gamma {u_2}(t)$直接输入和带时间延迟混合输入$(1 - \gamma ){u_2}(t - s)$两个部分，同时，${u_1}(t)$与${x_3}(t - r)$组成比例(proportion，P)控制器实现对反应器内部温度的控制。该过程最优控制的目标是在$\gamma = 0.9$下使反应器在0.2 h内保持在最优平衡状态，其TDOCP问题数学模型表述如下：

$\begin{gathered} \min \;{g_0}(u) = \int_0^{0.2} {(\left\| {\mathit{\boldsymbol{x}}(t)} \right\|_2^2 + 0.01\left\| {\mathit{\boldsymbol{u}}(t)} \right\|_2^2){\rm{d}}t} \\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;{{\dot x}_1}(t) = - {x_1}(t) - R(t), \;\;{{\dot x}_2}(t) = - {x_2}(t) + 0.9{u_2}(t - s) + 0.1{u_2}(t) \\ \;\;\;\;\;\;\;\;{{\dot x}_3}(t) = - 2{x_3}(t) + 0.25R(t) - 1.05{u_1}(t){x_3}(t - r) \\ \;\;\;\;\;\;\;\;R(t) = (1 + {x_1}(t))(1 + {x_2}(t)){\exp _{}}(25{x_3}(t)/(1 + {x_3}(t))) \\ \;\;\;\;\;\;\;\;{x_3}(t) = - 0.02, \;\; - r \leqslant t \leqslant 0, \;\;{u_2}(t) = 1, \;\; - s \leqslant t \leqslant 0, \;\;|{u_1}| \leqslant 5 \\ \;\;\;\;\;\;\;\;r = 0.015, \;\;s = 0.02, \;\;\mathit{\boldsymbol{x}}(0) = [0.49\; - 0.0002\; - 0.02] \\ \end{gathered} $

(16)

采用Extended-CVP方法对该问题进行求解，控制时域采用均匀时间网格划分，分别设置20和40两种分段数，初始控制参数取${\mathit{\boldsymbol{\sigma}} _0} = 0.2$，NLP求解器优化精度为${10^{ - 4}}$。表 1列出了Yu等^[7]提出的混合时间尺度变换HTS-CVP方法、Hirmajer等^[17]提出的DOTCVP优化软件、Betts等^[18]的MOS(method of steps)优化工具包和Campbell等^[19]实现的直接转换法(direct transcription，DT)与本文方法的结果对比。分析表 1可知，与HTS-CVP方法和DOTCVP优化软件相比，本文Extended-CVP方法在20和40分段数下均取得了与文献结果一致的最优目标函数值0.021 3，显示出本文方法在优化计算上的正确性。同时，在40个优化参数下本文方法耗时5.2 s取得了与MOS方法160个优化参数一样的优化结果，较少的优化参数表明本文扩展变分系统在TDOCP问题计算上的有效性。图 3、4分别给出了40分段数下Extended-CVP方法得到的最优控制曲线和状态曲线，与文献结果对比一致，说明本文方法的正确性。

4.2 状态延迟LQR最优控制

本例以Jajarmi等^[20]研究的带状态延迟线性二次型调节器(linear quadratic regulator，LQR)为对象进行测试，其TDOCP问题模型如下所示，模型中所采用的参数数值见表 2。表中a、c为状态变量x(t)输入矩阵的取值参数，b为延迟状态变量x(t-h)输入矩阵的取值参数，h为时间延迟量，t_f为最优控制终值时间，S为终值状态加权矩阵，Q为状态变量加权矩阵，R为控制变量加权矩阵。

$\begin{gathered} \min \;{g_0}(u) = \frac{1}{2}{x^{\rm{T}}}({t_f})\mathit{\boldsymbol{Sx}}({t_f}) + \frac{1}{2}\int_0^{{t_f}} {({x^{\rm{T}}}(t)\mathit{\boldsymbol{Qx}}(t){\rm{ + }}{\mathit{\boldsymbol{u}}^{\rm{T}}}(t)\mathit{\boldsymbol{Ru}}(t)){\rm{d}}t} \\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;{\mathit{\boldsymbol{\dot x}}}(t) = {A_1}(t)\mathit{\boldsymbol{x}}(t) + {A_2}(t)\mathit{\boldsymbol{x}}(t - h) + B(t)\mathit{\boldsymbol{u}}(t);\;\;\mathit{\boldsymbol{x}}(t) = {[1, 0]^{\rm{T}}}, \;\; - h \leqslant t \leqslant 0 \\ \;\;\;\;\;\;{A_1}(t) = \left[ {\begin{array}{*{20}{c}} 0&1 \\ { - 4{\pi ^2}(a + c\cos 2\pi t)}&0 \end{array}} \right], \;\;{A_2}(t) = \left[ {\begin{array}{*{20}{c}} 0&0 \\ { - 4{\pi ^2}b\cos 2\pi t}&0 \end{array}} \right], \;\;B(t) = {\left[ {\begin{array}{*{20}{c}} 0&1 \end{array}} \right]^{\rm{T}}} \\ \end{gathered} $

(17)

测试中，控制时域采用等间隔方式划分，分段数设置为21、40和150，初始控制参数取${\mathit{\boldsymbol{\sigma}} _0} = 1$，NLP求解器优化精度设置为${10^{ - 4}}$。对于上述LQR最优控制问题，文献[20]以CVP方法为基础推导了循环单射(recursive shooting，RS)RS-CVP方法，通过8次循环迭代得到了18.110 3的当前文献最优目标函数值；文献[9]提出了改进的有限差分(finite difference，FD)方法，在21个优化参数下得到了18.450 2的优化结果；文献[21]推导了迭代对偶谱方法(iterative symplectic pseudospectral method，ISPM)，在双重配点下计算得到的目标函数值为18.347 6。表 3列出了文献方法和本文方法的求解结果，比较可知，在21个优化参数下，本文方法求得的优化结果为18.446 4，优于FD方法在相同优化参数下的结果，虽然本文方法的优化结果相比于RS-CVP方法存在0.336 1的差距，但是本文方法只需要一次参数化迭代；同时，与文献[21]所求得结果进行对比发现，本文方法与ISPM方法的目标函数、系统状态曲线一致，显示出本文方法的实用性和正确性。图 5、6分别给出了40分段数下Extended-CVP方法求得的LQR问题最优控制曲线和状态曲线，与文献结果一致，同时也很好地满足了控制要求，进一步表明了本文方法的效果。

5 结论

本文提出了一种用于求解工业过程TDOCP问题的扩展控制变量参数化算法。该方法通过引入辅助状态变量推导了时间延迟状态变量扩展变分系统，拓展了传统CVP方法对于带时间延迟最优控制问题的求解能力；同时，给出了目标函数和约束函数对控制参数梯度信息的求解方法，为基于梯度信息的NLP求解器高效求解提供了支撑，保证了优化求解的精度和效率。所提出的方法针对连续搅拌釜反应器和LQR控制器两个典型TDOCP问题进行了测试，结果显示提出方法能在较少优化参数下取得与文献结果相近甚至更优的结果，表明了本文方法的有效性和正确性。考虑非均匀时间间隔控制变量参数化和不确定参数进一步提升本文方法对TDOCP问题的适用性是下一步工作计划开展的课题。