西南大学学报 (自然科学版)  2019, Vol. 41 Issue (9): 68-76.  DOI: 10.13718/j.cnki.xdzk.2019.09.009
0
Article Options
  • PDF
  • Abstract
  • Figures
  • References
  • 扩展功能
    Email Alert
    RSS
    本文作者相关文章
    曾亮
    欢迎关注西南大学期刊社
     

  • 改进的灰色多变量GM(1, N)模型及其应用    [PDF全文]
    曾亮     
    广东理工学院 基础教学部, 广东 肇庆 526100
    摘要:针对传统灰色多变量GM(1,N)模型的时间响应式不精确和建模精度不高等问题,利用各变量的一阶累加序列具有近似非齐次指数增长率的特性,对传统GM(1,N)模型的背景值进行了重构,提出了一种改进的GM(1,N)模型,并给出了该模型的参数估计式和近似时间响应式.最后通过算例模拟和应用实例表明了改进的GM(1,N)模型的有效性和实用性.
    关键词灰色系统    GM(1, N)模型    非齐次指数增长率    背景值    预测    

    自灰色系统理论[1]提出以来,灰色预测模型已广泛应用于工业、农业、气象、水利、交通、地质、经济和管理等众多领域[2-4].灰色预测模型可划分为单变量和多变量两类模型,其中单变量模型以GM(1,1)模型为代表,它针对的是单时间序列建模,并没有考虑系统相关因素的影响.

    多变量模型以GM(1,N)模型为代表,该模型包含了一个系统特征变量和N-1个影响因子变量,可以分析多个影响因子变量对系统行为的作用.但由于该模型时间响应式不精确,导致了模型精度不高,所以在应用中受限.为此,很多学者对GM(1,N)模型进行改进研究,得到了一些很有价值的成果[4-15].

    以上研究虽然取得了一些成果,但并未真正解决时间响应式不精确的问题.本文在现有研究的基础上,考虑到模型误差主要来源于从连续到离散的转化误差和时间响应函数的近似计算误差,利用一阶累加序列具有近似非齐次指数增长率的特性,对模型进行了两方面的改进:一是对背景值进行了重构,尽量减少转化误差;二是对模型时间响应式进行了近似计算,尽量减少因假定相关因素变量为灰常量而带来的模拟计算误差.数值仿真和应用实例表明,改进后的模型能显著提高建模精度.

    1 传统GM(1,N)模型

    定义1  若记X1(0)=(x1(0)(1),x1(0)(2),…,x1(0)(n))为系统的特征数据序列,Xi(0)=(xi(0)(1),xi(0)(2),…,xi(0)(n))(i=2,3,…,N)为系统的相关因素数据序列,Xi(1)=(xi(1)(1),xi(1)(2),…,xi(1)(n))(i=1,2,…,N)为Xi(0)(i=1,2,…,N)的一阶累加生成序列,其中

    $ x_i^{(1)}(k) = \sum\limits_{j = 1}^k {x_i^{(0)}} (j)\quad k = 1,2, \cdots ,n $ (1)

    Z1(1)=(z1(1)(2),z1(1)(3),…,z1(1)(n))为X1(1)的紧邻均值生成序列,其中

    $ z_1^{(1)}(k) = \frac{1}{2}\left( {x_1^{(1)}(k - 1) + x_1^{(1)}(k)} \right) $ (2)

    则称

    $ x_1^{(0)}(k) + az_1^{(1)}(k) = \sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(k) $ (3)

    为GM(1,N)模型.

    定理1  Xi(0)Xi(1)(i=1,2,…,N)和Z1(1)参见定义1,则GM(1,N)模型参数的最小二乘估计值$\mathit{\boldsymbol{\hat a}}$=(ab2b3,…,bn)T满足

    (1) 当nN+1,且|BTB|≠0时,$\mathit{\boldsymbol{\hat a}}$=BT(BTB)-1Y

    (2) 当n=N+1,且|B|≠0时,$\mathit{\boldsymbol{\hat a}}$=B-1Y

    (3) 当nN+1,且|BTB|≠0时,$\mathit{\boldsymbol{\hat a}}$=(BTB)-1BTY

    其中

    $ \mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} { - z_1^{(1)}(2)}&{x_2^{(1)}(2)}& \cdots &{x_N^{(1)}(2)}\\ { - z_1^{(1)}(3)}&{x_2^{(1)}(3)}& \cdots &{x_N^{(1)}(3)}\\ \vdots & \vdots & \ddots & \vdots \\ { - z_1^{(1)}(n)}&{x_2^{(1)}(n)}& \cdots &{x_N^{(1)}(n)} \end{array}} \right]\quad \mathit{\boldsymbol{Y}} = \left( {\begin{array}{*{20}{c}} {x_1^{(0)}(2)}\\ {x_1^{(0)}(3)}\\ \vdots \\ {x_1^{(0)}(n)} \end{array}} \right) $ (4)

        证明过程可参考文献[10],在此略.

    定义2  若$\mathit{\boldsymbol{\hat a}}$=(ab2b3,…,bn)T,则称

    $ \frac{{{\rm{d}}x_1^{(1)}(t)}}{{{\rm{d}}t}} + ax_1^{(1)}(t) = \sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(t) $ (5)

    为GM(1,N)模型的白化微分方程.

    定理2  Xi(0)Xi(1)(i=1,2,…,N),Z1(1)$\mathit{\boldsymbol{\hat a}}$参见定义1和定理1,初始条件$\hat x$1(1)(1)=x1(1)(1)=x1(0)(1),当Xi(1)(i=2,3,…,N)变化幅度不大时,GM(1,N)模型的时间响应式为:

    $ \hat x_1^{(1)}(k) = \left[ {x_1^{(0)}(1) - \frac{1}{a}\sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(1)} \right]{{\rm{e}}^{ - a(k - 1)}} + \frac{1}{a}\sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(k) $ (6)

    还原式为:

    $ \begin{array}{*{20}{c}} {\hat x_1^{(0)}(1) = x_1^{(1)}(1)}\\ {\hat x_1^{(0)}(k) = \hat x_1^{(1)}(k) - \hat x_1^{(1)}(k - 1)\quad k = 2,3, \cdots ,n} \end{array} $ (7)

      求解方程(5)可得

    $ x_1^{(1)}(t) = {{\rm{e}}^{ - at}}\left\{ {\int {{{\rm{e}}^{at}}} \left[ {\sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(t)} \right]{\rm{d}}t + C} \right\} $ (8)

    Xi(1)(i=2,3,…,N)变化幅度不大时,可以认定$\sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(t)$为灰常量,则GM(1,N)模型的时间响应式为:

    $ \hat x_1^{(1)}(k) = C{{\rm{e}}^{ - ak}} + \frac{1}{a}\sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(k) $ (9)

    若令$\hat x$1(1)(1)=x1(1)(1)=x1(0)(1),则可得

    $ C = {{\rm{e}}^a}\left[ {x_1^{(0)}(1) - \frac{1}{a}\sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(1)} \right] $

    于是

    $ \hat x_1^{(1)}(k) = \left[ {x_1^{(0)}(1) - \frac{1}{a}\sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(1)} \right]{{\rm{e}}^{ - a(k - 1)}} + \frac{1}{a}\sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(k) $ (10)
    2 改进的GM(1,N)模型

    对方程(5)两边同时在[k-1,k]上积分,有

    $ \int_{k - 1}^k {\rm{d}} x_1^{(1)}(t) + \int_{k - 1}^k a x_1^{(1)}(t){\rm{d}}t = \sum\limits_{i = 2}^N {\int_{k - 1}^k {{b_i}} } x_i^{(1)}(t){\rm{d}}t $ (11)

    $ x_1^{(1)}(k) - x_1^{(1)}(k - 1) + a\int_{k - 1}^k {x_1^{(1)}} (t){\rm{d}}t = \sum\limits_{i = 2}^N {{b_i}} \int_{k - 1}^k {x_i^{(1)}} (t){\rm{d}}t $ (12)

    $ x_1^{(0)}(k) + a\int_{k - 1}^k {x_1^{(1)}} (t){\rm{d}}t = \sum\limits_{i = 2}^N {{b_i}} \int_{k - 1}^k {x_i^{(1)}} (t){\rm{d}}t $ (13)

    若将z1(1)(k)=$\frac{1}{2}$[x1(1)(k-1)+x1(1)(k)]近似代替$\int_{k - 1}^k {x_1^{(1)}} (t){\rm{d}}t$,并视xi(1)(t)(i=2,3,…,N)为灰常量,则该模型为传统的GM(1,N)模型:

    $ x_1^{(0)}(k) + az_1^{(1)}(k) = \sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(k) $

    事实上,将xi(1)(t)(i=2,3,…,N)视为灰常量,即认为Xi(1)(i=2,3,…,N)变化幅度不大,而在实际中Xi(1)(i=2,3,…,N)变化幅度很大的情况并不少见,所以在离散化过程中易带来较大误差.为尽量减少误差,考虑到一阶累加后的序列Xi(1)(i=1,2,…,N)一般具有灰指数率,而且一般是近似非齐次指数增长率,于是可假设

    $ x_i^{(1)}(t) = {\alpha _i}{{\rm{e}}^{{\beta _i}t}} + {\gamma _i},i = 1,2, \cdots ,N $ (14)

    其中αiβiγi为待定常数,且满足

    $ x_i^{(1)}(k) = {\alpha _i}{{\rm{e}}^{{\beta _i}k}} + {\gamma _i}\quad k = 1,2, \cdots ,n\quad i = 1,2, \cdots ,N $ (15)

    在区间[k-1,k]上,令

    $ \begin{array}{l} D_i^{(1)}(k) = \int_{k - 1}^k {x_i^{(1)}} (t){\rm{d}}t = \int_{k - 1}^k {\left( {{\alpha _i}{{\rm{e}}^{{\beta _i}t}} + {\gamma _i}} \right)} {\rm{d}}t = \left. {\left( {\frac{{{\alpha _i}}}{{{\beta _i}}}{{\rm{e}}^{{\beta _i}t}} + {\gamma _i}t} \right)} \right|_{k - 1}^k = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{{{\beta _i}}}\left[ {{\alpha _i}{{\rm{e}}^{{\beta _i}k}} - {\alpha _i}{{\rm{e}}^{{\beta _i}(k - 1)}}} \right] + {\gamma _i} = \frac{1}{{{\beta _i}}}\left[ {x_i^{(1)}(k) - x_i^{(1)}(k - 1)} \right] + {\gamma _i} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{{{\beta _i}}}x_i^{(0)}(k) + {\gamma _i} \end{array} $ (16)

    为了求出αiβiγi的值,在区间[k-1,k]上充分利用xi(1)(k-1)和xi(1)(k)的值,并利用初始点xi(1)(1)的值建立方程组:

    $ \left\{ {\begin{array}{*{20}{l}} {x_i^{\left( 1 \right)}(k - 1) = {\alpha _i}{{\rm{e}}^{{\beta _i}(k - 1)}} + {\gamma _i}\;\;\;\;\;\;\;\left( {17} \right)}\\ {x_i^{\left( 1 \right)}(k) = {\alpha _i}{{\rm{e}}^{{\beta _i}k}} + {\gamma _i}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {18} \right)}\\ {x_i^{\left( 1 \right)}(1) = {\alpha _i}{{\rm{e}}^{{\beta _i}}} + {\gamma _i}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {19} \right)} \end{array}} \right. $

    由式(18)减式(17)可得

    $ x_i^{(1)}(k) - x_i^{(1)}(k - 1) = {\alpha _i}{{\rm{e}}^{{\beta _i}k}} - {\alpha _i}{{\rm{e}}^{{\beta _i}(k - 1)}} = {\alpha _i}{{\rm{e}}^{{\beta _i}(k - 1)}}\left( {{{\rm{e}}^{{\beta _i}}} - 1} \right) $

    $ x_i^{(0)}(k) = {\alpha _i}{{\rm{e}}^{{\beta _i}(k - 1)}}\left( {{{\rm{e}}^{{\beta _i}}} - 1} \right) $ (20)

    于是

    $ x_i^{(0)}(k - 1) = {\alpha _i}{{\rm{e}}^{{\beta _i}(k - 2)}}\left( {{{\rm{e}}^{{\beta _i}}} - 1} \right) $ (21)

    由式(20)和式(21)可得

    $ \frac{{x_i^{(0)}(k)}}{{x_i^{(0)}(k - 1)}} = \frac{{{\alpha _i}{{\rm{e}}^{{\beta _i}(k - 1)}}\left( {{{\rm{e}}^{{\beta _i}}} - 1} \right)}}{{{\alpha _i}{{\rm{e}}^{{\beta _i}(k - 2)}}\left( {{{\rm{e}}^{{\beta _i}}} - 1} \right)}} = {{\rm{e}}^{{\beta _i}}} $

    $ {\beta _i} = \ln \frac{{x_i^{(0)}(k)}}{{x_i^{(0)}(k - 1)}} = \ln x_i^{(0)}(k) - \ln x_i^{(0)}(k - 1) $ (22)

    由式(20)和式(22)可得

    $ \begin{array}{l} {\alpha _i} = \frac{{x_i^{(0)}(k)}}{{{{\rm{e}}^{{\beta _i}(k - 1)}}\left( {{{\rm{e}}^{{\beta _i}}} - 1} \right)}} = \frac{{x_i^{(0)}(k)}}{{{{\left[ {\frac{{x_i^{(0)}(k)}}{{x_i^{(0)}(k - 1)}}} \right]}^{k - 1}} \cdot \left[ {\frac{{x_i^{(0)}(k)}}{{x_i^{(0)}(k - 1)}} - 1} \right]}}{\rm{ = }}\\ \;\;\;\;\;\;\;\frac{{x_i^{(0)}(k) \cdot {{\left[ {\frac{{x_i^{(0)}(k - 1)}}{{x_i^{(0)}(k)}}} \right]}^{k - 1}}}}{{\frac{{x_i^{(0)}(k) - x_i^{(0)}(k - 1)}}{{x_i^{(0)}(k - 1)}}}} = \\ \;\;\;\;\;\;\;\frac{{{{\left[ {x_i^{(0)}(k - 1)} \right]}^k}}}{{{{\left[ {x_i^{(0)}(k)} \right]}^{k - 2}} \cdot \left[ {x_i^{(0)}(k) - x_i^{(0)}(k - 1)} \right]}} \end{array} $ (23)

    由式(19)、式(22)和式(23)可得

    $ \begin{array}{l} {\gamma _i} = x_i^{(1)}(1) - {\alpha _i}{{\rm{e}}^{{\beta _i}}} = x_i^{(0)}(1) - \frac{{{{\left[ {x_i^{(0)}(k - 1)} \right]}^k}}}{{{{\left[ {x_i^{(0)}(k)} \right]}^{k - 2}} \cdot \left[ {x_i^{(0)}(k) - x_i^{(0)}(k - 1)} \right]}} \cdot \frac{{x_i^{(0)}(k)}}{{x_i^{(0)}(k - 1)}} = \\ \;\;\;\;\;\;x_i^{(0)}(1) - \frac{{{{\left[ {x_i^{(0)}(k - 1)} \right]}^{k - 1}}}}{{{{\left[ {x_i^{(0)}(k)} \right]}^{k - 3}} \cdot \left[ {x_i^{(0)}(k) - x_i^{(0)}(k - 1)} \right]}} \end{array} $ (24)

    将式(22)、式(23)和式(24)代入Di(1)(k),可得

    $ D_i^{(1)}(k) = \frac{{x_i^{(0)}(k)}}{{\ln x_i^{(0)}(k) - \ln x_i^{(0)}(k - 1)}} + x_i^{(0)}(1) - \frac{{{{\left[ {x_i^{(0)}(k - 1)} \right]}^{k - 1}}}}{{{{\left[ {x_i^{(0)}(k)} \right]}^{k - 3}} \cdot \left[ {x_i^{(0)}(k) - x_i^{(0)}(k - 1)} \right]}} $ (25)

    Di(1)(k)代替$\int_{k - 1}^k {x_i^{(1)}} (t){\rm{d}}t$的值,可得到

    $ x_1^{(0)}(k) + aD_1^{(1)}(k) = \sum\limits_{i = 2}^N {{b_i}} D_i^{(1)}(k) $

    其中

    $ D_i^{\left( 1 \right)}(k) = \frac{{x_i^{(0)}(k)}}{{\ln x_i^{(0)}(k) - \ln x_i^{(0)}(k - 1)}} + x_i^{(0)}(1) - \frac{{{{\left[ {x_i^{(0)}(k - 1)} \right]}^{k - 1}}}}{{{{\left[ {x_i^{(0)}(k)} \right]}^{k - 3}} \cdot \left[ {x_i^{(0)}(k) - x_i^{(0)}(k - 1)} \right]}} $

    定义3  Xi(0)Xi(1)(i=1,2,…,N)见定义1,则称

    $ x_1^{(0)}(k) + aD_1^{(1)}(k) = \sum\limits_{i = 2}^N {{b_i}} D_i^{(1)}(k) $ (26)

    为改进的GM(1,N)模型,其中

    $ D_i^{(1)}(k) = \frac{{x_i^{(0)}(k)}}{{\ln x_i^{(0)}(k) - \ln x_i^{(0)}(k - 1)}} + x_i^{(0)}(1) - \frac{{{{\left[ {x_i^{(0)}(k - 1)} \right]}^{k - 1}}}}{{{{\left[ {x_i^{(0)}(k)} \right]}^{k - 3}} \cdot \left[ {x_i^{(0)}(k) - x_i^{(0)}(k - 1)} \right]}} $ (27)

    方程

    $ \frac{{{\rm{d}}x_1^{(1)}(t)}}{{{\rm{d}}t}} + ax_1^{(1)}(t) = \sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(t) $

    为改进的GM(1,N)模型的白化微分方程.

    定理3  设$\mathit{\boldsymbol{\hat a}}$=[ab2b3,…,bN]T

    $ \mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} { - D_1^{(1)}(2)}&{D_2^{(1)}(2)}& \cdots &{D_N^{(1)}(2)}\\ { - D_1^{(1)}(3)}&{D_2^{(1)}(3)}& \cdots &{D_N^{(1)}(3)}\\ \vdots & \vdots & \ddots & \vdots \\ { - D_1^{(1)}(n)}&{D_2^{(1)}(n)}& \cdots &{D_N^{(1)}(n)} \end{array}} \right]\quad \mathit{\boldsymbol{Y}} = \left( {\begin{array}{*{20}{c}} {x_1^{(0)}(2)}\\ {x_1^{(0)}(3)}\\ \vdots \\ {x_1^{(0)}(n)} \end{array}} \right) $

    1) 当nN+1且|BTB|≠0时,$\mathit{\boldsymbol{\hat a}}$=BT(BTB)-1Y

    2) 当n=N+1且|B|≠0时,$\mathit{\boldsymbol{\hat a}}$=B-1Y

    3) 当nN+1且|BTB|≠0时,$\mathit{\boldsymbol{\hat a}}$=(BTB)-1BTY.

        证明方法同定理1,此略.

    定理4  改进的GM(1,N)模型见定义3,则其时间响应函数为:

    1) 当k=1时,

    $ \hat x_1^{(1)}(1) = C{{\rm{e}}^{ - a}} $ (28)

    2) 当k≥2时,

    $ \hat x_1^{(1)}(k) = {{\rm{e}}^{ - ak}}\left\{ {\sum\limits_{p = 2}^k {\sum\limits_{i = 2}^N {\left[ {\frac{{{b_i}{\alpha _i}}}{{a + {\beta _i}}}{{\rm{e}}^{\left( {a + {\beta _i}} \right)(p - 1)}}\left( {{{\rm{e}}^{a + {\beta _i}}} - 1} \right) + \frac{{{b_i}{\gamma _i}}}{a}{{\rm{e}}^{a(p - 1)}}\left( {{{\rm{e}}^a} - 1} \right)} \right]} } + C} \right\} $ (29)

    其中

    $ {\alpha _i} = \frac{{{{\left[ {x_i^{(0)}(p - 1)} \right]}^p}}}{{{{\left[ {x_i^{(0)}(p)} \right]}^{p - 2}} \cdot \left[ {x_i^{(0)}(p) - x_i^{(0)}(p - 1)} \right]}} $
    $ {\beta _i} = \ln x_i^{(0)}(p) - \ln x_i^{(0)}(p - 1) $
    $ {\gamma _i} = x_i^{(0)}(1) - \frac{{{{\left[ {x_i^{(0)}(p - 1)} \right]}^{p - 1}}}}{{{{\left[ {x_i^{(0)}(p)} \right]}^{p - 3}} \cdot \left[ {x_i^{(0)}(p) - x_i^{(0)}(p - 1)} \right]}} $

    C为待定常数.

        由改进的GM(1,N)模型的白化微分方程可解得

    $ \begin{array}{l} x_1^{(1)}(t) = {{\rm{e}}^{ - at}}\left\{ {\int {{{\rm{e}}^{at}}} \left[ {\sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(t)} \right]{\rm{d}}t + {C_1}} \right\} = \\ \;\;\;\;\;\;\;\;\;\;\;\;{{\rm{e}}^{ - at}}\left\{ {\int_1^t {{{\rm{e}}^{a\omega }}} \left[ {\sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(\omega )} \right]{\rm{d}}\omega + C} \right\} \end{array} $

    $ \hat x_1^{(1)}(k) = {{\rm{e}}^{ - ak}}\left\{ {\int_1^k {{{\rm{e}}^{a\omega }}} \left[ {\sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(\omega )} \right]{\rm{d}}\omega + C} \right\} $

    k=1时,$\hat x$1(1)(1)=${{\rm{e}}^{ - a}}\left\{ {\int_1^1 {{{\rm{e}}^{a\omega }}} \left[ {\sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(\omega )} \right]{\rm{d}}\omega + C} \right\} = C{{\rm{e}}^{ - a}}$;当k≥2时,考虑到累加后序列只是近似服从指数增长,为了充分利用每个数据,尽量避免由于个别数据误差较大而带来的影响,本文将区间[1,k]分解为多个长度为1的区间之和,即

    $ \left[ {1,k} \right] = \left[ {1,2} \right] \cup \left[ {2,3} \right] \cup \cdots \cup \left[ {k - 1,k} \right] = \bigcup\limits_{p = 2}^k {\left[ {p - 1,p} \right]} $

    于是

    $ \hat x_1^{(1)}(k) = {{\rm{e}}^{ - ak}}\left\{ {\sum\limits_{p = 2}^k {\int_{p - 1}^p {{{\rm{e}}^{a\omega }}} } \left[ {\sum\limits_{i = 2}^N {{b_i}} x_i^{(1)}(\omega )} \right]{\rm{d}}\omega + C} \right\} $

    xi(1)(ω)=αieβiω+γi,其中

    $ {{\alpha _i} = \frac{{{{\left[ {x_i^{(0)}(p - 1)} \right]}^p}}}{{{{\left[ {x_i^{(0)}(p)} \right]}^{p - 2}} \cdot \left[ {x_i^{(0)}(p) - x_i^{(0)}(p - 1)} \right]}}} $
    $ {\beta _i} = \ln x_i^{(0)}(p) - \ln x_i^{(0)}(p - 1) $
    $ {\gamma _i} = x_i^{(0)}(1) - \frac{{{{\left[ {x_i^{(0)}(p - 1)} \right]}^{p - 1}}}}{{{{\left[ {x_i^{(0)}(p)} \right]}^{p - 3}} \cdot \left[ {x_i^{(0)}(p) - x_i^{(0)}(p - 1)} \right]}} $

    $ \begin{array}{l} \hat x_1^{\left( 1 \right)}(k) = {{\rm{e}}^{ - ak}}\left\{ {\sum\limits_{p = 2}^k {\int_{p - 1}^p {{{\rm{e}}^{a\omega }}} } \left[ {\sum\limits_{i = 2}^N {{b_i}} \left( {{\alpha _i}{{\rm{e}}^{{\beta _i}\omega }} + {\gamma _i}} \right)} \right]{\rm{d}}\omega + C} \right\} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;{{\rm{e}}^{ - ak}}\left\{ {\sum\limits_{p = 2}^k {\int_{p - 1}^p {\left[ {\sum\limits_{i = 2}^N {{b_i}} {\alpha _i}{{\rm{e}}^{\left( {a + {\beta _i}} \right)\omega }} + \sum\limits_{i = 2}^N {{b_i}} {\gamma _i}{{\rm{e}}^{a\omega }}} \right]} } {\rm{d}}\omega + C} \right\} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;{{\rm{e}}^{ - ak}}\left\{ {\sum\limits_{p = 2}^k {\left[ {\left. {\sum\limits_{i = 2}^N {\frac{{{b_i}{\alpha _i}}}{{a + {\beta _i}}}} {{\rm{e}}^{\left( {a + {\beta _i}} \right)\omega }}} \right|_{p - 1}^p + \left. {\sum\limits_{i = 2}^N {\frac{{{b_i}{\gamma _i}}}{a}} {{\rm{e}}^{a\omega }}} \right|_{p - 1}^p} \right]} + C} \right\} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;{{\rm{e}}^{ - ak}}\left\{ {\sum\limits_{p = 2}^k {\sum\limits_{i = 2}^N {\left[ {\frac{{{b_i}{\alpha _i}}}{{a + {\beta _i}}}{{\rm{e}}^{\left( {a + {\beta _i}} \right)(p - 1)}}\left( {{{\rm{e}}^{a + {\beta _i}}} - 1} \right) + \frac{{{b_i}{\gamma _i}}}{a}{{\rm{e}}^{a(p - 1)}}\left( {{{\rm{e}}^a} - 1} \right)} \right]} } + C} \right\} \end{array} $

    对于定理4中常数C的确定一般有两种方法:

    1) 取定$\hat x$1(1)(1)=x1(1)(1)=x1(0)(1),可得$\hat x$1(1)(1)=Ce-a,则

    $ C = \hat x_1^{(1)}(1){{\rm{e}}^a} $ (30)

    2) 以绝对误差的平方和最小为准则来确定,见定理5.

    定理5  若以绝对误差的平方和最小为准则,改进的GM(1,N)模型的时间响应函数中常数C的最优值为

    $ C = \frac{{\sum\limits_{k = 1}^n {{{\rm{e}}^{ - ak}}} x_1^{(1)}(k) - \sum\limits_{k = 2}^n {{{\rm{e}}^{ - 2ak}}} \left\{ {\sum\limits_{p = 2}^k {\sum\limits_{i = 2}^N {\left[ {\frac{{{b_i}{\alpha _i}}}{{a + {\beta _i}}}{{\rm{e}}^{\left( {a + {\beta _i}} \right)(p - 1)}}\left( {{{\rm{e}}^{a + {\beta _i}}} - 1} \right) + \frac{{{b_i}{\gamma _i}}}{a}{{\rm{e}}^{a(p - 1)}}\left( {{{\rm{e}}^a} - 1} \right)} \right]} } } \right\}}}{{\sum\limits_{k = 1}^n {{{\rm{e}}^{ - 2ak}}} }} $ (31)

      若记

    $ m(k) = \sum\limits_{p = 2}^k {\sum\limits_{i = 2}^N {\left[ {\frac{{{b_i}{\alpha _i}}}{{a + {\beta _i}}}{{\rm{e}}^{\left( {a + {\beta _i}} \right)(p - 1)}}\left( {{{\rm{e}}^{a + {\beta _i}}} - 1} \right) + \frac{{{b_i}{\gamma _i}}}{a}{{\rm{e}}^{a(p - 1)}}\left( {{{\rm{e}}^a} - 1} \right)} \right]} } $ (32)

    则模型的绝对误差平方和

    $ \begin{array}{l} f(C) = \sum\limits_{k = 1}^n {{{\left[ {\hat x_1^{(1)}(k) - x_1^{(1)}(k)} \right]}^2}} = {\left[ {\begin{array}{*{20}{c}} {\hat x_1^{(1)}(1) - x_1^{(1)}(1)} \end{array}} \right]^2} + \sum\limits_{k = 2}^n {{{\left[ {\hat x_1^{(1)}(k) - x_1^{(1)}(k)} \right]}^2}} = \\ \;\;\;\;\;\;\;\;\;\;\;{\left[ {C{{\rm{e}}^{ - a}} - x_1^{(1)}(1)} \right]^2} + \sum\limits_{k = 2}^n {{{\left\{ {{{\rm{e}}^{ - ak}}[m(k) + C] - x_1^{(1)}(k)} \right\}}^2}} \end{array} $

    求导可得

    $ \begin{array}{l} f'(C) = 2\left[ {C{{\rm{e}}^{ - a}} - x_1^{(1)}(1)} \right]{{\rm{e}}^{ - a}} + 2\sum\limits_{k = 2}^n {\left\{ {{{\rm{e}}^{ - ak}}[m(k) + C] - x_1^{(1)}(k)} \right\}} {{\rm{e}}^{ - ak}} = \\ \;\;\;\;\;\;\;\;\;\;\;\;2C{{\rm{e}}^{ - 2a}} - 2x_1^{(1)}(1){{\rm{e}}^{ - a}} + 2C\sum\limits_{k = 2}^n {{{\rm{e}}^{ - 2ak}}} + 2\sum\limits_{k = 2}^n {{{\rm{e}}^{ - 2ak}}} m(k) - 2\sum\limits_{k = 2}^n {{{\rm{e}}^{ - ak}}} x_1^{(1)}(k) = \\ \;\;\;\;\;\;\;\;\;\;\;\;2C\sum\limits_{k = 1}^n {{{\rm{e}}^{ - 2ak}}} + 2\sum\limits_{k = 2}^n {{{\rm{e}}^{ - 2ak}}} m(k) - 2\sum\limits_{k = 1}^n {{{\rm{e}}^{ - ak}}} x_1^{(1)}(k) \end{array} $

    f(C)=0,得

    $ C\sum\limits_{k = 1}^n {{{\rm{e}}^{ - 2ak}}} = \sum\limits_{k = 1}^n {{{\rm{e}}^{ - ak}}} x_1^{(1)}(k) - \sum\limits_{k = 2}^n {{{\rm{e}}^{ - 2ak}}} m(k) $

    即得驻点

    $ \begin{array}{*{20}{c}} {C = {C_0} = \frac{{\sum\limits_{k = 1}^n {{{\rm{e}}^{ - ak}}} x_1^{(1)}(k) - \sum\limits_{k = 2}^n {{{\rm{e}}^{ - 2ak}}} m(k)}}{{\sum\limits_{k = 1}^n {{{\rm{e}}^{ - 2ak}}} }} = }\\ {\frac{{\sum\limits_{k = 1}^n {{{\rm{e}}^{ - ak}}} x_1^{(1)}(k) - \sum\limits_{k = 2}^n {{{\rm{e}}^{ - 2ak}}} \left\{ {\sum\limits_{p = 2}^k {\sum\limits_{i = 2}^N {\left[ {\frac{{{b_i}{\alpha _i}}}{{a + {\beta _i}}}{{\rm{e}}^{\left( {a + {\beta _i}} \right)(p - 1)}}\left( {{{\rm{e}}^{a + {\beta _i}}} - 1} \right) + \frac{{{b_i}{\gamma _i}}}{a}{{\rm{e}}^{a(p - 1)}}\left( {{{\rm{e}}^a} - 1} \right)} \right]} } } \right\}}}{{\sum\limits_{k = 1}^n {{{\rm{e}}^{ - 2ak}}} }}} \end{array} $

    又由于

    $ f'(C) = 2\sum\limits_{k = 1}^n {{{\rm{e}}^{ - 2ak}}} > 0 $

    C=C0f(C)的极小值点.由于极值点唯一,则极小值点同时为最小值点,即f(C)在C=C0处取得最小值.

    3 模型的检验

    文献[10]和文献[11]通过借用经典GM(1,1)幂模型的建模机理,提出了灰色多变量GM(1,N)幂模型,并通过实际应用说明了模型具有较好的实用性,是对传统GM(1,N)模型的有效推广.为了验证本文提出的改进的GM(1,N)模型的模拟预测效果,下面通过两个实例对传统GM(1,N)模型、GM(1,N)幂模型和本文中两种不同初值条件下的改进的GM(1,N)模型加以比较.为方便起见,记利用点(1,x1(0)(1))确定初值的改进的GM(1,N)模型为改进的GM(1,N)模型1,记绝对误差平方和最小意义下的改进的GM(1,N)模型为改进的GM(1,N)模型2.

    3.1 数值模拟

    为了验证改进的GM(1,N)模型的模拟预测效果,现采用一个算例进行模拟分析.取x1(0)(k)=20e0.09kx2(0)(k)=30e0.06k,令k=1,2,…,6,可得系统的特征序列:

    $ \mathit{\boldsymbol{X}}_1^{(0)} = (21.8835,23.9443,26.1993,28.6666,31.3662,34.3201) $

    系统的相关因素序列:

    $ \mathit{\boldsymbol{X}}_2^{(0)} = (31.8551,33.8249,35.9165,38.1375,40.4958,42.9999) $

    对原始序列Xi(0)(i=1,2)分别建立传统GM(1,2)模型、GM(1,2)幂模型、改进的GM(1,2)模型1和改进的GM(1,2)模型2,得到的具体模拟结果和相对误差见表 1.

    表 1 3种模型的模拟结果和相对误差

    表 1可看出,传统的GM(1,2)模型和GM(1,2)幂模型的平均相对误差分别为6.50%和4.39%,而两种改进的GM(1,2)模型的平均相对误差分别为2.47%和2.48%,建模精度显著提高.

    3.2 应用实例

    表 2为2009-2014年广东省高新技术投入(R & D占GDP比例)和产出(产品产值)数据(来源:广东科技统计网,http://www.sts.gd.cn/).

    表 2 2009-2014年广东省高新技术投入产出数据

    以广东省高新技术产品产值序列为系统特征序列,即

    $ X_1^{(0)} = (25595.12,29113.16,35592.36,39120.42,45142.12,49256.22) $

    以R & D占GDP比例序列为系统相关因素序列,即

    $ X_2^{(0)} = (1.65,1.76,1.96,2.17,2.32,2.37) $

    下面对X1(0)X2(0)分别建立传统的GM(1,2)模型、GM(1,2)幂模型、改进的GM(1,2)模型1和改进的GM(1,2)模型2,得到的模拟结果见表 3.

    表 3 3种模型对广东省高新技术产业总产值的模拟结果

    表 3可以看出,传统GM(1,2)模型和GM(1,2)幂模型的平均相对误差分别为13.86%和7.29%,而本文提出的两种改进的GM(1,2)模型的平均相对误差分别为5.78%和6.13%,说明改进的GM(1,2)模型能更准确反映广东省高新技术产业总产值的变化规律,更具实用性.

    4 结论

    灰色多变量模型能够有效反映现实经济社会系统中普遍存在着的相互依赖、相互制约的关系.本文通过分析传统GM(1,N)模型的不足,重构了传统GM(1,N)模型的背景值,并对模型的时间响应式进行了近似计算.在初值的选取上,提供了两种不同的方式确定,在实际应用时可根据具体情况灵活选择.最后数值实例表明改进后的GM(1,N)模型能有效提高建模精度.本文提出的模型是在一阶累加序列具有近似非齐次指数增长率的前提下改进的,如果原始序列具有近似齐次指数率,那么改进效果显著;如果原始序列为振荡序列,那么改进效果有限,甚至可能出现精度反而降低的现象.因此,如何有效增强GM(1,N)模型的适应性,是值得进一步研究的问题.

    参考文献
    [1]
    DENG J L. Control Problems of Grey Systems[J]. Systems & Control Letters, 1982, 1(5): 288-294.
    [2]
    邓聚龙. 灰理论基础[M]. 武汉: 华中科技大学出版社, 2002.
    [3]
    刘帅, 袁汝华. 基于Leslie矩阵的水资源供需平衡灰色预测研究[J]. 贵州师范大学学报(自然科学版), 2016, 34(1): 22-27. DOI:10.3969/j.issn.1004-5570.2016.01.004
    [4]
    TIEN T L. A Research on the Grey Prediction Model GM(1, n)[J]. Applied Mathematics and Computation, 2012, 218(9): 4903-4916. DOI:10.1016/j.amc.2011.10.055
    [5]
    何满喜. 建立GM(1, N)预测模型的新方法[J]. 农业系统科学与综合研究, 1997, 13(4): 241-244.
    [6]
    GUO H, XIAO X P, FORREST J. A Research on a Comprehensive Adaptive Grey Prediction Model CAGM(1, N)[J]. Applied Mathematics and Computation, 2013, 225(1): 216-227.
    [7]
    SHEN Y, SUN H Y. Optimization of Background Values of GM(1, N) Model and Its Application[J]. International Journal of Information Processing and Management, 2013, 4(1): 58-64. DOI:10.4156/ijipm
    [8]
    WANG Z X, HAO P. An Improved Grey Multivariable Model for Predicting Industrial Energy Consumption in China[J]. Applied Mathematical Modelling, 2016, 40(11-12): 5745-5758. DOI:10.1016/j.apm.2016.01.012
    [9]
    ZENG B, LUO C M, LIU S F, et al. Development of an Optimization Method for the GM(1, N) Model[J]. Engineering Applications of Artificial Intelligence, 2016, 55: 353-362. DOI:10.1016/j.engappai.2016.08.007
    [10]
    王正新. 灰色多变量GM(1, N)幂模型及其应用[J]. 系统工程理论与实践, 2014, 34(9): 2357-2363.
    [11]
    WANG Z X, YE D J. Forecasting Chinese Carbon Emissions from Fossil Energy Consumption Using Non-Linear Grey Multivariable Models[J]. Journal of Cleaner Production, 2017, 142: 600-612. DOI:10.1016/j.jclepro.2016.08.067
    [12]
    王正新. 多变量时滞GM(1, N)模型及其应用[J]. 控制与决策, 2015, 30(12): 2298-2304.
    [13]
    毛树华, 高明运, 肖新平. 分数阶累加时滞GM(1, N, τ)模型及其应用[J]. 系统工程理论与实践, 2015, 35(2): 430-436.
    [14]
    王正新. 具有交互效应的多变量GM(1, N)模型[J]. 控制与决策, 2017, 32(3): 515-520.
    [15]
    HSU L C. Forecasting the Output of Integrated Circuit Industry Using Genetic Algorithm Based Multivariable Grey Optimization Models[J]. Expert Systems with Applications, 2009, 36(4): 7898-7903. DOI:10.1016/j.eswa.2008.11.004
    An Improved Grey Multivariable GM(1, N) Model and Its Application
    ZENG Liang     
    Department of Basic Courses, Guangdong Polytechnic College, Zhaoqing Guangdong 526100, China
    Abstract: For the traditional grey multivariable GM(1, N) model, the time response formula is inaccurate and the modeling accuracy is not high. Using the characteristics of the first-order cumulative sequence of each variable with an approximately non-homogeneous exponential increment rate, the background value of the traditional GM(1, N) model is reconstructed, then an improved GM(1, N) model is proposed, and the parameter estimation formula and approximate time response formula are given. Finally, a numerical simulation is made and an example is given to prove the effectiveness and practicability of the improved GM(1, N) model.
    Key words: grey system    GM(1, N) model    non-homogeneous exponential increment rate    background value    prediction    
    X