文章快速检索  
  高级检索
分数阶累加时滞GM(1,N,τ)模型及其应用
毛树华1,2, 高明运1, 肖新平1    
1. 武汉理工大学 理学院, 武汉 430070;
2. 武汉理工大学 可靠性研究中心, 武汉 430070
摘要:基于系统的时滞性, 本文建立了时滞灰色GM(1,N,τ)模型, 给出了模型的最小二乘参数估计公式以及模型的解析解.在引入分数阶累加生成算子后, 将原模型扩展为分数阶累加GM(1,N,τ) 模型, 当时滞值为非整数情况时, 采用相邻整数点加权构造法, 完善了模型; 通过粒子群算法确定模型最优的分数阶累加生成阶数. 最后本文结合武汉市1995-2008年14 年科技投入及经济增长的实际背景, 分别建立了经典时滞GM(1,N,τ)和分数阶累加时滞GM(1,N,τ) 模型对GDP数据做了预测, 比较了两个模型预测结果, 发现分数阶累加时滞GM(1,N,τ) 模型具有更高的建模精度.
关键词灰色模型     GM(1,N,τ)模型     分数阶累加生成     时滞     粒子群算法    
Fractional order accumulation time-lag GM(1,N,τ) model and its application
MAO Shu-hua1,2, GAO Ming-yun1, XIAO Xin-ping1     
1. College of Science, Wuhan University of Technology, Wuhan 430070, China;
2. Reliability Engineering Institute, Wuhan University of Technology, Wuhan 430070, China
Abstract:Based on the time lag effect of system, this paper constructs the time-lag GM(1,N,τ) model, and provides its least squares parameter estimation formula and analytical solution. Then introducing fractional order accumulation generation operator, GM(1,N,τ) is transformed into fractional order accumulation time-lag GM(1,N,τ), offering adjacent integer points weighted method when time-lag value is not integer, determining the time-lag value via particle swarm optimization. Finally, we test the original GM(1,N,τ) and fractional order accumulation time-lag GM(1,N,τ) by using the real data of R&D investment and GDP in Wuhan from 1995 to 2008. The example indicates that the accuracy of fractional order accumulation time-lag GM(1,N,τ) is more satisfactory.
Key words: grey model     GM(1,N,τ) model     fractional order accumulation     time-lag     particle swarm optimization    

0 引言

自邓聚龙教授提出灰色系统理论以来, 灰色预测模型已广泛应用于工业、农业、水利、地质、科教、军事等众多领域[1]. 灰色模型是灰色系统理论的核心和重要组成部分, 它是在时间序列的基础上,建立近似微分方程模型,简称灰建模. 灰建模是少数据的建模,其目的是在数据有限的条件下, 模仿微分方程建立具有部分微分方程性质的模型,即灰色模型, 其中GM(1,1) 模型是最基础的也是应用最广泛的模型, 它主要是针对单变量系统,建立关于时间(或广义时间)的预测模型. 对于多变量系统,类比GM(1,1)模型,学者们提出了GM$(1,N)$模型. GM$(1,N)$ 模型中包括一个行为变量,及$N-1$个因子变量. GM$(1,N)$为分析模型、因子模型,它不具有全信息. 然而, 当有必要对多因子的系统作整体的、全局的、动态的分析时, 就需要使用GM$(1,N)$ 模型[1, 2].

GM$(1,N)$模型,可以看成是GM(1,1)模型的推广形式, 文献[3]假定因子项为常数或者变化很小, 进而退化为一个特殊的GM(1,1)模型, 再运用GM(1,1)模型求解公式进行拟合和预测. 文献[4]基于多因子灰色模型的精确级差分析,建立了MGM$(1,N)$模型; 文献[5]在此基础上,将误差融入级差格式, 基于理想状态时的相对误差提出了MGMp$(1,N)$模型; 文献[6] 结合自记忆方程,建立了GM(1,2)自记忆模型; 文献[7]在白化GM$(1,N)$方程中考虑因素间线性关系派生出了GM$(1,N,x_{i}^{(0)})$和GM$(1,N,x_{i}^{(1)})$模型.

文献[8]首次考虑将时滞因子加入到GM(1,2)模型中, 并通过比较模型的精度来确定时滞值, 说明了带有时滞的GM(1,2)模型能更好地反映两变量间的动态关系. 文献[9]在研究娘子关泉水流量与降水量关系时, 通过时滞GM(1,2)模型对娘子关泉水流量进行了预测,结果与实际相接近. 这些模型对于变量间存在整数时滞值,具有很好的预测能力. 文献[10]采用灰色关联度的方法计算模型的时滞值, 进一步完善带时滞的GM$(1,N)$模型.

