时滞微分方程是具有时间滞后的微分方程, 它用于描述既依赖于当前状态也依赖于过去状态的发展系统.对于一些事物的发展趋势, 若不考虑时滞, 则这样的模型是毫无意义的.因此对时滞系统模型的研究是有十分重要意义的课题.当对参数变化情况掌握之后, 其对现实生活有着指导性的作用, 它能更准确的刻画实际问题.在流体力学、电子学、种群生态学、核反应堆动力学、经济学及现代控制理论等领域中存在着大量的时滞微分方程, 例如, 经济学中价值法则的作用是由于生产与消费之间的时滞造成的; 在生态系统中, 为了更真实地反映自然, 时滞常是一种不应忽略的因素等等^[1].

自从Mackey和Glass^[2]发现时滞系统中的混沌现象以来, 时滞混沌系统引起了人们浓厚的兴趣. Mackey-Glass混沌时滞微分方程具有如下的形式:

$ \frac{{dy(t)}}{{dt}} = - {p_1}y(t) + \frac{{{p_2}y(t - \tau )}}{{1 + {y^{{p_3}}}(t - \tau )}}, $

(1)

初始条件为

$ y({t_0}) = {y_0}. $

此方程最初是用于描述白血球繁殖的模型, 后来成为了混沌理论中超混沌系统的典型代表, 其中$y$代表循环白血球的浓度, ${p_1}, {p_2}, {p_3}$是参数, $\tau $为时滞参数, ${y_0}$为初始值, 通常取${p_1}=0.1, {p_2}=0.2, {p_3}=10$. Famer^[3]对Mackey-Glass方程的行为特性做了深入的研究, 当$\tau > 17$时表现出混沌性, $\tau $的值越大, 其混沌程度就越高.估计时滞混沌系统的未知参数是混沌控制同步中必须解决的关键问题.近几十年来, 关于参数反演的研究十分广泛, 众多学者相继提出了许多有效的方法, 这些方法^{[4, 6]}解反演问题时大都是把反演转化为优化问题, 而微分进化算法是求解最优化问题的一种非常有效的方法, 它属于全局最优化算法, 具有鲁棒性强, 收敛速度快, 计算精度高的优点.据文献报告, 微分进化算法在收敛速度和稳定性方面都超过了其他几种知名的随机算法, 如退火单纯形算法(annealed nelder and mend strategy, 简称ANM)、自适应模拟退火算法(adaptive simulate annealing, 简称ASA)、演化策略(evolution strategies, 简称ES)和随机微分算法(stochastic differential equations).本文以微分进化算法为基础, 讨论了Mackey-Glass混沌时滞微分方程的参数反演问题, 并对参数的灵敏度进行了详细的分析, 最后给出了数值模拟, 结果表明了该算法的可行性与有效性.

考虑如下形式的Mackey-Glass混沌时滞微分方程

$ \frac{{dy(t)}}{{dt}} = - {p_1}y(t) + \frac{{{p_2}y(t - 17)}}{{1 + {y^{10}}(t - 17)}}, $

(2)

初始条件为

$ y({t_0}) = {p_3}, $

(3)

其中${p_1}, {p_2}, {p_3}$为参数.如果${p_1}, {p_2}, {p_3}$已知, 通过求解(2) 便可得到该方程的变化规律.而在实际问题中, ${p_1}, {p_2}, {p_3}$常常是未知的, 如何根据上述方程及其它已知信息确定未知参数, 便构成了一类参数识别反演问题.

为了反演参数${p_1}, {p_2}, {p_3}, $一般要附加先验条件, 不妨取解在有限个时间点上的值是已知的, 例如, $y({t_l})$为已知, $l = 1, 2, \cdots, r$.这时, 参数反演问题就可以表示为求下面的最小值问题

$ \mathop {\min }\limits_{({p_1}, {p_2}, {p_3}) \in {R^3}} f({p_1}, {p_2}, {p_3}) = \sum\limits_{l = 1}^r {{{[y({t_l}, {p_1}, {p_2}, {p_3})-y({t_l})]}^2}}, $

(4)

其中$y({t_l}, {p_1}, {p_2}, {p_3})$是以${p_1}, {p_2}, {p_3}$为参数求解原问题(2) 获得的相应的计算值.微分进化算法是通过种群内个体间的合作与竞争来实现对优化问题的求解, 其本质上是一种基于实数编码的具有保优思想的进化算法.该算法得到了广泛关注, 因为它的实现技术简单, 在对各种测试问题的实验中表现优异.对于连续变量的函数优化, 微分进化算法能更快、更稳定地收敛到问题的全局最优解^{[7, 8]}.

