自上个世纪80年代邓聚龙教授提出灰色系统理论^[1]以来, 灰色系统理论在许多领域得到了广泛的应用和发展^[2-6]. GM(1, 1) 模型是灰色系统理论的核心内容之一, 它主要针对“小样本、贫信息”的不确定系统.但是, 传统的GM(1, 1) 模型具有预测精度不高的问题, 如何提高GM(1, 1) 模型的预测精度, 已经成为广大研究人员最关注的问题.而影响GM(1, 1) 模型预测精度的主要因素有初始条件的选取、背景值的重构和参数估计方法的改进, 其中背景值的重构, 具有决定性的作用.传统方法实质上是使用紧邻均值构造背景值, 误差通常较大, 从而导致模型预测的偏差也较大, 预测精度自然达不到要求.近年来不少学者提出了提高灰色GM(1, 1) 模型预测精度的方法, 文献[7-9]从背景值的几何意义出发, 进行了系列研究, 提出了若干种背景值的构造方式; 文献[10]利用拉格朗日插值公式对背景值进行重构, 对传统模型的背景值进行了改进, 并利用最小二乘法对初始值进行了优化; 李俊蜂等^[11]提出一种基于数值分析中的插值法和Newton-Cotes公式的背景值构造方法.但是, 众所周知, 当$n$较大时, 高次插值会出现Runge现象, 造成的误差可能会很大, 故不能通过提高阶的方法来提高求积精度; 文献[12]提出一种利用梯形面积和矩形面积相结合的正负误差补偿方法进行背景值的重构方法, 取得一定的效果, 但是$n_{1}$值的选取很难把握; 文献[13]提出用Simpson公式构造模型背景值的方法, 拟合精度在一定程度上得到了改进; 王晓佳等^[14]提出利用分段插值函数与Newton插值函数相结合的组合插值方法, 但分段过程比较繁琐; 蒋诗泉等^[15]提出基于复化梯形公式的背景值优化方法, 但是仅选取区间等分数$n=4$的情形, 而对于$n$值取更大的情况未予以考虑.

文中构造出带形状参数的三次Bézier基函数, 然后给出插值函数, 提出了基于复化梯形公式的背景值构造方法, 给定误差限, 结合提出的优化算法, 计算出对应的背景值, 提高了GM(1, 1) 模型的预测精度.以工程实际应用中无镍铸造装甲钢的断裂韧度值为例, 实例研究表明, 提出的方法不仅提高了GM(1, 1) 模型的预测精度, 并且扩展了GM(1, 1) 模型的应用范围.

给定x⁽⁰⁾(t) = {x⁽⁰⁾(1), x⁽⁰⁾(2), …, x⁽⁰⁾(n), x⁽⁰⁾(t)＞0, i = 1, 2, …, n}为原始序列, 可以建立基本的GM(1, 1) 模型, 将序列x⁽⁰⁾(t)进行一次累加(1-AGO), 得到累加序列x⁽¹⁾(t)为x⁽¹⁾(t) ={x⁽¹⁾(1), x⁽¹⁾(2), …, x⁽¹⁾(n), …, x⁽¹⁾(n)}, 其中

$x^{(1)}(k)=\sum\limits_{t=1}^{k}x^{(0)}(t).$

序列的白化微分方程为

$\begin{aligned} \frac{dx^{(1)}t}{dt}+ax^{(1)}(t)=b. \end{aligned}$

(2.1)

灰色GM(1, 1) 模型的基本形式如下

$\begin{aligned} x^{(0)}(k)+az^{(1)}(k)=b, \end{aligned}$

(2.2)

其中$z^{(1)}(k)=\frac{1}{2}(x^{(1)}(k)+x^{(1)}(k-1)), k=2, 3, \dots, n$为背景值.

若$\hat{a}=[a, b]^{T}$为参数列, 且