20世纪70年代以来, 研究人员发现分数阶微积分可以作为一种很好的描述与刻画手段, 应用于分形几何、幂律现象与记忆过程等相关现象或过程. 文献[11] 将分数阶累加生成应用到灰色模型中, 建立了基于分数阶累加生成技术的灰色系统模型.

上述学者的研究成果极大地推动了灰色模型的发展, 然而随着灰色预测模型应用范围不断扩大,各式各样的新问题层出不穷, 这就需要对GM$(1,N)$ 模型做相应改进. 本文基于现有的研究基础. 在第1节提供了时滞GM$(1,N,\tau)$模型的最小二乘参数估计方法以及模型的解. 第2节提供了模型的评价标准. 在第3节,首先介绍了分数阶累加生成矩阵; 随后对该模型的累加生成方式做了推广, 将分数阶累加生成引入时滞GM$(1,N,\tau)$模型, 建立了分数阶GM$(1,N,\tau)$ 模型; 然后对非整数时滞值情况下, 对模型进行了完善,最后基于粒子群算法给出了确定模型阶数的方法. 在第4节,结合实际数据, 分别建立了传统灰色模型和分数阶累加生成灰数模型, 并对预测结果进行分析比较.\vspace{-1.5mm}

1 时滞GM$(1,N,\tau)$模型

GM$(1,N,\tau)$是多变量时滞灰色模型的符号. GM$(1,N,\tau)$ 中包括一个行为变量,记为$y$; $N-1$个因子变量,记为${{x}_{i}}$, $i=1,2,\cdots ,N-1$; 以及因子变量与行为变量的时间差, 即时滞值$\tau$.

定义 1 设数据序列$\displaystyle{{y}^{\left( 0 \right)}}=\{{{y}^{(0)}}(1),\ {{y}^{(0)}}(2),\ \cdots ,{{y}^{(0)}}(n)\}$ 为系统行为序列,序列$\displaystyle x_{i}^{(0)}=\{x_{i}^{(0)}(1),\ x_{i}^{(0)}(2), \cdots ,x_{i}^{(0)}(n)\}$为相关因素序列,$i=1,2,\cdots ,N-1$. ${{y}^{(1)}}$ 为${{y}^{(0)}}$ 的一次累加生成序列, 即${{y}^{(1)}}(k)=\sum_{j=1}^{k}{{{y}^{(0)}}(j)}$; $x_{i}^{(1)}$ 为$x_{i}^{(0)}$的一次累加生成序列, 即$x_{i}^{(1)}(k)=\sum_{j=1}^{k}{x_{i}^{(0)}(j)}$; ${{z}_{y}}^{(1)}$ 为${{y}^{(1)}}$ 的紧邻均值生成序列, 即${{z}_{y}}^{(1)}(k)=0.5\left( {{y}^{(1)}}(k)+{{y}^{(1)}}(k-1) \right)$,则称

\begin{equation} {{y}^{(0)}}(k)+a{{z}_{y}}^{(1)}(k)=\sum\limits_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(1)}\left( k-\tau \right)} \end{equation} (1)
为GM$(1,N,\tau)$模型.

定理 1 设${{y}^{(0)}}$为系统行为序列, $x_{i}^{(0)}$$(i=1,2,\cdots ,N-1)$为系统因子数据序列, ${{y}^{(1)}}$,$x_{i}^{(1)}$分别为${{y}^{(0)}}$, $x_{i}^{(0)}$的一次累加生成序列, $z_{y}^{(1)}$为${{y}^{(1)}}$的紧邻均值生成序列,

$${B}=\left[\begin{matrix} {{z}_{y}}^{(1)}(\tau +2) & x_{1}^{(1)}(2) & \cdots & x_{N-1}^{(1)}(2) \\ {{z}_{y}}^{(1)}(\tau +3) & x_{1}^{(1)}(3) & \cdots & x_{N-1}^{(1)}(3) \\ \vdots & \vdots & \ddots & \vdots \\ {{z}_{y}}^{(1)}(n) & x_{1}^{(1)}(n-\tau ) & \cdots & x_{N-1}^{(1)}(n-\tau ) \\ \end{matrix} \right],\qquad {Y}=\left[\begin{matrix} {{y}^{(0)}}(\tau +2) \\ {{y}^{(0)}}(\tau +3) \\ \vdots \\ {{y}^{(0)}}(n) \\ \end{matrix} \right],$$
则参数列$\hat{{P}}=\left[a,b_{1},b_{2},\cdots ,{b}_{N-1} \right]^{\rm T}$的最小二乘估计 满足
\begin{equation} \hat{{P}}={{\left( {{{B}}^{\rm T}}{B} \right)}^{-1}}{{{B}}^{\rm T}}{Y} \end{equation} (2)