微分进化算法求解非线性最小二乘问题$\mathop {\min }\limits_{({p_1}, {p_2}, {p_3}) \in {R^3}} f({p_1}, {p_2}, {p_3})$的算法流程如下:

步骤 1  种群初始化

首先, 输入演化参数:种群规模N, 最大演化代数${\rm sub}\_{\rm gen}$, 交叉概率${P_c}$, 交叉因子$xf \in (0, 1)$, 自变量的下界$lb$和上界$ub$; 并令演化代数$t = 0$, 随机生成初始种群$Z(0) = \{ {P_1}, {P_2}, \cdots, {P_N}\} $, 其中${P_i} = \{ {p_{i1}}, {p_{i2}}, {p_{i3}}\} $.

随机生成初始种群的具体方法如下:

$ {p_{ij}}(t) = {p_{\min, j}} + {\rm rand}_{j}(0, 1)({p_{\max, j}} - {p_{\min, j}}), ~~~ i = 1, 2, \cdots, N;~~~j = 1, 2, 3, $

(5)

其中$3$为变量空间的维数, 也就是所需要反演参数的个数; ${\rm rand}_{j}(0, 1)$是[0, 1]区间满足正态分布的随机数; $[{p_{\min, j}}, {p_{\max, j}}]$是搜索空间第$j$个变量的下界与上界.

步骤 2  个体评价

计算每个个体${P_i}(t)$的适应值, 其中评价函数$f$就是优化的目标函数$f({p_1}, {p_2}, {p_3})$.

步骤 3  变异

对于每一个目标个体向量${P_i}(t), i = 1, 2, \cdots, N$, 按照下面的方式产生变异个体${V_i}(t)$, 其中${V_i}(t) = ({v_{i1}}, {v_{i2}}, {v_{i3}})$, 具体方式如下:

$ {v_{ij}}(t) = {p_{{r_1}j}}(t) + xf \cdot ({p_{{r_2}j}}(t) - {p_{{r_3}j}}(t)), $

(6)

其中${r_1}, {r_2}, {r_3} \in \left\{ {1, 2, \cdots, N} \right\}$并且$i \ne {r_1} \ne {r_2} \ne {r_3}$, 所以种群的大小$N$是大于或者等于4的整数.参数$xf$是指控制差分变量$({p_{{r_2}j}}(t) - {p_{{r_3}j}}(t))$的缩放因子, 它通常在[0, 2]上取值.

步骤 4  杂交

为了增加群体多样性, 微分进化算法引入了杂交操作, 最终来产生新的个体${P'_i}(t) = \left( {{{p'}_{i1}}(t), {{p'}_{i2}}(t), {{p'}_{i3}}(t)} \right)$; 具体方法如下所示:

$ {p'_{ij}}(t) = \left\{ \begin{array}{l} {v_{ij}}(t), ~~~~~{\hbox {如果}}~~~{\rm rand}\left[{0, 1} \right] < {P_c}~~~ {\hbox{ 或者}}~~~ j = {j_{\rm rand}}, \\ {p_{ij}}(t), ~~~~~{\hbox{ 其他}}, \end{array} \right. $

(7)

其中$i = 1, 2, \cdots, N$且$j = 1, 2, 3 $; ${P_c}$是交叉概率, ${j_{\rm rand}}$是[1, 3]之间的一个随机整数, 它保证了新个体${P'_i}(t)$至少有一个分量来自变异个体${V_i}(t)$, 从而避免与父个体${P_i}(t)$相同.微分进化算法常用的杂交方式有:指数杂交、二进制杂交、二项式杂交.

步骤 5  选择

微分进化算法与其它智能优化算法一样, 也采用的是达尔文的适者生存定律, 即是适应值好的个体进入下一代, 适应值差的个体被淘汰.微分进化算法采用的是贪婪选择策略.选择过程如下:

$ {P_i}(t + 1) = \left\{ {\begin{array}{l} {{{P'}_i}(t), ~~~~~{\hbox{如果}}~~~f({{P'}_i}(t)) \le f({P_i}(t))}, \\ {{P_i}(t), ~~~~~~~{\hbox{其他}}}, \end{array}} \right. $

(8)

其中$f({P_i}(t))$为个体${P_i}(t)$的适应值, 由于微分进化算法采用一对一选择操作, 与排序或竞标赛选择相比更能维持群体的多样性, 接下来令$t = t + 1$.

步骤 6  终止检验

如果演化代数$t ={\rm sub}\_{\rm gen}$, 则满足终止条件并输出种群${Z}(t)$中具有最小适应值的个体为最优解.否则转到步骤2.

灵敏度分析^{[9, 11]}是研究与分析一个系统(或模型)的状态或输出变化对系统参数或周围条件变化的敏感程度的方法.在最优化方法中经常利用灵敏度分析来研究原始数据不准确或发生变化时最优解的稳定性.通过灵敏度分析还可以决定哪些参数对系统或模型有较大的影响.因此, 灵敏度分析几乎在对各种方案进行评价时都是很重要的^[12].

考虑如下Mackey-Glass混沌时滞微分方程

$ \frac{{dy(t)}}{{dt}} = - {p_1}y(t) + \frac{{{p_2}y(t - 17)}}{{1 + {y^{10}}(t - 17)}}, y({t_0}) = {p_3}, y \in R, P \in {R^3}, $

(9)

也可写成如下的形式

$ y'(t) = f(t, y, P), y({t_0}) = {p_3}, y \in R, P \in {R^3}, $

(9')

其中$P = ({p_1}, {p_2}, {p_3})$为参数, 变量$y$对${p_i}$的绝对灵敏度定义为

$ {S_j}(t) = \frac{{\partial y(t, p)}}{{\partial {p_j}}}, j = 1, 2, 3, $

(10)

为了消除量纲的影响, 也可定义相对灵敏度.其表达式为

$ {S_j}(t) = \frac{{\partial y(t, p)}}{{\partial {p_j}}}\frac{{{p_j}}}{y}, j = 1, 2, 3. $

(11)

本文对这两种定义都给出了计算.为了得到${S_j}(t)$满足的关系式, 方程$(9)$两边对${p_j}$求导

$ \frac{{\partial y'(t, y, p)}}{{\partial {p_j}}} = \frac{d}{{dt}}(\frac{{\partial y(t, p)}}{{\partial {p_j}}}) = \frac{d}{{dt}}{S_j}(t) = \frac{{\partial f}}{{\partial y}} \bullet \frac{{\partial y}}{{\partial {p_j}}} + \frac{{\partial f}}{{\partial {p_j}}}, $

其中$j = 1, 2, 3$.上式中$f$是已知的, $\frac{{\partial f}}{{\partial y}}$和$\frac{{\partial f}}{{\partial {p_j}}}$可通过复合函数求偏导获得.与方程$(9)'$联合, 可得到求灵敏度${S_j}(t)$的微分方程组, 写成矩阵如下所示

$ Y'(t) = F(t, Y, P), Y({t_0}) = {Y_0}, Y \in {R^4}, P \in {R^3}, $

(12)

其中

$ \begin{eqnarray*} Y(t) = \left( {\begin{array}{*{20}{c}} {y(t)}\\ {{S_1}(t)}\\ {{S_2}(t)}\\ {{S_3}(t)} \end{array}} \right), ~~~~F(t, Y, P) = \left( {\begin{array}{*{20}{c}} {f(t, y, p)}\\ {\frac{{\partial f}}{{\partial y}}{S_1}(t) + \frac{{\partial f}}{{\partial {p_1}}}}\\ {\frac{{\partial f}}{{\partial y}}{S_2}(t) + \frac{{\partial f}}{{\partial {p_2}}}}\\ {\frac{{\partial f}}{{\partial y}}{S_3}(t) + \frac{{\partial f}}{{\partial {p_3}}}} \end{array}} \right), \\ ({S_i}(t) = \frac{{\partial y(t, p)}}{{\partial {p_i}}}, ~~~~{S'_i}(t) = \frac{{\partial f}}{{\partial y}}{S_i}(t) + \frac{{\partial f}}{{\partial {p_i}}}, i = 1, 2, 3, \\ Y({t_0}) = {(y({t_0}), {S_1}({t_0}), {S_2}({t_0}), {S_3}({t_0}))^{\rm T}}). \end{eqnarray*} $

由于方程$(9)'$的初始条件为$y({t_0}) = {p_3}$, 只与${p_3}$有关, 与参数${p_1}, {p_2}$无关, 便可得${S_1}({t_0}) = {S_2}({t_0}) = 0, {S_3}({t_0}) = 1$. (12) 式是一个微分方程组, 可通过龙格-库塔法求出其数值解, 从而可得到变量$y$对参数${p_j}$在任意时刻$t$处的绝对灵敏度${S_j}(t)$, 其相对灵敏度也可用类似的方法求得.

以上述的Mackey-Glass混沌时滞微分方程为例进行数值模拟.为了反演参数${p_1}, {p_2}, {p_3}$需要附加先验条件, 模拟时先给${p_1}, {p_2}, {p_3}$赋真值, 取${p_1} = 0.1$, ${p_2} = 0.2, {p_3} = 1.2$, 可通过龙格-库塔法求解出$y$在任意时刻$t$处的数值解, 并把$y$的值作为已知条件来反求参数.利用微分进化算法反演参数时, 控制参数为交叉因子$xf$=0.9, 交叉概率${P_c}$=0.3, 种群规模N=150.其反演结果见表4-1.

从表4-1中可以看出, 经过1000次迭代的反演值是${p_1} = 0.0999999999, $ ${p_2} = 0.2000000000, $ ${p_3}= 1.2000000001, $与真值${p_1} = 0.1, {p_2} = 0.2, {p_3} = 1.2$相比精度很高, 误差只有3.115144e-019, 而经过2000次迭代的反演值几乎与真值相差无几, 误差只有4.972597e-029;另外一方面, 本文算法与其他智能算法相比, 其运行时间也相对较短.这充分说明了本文算法在解决此类Mackey-Glass混沌时滞微分方程的参数反演问题时是可行的且十分有效的.

参数对每个变量的影响程度可以由计算灵敏度得到, 以上述的Mackey-Glass方程为例, 求其绝对灵敏度的方程组如下所示:

$ \left\{ \begin{array}{l} y' = \frac{{{p_2}y(t - 17)}}{{1 + {y^{10}}(t - 17)}} - {p_1}y(t), \\ {{S'}_1} = \frac{{{p_2}{S_1}(t - 17)[1-9{y^{10}}(t-17)]}}{{{{[1 + {y^{10}}(t-17)]}^2}}} - y(t) - {p_1}{S_1}, \\ {{S'}_2}= \frac{{y(t - 17) + {y^{11}}(t - 17) + {p_2}{S_2}(t - 17) - 9{p_2}{S_2}(t - 17){y^{10}}(t - 17)}}{{{{[1 + {y^{10}}(t-17)]}^2}}} - {p_1}{S_2}, \\ {{S'}_3} = \frac{{{p_2}{S_3}(t - 17)[1-9{y^{10}}(t- 17)]}}{{{{[1 + {y^{10}}(t-17)]}^2}}} - {p_1}{S_3}, \end{array} \right.~~~~~ \left\{ \begin{array}{l} y(0) = {p_3}, \\ {S_1}(0) = 0, \\ {S_2}(0) = 0, \\ {S_3}(0) = 1, \end{array} \right. $

(13)

其中${S_i} = \frac{{\partial y(t, p)}}{{\partial {p_i}}}, {S'_i} = \frac{{\partial {S_i}}}{{\partial t}}, i = 1, 2, 3$. (13) 式是一个时滞微分方程组, 通过龙格-库塔法求出其数值解, 从而可得到变量$y$对参数${p_j}$在任意时刻\(t\)处的绝对灵敏度${S_j}(t)$, 其相对灵敏度也可用类似的方法求得, 如图 1所示.

由图a和图d可以看出, $y$对${p_1}$的绝对灵敏度及相对灵敏度均小于零, 这说明${p_1}$的增加将会引起$y$的减小.例如, 当$t=0.5$时, 通过计算可知$\frac{{\partial y}}{{\partial {p_1}}}{|_{t = 0.5}} = {\rm{ - 0}}{\rm{.5747}}$, 也就是说$t=0.5$时${p_1}$增加$10\% $, 将会导致$y$减少$5.747\% $.由图b和图e可以看出$y$对${p_2}$的绝对灵敏度及相对灵敏度均大于零, 这说明${p_2}$的增加将会引起$y$的增加.通过计算可知当$t=19$时, $\frac{{\partial y}}{{\partial {p_2}}}{|_{t = 19}} = 1.5314$, 也就是说此时${p_2}$增加$10\% $, 将会导致$y$增加$15.341\% $.关于$y$对${p_3}$的灵敏度也可以类似的分析.

本文应用微分进化算法对Mackey-Glass混沌时滞微分方程的未知参数进行了反演, 同时对此类问题的灵敏度给出了具体的分析, 数值模拟的结果表明了该算法对解决此类时滞微分方程参数反演的问题是可行而且十分有效的.文中选择的数值算例是一个非线性方程和三个参数, 其实对多个非线性方程组成的方程组及多个参数本文方法也是可行的.在微分进化算法中, 交叉因子$xf$, 交叉概率${P_c}$, 种群规模$N$的选择对算法的收敛性有一定的影响, 经多次计算可知, 选取$xf=0.9, {P_c}=0.3$, $N=150$时反演效果最佳.