$Y = \left[ {\begin{array}{*{20}{c}} {{x^0}(2)}\\ {{x^0}(3)}\\ \vdots \\ {{x^0}(n)} \end{array}} \right],B = \left[ {\begin{array}{*{20}{c}} { - {z^1}(2)1}\\ { - {z^1}(3)1}\\ { \vdots \vdots }\\ { - {z^1}(n)1} \end{array}} \right],$

(2.3)

利用最小二乘法求解参数$a$和$b$,

$\begin{aligned} (B^{T}B)^{-1}B^{T}Y=\hat{a}. \end{aligned}$

(2.4)

$x^{(1)}(k)$的GM(1, 1) 模型为

$\begin{aligned} \hat{x}^{(1)}(k+1)=(x^{(0)}(1)-\frac{b}{a})e^{-ak}+\frac{b}{a}, k=1, 2, \dots, n. \end{aligned}$

(2.5)

由于预测方程是对累加数据序列$x^{(1)}(t)$的预测, 进行累减还原, 则可以得到原始数据序列的预测公式

$\begin{aligned} \hat{x}^{(0)}(k+1)=\hat{x}^{(1)}(k+1)-\hat{x}^{(1)}(k)=(1-e^{a})(x^{(0)}(1)-\frac{b}{a})e^{-ak}. \end{aligned}$

(2.6)

从公式(2.6) 可以看出, 模型的拟合预测精度取决于参数$a, b$及初始值$x^{(0)}(1)$的值.而参数$a, b$的值又取决于原始数据序列和背景值$z^{(1)}(k)$.因此, 背景值构造公式的合理性直接影响模型预测精度.传统的背景值计算实质上是梯形的面积, 而实际值正如图 1所示, 应该等于曲线$x^{(1)}(t)$在区间$[k-1, k]$上与$t$轴围成的面积, 这也是传统背景值计算公式的误差产生的根源所在.将积分区间分成若干小区间, 在每个小区间上用低阶求积公式, 然后将这些小区间上的积分值相加作为函数在整个区间上的积分的近似值, 这就是复化求积的思想.下述内容将给出利用三次Bézier基函数构造插值函数, 结合复化梯形公式的GM(1, 1) 模型的背景值优化方法及算法.

1962年, 法国雷诺(Renault)汽车公司的工程师Bézier提出了一种新的方法用来构造著名的Bézier曲线^[16], 设计人员只需要移动控制多边形的控制节点就可以方便地修改曲线的形状, 而且形状的变化完全在预料之中, 但是它无法插值控制节点, 并且需要移动控制节点的位置来改变曲线的形状.下述将构造出一种带形状参数$\lambda$的三次Bézier基函数, 不用调节控制多边形的节点位置, 只须通过调节形状参数$\lambda$的值, 即可调节设计曲线的形状, 而且当参数值$\lambda=2$时, 曲线插值中间控制节点.文中将利用参数$\lambda=2$时, 三次Bézier基函数的插值性质, 来构造插值函数, 从而计算GM(1, 1) 模型的背景值.

定义3.1^[1]    这是定义内容.对$t\in[0,1], \lambda\in R$, 则称关于$t$的多项式