{\HT 证明} 将数据代入模型(1)中得到

$$\begin{matrix} \displaystyle{y^{\left( 0 \right)}}\left( \tau +2 \right)=\sum\limits_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(1)}\left( 2 \right)}- a{{z}_{y}}^{\left( 1 \right)}\left( \tau +2 \right),\\ \displaystyle{y^{\left( 0 \right)}}\left( \tau +3 \right)=\sum\limits_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(1)} \left( 3 \right)}-a{{z}_{y}}^{\left( 1 \right)}\left( \tau +3 \right) ,\\ \vdots \\ \displaystyle{{y}^{\left( 0 \right)}}\left( n \right)=\sum\limits_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(1)} \left( n-\tau \right)}-a{{z}_{y}}^{\left( 1 \right)}\left( n \right). \\ \end{matrix}$$

上述方程组写成矩阵形式为${Y}\!=\!{B}\hat{{P}}$, 对于$\hat{{P}}={{[a,{{b}_{1}},{{b}_{2}},\cdots ,{{b}_{N-1}} ]}^{\rm T}}$ 的一组参数估计列,以 $\sum_{i=1}^{N-1}{{b}_{i}}x_{i}^{(1)}( k\\-\tau)$$-a{{z}_{y}}^{\left( 1 \right)}\left( k \right)$代替${{y}^{\left( 0 \right)}}\left( k \right)\left( k=\tau +2,\tau +3,\cdots ,n \right)$可得误差序列$\boldsymbol{\varepsilon} ={Y}-{B}\hat{{P}}. $

设$\displaystyle{s={{\boldsymbol{\varepsilon }}^{\rm T}}\boldsymbol{\varepsilon} ={{( {Y}-{B}\hat{{P}} )}^{\rm T}}( {Y}-{B}\hat{{P}})=\sum\limits_{k=\tau +2}^{n}{{{\left( {{y}^{(0)}}\left( k \right)-\sum\limits_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(1)}\left( k-\tau \right)}+a{{z}_{y}}^{(1)}\left( k \right) \right)}^{2}}}}$.

使$s$最小的$a,{{b}_{1}},{{b}_{2}},\cdots , {{b}_{N-1}}$ 应满足

\[\left\{ \begin{array}{l} \displaystyle\frac{{\partial s}}{{\partial a}} = 2\sum\limits_{k = \tau + 2}^n {\left( {{y^{(0)}}(k) - \sum\limits_{i = 1}^{N - 1} {{b_i} x_i^{(1)}(k - \tau )} + az_y^{(1)}(k)} \right)z_y^{(1)}(k) = 0}, \\ \displaystyle\frac{{\partial s}}{{\partial {b_1}}} = - 2\sum\limits_{k = \tau + 2}^n {\left( {{y^{(0)}}(k) - \sum\limits_{i = 1}^{N - 1} {{b_i}x_i^{(1)} (k - \tau )} + az_y^{(1)}(k)} \right)x_1^{(1)}(k - \tau ) = 0}, \\ \displaystyle\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \vdots \\ \displaystyle\frac{{\partial s}}{{\partial {b_{N - 1}}}} = - 2\sum\limits_{k = \tau + 2}^n {\left( {{y^{(0)}}(k) - \sum\limits_{i = 1}^{N - 1} {{b_i}x_i^{(1)}(k - \tau )} + az_y^{(1)}(k)} \right)x_{N - 1}^{(1)}(k - \tau ) = 0}. \end{array} \right.\]

$\sum\limits_{k = \tau + 2}^n {\left( {{y^{(0)}}(k) - \sum\limits_{i = 1}^{N - 1} {{b_i}x_i^{(1)}(k - \tau )} + az_y^{(1)}(k)} \right)} z_y^{(1)}(k) = 0,\\ \displaystyle\sum\limits_{k = \tau + 2}^n {\left( {{y^{(0)}}(k) - \sum\limits_{i = 1}^{N - 1} {{b_i}x_i^{(1)}(k - \tau )} + az_y^{(1)}(k)} \right)} x_1^{(1)}(k - \tau ) = 0,\\ \begin{array}{*{20}{c}} {}&{}&{}&{\begin{array}{*{20}{c}} {}&{}&{}& \qquad \vdots \end{array}} \end{array}\\ \displaystyle\sum\limits_{k = \tau + 2}^n {\left( {{y^{(0)}}(k) - \sum\limits_{i = 1}^{N - 1} {{b_i}x_i^{(1)}(k - \tau )} + az_y^{(1)}(k)} \right)} x_{N - 1}^{(1)}(k - \tau ) = 0.$

即 ${B}^{\text{T}}\boldsymbol{\varepsilon }=\mathbf{0}$, 于是${B}^{\text{T}}({Y}-{B}\hat{{P}} )=\mathbf{0}$, 或${{B}^{\text{T}}}{Y}-{{B}^{\text{T}}}{B}\hat{{P}}=\mathbf{0}$, 即${{B}^{\text{T}}}{B}\hat{{P}}={{B}^{\text{T}}}{Y}$.

所以

$$\hat{{P}}={{\left( {{B}^{\text{T}}}\mathbf{B} \right)}^{-1}}{{B}^{\text{T}}}{Y}.$$

由(1)可得到GM$(1,N,\tau)$模型的白化形式如下

\begin{equation} \frac{{d}{{y}^{(1)}}}{{ d}t}+a{{y}^{(1)}}\left(t\right)=\sum\limits_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(1)}\left( t-\tau \right)} \end{equation} (3)
定理 2如(3)所示微分方程的解为
\begin{equation} {{y}^{(1)}}\left( t+1 \right)=\left( {{y}^{(1)}}(\tau +1)-\sum\limits_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(1)}(t+1)} \right){{\rm e}^{-a(t-\tau )}}+\frac{1}{a}\sum\limits_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(1)}(t+1)},\quad t\ge \tau +1 \end{equation} (4)
证明
\begin{eqnarray*} & \displaystyle\frac{{\rm d}{y^{(1)}}}{{\rm d}t}+a{y^{(1)}}\left(t\right)=\sum\limits_{i=1}^{N-1}{{b_{i}}x_{i}^{(1)}\left(t-\tau \right)},\\ & \displaystyle\frac{{\rm d}{{y}^{(1)}}}{a{{y}^{(1)}}+\sum_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(1)}\left( t-\tau \right)}}={\rm d}t,\quad \text{即} \int{\frac{{\rm d}{{y}^{(1)}}}{a{{y}^{(1)}}+\sum_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(1)}\left( t-\tau \right)}}}=\int{{\rm d}t}. \end{eqnarray*}

因此$\ln ( a{{y}^{(1)}}+\frac{1}{a}\sum_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(1)}( t-\tau )} )=a( t+C )$,其中$C$ 为未知常数.

于是

$$\displaystyle a{{y}^{(1)}}+\frac{1}{a}\sum\limits_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(1)}\left( t-\tau \right)}={{C}_{1}}\exp \left( at \right),$$

其中$\displaystyle{{C}_{1}}=\exp \left( aC \right)$.

所以(4)式得证,故模型预测值为

\[\displaystyle{{\hat{y}}^{(1)}}\left( t+1 \right)=\left( {{y}^{(1)}}(\tau +1)-\sum\limits_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(1)}(t+1)} \right){{\rm e}^{-a(t-\tau )}}+\frac{1}{a}\sum\limits_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(1)}(t+1)}.\]

模型还原值

\begin{equation} {{\hat{y}}^{\left( 0 \right)}}\left( t+1 \right)={{\hat{y}}^{\left( 1 \right)}}\left( t+1 \right)-{{\hat{y}}^{\left( 1 \right)}}\left( t \right),\quad t\ge \tau +1 \end{equation} (5)

2 模型的误差计算

模型的误差是评价一个模型的可靠度与实用度的重要指标. 本文中, 主要以平均绝对误差百分比(MAPE,即mean absolute percentage error) 来表示模型的误差,计算公式如下

\begin{equation} \text{MAPE}=\frac{1}{n}\sum\limits_{k=1}^{n}{( {| {{y}^{\left( 0 \right)}}\left( k \right)-{{{\hat{y}}}^{\left( 0 \right)}}\left( k \right)|}/{{{y}^{\left( 0 \right)}}\left( k \right)} )\times 100\%} \end{equation} (6)

显然模型的MAPE值越大,模型的可靠度和实用价值越低.

3 分数阶累加GM$(1,N,\tau)$模型 3.1 分数阶累加生成矩阵