$ \begin{array}{rcl} \left\{\begin{array}{l} b_1(\lambda, t)=(1-\lambda t)(1-t)^{2}, \\ b_2(\lambda, t)=(2+\lambda)(1-t)t, \quad \quad t\in [0,1]\\ b_3(\lambda, t)=(1-\lambda+\lambda t)t^2, \end{array} \right. \end{array} $

(3.1)

为带形状参数$\lambda$的三次Bézier基函数.当$\lambda=0$时, 基函数退化为二次Bézier基函数.因此它是二次Bézier基函数的扩展.当$\lambda=2$时, 公式(3.1) 将简化为

$ \begin{array}{rcl} \left\{\begin{array}{l} b_1(t)=(1-2 t)(1-t)^{2}, \\ b_2(t)=4(1-t)t, \quad \quad \quad \quad \quad t\in [0,1].\\ b_3(t)=t^2(2t-1), \end{array} \right. \end{array} $

(3.2)

若任意给定三个插值节点$x^{(1)}(1), x^{(1)}(2), x^{(1)}(3) $, 插值函数为

$\begin{aligned} f(t)=b_{1}(t)x^{(1)}(1)+b_{2}(t)x^{(1)}(2)+b_{3}(t)x^{(1)}(3), t\in [0,1]. \end{aligned}$

(3.3)

当$t=0.5$时, 基函数$b_{1}(t), b_{3}(t)$的值均为0, 而$b_{2}(t)$的值为1, 则公式(3.3) 为$f(0.5)=x^{(1)}(2)$, 说明函数$f(t)$是插值于中间控制节点的. 图 2为给定三个插值节点(-2, 0), (0, 1), (2, 0), 所得到的插值函数图形, 从图中可以看出, 插值函数是插值中间控制点的.

定义3.2     这是定义内容.设在区间$[a, b]$上, 给定三个插值节点$(k, x^{(1)}(k))$, $(k+1, x^{(1)}(k+1))$, $(k+2, x^{(1)}(k+2))$, 再利用公式(3.2) 所定义的基函数, 则可得区间$[a, b]$上的插值函数

$\begin{aligned} f(t)=b_{1}(t-a)x^{(1)}(k)+b_{2}(t-a)x^{(1)}(k+1)+b_{3}(t-a)x^{(1)}(k+2), t\in [0,1]. \end{aligned}$

(3.4)

定义3.3 ^[17]    将区间$[a, b]$划分为$n$等份, 分点$x_{k}=a+kh$, $h=\frac{b-a}{n}$, $k=0, 1, \cdots, n$, 在每个子区间$[x_{k}, x_{k+1}], (k=0, 1, \cdots, n-1)$上采用梯形公式

$\displaystyle\int_{a}^{b}f(x)dx\approx\frac{b-a}{2}[f(a)+f(b)], $

可得

$ \begin{aligned} I=\int_{a}^{b}f(x)dx=\sum\limits_{k=0}^{n-1}\int_{x_{k}}^{x_{k+1}}f(x)dx=\frac{h}{2}\sum\limits_{k=0}^{n-1}[f(x_{k})+f(x_{k+1})]+R_{n}(f). \end{aligned}$

(3.5)

记

$ \begin{aligned} T_{n}=\frac{h}{2}\sum\limits_{k=0}^{n-1}[f(x_{k})+f(x_{k+1})]=\frac{h}{2}[f(a)+2\sum\limits_{k=1}^{n-1}f(x_{k})+f(b)], \end{aligned}$

(3.6)

称为复化梯形公式, 其中

$ \begin{aligned} R_{n}(f)=I-T_{n}=\sum\limits_{k=0}^{n-1}[-\frac{h^{3}}{12}f^{\prime\prime}(\eta_{k})], \eta_{k}\in(x_{k}, x_{k+1}). \end{aligned}$

(3.7)

定理3.1^[17]  这是定理内容.当$f(x)\in C^{2}[a, b]$时, 则

$ \begin{aligned} \lim\limits_{n\rightarrow\infty}T_{n}=\int_{a}^{b}f(x)dx, \end{aligned}$

(3.8)

即当$n\rightarrow\infty$时, 复化梯形公式是收敛的.

定理3.2     这是定理内容.对复化梯形公式$T_{n}$, 若将区间$[a, b]$进行$2n$等分, 则得到$T_{2n}$:

$ \begin{aligned} T_{2n}=\frac{h}{4}\sum\limits_{k=1}^{n-1}[f(x_{k})+2f(x_{k+\frac{1}{2}})+f(x_{k+1})]=\frac{1}{2}T_{n}+\frac{h}{2}\sum\limits_{k=1}^{n-1}f(x_{k+\frac{1}{2}}). \end{aligned}$

(3.9)

由定义3.1和定理3.2可得

$T_{2n}-T_{n}=\frac{1}{2}(h\sum\limits_{k=1}^{n-1}f(x_{k+\frac{1}{2}})-T_{n}).$

在使用复化梯形公式进行求积计算时, 必须先给出步长$h$或者等分数$n$的值, 步长$h$取得太大, 则精度无法满足要求; 步长取得太小, 则导致计算量过大.在使用复化梯形进行实际计算时, 给定误差值$\varepsilon$, 用$|T_{2n}-T_{n}|＜\varepsilon$成立与否, 作为计算的终止条件.如果满足给定的误差限, 则取$T_{2n}$作为所求的值; 否则, 再将子区间对分, 进行计算, 直到满足给定的误差要求为止.下一部分将给出利用复化梯形公式(3.6) 进行背景值优化的算法.

步骤1    输入原始序列${x^0}\left( i \right),{x^0}\left( i \right) ＞ 0,i = 1,2, \cdots ,n,\varepsilon $以及端点$a, b$的值, 置

$n\leftarrow1, h\leftarrow b-a;$

步骤2    计算利用公式(3.4), 计算插值函数, 并计算

$T\leftarrow\frac{1}{2}[f(a)+f(b)];$

步骤3    计算

$S=\frac{1}{2}[T+h\sum\limits_{j=0}^{n-1}f(a+\frac{1}{2}h+jh)];$

步骤4    如果$|S-T|＞\varepsilon, $那么$T\leftarrow S, n\leftarrow2n, h\leftarrow\frac{1}{2}h$, 转步骤3;

其他, 转步骤5;

步骤5    转向$T$; 结束.

传统灰色GM(1, 1) 模型的背景值$z^{(1)}(k)$的计算是与规模问题$n$无关的常数, 所以时间复杂度为$O(1)$.但对于使用复化梯形公式进行优化计算的背景值来说, 时间复杂度与区间等分数$n$相关, 时间复杂度为$O(n)$.

以无镍铸造装甲钢钢的断裂韧度值为例^[18], 分别计算传统灰色GM(1, 1) 模型和所提出的优化GM(1, 1) 模型在回火温度为600 ℃时试样断裂韧度的预测值, 以温度为200 ℃, 300 ℃, 400 ℃, 500 ℃时的断裂韧度值数据作为原始数据, 对温度为600 ℃时的韧度值进行预测, 并以600 ℃时的实际数据作为检验数据.

给定误差限$\varepsilon$的值, 即可以利用提出的算法, 进行背景值计算.此处取$\varepsilon=0.01$, 选取插值节点$(1, x^{(1)}(1))$, $(2, x^{(1)}(2))$, $(3, x^{(1)}(3))$为例, 利用上述给出的插值函数公式及背景值优化算法, 得到$z^{1}(2)=149.226$.同理可以得到$z^{1}(3)=220.098$, $z^{1}(4)=323.240$.

将利用复化梯形公式优化算法计算得到的背景值, 代入优化的灰色GM(1, 1) 模型, 即可得到预测值, 并计算出它的相对误差值.可以看出, 与传统的GM(1, 1) 模型比较, 优化后的灰色GM(1, 1) 模型的平均相对误差为$4.5457\%$, 比传统的GM(1, 1) 模型降低了$0.6235\%$, 有明显改进.预测误差也由$16.8118\%$降低至$13.1823\%$, 精度提高了$3.6295\%$.可见, 文中所提出的方法与传统模型方法相比, 具有一定的优越性.结果见表 1.

文中给出了基于带三次Bézier插值基函数插值, 结合复化梯形公式的背景值的优化方法, 给定误差限, 利用给出的算法, 即可计算出灰色GM(1, 1) 模型的背景值.与传统的背景值构造方法采用Lagrange ^[10]、Newton ^[11]插值或者组合插值^[13]相比, 文中构造出一组插值基函数, 方法简便, 只须将累加序列节点, 代入即可求出插值函数, 易于操作.实例研究表明, 文中提出的优化模型及算法的有效性, 比传统的灰色GM(1, 1) 模型误差有明显的改进, 扩大了GM(1, 1) 模型的适应性, 优化后的GM(1, 1) 模型更合理, 有助于提高GM(1, 1) 模型的预测精度, 扩展了GM(1, 1) 模型的应用范围.