对于序列$x^{(0)}$的一次累加生成,可以写为${{x}^{(1)}}\left( k \right)=\sum_{i=1}^{k}{{{x}^{(0)}}\left( i \right)}$, 采用矩阵表述可以写为${{x}^{(1)}}={A}{{x}^{(0)}}$, 其中${A}=\begin{smallmatrix} 1 & 0 & \cdots & 0 \\ 1 & 1 & \cdots & 0 \\ \scriptstyle \vdots & \vdots & \ddots & \vdots \\ 1 & 1 & \cdots & 1 \\ \end{smallmatrix} $,称${A}$为一次累加生成矩阵. 如果对原始序列进行$r(r\ge 2)$ 次累加生成, 那么可以写为 ${{x}^{(r)}}={{{A}}^{r}}{{x}^{(0)}}$, 称${{{A}}^{r}}$ 为$r$ 阶累加生成矩阵,当$r$为整数时, 易得\({{{A}}^{r}}={{(a_{ij}^{r})}_{n\times n}}\).

其中$\displaystyle{{(a_{ij}^{r})}_{n\times n}=\left\{\begin{array}{c@{\quad \quad}l}C_{i-j+r-1}^{r-1}=\frac{(i-j+r-1)!}{(r-1)!(i-j)!} & i\ge j\\ 0 & i<j \end{array}\right.} $,

这样可以得到

$$\displaystyle{{{A}}^{r}}=\left[\begin{matrix} 1 & 0 & 0 & \cdots & 0 \\ r & 1 & 0 & \cdots & 0 \\ \frac{r(r+1)}{2!} & r & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \frac{r(r+1)\cdots (r+n-2)}{(n-1)!} & \frac{r(r+1)\cdots (r+n-3)}{(n-2)!} & \frac{r(r+1)\cdots (r+n-4)}{(n-3)!} & \cdots & 1 \\ \end{matrix} \right].$$

根据这个公式,结合组合数的广义定义[12],可对${{{A}}^{r}}$ 进行推广,称为分数阶累加生成矩阵.

{\HT 注} 当$r$为正整数时,${{{A}}^{r}}$等同于矩阵运算, 是矩阵${A}$的$r$次方; 当$r$不是正整数时, ${{{A}}^{r}}$只是一个记号. 例如,当$n=4$时, 阶数为$\frac{1}{2},\frac{5}{4}$ 时,累加生成矩阵如下

$${{{A}}^{\frac{1}{2}}}=\left[\begin{matrix} 1 & 0 & 0 & 0 \\ 0.5 & 1 & 0 & 0\\ 0.375 & 0.5 & 1 & 0 \\ 0.3125 & 0.375 & 0.5 & 1 \end{matrix} \right],\qquad {{{A}}^{\frac{5}{4}}}=\left[\begin{matrix} 1 & 0 & 0 & 0\\ 1.25 & 1 & 0 &0 \\ 1.4063 & 1.25 & 1 & 0\\ 1.5234 & 1.4063 & 1.25 &1 \end{matrix} \right].$$

用$-r$代替上式中$r$,可以得到 对应的分数阶累减生成矩阵,形式如下

结合矩阵运算公式,易证${{{A}}^{r}}{{{A}}^{-r}}={E}$, 其中${E}$ 为$n$阶单位阵. 于是,上述两例的逆矩阵如下

$${{{A}}^{-\frac{1}{2}}}=\left[ \begin{matrix} 1 & 0 & 0 & 0 \\ -0.5 & 1 & 0 & 0\\ -0.125 & -0.5 & 1 & 0 \\ -0.0625 & -0.125 & -0.5 & 1 \end{matrix} \right],\qquad {{{A}}^{-\frac{5}{4}}}=\left[ \begin{matrix} 1 & 0 & 0 & 0\\ -1.25 & 1 & 0 &0 \\ 0.1563 & -1.25 & 1 & 0\\ 0.0391 & 0.1563 & -1.25 &1 \end{matrix} \right].$$

3.2 分数阶累加GM$(1,N,\tau)$模型

当我们采用分数阶累加生成矩阵${{{A}}^{r}}$代替原模型中的一次累加生成, 我们称所得到的模型为分数 阶累加GM$(1,N,\tau)$,此时模型白化方程

\begin{equation} \frac{{\rm d}{{y}^{(r)}}}{{\rm d}t}+a{{y}^{(r)}}= \sum\limits_{i=1}^{N-1}{{{b}_{i}}}x_{i}^{(r)}\left( t-\tau \right) \end{equation} (7)

灰色方程

\begin{equation} {{y}^{\left( r-1 \right)}}\left( k \right)+a{{z}_{y}}^{\left( r \right)}\left( k \right)=\sum\limits_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(r)}\left( k-\tau \right)},\quad k=\tau +2,\tau +3,\cdots ,n \end{equation} (8)

参数列$\hat{{P}}={{\left[a,{{b}_{1}},{{b}_{2}},\cdots ,{{b}_{N-1}} \right]}^{\text{T}}}$的最小二乘估计满足

\begin{equation} \hat{{P}}={{( {{B}^{*}}^{\rm T}{{B}^{*}} )}^{\mathbf{-1}}}{{B}^{*}}^{\rm T}{{{Y}}^{*}} \end{equation} (9)

$\text{其中}\,{{B}^{*}}=\left[\begin{matrix} {{z}_{y}}^{(r)}\left( \tau +2 \right) & x_{1}^{(r)}\left( 2 \right) & \cdots & x_{N-1}^{(r)}\left( 2 \right) \\ {{z}_{y}}^{(r)}\left( \tau +3 \right) & x_{1}^{(r)}\left( 3 \right) & \cdots & x_{N-1}^{(r)}\left( 3 \right) \\ \vdots & \vdots & \ddots & \vdots \\ {{z}_{y}}^{(r)}\left( n \right) & x_{1}^{(r)}\left( n-\tau \right) & \ldots & x_{N-1}^{(r)}\left( n-\tau \right) \\ \end{matrix} \right]$,${{{Y}}^{*}}=\left[\begin{matrix} {{y}^{(r)}}(\tau +2)-{{y}^{(r)}}(\tau +1) \\ {{y}^{(r)}}(\tau +3)-{{y}^{(r)}}(\tau +2) \\ \vdots \\ {{y}^{(r)}}(n)-{{y}^{(r)}}(n-1) \\ \end{matrix} \right].$

类似(5),模型的$r$次累加生成序列的预测值为

\begin{equation} {{\hat{y}}^{(r)}}\left( t+1 \right)=\left( {{y}^{(r)}}(\tau +1)-\sum\limits_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(r)}(t+1)} \right){{\rm e}^{-a(t-\tau )}}+\frac{1}{a}\sum\limits_{i=1}^{N-1}{{{b}_{i}}x_{i}^{(r)}(t+1)} \end{equation} (10)

同时我们可以得到序列的还原值

\begin{equation} {{\hat{y}}^{(0)}}\text{=}{{{A}}^{-r}}{{\hat{y}}^{(r)}} \end{equation} (11)

3.3 非整数时滞值下模型完善

上述模型主要是针对整数时滞值情景. 由于现实中的时滞值往往是非整数的,此时${{x}^{(r)}}(t-\tau )$数据不存在. 为解决这一种问题,本文提出一种相邻整数点加权构造法进行完善.

假设$x=\{x(1),x(2),\cdots ,x(n)\}$为一个原始数据序列, $k\text{=}1,2,\cdots ,n$为序列序号,显然针对该序列, 序号和序列数据是一一映射关系,即针对每一个$k$值, 有且仅有一个数据$x(k)$与之对应. 设实数$p\in \left[1,\ n \right]$, 函数$q=\left[p \right]$ 表示向下取整运算. 那么序列中, 序号$p$对应数据$x(p)$的计算定义如下

\begin{equation} x\left( p \right)=\left( q+1-p \right)x\left( q \right)+(p-q)x\left( q+1 \right) \end{equation} (12)

若当时滞值为1.8时,则数据$x(1.8)$的线性加权构造公式为

$$\displaystyle x(1.8)=(2-1.8)x(1)+(1.8-1)x(2)=0.2x(1)+0.8x(2).$$

结合(11),可将预测模型(10)改写为(13)形式.

\begin{equation} \begin{array}{r@{~}l} {{\hat{y}}^{(r)}}\left( t+1 \right)=& \left( {{y}^{(r)}}(\tau +1)-\sum\limits_{i=1}^{N-1}{{{b}_{i}}\left( {{k}_{1}}x_{i}^{(r)}\left( t-\left[ \tau \right] \right)+{{k}_{2}}x_{i}^{(r)}\left( t-\left[\tau \right]+1 \right) \right)} \right){{\rm e}^{-a(t-\tau )}}+ \\ & \displaystyle{+\frac{1}{a}\sum\limits_{i=1}^{N-1}{{{b}_{i}}\left( {{k}_{1}}x_{i}^{(r)}\left( t-\left[\tau \right] \right)+{{k}_{2}}x_{i}^{(r)}\left( t-\left[\tau \right]+1 \right) \right)}} \end{array} \end{equation} (13)
其中${{k}_{1}}=\tau -\left[\tau \right]$,${{k}_{2}}=\left[\tau \right]+1-\tau $,$t\ge \left[\tau \right]+1$.

3.4 模型阶数的确定

分数阶累加GM$(1,N,\tau)$模型实际上是GM$(1,N,\tau)$模型的一种拓展和补充, 当$r=1$时,该模型退化为\linebreak GM$(1,N,\tau)$模型,当$\tau=0$ 时,模型即为GM$(1,N)$ 模型. 对于模型阶数的确定, 本文选取使模型MAPE值最小的$r$,作为建模阶数. 这里以MAPE值最小为优化目标,建立关于$r$的优化模型,表述如下

$$\begin{array}{r@{~}l} & \min \ \operatorname{MAPE}(r) \\ & {\rm s.t.}\ \quad r>0 \\ \end{array}$$

由于MAPE的计算过程中,涉及到绝对值符号运算,MAPE关于$r$ 是不可导的,难以给出模型的解析式. 所以采用群体智能优化算法, 如遗传算法,粒子群算法等搜索出模型最优解.

事实上, 关于分数阶累加GM$(1,N,\tau)$模型其实可以理解为,将原始序列, 如${{x}^{(0)}},{{y}^{(0)}}$,通过广义累加生成矩阵, 如${{{A}}^{r-1}}$ 作用,得到序列${{x}^{(r-1)}}$, ${{y}^{(r-1)}}$,然后对作用后的序列, 进行传统的GM$(1,N,\tau)$模型建模. 另一方面, 分数阶累加是对传统的灰色建模的推广: 对累加生成方式的推广, 由原来的一阶累加推广到分数阶累加.

4 实例分析

本文收集武汉市1995-2008年14年间的科技投入(R\&D)和经济发展(GDP)数据[13, 14]表 1. 为了避免受通货膨胀因素影响,GDP和R\&D 资本额均以1978年不变价格进行调整,并分别进行了对数化运算, 得到LNPGDP和LNR\&D.

表 1 武汉市1995-2008年14年间的科技投入和经济发展数据

4.1 时滞值的估计

对于科技投入与经济增长之间时滞值的计算; 以LNPR\&D序列数据为科技投入数据,以LNPGDP的数据为经济增长数据, 利用文献[15]中所述时滞关联分析,可以得到时滞系数进行下面分析.

表 2 1996-2008年经济增长与科技投入灰关联度表

鉴于科技投入效益"长效性"的特点, 设定从科技投入周期开始四年后的灰关联度不作为判断对象, 其数值可能是投入的长期效益,也可能是以后时段的投入影响, 所以在上表中不列出. 计算取各极值点的年限数, 并以其算术平均值为时滞值, 得出近年来武汉市经济增长与科技投入的时滞值为\(\tau =\frac{1}{8}( 2+1+2+2+0+2\\+2+1 )=1.5\). 由此我们可以得到, 武汉市1996-2008年科技投入与经济发展的平均滞后期为1.5年,表明: 科技投入大约在1.5年之后开始产生收益.

4.2 $r$阶GM$(1,N,\tau)$结果比较

基于模型(13),我们利用粒子群算法,以MAPE值最小为目标, 搜索得到模型最优阶数为$r=0.039659$,模型参数为$a=-0.53572, {{b}_{1}}=0.24256$,根据(13)得到的预测模型为

\begin{eqnarray} {{\hat{y}}^{(r)}}(k\!+\!1)=(0.4308\!-\!0.2264( {{x}^{(r)}}(k\!-\!1 )+{{x}^{(r)}}(k\!-\!2 )))\\ {\operatorname{e}}^{-0.5357(k-2)}+\text{0.2264}( {{x}^{(r)}}(k\!-\!1 )+{{x}^{(r)}}(k\!-\!2)) \end{eqnarray} (14)
\begin{eqnarray} {{\hat{y}}^{(0)}}={{{A}}^{-0.039659}}{{\hat{y}}^{(r)}} \end{eqnarray} (15)

为了便于比较,本文同时计算了时滞值为1年和2年时模型的预测值, 如表 3所示.

表 3 GM$(1,N,\tau)$及$r$阶GM$(1,N,\tau)$在不同时滞下预测结果

4.3 结论及分析

通过表 3可以发现,针对该数据采用经典的GM$(1,N,\tau)$模型完全不可行, 该模型的预测效果极差,甚至预测 出负数, 这种预测显然完全背离了社会实际. 相反, 如果采用分数阶累加GM$(1,N,\tau)$模型,则该模型可以取得 很好的精度, 并且模型的最终预测结果相对满意, 充分说明了分数阶累加GM$(1,N,\tau)$模型在实际应用中具有 更高的适应度. 同时也可以发现,时滞值取$\tau=1.5$时,模型较为合理.

参考文献
[1] 邓聚龙.灰理论基础[M].武汉:华中科技大学出版社, 2008: 282-300.Deng Julong. Grey theory foundation[M]. Wuhan: Huazhong University of Science & Technology Press, 2008: 282-300.
[2] 肖新平,宋中民,李峰.灰技术基础及其应用[M].北京:科学出版社, 2005: 117-123.Xiao Xinping, Song Zhongmin, Li Feng. Grey technology foundation and its application[M]. Beijing: Science Press, 2005: 117-123.
[3] 刘思峰,党耀国,方志耕,等.灰色系统理论及其应用[M]. 5版.北京:科学出版社, 2010: 169-171.Liu Sifeng, Dang Yaoguo, Fang Zhigeng, et al. Grey system theory and its application[M]. 5th ed. Beijing: Science Press, 2010: 169-171.
[4] 翟军,盛建明,冯英浚. MGM(1,N)灰色模型及应用[J].系统工程理论与实践, 1997, 17(5): 109-113.Zhai Jun, Sheng Jianming, Feng Yingjun. The grey model MGM(1,N) and its application[J]. Systems Engineering —— Theory & Practice, 1997, 17(5): 109-113.
[5] 李小霞,同小军,陈绵云. 多因子灰色MGMp(1,N)优化模型[J]. 系统工程理论与实践, 2003, 23(4): 47-51.Li Xiaoxia, Tong Xiaojun, Chen Mianyun. MGMp(1,N) optimization model[J]. Systems Engineering —— Theory & Practice, 2003, 23(4): 47-51.
[6] 张高锋,姜国辉,付永锋,等. 灰色 GM(1,N)自记忆模型研究[J]. 沈阳农业大学学报, 2009, 40(2): 210-214.Zhang Gaofeng, Jiang Guohui, Fu Yongfeng, et al. Establishment of the grey GM(1,N) self-memory model[J]. Journal of Shenyang Agricultural University, 2009, 40(2): 210-214.
[7] 周伟,方志耕.非线性优化GM(1,N)模型及其应用研究[J].系统工程与电子技术, 2010, 32(2): 317-320.Zhou Wei, Fang Zhigeng. Nonlinear optimization method of gray GM(1,N) model and application[J]. Systems Engineering and Electronics, 2010, 32(2): 317-320.
[8] 翟军,冯英俊,盛建明.带有时滞的GM(1,2)模型及应用[J].系统工程, 1996, 14(6): 66-68.Zhai Jun, Feng Yingjun, Sheng Jianming. A time lag existing GM(1,2) model and its application[J]. Systems Engineering, 1996, 14(6): 66-68.
[9] 郝永红,黄登宇,高红波,等. 娘子关泉水流量的GM(1,2)时滞预测模型[J].中国岩溶, 2004(3): 43-47.Hao Yonghong, Huang Dengyu, Gao Hongbo, et al. GM(1,2) time lag model for discharge of Niangziguan Spring[J]. Carsologica Sinica, 2004(3): 43-47.
[10] Hao Y H, Wang Y J, Zhao J J, et al. Grey system model with time lag and application to simulation of karst spring discharge[J]. Grey Systems: Theory and Application, 2011, 1(1): 47-56.
[11] Wu L F, Liu S F, Yao L G, et al. Grey system model with the fractional order accumulation[J]. Communications in Nonlinear Science and Numerical Simulation, 2013, 18(7): 1775-1785.
[12] 屈婉玲.组合数学[M].北京:北京大学出版社, 1989, 11: 48-51.Qu Wanling. Combinatorial mathematics[M]. Beijing: Beijing University Press, 1989, 11: 48-51.
[13] 武汉市统计局.武汉市统计年鉴[R].北京:中国统计出版社, 2010: 3-20, 56-70.Wuhan Statistics Bureau. Wuhan statistical yearbook[R]. Beijing: China Statistics Press, 2010: 3-20, 56-70.
[14] 武汉市科学技术局,武汉市统计局.武汉科技统计年鉴[R].武汉:武汉市统计局, 2008: 29-36.Wuhan Science and Technology Bureau, Wuhan Statistics Bureau. Wuhan statistical yearbook of science and technology[R]. Wuhan: Wuhan Statistics Bureau, 2008: 29-36.
[15] 肖新平, 毛树华. 灰预测与决策方法[M].北京:科学出版社, 2013: 86-88.Xiao Xinping, Mao Shuhua. The method of grey prediction and decision[M]. Beijing: Science Press, 2013: 86-88.
0

文章信息

毛树华, 高明运, 肖新平
MAO Shu-hua, GAO Ming-yun, XIAO Xin-ping
分数阶累加时滞GM(1,N,τ)模型及其应用
Fractional order accumulation time-lag GM(1,N,τ) model and its application
系统工程理论与实践, 2015, 35(2): 430-436
Systems Engineering - Theory & practice, 2015, 35(2): 430-436.

文章历史

收稿日期:2013-07-16

相关文章

工作空间