灰色系统理论

灰色系统理论

-

灰色系统理论及其应用

第一章灰色系统的概念与基本原理

1.1灰色系统理论的产生和发展动态

1982年，北荷兰出版公司出版的《系统与控制通讯》杂志刊载了我国学者邓聚龙

教授的第一篇灰色系统理论论文”灰色系统的控制问题”，同年，《华中工学院学

报》发表邓聚龙教授的第一篇中文论文《灰色控制系统》，这两篇论文的发表标

志着灰色系统这一学科诞生。

1985灰色系统研究会成立，灰色系统相关研究发展迅速。1989海洋出版社出版

英文版《灰色系统论文集》，同年，英文版国际刊物《灰色系统》杂志正式创

刊。目前，国际、国内300多种期刊发表灰色系统论文，许多国际会议把灰色系

统列为讨论专题。国际著名检索已检索我国学者的灰色系统论著3000多次。灰

色系统理论已应用范围已拓展到工业、农业、社会、经济、能源、地质、石油等

众多科学领域，成功地解决了生产、生活和科学研究中的大量实际问题，取得了

显著成果。

1.2几种不确定方法的比较

（系统科学---系统理论）概率统计，模糊数学和灰色系统理论是三种最常用的不

确定系统研究方法。其研究对象都具有某种不确定性，是它们共同的特点。也正

是研究对象在不确定性上

的区别，才派生了这三种各具特色的不确定学科。

模糊数学着重研究“认识不确定”问题，其研究对象具有“内涵明确，外延不明确”

的特点。比如“年轻人”内涵明确，但要你划定一个确定的范围，在这个范围内是

年轻人，范围外不是年轻人，则很难办到了。

概率统计研究的是“随机不确定”现象，考察具有多种可能发生的结果之“随机不确

定”现象中每一种结果发生的可能性大小。要求大样本，并服从某种典型分布。

灰色系统理论着重研究概率统计，模糊数学难以解决的“小样本，贫信息”不确定

性问题，着重研究“外延明确，内涵不明确”的对象。如到2050年，中国要将总

人口控制在15亿到16亿之间，这“15亿到16亿之间“是一个灰概念，其外延很

清楚，但要知道具体数值，则不清楚。

三种不确定性系统研究方法的比较分析

1.3灰色系统理论的基本概念定义1.3.1信息完全明确的系统称为白色系统。

定义1.3.2信息未知的系统称为黑色系统。

定义1.3.3部分信息明确，部分不明确的系统称为灰色系统。

在工程技术、社会、经济、农业、生态、环境等各种系统中经常会遇到信息不完

全的情况。比如：农业方面，农田耕作面积往往因许多非农业的因素而改变，因

此很难准确计算农田产量、产值，这是缺乏耕地面积信息；生物防治方面，害虫

与天敌间的关系即使是明确的，但天敌与饵料、害虫与害虫间的许多关系却不明

确，这是缺乏生物间的关联信息；一项土建工程，尽管材料、设备、施工计划、

图纸是齐备的，可是还很难估计施工进度与质量，这是缺乏劳动力及技术水平的

信息；一般社会经济系统，除了输出的时间数据列（比如产值、产量、总收入、

总支出等）外，其输入数据列不明确或者缺乏，因而难以建立确定的完整的模

型，这是缺乏系统信息；工程系统是客观实体，有明确的“内”、“外”关系（即系

统内部与系统外部，或系统本体与系统环境），可以较清楚地明确输入与输出，

因此可以较方便地分析输入对输

出的影响，可是社会、经济系统是抽象的对象，没有明确的“内”、“外”关系，不

是客观实体，因此就难以分析输入（投入）对输出（产出）的影响，这是缺乏

“模型信息”（即用什么模型，用什么量进行观测控制等信息）。信息不完全的情

况归纳起来有：元素（参数）信息不完全；结构信息不完全；关系信息（特指

“内”、“外”关系）不完全；运行的行为信息不完全。

一个商店可看作是一个系统，在人员、资金、损耗、销售信息完全明确的情况

下，可算出该店的盈利大小、库存多少，可以判断商店的销售态势、资金的周转

速度等，这样的系统是白色系统。

遥远的某个星球，也可以看作一个系统，虽然知道其存在，但体积多大，质量多

少，距离地球多远，这些信息完全不知道，这样的系统是黑色系统。

人体是一个系统，人体的一些外部参数（如身高、体温、脉搏等）是已知的，而

其他一些参数，如人体的穴位有多少，穴位的生物、化学、物理性能，生物的信

息传递等尚未知道透彻，这样的系统是灰色系统。

显然，黑色、灰色、白色都是一种相对的概念。世界上没有绝对的白色系统，因

为任何系统总有未确知的部分，也没有绝对的黑色系统，因为既然一无所知，也

就无所谓该系统的存在了。

1.4灰色系统理论的基本原理

公理1(差异信息原理)“差异“即信息，凡信息必有差异。公理2（解的非唯一性原

理）信息不完全，不确定的解是非唯一的。

公理3（最少信息原理）灰色系统理论的特点是充分开发利用已占有的“最少信息

“。

公理4(认知根据原理)信息是认知的根据。

公理5（新信息优先原理）新信息对认知的作用大于老信息。公理6（灰性不灭

原理）：

“信息完全”是相对的,“信息不完全”是绝对的。

1.5灰色系统理论的主要内容

灰色系统理论经过20多年的发展，现在已经基本建立起一门新兴学科的结构体

系。其主要内容包括以灰色代数系统，灰色方程、灰色矩阵等为基础的理论体

系。以灰色序列生成为基础的方法体系，以灰色关联空间为依托的分析体系。以

灰色模型（GM）为核心的模型体系，以系统分析，评估，建模，预测，决策，

控制，优化为主体的技术体系。

灰色系统的特点

灰色系统理论以“部分信息已知、部分信息未知”的“小样本”、“贫信息”不确定型

系统的研究对象。

（1）用灰色数学来处理不确定量，使之量化。

在数学发展史上，最早研究的是确定型的微分方程，即在拉普拉斯决定论框架内

的数学。他认为一旦有了描写事物的微分方程及初值，就能确知事物任何时候的

运动。随后发展了概率论与数理统计，用随机变量和随机过程来研究事物的状态

和运动。模

糊数学则研究没有清晰界限的事物，如儿童和少年之间没有确定的年龄界限加以

截然划分等，它通过隶属函数来使模糊概念量化，因此能用模糊数学来描述如语

言、不精确推理以及若干人文科学。灰色系统理论则认为不确定量是灰数，用灰

色数学来处理不确定量，同样能使不确定量予以量化。

的方法研究）

1、概率论与数理统计；2、模糊数学；3、灰色数学（灰色系统理论）

（2）充分利用已知信息寻求系统的运动规律。

研究灰色系统的关键是如何使灰色系统白化、模型化、优化。灰色系统视不确定

量为灰色量。提出了灰色系统建模的具体数学方法，它能利用时间序列来确定微

分方程的参数。灰色预测不是把观测到的数据序列视为一个随机过程，而是看作

随时间变化的灰色量或灰色过程，通过累加生成和累减生成逐步使灰色量白化，

从而建立相应于微分方程解的模型并做出预报。这样，对某些大系统和长期预测

问题，就可以发挥作用。

（3）灰色系统理论能处理贫信息系统。

灰色预测模型只要求较短的观测资料即可，这和时间序列分析，多元分析等概率

统计模型要求较长资料很不一样。因此，对于某些只有少量观测数据的项目来

说，灰色预测是一种有用的工具。

1.6灰数

灰数是灰色系统理论的基本“单元“或”细胞“。我们把只知道大概范围而不知道其

确切值的数称为灰数。在应用中，灰1，2，3量化（用确定量

数实际上指在某一个区间或某个一般的数集内取值的不确定数。通常用记号“⊗”

表示灰数。

灰数有以下几类：

1.仅有下界的灰数。有下界而无上界的灰数记为⊗∈

[a,∞]，其中-是灰数⊗的下确界，是确定的数，我

,∞]为⊗的取数域，简称⊗的灰域。们称[a-

2.仅有上界的灰数。有上界而无下界的灰数记为⊗∈

[-∞,a]，其中a---是灰数⊗的上确界，是确定的数。

3.区间灰数。既有下界又有上界的灰数称为区间灰数，

记为⊗∈[--a,a]--

4.

5.连续灰数与离散灰数。黑数与白数。当⊗∈[-∞,+∞]，称⊗为黑数；当⊗∈

[a,a]且a=a--------时，称⊗为白数。

6.本征灰数与非本征灰数。

本征灰数是指不能或暂时还不能找到一个白数作为其“代表”的灰数，比如一般的

事前预测值，宇宙的总能量等。

非本征灰数是指凭先验信息或某种手段，可以找到一个白数作为其代表的灰数。

我们称此白数为相应灰数的白化值。

第二章序列算子与灰色序列生成

灰色系统理论的主要任务之一，是根据社会，经济，生态等系统的行为特征数

据，寻找不同系统变量之间或某些系统变

量自身的数学关系和变化规律。灰色系统理论认为任何随机过程都是在一定幅值

范围和一定时区内变化的灰色量，并把随机过程看成灰色过程。

灰色系统理论是通过对原始数据的挖掘，整理来寻求其变化规律的，这是一种就

数据寻找数据的现实规律的途径，我们称为灰色序列生成。灰色系统理论认为，

尽管客观系统表象复杂，数据离乱，但它总是有整体功能的，因此必然蕴含某种

内在规律。关键在于如何选择适当的方式去挖掘它和利用它。一切灰色序列都能

通过某种生成弱化其随机性，显现其规律性。

例如考虑4个数据，记为X(0)(1),X(0)(2),X(0)(3),X(0)(4)，其数据见下表：

将上表数据作图得

上图表明原始数据X(0)没有明显的规律性，其发展态势是

摆动的。如果将原始数据作累加生成，记第K个累加生成为X(1)(K),并且

X(1)(1)=X(0)(1)=1

X(1)(2)=X(0)(1)+X(0)(2)=1+2=3

X(1)(3)=X(0)(1)+X(0)(2)+X(0)(3)=1+2+1.5=4.5

X(1)(4)=X(0)(1)+X(0)(2)+X(0)(3)+X(0)(4)=1+2+1.5+3=7.5

得到数据如下表所示

上图表明生成数列X(1)是单调递增数列。

2.1冲击扰动系统与序列算子

定义2.1.1设

X0=(x0(1),x0(2),,x0(n))为系统真实行为序列，而观察到的系统行为数据序列为

X=(x(1),x(2),,x(n))=(x0(1)+ε1,x0(2)+ε2,,x0(n)+εn)=X0+ε

其中，ε=(ε1,ε2εn)为冲击扰动项（干扰项）。X称为动序列。

所以本章我们的讨论围绕：由X→X0展开（扰动还原真实）

2.2缓冲算子公理

定义2.2.1设系统行为数据序列为X

1.若∀k=2,3,n,x(k)-x(k-1)>0，则称=(x(1),x(2),,x(n))，X为单调增长序列；

2.若1中不等号反过来成立，则称X为单调衰减序列；

3.若∃k,k'∈{2,3,n},有x(k)-x(k-1)>0,x(k')-x(k'-1)<0，则称X为随机振荡序列。

4.设M=max{x(k)|k=12,3,，,n},m={x(k)|k=12,3,，,n}，则称M-m为序列X的振幅

定义2.2.2设X=(x(1),x(2),,x(n))为系统行为数据系列，D为作用于X的算子，X

经过算子D作用后所得序列记为

XD=(x(1)d,x(2)d,,x(n)d)

称D为序列算子，称XD为一阶算子作用序列。

序列算子的作用可以多次，相应的，若D1,D2,D3都是序列算子，

我们称D1D2为二阶算子，并称

XD1D2=(x(1)d1d2,x(2)d1d2,,x(n)d1d2)

为二阶算子作用序列，同理，D1D2D3为三阶序列算子„„

定义2.2.3称下述三公理为缓冲算子三公理，满足缓冲算子三公理的序列算子D

称为缓冲算子，一阶，二阶，三阶„„缓冲算子作用序列称为一阶，二阶，三阶„„

缓冲序列。

公理1（不动点公理）设X=(x(1),x(2),,x(n))为系统行为

数据系列，D为序列算子，则D满足x(n)d=x(n)。

不动点公理限定在序列算子作用下，系统行为数据序列的数据x(n)保持不变。

根据定性分析的结论，亦可使x(n)以后的若干个数据在序列算子作用下保持不

变。例如，令

x(j)d≠x(j)且x(i)d=x(i)

其中,j=1,2,k-1i=k,k+1,,n.

公理2.（信息充分利用公理）系统行为数据序列X中的每一个数据x(k),k=1,2,，

都要充分地参与算子的作用全过程

公理3（解析化、规范化公理）任意的

皆可由一个统一的x(1),x(2),x(k)d,(=k1,2，,,x(n)的初等解析式表达。

定义2.2.4设X为原始数据序列，D为缓冲算子，当X分别为增长序列，衰减序

列或振荡序列时：

1.若缓冲序列XD比原始序列X的增长速度（或衰减速度）减缓或振幅减小，则

称缓冲算子D为弱化算子。

2.若缓冲序列XD比原始序列X的增长速度（或衰减速度）加快或振幅增大，则

称缓冲算子D为强化算子。

2.3实用缓冲算子的构造

定理2.3.1设原始数据序列X=(x(1),x(2),,x(n))令缓

(1)d,x(2d),x,n(d冲序列XD=(x

其中x(k)d=1[x(k)+x(k+1)+n-k+1+x(n)]；k=1,2,„„，n,则当X为增长序列，衰减序

列或振荡序列时，D为弱化算子，并称为平均

弱化缓冲算子（AWBO）

证明：直接利用x(k)d,(k=1,2,)的定义，可知定理成立。推论2.3.1对于定理1中

定义的弱化算子D，

令XD2=XDD=(x(1)d2,x(2)d2,

x(k)d2=1[x(k)d+x(k+1)d+n-k+1,x(n)d2)+x(n)d],k=1,2n，

则D2对于增长序列，衰减序列或振荡序列时，皆为二阶弱化算子。

定理2.3.2设原始序列和其缓冲算子序列分别为

X=(x(1),x(2),,x(n))

XD=(x(1)d,x(2)d,

其中x(k)d=

x(1)+x(2)+,x(n)d),k=1,2,n-1+x(k-1)+kx(k)2k-1x(n)d=x(n)

则当X为增长序列（越来越大），衰减序列或振荡序列时，D为强化算子。

推论2.3.2设D为定理2中定义的强化算子，令

2=XD=D((1x)2d,XD2(x2)d,2,x，其中(n)d)

x(n)d2=x(n)d=x(n)

x(k)d2=x(1)d+x(2)d++x(k-1)d+kx(k)d

2k-1,k=1,2,n-1，

则D2对于增长序列，衰减序列或振荡序列皆为二阶强化算子。

定理2.3.3原始数据序列和其缓冲算子序列分别为

X=(x(1),x(2),,x(n))

XD=(x(1)d,x(2)d,,x(n)d)

其中x(k)d=kx(k)+(k+1)x(k+1)++nx(n)

(n+k)(n-k+1)/2,k=1,2则当,n，X为增长

序列，衰减序列或振荡序列时，D为弱化算子，并称D为加权平均弱化缓冲算子

（WAWBO）

定理2.3.4设X=(x(1),x(2),,x(n))为非负的系统行为数据序列，令XD=(x(1)d,x(2)d,

其中x(k)d=[x(k)x(k+1)x(n)]1

n-k+1n,x(n)d)1n-k+1=[∏x(i)]

i=k,k=1,2,n。

则当X为增长序列，衰减序列或振荡序列时，D为弱化缓冲算子，并称D为几何

平均弱化缓冲算子（GAWBO）

定理2.3.5设X=(x(1),x(2),,x(n))为系统行为数据序列，各时点的权重向量为

ω=(ω1,ω2

XD=(x(1)d,x(2)d,,x(n)d)ωn)，则

其中x(k)d=ωkx(k)+ωk+1x(k+1)++ωnx(n)

ωk+ωk+1++ωn,k=1,2,n。则当XD皆

为弱化缓冲算子，并称D为加权平均弱化缓冲算子（WAWBO）。

定理2.3.6设X

量为ω=(ω1,ω2=(x(1),x(2),,x(n))，各时点的权重向ωn)>0，

,x(n)d)令XD=(x(1)d,x(2)d,

其中

1

x(k)d=[x(k)xωkωk+1(k+1)x(n)]ωnk+k+1++n=[∏x(i)]

i=kn1k+k+1++n,k=1,2,n

则当XD为弱缓冲算子，并称D为加权几何平均弱化缓冲算子（WGAWBO）。

定理2.3.7设X=(x(1),x(2),,x(n))为系统行为数据序

列，令XD=(x(1)d,x(2)d,,x(n)d),k=1,2,n。(n-k+1)x2(k)其中

x(k)d=x(k)+x(k+1)++x(n)

则当X为增长序列，衰减序列或振荡序列时，D为强化缓冲算子，并称D为平均

强化缓冲算子（ASBO）

定理2.3.8设X=(x(1),x(2),,x(n))为非负的系统行为数据序列，令XD=(x(1)d,x(2)d,

其中x(k)d=x2(k)

[x(k)x(k+1)x(n)]1

n-k+1,x(n)d)x2(k)n1n-k+1=,k=1,2,n。[∏x(i)]

i=k

则当X为增长序列，衰减序列或振荡序列时，D为强化缓冲算子，并称D为几何

平均强化缓冲算子（GASBO）

以上列举了部分缓冲算子，当然，我们还可以考虑构造其它形式的实用缓冲算

子，缓冲算子不仅可以用于灰色系统建模，而且还可以用于其它各种模型建模。

通常在建模之前根据定性分析结论对原始数据序列施以缓冲算子，淡化或消除冲

击扰动对系统行为数据序列的影响，往往会收到预期的效果。

[例2.3.1]河南省长葛县乡镇企业产值数据（1983年-1986年）为

X=(10155,12588,23480,35388)

其增长势头很猛，1983-1986年每年平均递增51.6%，尤其是1984-1986年，每年

平均递增67.7%。因此普遍认为今后不可能一直保持这么高的发展速度。经过认

真分析，大家认识到增长速度高主要是基数低，而基数低的原因是过去对有利与

乡镇企业发展的政策没有用足，用活，用好。要弱化序列增长趋势，就要将

乡镇企业发展比较有利的现行政策因素附加到过去的年份中去，为此，引进推论

1所示的二阶弱化算子，得二阶缓冲序列XD2=(27260,29547,32411,35388)

用XD2建模预测得，1986-2000年该县乡镇企业每年平均递增

9.4%，这一结果是1987年得到的，与“八五”后半期和“九五”期间该县乡镇企业发

展实际基本吻合。

2.4均值生成算子

在收集数据时，常常由于一些不易克服的困难导致数据序列出现空缺（也称空

穴），有些数据序列虽然完整，但由于系统行为在某个时点上发生突变而形成异

常数据，剔除异常数据就会留下空穴，如何填补空穴，自然成为数据处理过程中

首先遇到的问题，均值生成是常用的构造新数据，填补原序列空穴，生成新序列

的方法。

定义2.4.1设序列X在k出现有空穴，记为∅(k)，即

X=(x(1),x(2),,x(k-1),∅(k),x(k+1),,x(n))则称x(k-1)和x(k+1)为∅(k)的界值，x(k-1)为

前界，x(k+1)为后界

当∅(k)是由x(k-1)和x(k+1)生成时，称生成值x(k)为[x(k-1)，x(k+1)]的内点

定义2.4.2设序列

X=(x(1),x(2),,x(k-1),∅(k),x(k+1),,x(n))为k处有空穴∅(k)的序列，而∅(k)

=x*(k)=0.5x(k-1)+0.5x(k)称为非紧邻均值生成数，所得序列称为非紧邻生成序

列。

定义2.4.3设序列X=(x(1),x(2),,x(n))，若

x*(k)=0.5x(k-1)+0.5x(k)，则称x*(k)为紧邻生成数，由紧邻生

成数构成的序列称为紧邻均值生成序列。

2.5序列的光滑性

定义2.5.1设序列X

是X的均值生成序列：

Z=(z(1),z(2),=(x(1),x(2),,x(n),x(n+1))，Z,z(n))，其中z(k)=0.5x(k-1)+0.5x(k)，X*是

某一可导函数的代表序列，d为n维空间的距离函数，我们将X删去x(n+1)后所

得的序列仍记X，若X满足

1.当k充分大时，x(k)<∑x(i)

i=1k-1

x*(k)-x(k)≥maxx*(k)-z(k)1≤k≤n1≤k≤n

则称Ｘ为光滑序列，１，２为序列光滑条件。

定义２.５.２称ρ(k)=x(k)

∑x(i)

i=1k-1;k=2,3,n

为序列Ｘ的光滑比。

定义２.５.３若序列Ｘ满足

１.ρ(k+1)<1;ρ(k)k=2,3,n-1

２.ρ(k)∈[0,ε];

３.ε<0.5

则称Ｘ为准光滑序列。

２.６级比生成算子k=3,4,n

定义２.６.１设序列X=(x(1),x(2),,x(n))，则称

σ(k)=x(k);x(k-1)k=2,3,n

为序列Ｘ的级比。

２.７累计生成算子与累减生成算子

累加生成是使灰色过程由灰变白的一种方法，它在灰色系统理论中占有极其重要

的地位。通过累加可以看出灰量积累过程的发展态势，使离乱的原始数据中蕴含

的积分特性或规律充分显露出来。

定义２.７.１设X0=(x0(1),x0(2),,x0(n))，Ｄ为序列算子

X0D=(x0(1)d,x0(2)d,

x(k)d=∑x0(i);0

i=1k,x0(n)d)，其中,n。k=1,2,3,

则称Ｄ为X0的一次累加生成算子，记为１-ＡＧＯ（AccumulatingGeneration

Operator）,称r阶算子Dr为X0的r次累加生成算子，记为r-AGO,习惯上，我们

记

X0D=X1=(x1(1),x1(2),,x1(n))

X0Dr=Xr=(xr(1),xr(2),,xr(n))

其中x(k)d=∑xr-1(i);r

i=1kk=1,2,3,,n

定义2.7.2（累减）设X0=(x0(1),x0(2),

算X0D=(x0(1)d,x0(2)d,,x0(n)d)，其中

x0(k)d=x0(k)-x0(k-1)k=1,2,3,,n,x0(n))，D为序列

2.8灰指数律

定义2.8.1设序列X

1.x(k)=ceak;

2.x(k)=ceak;=(x(1),x(2),,x(n))，若对于c,a≠0;k=1,2c,a,b≠0;k=1,2n则称X为齐次指

数序列。X为齐次指数序列。

,x(n))若n，称定义2.8.2设序列X

1.∀k,σ(k)=

律。2.

律。3.∀k,σ(k)=∀k,σ(k)==(x(1),x(2),x(k)∈(0,1]，则称序列x(k-1)X具有负的灰指

数规x(k)∈(1,b]，则称序列x(k-1)X具有正的灰指数规x(k)∈[a,b],b-a=δx(k-1)则称

序列X具有绝对灰度为δ的灰指数规律。

4.δ

定理2.8.1设序列X0=(x0(1),x0(2),,x0(n))为非负准光滑序<0.5时，称X具有准指

数规律。列，则X0的一次累加生成序列X1具有准指数规律。注：定理2.8.1是

灰色系统建模的理论基础

第三章灰色关联分析

对两个系统或两个因素之间关联性大小的量度，称为关联度。它描述系统发展过

程中因素间相对变化的情况，也就是变化大小、方向及速度等指标的相对性。如

果两者在系统发展过

程中相对变化基本一致，则认为两者关联度大；反之，两者关联度就小。灰色系

统理论的关联度分析与数理统计学的相关分析是不同的，两者的区别在于第一，

它们的理论基础不同。关联度分析基于灰色系统的灰色过程，而相关分析则基于

概率论的随机过程；第二，分析方法不同。关联分析是进行因素间时间序列的比

较，而相关分析是因素间数组的比较；第三，数据量要求不同。关联分析不要求

数据太多，而相关分析则需有足够的数据量；第四，研究重点不同。关联度分析

主要研究动态过程，而相关分析则以静态研究为主。因此，关联度分析适应性更

广，在用于社会经济系统中的应用更有其独到之处。

一般的抽象系统，如社会系统，经济系统，农业系统，生态系统等都包含有许多

种因素，多种因素共同作用的结果决定了该系统的发展态势。我们常常希望知道

众多的因素中，哪些是主要因素，哪些是次要因素，哪些因素对系统发展影响

大，哪些因素对系统发展影响小，哪些因素对系统发展起推动作用需加强，哪些

因素对系统发展起阻碍作用需抑制„„

数理统计中的回归分析，方差分析，主成分分析等都是用来进行系统特征分析的

方法。但数理统计中的分析方法往往需要大量数据样本，且服从某个典型分布。

灰色关联分析方法弥补了采用数理统计方法作系统分析所导致的缺憾.它对样本量

的多少和

样本有无规律都同样适用,而且计算量小,十分方便,更不会出现量化结果与定性分

析结果不符的情况。

灰色关联分析的基本思想是根据序列曲线几何形状的相似程度来判断其联系是否

紧密。曲线越接近，相应序列之间关联度就越大，反之就越小。例如某地区农业

总产值X0,种植业总产值

畜牧业总产值X2和林业总产值X3，从X1，

统计数据如下：

X0=（18，20，22，35，41，46）

X1=（8，11，12，17，24，29）

X2=（3，2，7，4，11，6）

X0=（5，7，7，11，5，10）

1997-2002年共6年的

从直观上看，与农业总产值曲线最相似的是种植业总产值曲线，而畜牧业总产值

曲线和林果业总产值去与农业总产值曲线在几何形状上差别较大。因此我们可以

说该地区的农业仍然是以种植业为主的农业，畜牧业和林果业还不够发达。

3.1灰色关联因素和关联算子集

进行系统分析，选准系统行为特征的映射量后，还需进一步明确影响系统行为的

有效因素。如要作量化研究分析，则需要对系统行为特征映射量和各有效因素进

行处理，通过算子作用，使之化为数量级大体相近的无量纲数据，并将负相关因

素转化为正相关因素。

定义3.1.1设Xi=(xi(1),xi(2),序列，D1为序列算子，且

XiD1)d1,xi(2)d1,1=(xi(

,xi(n)d1)

n

,xi(n))为因素Xi的行为

其中xi(k)d1=

xi(k)

xi(1)

xi(1)≠0;k=1,2

，则称D1为

初值化算子。XiD1为Xi在初值化算子D1的象，简称初值象。

定义3.1.2设Xi=(xi(1),xi(2),序列，D2为序列算子，且

XiD2=(xi(1)d2,xi(2)d2,

,xi(n)d2)

n

,xi(n))为因素Xi的行为

其中xi(k)d2

=

xi(k)1n

∑x(k)

i

k=1

n

;k=1,2

，则称D2为均

值化算子。XiD2为Xi在均值化算子D2的象，简称均值象。

定义3.1.3设Xi=(xi(1),xi(2),序列，D3为序列算子，且

XiD3=(xi(1)d3,xi(2)d3,

,xi(n)d3)

;k=1,2

n

,xi(n))为因素Xi的行为

其中xi(k)d3

=

xi(k)-minxi(k)

k

maxxi(k)-minxi(k)

k

k

，

则称D3为区间化算子。XiD3为Xi在区间化算子D3的象，简称区间

值象。

命题4.1.1初值化算子、均值化算子和区间值化算子皆可以使系

统行为序列无量纲化，且在数量上规一。一般地，不宜混合、重叠使用。

定义3.1.4设Xi=(xi(1),xi(2),

素Xi的行为序列，D4为序列算子，且

XiD4=(xi(1)d4,xi(2)d4,,xi(n)d4)

n,xi(n))，x(k)∈[0,1]为因其中xi(k)d4=1-xi(k);k=1,2，则称D4为逆化

算子。XiD4为Xi在逆化算子D4的象，简称逆化象。

定义3.1.5设Xi=(xi(1),xi(2),

素Xi的行为序列，D5为序列算子，且

XiD5=(xi(1)d5,xi(2)d5,,xi(n)d5),xi(n))，x(k)∈[0,1]为因n其中xi(k)d5=1

xi(k)xi(k)≠0;k=1,2，则称

D5为倒数算子。XiD5为Xi在倒数化算子D5的象，简称倒数化象。

命题4.1.3若系统因素Xi=(xi(1),xi(2),

为呈负相关关系，则Xi=(xi(1),xi(2),,xi(n))与系统主行,xi(n))的逆化算子作

用像和倒数化作用像与X0具有正相关关系。

定义3.1.6称D={Di|i=1,2,3,4,5}为灰色关联算子集（以上五个）。

定义3.1.7设X为系统因素集合，D为灰色关联算子集，称（X,D）为灰色关联因

子空间。

3.2灰色关联公理与灰色关联度

灰色系统理论提出了一种新的分析方法—关联度分析方法，即根据因素之间发展

态势的相似或相异程度来衡量因素间关联的程度，它揭示了事物动态关联的特征

与程度。由于以发展态势为立足点，因此对样本量的多少没有过分的要求，也不

需要典型的分布规律，计算量少到甚至可用手算，且不致出现关联度的量化结果

与定性分析不一致的情况。这种方法已应用到农业经济、水利、宏观经济等各方

面，都取得了较好的效果。

灰色系统理论建模的主要任务是根据具体灰色系统的行为特征数据，充分开发并

利用不多的数据中的显信息和隐信息，寻找因素间或因素本身的数学关系。通常

的办法是采用离散模型，建立一个按时间作逐段分析的模型。但是，离散模型只

能对客观系统的发展做短期分析，适应不了从现在起做较长远的分析、规划、决

策的要求。尽管连续系统的离散近似模型对许多工程应用来讲是有用的，但在某

些研究领域中，人们却常常希望使用微分方程模型。事实上，微分方程的系统描

述了我们所希望辨识的系

统内部的物理或化学过程的本质。

大千世界里的客观事物往往现象复杂，因素繁多。我们往往需要对系统进行因素

分析，这些因素中哪些对系统来讲是主要的，哪些是次要的，哪些需要发展，哪

些需要抑制，哪些是潜在的，哪些是明显的。一般来讲，这些都是我们极为关心

的问题。事实上，因素间关联性如何、关联程度如何量化等问题是系统分析的关

键和起点。

对于两系统之间的因素，其随时间或不同对象而变化的关联性的大小的量度，称

为关联度。在系统发展过程中，若两个因素变化的趋势具有一致性，即变化程度

较高，即可谓二者的关联度较高；反之，则较低。因此，灰色关联度分析方法，

是根据因素之间发展趋势的相似或相异程度，即“灰色关联度”作为衡量因素之间

关联程度的一种方法。灰色系统理论提出了对各子系统进行灰色关联度分析的概

念，意图透过一定方法，去寻求系统各子系统（或因素）之间数值的关系。因

此，灰色关联度分析对于一个系统的发展变化态势提供了量化的度量，非常适合

动态历程分析。

关联分析实际上是动态过程发展态势的量化比较分析。

所谓发展态势比较，也就是系统各时期有关统计数据的几何关系的比较。

例如，某地区1977～1983年总收入与养猪、养兔收入资料见表1。

定义3.2.1（灰色关联公理）

设X0=(x0(1),x0(2),

X1=(x1(1),x1(2),,x0(n))为，且,x1(n))

„„„„„„„„„„„„„„

Xi=(xi(1),xi(2),,xi(n))

„„„„„„„„„„„„„„

Xm=(xm(1),xm(2),,xm(n))

为相关因素序列，给定实数r(x0(k),xi(k))，若实数

1n

r(X0,Xi)=∑r(x0(k),xi(k))，满足nk=1

1.规范性0

2.整体性对于Xi,Xj∈X={XS|s=0,1,,m;m≥2}有

r(Xi,Xj)≠r(Xj,Xi)i≠j

3.偶对对称性Xi,Xj∈X，有

r(Xi,Xj)=r(Xj,Xi)⇔X={Xi,Xj}

4.接近性|x0(k-)ix(k越小，)|r(x0(k),xi(k))越大。1n则称r(X0,Xi)=∑r(x0(k),xi(k))为

Xi,Xj∈X的灰色关联度，其中nk=1

r(x0(k),xi(k))为Xi和Xj在k点的关联系数，并称条件1.2.3.4为灰色关联四公

理。

在灰色关联公理中，规范性0

严格无关联。整体性则体现了环境对灰色关联比较的影响，环境不同，灰色关联

度也随之变化，因此对称性不一定满足。

偶对对称性表明，当灰色关联因子集中只有两个序列时，满足对称性。

接近性是对关联量化的约束。

定义3.2.2设系统行为序列

X0=(x0(1),x0(2),,x0(n))

X1=(x1(1),x1(2),,x1(n))

„„„„„„„„„„„„„„

Xi=(xi(1),xi(2),,xi(n))

„„„„„„„„„„„„„„

Xm=(xm(1),xm(2),,xm(n))

对于ξ∈(0,1)令

r(x0(k),xi(k))=minminx0(k)-xi(k)+ξ⋅maxmaxx0(k)-xi(k)ikik

x0(k)-xi(k)+ξ⋅maxmaxx0(k)-xi(k)ik

记r(x0(k),xi(k))为r0i(k)，1n1n

r(X0,Xi)=∑r(x0(k),xi(k))=∑r0i(k)nk=1nk=1，则

1n

r(X0,Xi)=∑r(x0(k),xi(k))满足灰色关联公理，其中ξnk=1称为分辨系

数。r(X0,Xi)称为X0,Xi的灰色关联度，记为r0i。根据关联度的定义，可得关联

度的计算步骤如下：

1．根据评价目的确定评价指标体系，收集评价数据

设m个数据序列形成如下矩阵：

(X0,X1⎛x0(1)x1(1)x0(2)x1(2),Xm)=x(n)x(n)1⎝0xm(1)⎫⎪xm(2)⎪⎪⎪xm(n)⎪⎭其

中n为指标的个数，

Xi=(xi(1),xi(2),,xi(n)),i=1,2,T,m．

m是比较数列，n是数据指标,X0是参考数据列

2．确定参考数据列X0

参考数据列应该是一个理想的比较标准，可以以各指标的最优值（或最劣值）构

成参考数据列，也可根据评价目的选择其它参照值．记作

X0=(x0(1),x0(2),,x0(m))

3．对指标数据序列用关联算子进行无量纲化（目的是消除数量级大小不同的影

响，以便于进行计算和比较分析），无量纲化后的数据序列形成如下矩阵：

(X0',X1',

⎛x'(1)x'(1)1

x'(2)x'(2)1

')=0,Xm

x'(n)x'(n)1⎝0xm'(1)⎫

⎪xm'(2)⎪

⎪⎪⎪'xm(n)⎪⎭

常用的无量纲化方法有均值化像法、初值化像法等．

1n

xi(k)∑nk=1

i=0,1,,m；k=1,2,

xi'(k)=

xi(k)

,xi'(k)=

,n.

xi(k)xi1

由于系统中各因素的量纲(或单位)不一定相同，如劳动力为人，

产值为万元，产量为吨等，且有时数值的数量级相差悬殊，如人均收入为几百

元，粮食每公顷产量为几千公斤，费用为几十万元，有些产业产值达百亿元，有

些产业才几万元，等等，这样的数据很难直接进行比较，且它们的几何曲线比例

也不同。因此，对原始数据需要消除量纲(或单位)，转换为可比较的数据序列。

目前，原始数据的变换有以下几种常用方法：

a)均值化变换。先分别求出各个序列的平均值，再用平均值去除对应序列中的各

个原始数据，所得到新的数据列，即为均值化序列。其特点是量纲为一，其值大

于0，并且大部分近于1，数列曲线互相相交。

b)初值化变换。分别用同一序列的第一个数据去除后面的各个原始数据，得到新

的倍数数列，即为初值化数列。量纲为一，各值均大于0，且数列有共同的起

点。

c)标准化变换。先分别求出各个序列的平均值和标准差，然后将各个原始数据减

去平均值后再除以标准差，这样得到的新数据序列即为标准化序列。量纲为一，

其均值为0，方差为1。

一般情况下，对于较稳定的社会经济系统数列作动态序列的关联度分析时，多采

用初值化变换，因为这样的数列多数是增长的趋势。若对原始数列只作数值间的

关联比较，可用均值化变换，譬如进行产业结构变化的关联分析，自然因素周期

性变化的关联分析等。

补充

4．逐个计算每个被评价对象指标序列与参考序列对应元素的绝对差值

即∆i(k)=x0'(k)-xi'(k)；k

nm=1,,ni=1,,mminx0'(k)-xi'(k)与5．确定M=mini=1k=1

m=maxmaxx0'(k)-xi'(k)i=1k=1nm

6．计算关联系数

分别计算每个比较序列与参考序列对应元素的关联系数r(x0'(k),xi'(k))=

k=1,m+ξ⋅M∆i(k)+ξ⋅M,n

式中ξ为分辨系数，在（0，1）内取值，ξ越小，关联系数间的差异越大，区分

能力越强．通常ξ取0.5．

ρ称为分辨系数，其意义是削弱最大绝对差数值太大引起的失真，提高关联系数

之间的差异显著性，ρ∈(0,1)，一般情况下可取0.1~0.5。

关联系数反映两个被比较序列在某一时刻的紧密(靠近)程度。如在∆min的时刻,

Lio＝1，而在∆max的时刻则关联系数为最小值。因此，关联系数的范围为0

≤1。

7.计算关联度

1n

r(X0,Xi)==∑r0i(k)nk=1

由以上所述可知，关联度分析实质上是对时间序列数据进行几何关系比较，若两

序列在各个时刻点都重合在一起，即关联系数均等于1，则两序列的关联度也必

等于1。另一方面，两比较序列在任何时刻也不可垂直，所以关联系数均大于

0，故关联度也都大于0。因此，两序列的关联度便以两比较序列各个时刻的关联

1n系数之平均值计算，即：r(X0,Xi)==∑r0i(k)nk=1

式中r0i为子序列i与母序列0的关联度，N为比较序列的长度(即数据个数)。

用几何坐标表示，即在横坐标为时间t、纵坐标为关联系数L的坐标图中，绘出

关联系数曲线(虚线)。该折线与横坐标间围成的面积，称为关联面积，记作S0i，

而母序列自身的关联系数处处为1。所以，取纵坐标L=1，作水平线与横坐标间

围成的面积为重合面积，记为S00，则关联度的几何意义为两面积之比，即

r0i=S0i

S00。因关联系数曲线为等时距，且S00＝1，故

1r0i=S0i=Nk=1∑Loi(k)N

。

8．依据各观察对象的关联

序，得出综合评价结果．

因素有关：

1)母序列X0不同，则关联度不同；

2)子序列Xi不同，则关联度不同；

3)参考点0(或数据变换)不同，关联度不同；

4)数据序列长度N不同，关联度不同；

5)分辨系数ρ不同，关联度不同。

一般来说，关联度也满足等价“关系”三公理，即：1)自反性:r00＝1；2)对称性:

r0i＝ri0；3)传递性:r0a>r0b,r0b>r0c，则r0a>r0c。

将m个子序列对同一母序列的关联度按大小顺序排列起来，便组成关联序，记为

{X}。它直接反映各个子序列对于母序列的“优劣”关系。若r0a>r0b，则称{Xa}

对于相同母序列{X0}有优于{Xb}的特点，记为{Xa不难看出，关联度与下列|X0}

{Xb|X0}；若

ar0a

{Xb|X0}；若

{Xa}对于母序列{Xo}等价于（或等于）{Xb}，记为{Xa│X0}~{Xb│X0}；若有

r0a≥r0b，称{Xa}对于母序列{Xo}优于或等于{Xb}，

{Xb|X0}记为{Xa|X0}；若有r0a≤r0b，则称{Xa}对于母序列{Xo}

劣~

{Xb|X0}

于或等于{Xb}，记为{Xa|X0}。~

根据上述几种关系，可定义两种有代表性的关联序，即“有

序”与“偏序”。若关联序{X}为有序，那么所有元素之间必存在以下几种关系之

一：“优于”(),“劣于”()，或“等价于”(~)。若关联序{X}为偏序，则不是所有元素

都可比较的。

一般而言，各因素只要能构成关系，算出关联度，则总是“有序”的。只有在无“参

考点”或无参考母序列的情况下，才可能出现“偏序”现象。

3.3灰色关联分析的应用举例

1、利用灰色关联分析对6位教师工作状况进行综合评价1．评价指标包括：专业

素质、外语水平、教学工作量、科研成果、论文、著作与出勤．

2．对原始数据经处理后得到以下数值，见下表

教师考评数据表

3．确定参考数据列：

X0={9,

9,9,9,8,9,

9}

''4．计算∆i(k)=x0(k)-xi(k)，见下表

n

m

minx0(k)-xi(k)=min(0,1,0,1,0,0)=05．mini=1k=1

maxmaxx0(k)-xi(k)=max(7,6,5,6,6,5)=7

i=1

k=1

nm

''6．依据式r(x0(k),xi(k))=

m+ξ⋅M

∆i(k)+ξ⋅M

，取ξ=0.5计算，

得

r01(1)=0+0.5⨯7=0.7781+0.5⨯7r01(2)=0+0.5⨯7=1.0000+0.5⨯7r01(3)=0.778

r01(4)=0.636r01(5)=0.467r01(6)=0.333r01(7)=1.000

同理得出其它各值，见表（12—5）

7．分别计算每个人各指标关联系数的均值（关联度）：

r01=0.778+1.000+0.778+0.636+0.467+0.333+1.000=0.7137

同理

r02=0.614，r03=0.680，r04=0.599，r05=0.683，r06=0.658，

8．如果不考虑各指标权重（认为各指标同等重要），六个被评价对象由好到劣

依次为1号（r01=0.713），5号（r05=0.683），

3号（r03=0.680），6号（r06=0.658），2号（r02=0.614），4号

（r04=0.599）．即

r01>r05>r03>r06>r02>r04Matlab程序:

clc;clear;

formatshortg

A=[8987529

7875738

9796647

6888436

8669838

8957648

];%已经处理后的数据x0=[9999899];

fori=1:size(A,1)

forj=1:size(A,2)

t(i,j)=abs(A(i,j)-x0(j));end

end

Min=min(min(t));Max=max(max(t));rho=0.5;

r=(Min+rho*Max)./(t+rho*Max);rend=sum(r')/size(r,2);

[rs,rind]=sort(rend,'descend')

2、补例

3、补例

建模步一：数据处理(可以不用无量纲化)

x0(k)={135048,143875,...,274618}

x1(k)={98855.136,105028.75,...,212279.714}

x2(k)={23228.256,23451.625,...,27187.182}

x3(k)={3646.296,4028.5,...,11259.338}

x4(k)={9318.312,11366.125,...,23891.766}

步二：计算关联度。经计算

minminx0(k)-xi(k)=36192.864

maxmaxx0(k)-xi(k)=263358.662

εi(k)=minminx0(k)-xi(k)+ρmaxmaxx0(k)-xi(k)

x0(k)-xi(k)+ρmaxmaxx0(k)-xi(k)=

=36192.864-131679.331x0(k)-xi(k)+131679.331

设分辨系数为0.5

将相应x0(k)与xi(k)的数值代入式εi(k)=

ε1=∆min+∆max∆+ρ∆max中，得，(1,09844,0.9783

0.973,0.95,0.932,0.916,0.899,0.874,0.865)

0.666,0.649,0.601,0.554,0.521,0.497,0.477,ε2=(0.689,

0.46,0.443)

0.618,0.604,0.562,0.52,0.492,0.472,0.454,ε3=(0.638,

0.44,0.425)

0.635,0.62,0.576,0.535,0.506,0.485,0.467,ε4=(0.652,

0.454,0.439)

110步三：算出关联度R(x0,xi)，由公式Ri=∑εi(k)分别计算出原煤、ni=1

原油、天然气及水核风电关于能源总产量的关联度R1,R2,R3,R4。

110

R1=R(x0,xi)=∑ε1(k)=0.937110i=1

110

R2=R(x0,xi)=∑ε2(k)=0.555710i=1

110

R3=R(x0,xi)=∑ε3(k)=0.522410i=1

110

R4=R(x0,xi)=∑ε4(k)=0.5369ni=1

步四：比较关联度大小得出结论。由R1>R2>R4>R3说明原煤与能

源总量关系最密切，而原油、水核风电和天然气的密切程度依次

较差。4、

5、

Matlab程序

functionoutput=grayrela(x0)

%参考因子与比较因子共同存储在一个矩阵x0中,参考因子位于第一列

%斜率序列

fori=2:length(x0(:,1))

x1(i,:)=x0(i,:)-x0(i-1,:);

end

%标准化

m=length(x1(1,:));

fori=1:m

x2(:,i)=x1(:,i)/std(x1(:,i));

end

%排序

[y,pos]=sort(x2(:,1));

x2_sorted=x2(pos,:);

%判定关联性质

n=length(x1(:,1));

k=[1:n]';

forj=1:m

sig_j(j)=qiuhe(k.*x2_sorted(:,j))-qiuhe(x2_sorted(:,j))*qiuhe(k)/n;

end

%caculationofdistantion

forj=2:m

dist_0i(:,j)=abs(sign(sig_j(:,j)./sig_j(:,1)).*x2_sorted(:,j)-x2_sorted(:,

1));

end

%计算关联系数

fori=1:n

forj=1:m

coef_rela(i,j)=(min(dist_0i)+0.5*max(dist_0i))/(dist_0i(i,j)+0.5*max(dist_0i));

end

end

forj=1:m

output(j)=qiuhe(coef_rela(:,j))/n;

end

其中：

functionoutput=qiuhe(input)

output=0;

fori=1:length(input)

output=output+input(i);end

3.4广义灰色关联度

命题3.4.1设X0=（x0(1),x0(2),,x0(n)）

，而Xi=（xi(1),xi(2),,xi(n)）

X00=（x00(1),x00(2),,x00(n)）和Xi0=（xi0(1),xi0(2),,xi0(n)）

分别为X0与Xi的始点化像，即x00(k)=x0(k)-x0(1)，xi0(k)=xi(k)-xi(1)，则记

1

|s0|=|∑x0(k)+00(n)|，

2k=2

0n-1n-1

|si|=|∑xi0(k)+1i0(n)|及

k=2

n-1

2

1

|si-s0|=|∑(xi0(k)-x00(k))+i0(n)-x00(n))|

2k=2

定义3.4.1设序列，s0，si如命题3.4.1中所示，则称

ε0i=

1+|s0|+|si|

1+|s0|+|si|+|si-s0|

为X0与Xi的灰色绝对关联度，简称绝对关联度。绝对关联度满足灰色关联公理

中规范性，偶对对称性与接近性，但不满足整体性。

定理3.4.1灰色绝对关联度ε0i具有如下的性质：1.0<ε0i≤1

2.ε0i只与X0与Xi几何形状有关，而与其它无关，或者说，平移不改变绝对关

联度的值。

3.任何两个序列都不是绝对无关的，即ε0i恒不为零。4.X0与Xi几何形状相似

程度越大，ε0i越大。

5.当X0或Xi中任一观测数据变化了，ε0i将随之变化。6.X0与Xi长度变化，

ε0i也变化。7.ε00=εii=18.ε0i=εi03.5灰色相对关联度

定义3.5.1设序列X0与Xi长度相同，且初值皆不等于零，X0'与Xi'分别为X0与

Xi的初值像，则称X0'与Xi'的灰色绝对关联度为

X0与X1的灰色相对关联度，简称为相对关联度，记为R0i。

相对关联度表征了序列X0与Xi相对与始点的变化速率之间的关系，X0与Xi的

变化速率越接近，R0i越大，反之越小。

定理3.5.1设序列X0与Xi长度相同，且初值皆不等于零，若

X0=cXi，其中c>0为常数，则R0i=1

定理3.5.2灰色相对关联度R0i具有如下的性质：1.0

2.R0i只与序列X0与Xi的相对于始点的变化率有关，而与各观测值的大小无

关，或者说，数乘不改变相对关联度的值。

3.任何两个序列的变化率都不是毫无联系的，即R0i恒不为零。

4.X0与Xi相当于始点的变化速度越接近，R0i越大。5.当X0或Xi中任一观测

数据变化了，R0i将随之变化。6.X0与Xi长度变化，R0i也变化。

7.R00=Rii=1

8.R0i=Ri0

3.6灰色综合关联度

定义3.6.1设序列设序列X0与Xi长度相同，且初值皆不等于

零，ε0i和R0i分别为X0与Xi的灰色绝对关联度和相对关联度，

θ∈[0,1]则称

ρ0i=θε0i+(1-θ)R0i

为X0与Xi的灰色综合关联度，简称综合关联度。

综合关联度既体现了折线X0与Xi的相似程度，又反映了X0与Xi相对于始点的

变化速率的接近程度，是较为全面的表征序列之间联系是否紧密的一个数量指

标。

例3.6.1河南省长葛县乡镇企业经济的灰色关联分析

20世纪80年代中期，长葛县乡镇企业发展比较快，1983年到1986年，平均每年

增长51.6%，乡镇企业经济在全县经济发展中占有重要地位，1986全县乡镇企业

产值35388万元，占工农业总产值的60%。如何有效的加快乡镇企业的发展，是

大家普遍关心的问题。据分析，乡镇企业产值主要与固定资产，流动资产，劳动

力，企业留利四个因素有关。该县乡镇企业产值及相关因素行为数据如下表：

单位：万元

1.求绝对关联度。令Xi0=(xi(1)-xi(1),xi(2)-xi(1),xi(3)-xi(1),xi(4)-

xi(1))=(x(1),x(2),x(3),x(4));i=0,1,2,3,40

i0i0i0i即始点零化

像，则

0X0=(0,2433,13325,25233)

X10=(0,2433,13325,25233)

0X2=(0,-194,1661,3188)

X=(0,408,461,3001)0

3

0X4=(0,624,2030,3314)

由|si|=|∑xi0(k)+

k=2310i(4)|;i=0,1,2,3,4得2

s0=28374.5,s1=3058.5,s2=2369.5,s3=85580,s4=4311

由|si-s0|=|∑(xi0(k)-x00(k))+1i0(4)-x00(4))|;i=1,2,3,4得

k=2

3

2

|s1-s0|=25316,|s2-s0|=26005,|s3-s0|=57205,|s4-s0|=24063.5

由ε0i=

1+|s0|+|si|

；i=1,2,3,4得

1+|s0|+|si|+|si-s0|

ε01=0.554,ε02=0.542,ε03=0.666,ε04=0.576

2.求相对关联度。先求出Xi的初值像，由

⎛xi(1)xi(2)xi(3)xi(4)⎫

Xi'=xi'(1),xi'(2)xi'(3)xi'(4)=x1,x1,x1,x1⎪⎪iii⎝i⎭i=0,1,2,3,4

()

得

X0'=(1,1.2396,2.305)1,3.4848

X1'=(1,0.9489,1.4372,1.8379)X2'=(1,1.2329,1.2631,2.7129)X3'=(1,1.8550,2.3851,3.53

68)X4'=(1,1.5361,2.6924,3.8471)

诸Xi'的始点零化像为

Xi'0=xi'(1)-xi'(1),xi'(2)-xi'(1),xi'(3)-xi'(1),xi'(4)-

xi'(1)=(xi'0(1),xi'0(2),xi'0(3),xi'0(4));i=0,1,2,3,4

(

)

从而有

'0=(0,0.2396,1.3051,2.4848)X0

X1'0=(0,-0.0511,0.4372,0.8379)'0=(0,0.2329,0.2631,1.7129)X2

'0=(0,0.8850,1.3851,2.5368)X3

'0=(0,0.5361,1.6924,2.8471)X4

'0(k)+由|si'|=|∑xi

k=23

1'0

i(4)|;i=0,1,2,3,42

s0'=2.7871,s1'=0.80505,s2'=1.35245,s3'=3.5385,s4'=3.65205由|si'-s0'|=|∑(xi'0(k)-

x0'0(k))+1i'0(4)-x0'0(4))|;i=1,2,3,4得2k=2

'-s0'|=1.98205,|s2'-s0'|=0.43465,|s3'-s0'|=0.7514,|s4'-s0'|=0.86495|s13

由r0i='|+|si'|1+|s0；i=1,2,3,4得'|+|si'|+|si'-s0'|1+|s0

r01=0.6985,r02=0.7818,r03=0.9070,r04=0.8958

3.求综合关联度。取θ=0.5，由综合关联度ρ0i=θε0i+(1-θ)R0i得

ρ01=0.6263,ρ02=0.6619,ρ03=0.7865,ρ04=0.7359

4.结果分析。由ρ03>ρ04>ρ02>ρ01得知

相对X0来说，X3为最优因素，X4次之，X2又次之，X1最差。也

就是说，劳动力对乡镇企业的产值影响最大，企业留利对产值的影响仅次于劳动

力，固定资产对产值的影响最小。这一结果与该县的实际情况吻合。这个县的乡

镇企业主要是劳动密集型产业，产值的增长在很大程度上是靠增加劳动力来实现

的。

第四章灰色系统模型

研究一个系统，一般应首先建立系统的数学模型，进而对系统的整体功能，协调

功能以及系统各因素之间的关联关系，因果关系进行具体的量化研究。这种研究

必须以定性分析为先导，定量与定性紧密结合，系统模型的建立，一般要经过思

想开发，因素分析，量化，动态化，优化五个步骤。即语言模型，网络模型，量

化模型，动态模型，优化模型。

在建模过程中，要不断的将下一阶段中所得的结果回馈，经过多次循环往返，使

整个模型逐步趋于完善。

灰色预测方法的特点表现在：首先是它把离散数据视为连续变量在其变化过程中

所取的离散值，从而可利用微分方程式处理数据；而不直接使用原始数据而是由

它产生累加生成数，对生成数列使用微分方程模型。这样，可以抵消大部分随机

误差，显示出规律性。

灰色系统理论的建模思想

下面举一个例子，说明灰色理论的建模思想。考虑4个数据，记为(

0)(1),(0)(2),(0)(3),(0)(4)

上图表明原始数据X(0)没有明显的规律性，其发展态势是摆动的。如果将原始数

据作累加生成1－AGO，记第K个累加生成为X(1)(K),并且

X(1)(1)=X(0)(1)=1

X(1)(2)=X(0)(1)+X(0)(2)=1+2=3

X(1)(3)=X(0)(1)+X(0)(2)+X(0)(3)=1+2+1.5=4.5

X(1)(4)=X(0)(1)+X(0)(2)+X

(0)(3)+X(0)(4)=1+2+1.5+3=7.5

上图表明生成数列X是单调递增数列。

4.1GM(1,1)模型

灰色系统理论的微分方程成为Gm模型，G表示gray（灰色），m表示model

（模型），Gm（1，1）表示1阶的、1个变量的微分方程模型。

定义4.1.1设X0=(x0(1),x0(2),

X1=(x1(1),x1(2),

x0(k)+ax1(k)=b,x0(n)),x1(n))，则称

为GM(1,1)模型的原始形式。

定义4.1.2设X0=(x0(1),x0(2),

X1=(x1(1),x1(2),

Z1=(z1(2),z1(3),,x0(n))，,x1(n))，,z1(n))

k=1,2,n其中z1(k)=1(x1(k)+x1(k-1))2

则称x0(k)+az1(k)=b为GM(1,1)模型的基本形式。

定义4.1.3设X0为非负序列：

X0=(x0(1),x0(2),

X1为X0的,x0(n))1-AGO（即一次累加）序列：

,x1(n))，

nX1=(x1(1),x1(2),其中x1(k)=∑x0(i)；k=1,2,

i=1k；Z1为X1的紧邻均值生成序列

Z1=(z1(2),z1(3),,z1(n))，k=1,2,n其中z1(k)=1(x1(k)+x1(k-1))2

⎡x0(2)⎤⎢x(3)⎥∧T若a=[a,b]为参数列，且Y=⎢0⎥⎢⎥⎢⎥x(n)0⎣⎦⎡-Z1(2)1⎤⎢-

Z(3)1⎥⎥,B=⎢1⎢⎥⎢⎥-Z(n)1⎣1⎦

则GM(1,1)模型x0(k)+az1(k)=b的最小二乘估计参数列满足

a=[a,Tb=]∧T-1(BTB)BY定义4.1.４设X0为非负序列，X1为X0的1-AGO（即

一次累加）

dx1+ax1=bdt序列，Z1为X1的紧邻均值生成序列，则称

为ＧＭ（１，１）模型x0(k)+az1(k)=b的白化方程，也叫影子方程。

定理４.１.１设B,Y,a如定理４.１３所述，a=[a,b]T=(BTB)-1BTY，则∧∧

１.白化方程

bbx1(t)=(x1(1)-)e-at+aadx1+ax1=bdt的解（也称时间响应函数）为

２.GM(1,1)模型x0(k)+az1(k)=b的时间响应函数序列为b-akbx1(k+1=)x0(()-

1e+aa∧k=n1,2,,

３.还原值x0(k+1=)α

a∧1∧1x(k+-1=)x1eak∧k+(∧-1)x1b⎫⎛=(1-e)0x()1+⎪a⎭⎝b⎫⎛=(1-ea)

0x()1+⎪a⎭⎝;=k1,2,n,=k1,2,n,k()-aek;

定义4.1.５称GM(1,1)模型中的参数-a为发展系数，b为灰色作用量。

-a反映了x1∧与x的发展态势。一般情况下，系统作用量应0∧

是外生的或者前定的，而GM(1,1)是单序列建模，只用到系统的行为序列（或称

输出序列，背景值），而无外生作用序列（或称输入序列，驱动量）。GM(1,1)

模型中的灰色作用量b是从背景值挖掘出来的数据，它反映数据变化的关系，其

确定内涵是灰色的，灰色作用量是内涵外延化的具体表现，它的存在，是区别灰

色建模与一般输入输出建模的分水岭，也是区别灰色系统观点与灰箱观点的重要

标志。

定理４.１.２GM(1,1)模型

x0(k)+az1(k)=b

可转化为x0(k)=β-αx1(k-1)，其中β=

ba,α=1+0.5a1+0.5a

４.２残差ＧＭ（１，１）模型

当ＧＭ（１，１）模型的精度不符合要求时，可用残差序列建立ＧＭ（１，１）

模型，对原来的模型进行修正，以提高精度。

定义４.２.１设X0为原始序列，X1为X0的１-ＡＧＯ序列，ＧＭ（１，１）模

型的时间响应式

bb

x1(k+1)=(x0(1)-)e-ak+

aa

∧

k=1,2,

,n

则称

b

dx1(k+1)=-a(x0(1)-)e-ak

a

∧

k=1,2,

,n

为导数还原值。

命题４.２.１设dx

∧

∧

b-ak

(k+1)=-a(x1-)e()10

a

∧

k=1,2,

,n为

b⎫-aka⎛x0(k+1)=x1(k+1)-x1(k)=(-1ex1-())0⎪e导数还原值，a⎭⎝

k=1,2,,n

∧

为

累减还原值，则dx

∧

1

(k+1)≠x0(k+1)

∧

k=1,2,

,n

由此命题知，ＧＭ（１，１）模型既不是微分方程，也不是差分方程。但当a充

分小时，1-e

a

≈-a，有dx1(k+1)≈x0(k+1)。

∧∧

这说明了微分与差分的结果十分接近，因此ＧＭ（１，１）模型既可以看成微分

方程，也可以看成差分方程。

鉴于导数还原值与累减还原值不一致，为减少往返运算造成的误差，往往用X1

的残差修正X1的模拟值x

定义４.２.２设ε0

∧

∧

1

(k+1)。

=(ε0(1),ε0(2),ε0(n))

其中ε(k)=x1(k)-x1(k)为X1的残差序列。若存在k0，满足

１.∀k≥k0,ε0(k)的符号一致。

２.n-k0≥4，则称(ε0(k0),ε0(k0+1),为可建模残差尾段，仍记为ε0

命题４.２.2设ε0

ε0(n)

ε0(n))

)

=(ε0(k0),ε0(k0+1),

=(ε0(k0),ε0(k0+1),ε0(n))为可建模残差

尾段，其１-ＡＧＯ序列为

ε1=(ε1(k0),ε1(k0+1),

ε1(n))

的ＧＭ（１，１）的时间响应式为

ε(k+1)=(ε0(k0)-

1

∧

bε-aε(k-k0)bε

)e+aεaε

k≥k0

则残差尾段ε0的模拟序列为

∧

⎛∧0

ε=ε(k0),ε0(k0+1),

⎝∧0

ε

∧0

(n)⎪

⎭

⎫

其中ε(k+1)=-aε(ε0(k0)-

∧

bε-aε(k-k0)

)eaε

∧

∧

k≥k0

k=1,2,

,n，则相应的

定义４.２.3若x

∧

(k)=x1(k)-x1(k-1)

残差修正时间响应式

b-ak⎧a

1-e(x1-)e()()0⎪∧a⎪

x0(k+1)=⎨

⎪(1-ea)(x0(1)-b)e-ak±aε(ε0(k0)-bε)e-aε(k-k0)

aaε⎪⎩

k

称为累减还原式的残差修正模型。

定义４.２.4若x时间响应式

b-ak⎧

-a(x1-)e()()0⎪∧a⎪

x0(k+1)=⎨

⎪(-a)(x0(1)-b)e-ak±aε(ε0(k0)-bε)e-aε(k-k0)

aaε⎪⎩

k

∧

b

(k+1)=(-a)(x0(1)-)e-ak，则相应的残差修正

a

称为导数还原式的残差修正模型。

例4.2.1湖北省云梦县油菜发病率数据为

X0=(6,20,40,25,40,45,35,21,14,18,15.5,17,15)

建立GM（1，1）模型，得时间响应式为

x1(k+1)=-567.999e-0.06486k+573.999

作累减还原，得

⎛35.6704,33.4303,31.3308,29.3682,27.5192,25.7900,⎫X0=⎪

⎝24.1719,22.6534,21.2307,19.8974,18.6478,17.4768⎭

检验其精度：列出误差检验表如下。

由表可以看出，模拟误差较大，进一步计算残差平方和

s=εεT=57.18

平均相对误差

113

∆=∑∆k=30.11%12k=2

残差平方和很大，相对精度不到70%，需采用残差模型进行修正。取k0=9得残

差尾段

ε0=(ε0(9),ε0(10),ε0(11),ε0(12),ε0(13))

=(-8.6534,-3.2307,-43974,-1.6478,-2.4786)

用此尾段可建立残差尾段模型，取绝对值，得

ε0=(8.6534,3.2307,43974,1.6478,2.4786)

建立GM（1，1）模型，得ε0的1-AGO序列ε1的时间响应式

ε(k+1)=-24e-0.16855(k-9)+32.71

其导数还原值为

ε0(k+1)=(-0.16885)(-24)e-0.16855(k-9)=4.0452e-0.16855(k-9)

∧∧由x∧b⎫-aka⎛-0.06486k(k+1)=x(k+1)-x(k)=1-ex1-()()0110⎪e=38.0614ea⎭⎝

得累减还原式的残差修正模型为

-0.06486k⎧,k<9⎪38.0614ex0(k+1)=⎨-0.16855(k-9)-0.06486k-

4.0452e,k≥9⎪⎩38.0614e∧

其中，ε(k+1)的符合与原始残差序列的符合一致。0

按此模型，可对k=9,10,11,12,13四个模拟值进行修正，修正后

的精度有明显的提高：113计算可得平均相对误差降为∆=∑∆k=4.595%4k=10

灰色系统模型的检验

定义1.

设原始序列

X(0)=x(0)(1),x(0)(2),,x(0)(n){}

相应的模型模拟序列为

ˆ(0)=xˆ(0)(1),xˆ(0)(2),,xˆ(0)(n)X{}

残差序列

ˆ(0)(1),x(0)(2)-xˆ(0)(2),,x(0)(n)-xˆ(0)(n)}={x(0)(1)-x

相对误差序列

⎧ε(1)ε(2)ε(n)⎫∆=⎨(0),(0),,(0)⎬x(n)⎭⎩x(1)x(2)ε(0)={ε(1),ε(2),ε(n)}

={∆k}1n

(k)1.对于k＜n,称∆k=ε为k点模拟相对误差，称(0)x(k)

1n为滤波相对误差，称=∑∆k∆n=(0)nk=1x(n)ε(n)为平均模拟相对误差；

2.称1-为平均相对精度，1-∆n为滤波精度；

3.给定α，当<α，且∆n<α成立时，称模型为残差合格模型。

定义2

ˆ(0)为相应的模拟误差序列，ε为X(0)与设X(0)为原始序列，X

ˆ(0)的绝对关联度，X若对于给定的ε0>0,ε>ε0，则称模型为关联合格模型。

定义3

ˆ(0)为相应的模拟误差序列，ε(0)为残差设X(0)为原始序列，X

序列。

1n(0)=∑x(k)为X(0)的均值，nk=1

1n(0)s=∑(x(k)-)2为x(0)的方差，nk=1

1n=∑ε(k)为残差均值，nk=1

1n2s2=∑(ε(k)-)2为残差方差，nk=1

s1.称c=2为均方差比值；对于给定的c0>0，当c

型为均方差比合格模型。

2.称p=pε(k)-0,当p>p0时，称模型为小

误差概率合格模型。

应用举例

例6-1

设原始序列

X(0)=x(0)(1),x(0)(2),x(0)(3),x(0)(4),x(0)(5){}

=(2.874,3.278,3.337,3.390,3.679)

建立Gm(1,1)模型，并进行检验。

解：1）对X(0)作1-AGO，得

[D为X(0)的一次累加生成算子，记为1-AGO，AcumulatedGeneratingOperator]

X(1)={x(1)(1),x(1)(2),x(1)(3),x(1)(4),x(1)(5)}

=(2.874,6.152,9.489,12.579,16.558)

Matlabx0=[2.8743.2783.3373.3903.679];x1=zeros(size(x0));x1(1)=x0(1);fori=2:5

forj=1:ix1(i)=x1(i)+x0(j);endend

2)对X(1)作紧邻均值生成，令

Z(1)(k)=0.5x(1)(k)+0.5x(1)(k-1)

Z(1)={z(1)(1),z(1)(2),z(1)(3),z(1)(4),z(1)(5)}

=(2.874,4.513,7.820,11.84,14.718)

于是，

(0)⎡x(2)⎤1⎤⎡3.278⎤(0)⎢⎥1⎥，Y=⎢x(3)⎥=⎢3.337⎥

(0)1⎥x(4)⎥⎢3.390⎥⎢1⎥⎢(0)⎦⎣3.679⎥⎦⎢⎥x(5)⎣⎦

⎡-4.5131⎤.513-7.820-11.184-14.718⎤∙⎢-7.8201⎥BTB=⎡-41111⎥⎢⎣⎦⎢-

11.841⎥⎢⎣-14.7181⎥⎦

423.221-38.235⎤=⎡4⎢⎥⎣-38.235⎦⎡-z(1)(2)⎢-z(1)(3)B=⎢(1)⎢-z(4)(1)⎢⎣-

z(5)1⎤⎡-4.513⎥1⎢-7.820⎥=1⎥⎢-11.84⎢⎣-14.7181⎥⎦

423.221-38.235⎤==⎡0.0173180.165542⎤(BB)=⎡-

4⎥⎢⎢⎥⎣38.235⎦⎣0.16655421.832371⎦

1⎧38.235⎤⎫⎡4=⎪423.221⨯4-38.2352⎢⎥⎣38.235423.221⎦⎪⎪⎪⎨⎬

138.235⎤⎡4⎪=⎪38.235423.221⎥⎪230.969⎢⎪⎣⎦⎩⎭

ˆ=(BTB)-1BTYaT-1-1

⎡3.278⎤0.0173180.165542⎤∙⎡-4.513-7.820-11.184-

14.718⎤∙⎢3.337⎥=⎡0111⎥⎢⎥⎣.16655421.832371⎦⎢⎣1⎦⎢3.390⎥⎢⎣3.679⎥⎦

⎡3.278⎤0.0873860.030115-0.028143-0.089344⎤∙⎢3.337⎥=⎡1⎢⎣.0852800.537833-

0.01905110.604076⎥⎦⎢3.390⎥⎢⎣3.679⎥⎦

.037156⎤=⎡-30

⎣.065318⎥⎦⎢

3）确定模型

dx(1)

-0.037156x(1)=3.065318dt

及时间响应式

bb

ˆ(1)(k+1)=(x(0)(1)-)e-ak+x

aa

=85.3728e0.037156k-82.4986

4)求X(1)的模拟值

ˆ(1)={xˆ(1)(1),(1)(2),xˆ(1)(3),xˆ(1)(4),xˆ(1)(5)}Xx

=(2.8740,6.1058,9.4599,12.9410,16.5538)

(0)

5)还原出X的模拟值，由

ˆ(0)(k+1)=xˆ(1)(k+1)-xˆ(1)(k)x

ˆ(0)={xˆ(0)(1),xˆ(0)(2),xˆ(0)(3),xˆ(0)(4),xˆ(0)(5)}得X

=（2.8740，3.2318，3.3541，3.4811，3.6128）

6）误差检验

⎡ε(2)⎤⎢ε(3)⎥

s=εTε=[ε(2)ε(3)ε(4)ε(5)]∙⎢

ε(4)⎥⎢ε(5)⎥⎣⎦

⎡0.0462⎤

0.0171⎥=[0.0462-0.0171-0.09110.0662]∙⎢-

⎢-0.0911⎥0.0662⎢⎥⎣⎦

=0.0151085

平均相对误差

151

∆=∑∆k=(1.41%+0.51%+2.69%+1.80%)4k=14

=1.0625%

ˆ的灰色关联度计算X与X

S=1(x(k)-x(1)+(x(5)-x(1))∑2k=2

24=(3.278-2.874)+(3.337-2.874)+(3.390-2.874)+1(3.679-2.874)

=0.404+0.463+0.516+=1.7855

ˆ=S1ˆˆˆ(5)-xˆ(1)(x(k)-x(1)+(x∑2k=24

1=(3.2318-2.874)+(3.3541-2.874)+(3.4811-2.874)+(3.6128-2.874)2

=0.3578+0.4801+0.6071+0.=1.8144

ˆ-S=S1ˆˆˆ(5)-xˆ(1))[][(x(5)-x(1))-(x(x(k)-x(1))-(x(k)-x(1))+∑k=242

1=(0.3578-0.404)+(0.4801-0.463)+(0.6071-0.516)+(0.3694-0.4025)2

=-0.0462+0.0171+0.091-0.=0.04535

ε=ˆ1+S+S

ˆ+Sˆ-S1+S+S=1+1.7855+1.81444.5999=1+1.7855+1.8144+0.045354.64525

=0.9902＞0.90

精度为一级，可以用

ˆ(1)(k+1)=85.3728xe0.037156k-82.4986

ˆ(0)(k+1)=xˆ(1)(k+1)-xˆ(1)(k)预测。x

Matlabclc;clear;

x0=[2.8743.2783.3373.3903.679];

x1=zeros(size(x0));

x1(1)=x0(1);

fori=2:5

forj=1:i

x1(i)=x1(i)+x0(j);

end

end

%对作紧邻均值生成z1

z1=zeros(size(x1));

z1(1)=x1(1);

fori=2:5

z1(i)=0.5*x1(i)+0.5*x1(i-1);

end

%求参数a

B=[-z1(2:end)',ones(4,1)];

Y=x0(2:end)';

a=inv(B'*B)*B'*Y;

%系数求出，模型确定

%求出时间响应式Gx1(k+1)=85.3728*exp(0.03715*k)-82.4986

f=inline('85.3728*exp(0.037156*x)-82.4986');x=0:4;

Gx1=f(x);

%还原出x0的模拟值

Gx0=zeros(1,5);

Gx0(1)=Gx1(1);

fori=2:5

Gx0(i)=Gx1(i)-Gx1(i-1);

end

%误差检验

e=x0-Gx0;%残差

e=e(2:end);

e1=abs(e./x0(2:end))%相对误差Ae1=sum(e1)/4%平均相对误差

s=e*e'%残差平方和为

%计算X0与Gx0的灰色关联度ee

s=0;

fori=2:4

s=s+abs(x0(i)-x0(1));

end

s=s+abs(0.5*(x0(5)-x0(1)));

Gs=0;

fori=2:4

Gs=Gs+abs(Gx0(i)-Gx0(1));

end

Gs=Gs+abs(0.5*(Gx0(5)-Gx0(1)));

CHA=abs(Gs-s)

ee=(1+abs(s)+abs(Gs))/(1+abs(s)+abs(Gs)+CHA)

例6-2

试建立Gm(1,1)模型的白化方程及时间响应式，并对Gm(1,1)模型进行检验，预测

该企业2001-2005年产值。解：设时间序列为

=（27260，29547，62411，35388）

X(1)={x(1)(1),x(1)(2),x(1)(3),x(1)(4)}

)=(27260,56807,89218,124606

2)对X(1)作紧邻均值生成，令

Z(1)(k)=0.5x(1)(k)+0.5x(1)(k-1)

Z(1)={z(1)(1),z(1)(2),z(1)(3),z(1)(4)}

)=(27260,42033.5,73012.5,106912X(0)=x(0)(1),x(0)(2),x(0)(3),x(0)(4){}

于是，

(0)⎡-z(1)(2)1⎤⎡-42033⎡x(2)⎤⎡29547.51⎤⎤(1)(0)⎢⎥⎢⎥B=-z(3)1=⎢-73012.51，

Y=x(3)=⎢32411⎥⎢(1)⎥⎢1069121⎥⎢⎥(0)⎥⎢35388⎥⎦⎦⎢⎢⎣-

z(4)1⎥⎦⎣⎣x(4)⎥⎦⎣

ˆ=[a,b]T作最小二乘估计，得对参数列α

ˆ=(BTB)-1BTY=⎡-0.089995⎤a⎢⎣25790.28⎥⎦

dx(1)

-ax(1)=b设dt

由于

可得Gm(1,1)模型的白化方程

dx(1)

--0.089995x(1)=25790.28dta=-0.089995,b=25790.28

其时间响应式为

bb⎧(1)ˆ(k+1)=(x(0)(1)-)e-ak+=313834e0.089995k-286574⎪x

⎨aa⎪ˆ(0)(k+1)=xˆ(1)(k+1)-xˆ(1)(k)x⎩

由此得模拟序列

ˆ(0)={xˆ(0)(1),xˆ(0)(2),xˆ(0)(3),xˆ(0)(4)}X

=（27260，29553，32336，35381）检验：

残差序列为

ε(0)=(ε(0)(1),ε(0)(2),ε(0)(3),ε(0)(4))

=（0，-6，75，7）

⎧ε(0)(1)ε(0)(2)ε(0)(3)ε(0)(4)⎫∆=⎨(0),(0),(0),(0)⎬x(1)x(2)x(3)x(4)⎩⎭

=(0,0.0002,0.00231,0.0002)∆=(∆1,∆2,∆3,∆4)

平均相对误差

14

∆=∑∆k=0.00068=0.068%<0.014k=1

模拟误差∆4=0.0002=0.02%<0.01，精度一级

ˆ(0)的灰色关联度ε计算X(0)与X

S=

ˆ=S1(0)(0)(0)(x(k)-x(1)+(x(4)-x(0)(1))=11502∑2k=21ˆˆˆ(4)-xˆ(1)(x(k)-

x(1)+(x∑2k=233=11429.5

ˆ-S=Sˆ∑[(x

k=23(0)ˆ(0)(1))-(x(0)(k)-x(0)(1))+(k)-xˆ]1[(x2(0)ˆ(0)(1))-(x(0)(4)-x(0)(1))(4)-x=72.5

ε=ˆ1+S+S

ˆ+Sˆ-S1+S+S=1+11502+11429.5=0.997>0.901+11502+11429.5+72.5

精度为一级

计算均方差比

14(0)=∑x(k)=31151.54k=1

14(0)2S1=∑(x(k)-)2=9313116.25,S1=3051.744k=1

14=∑ε(K)=194K=1

142S2=∑(ε(k)-)2=1066.5,s2=32.664k=1

S32.66所以，c=2==0.011<0.35，均方差比值为一级S13051.74

计算小误差概率

0.6745S1=2058.40

ε(1)-=19,(2)-=25ε(3)-=56,ε(4)-=12

所以，P=P(ε(k)-0.95

小误差概率为一级，故可用

⎧xˆ(1)(k+1)=313834e0.089995k-286574⎨ˆ(0)(1)(1)ˆˆ⎩x(k+1)=x(k+1)-x(k)

进行预测，2001-2005年预测值为

ˆ(0)={xˆ(0)(5),xˆ(0)(6),xˆ(0)(7),xˆ(0)(8),xˆ(0)(9)}X

=（38713，42359，46318，50712，55488）

例6-3

解：X(0)={x(0)(1),x(0)(2),x(0)(3),x(0)(4),x(0)(5)}

=(1.67,1.51,1.03,2.14,1.99)

X(1)=x(1)(1),x(1)(2),x(1)(3),x(1)(4),x(1)(5){}

=(1.67,3.18,4.21,6.35,8.43)

对X(1)作紧邻均值生成，令

Z(1)(k)=0.5x(1)(k)+0.5x(1)(k-1)

Z(1)=z(1)(1),z(1)(2),z(1)(3),z(1)(4),z(1)(5){}

=(1.67,2.425,3.695,5.28,7.345)

于是，

(0)⎡x(2)⎤1⎤⎡1.51⎤(0)⎢⎥1⎥，Y=⎢x(3)⎥=⎢1.03⎥

(0)1⎥x(4)⎥⎢2.14⎥⎢1⎦⎥⎢1.99⎦⎥(0)⎣⎢⎥x(5)⎣⎦

⎡-2.4251⎤.425-3.695-5.28-7.345⎤∙⎢-3.6951⎥BTB=⎡-21111⎥⎢⎣⎦⎢-5.281⎥⎢⎣-

7.3451⎥⎦

101.361-18.745⎤=⎡4⎥⎢⎣-18.745⎦

0.073980.34669⎤(BTB)-1=⎡0⎢⎣.346691.8747⎥⎦

⎡1.51⎤.425-3.695-5.28-7.345⎤∙⎢1.03⎥BTY=⎡-21111⎥⎢⎣⎦⎢2.14⎥⎢⎣1.99⎥⎦

.38335⎤=⎡-33⎢⎣6.67⎥⎦

ˆ=⎡a⎤=(BTB)-1BTY=⎡0.0379800.34669⎤∙⎡-33.38335⎤

a⎢⎢⎣b⎥⎦⎣0.346691.8747⎥⎦⎢⎣6.67⎥⎦

.157⎤=⎡-00⎢⎣.931⎥⎦

dx(1)

-ax(1)=b方程为dt(1)dx-0.157x(1)=0.931dt⎡-z(1)(2)⎢-z(1)(3)B=⎢(1)⎢-z(4)(1)⎢⎣-

z(5)1⎤⎡-2.425⎥1⎢-3.695⎥=1⎥⎢-5.28⎢-7.345⎣1⎥⎦

及时间响应式

bbˆ(1)(k+1)=(x(0)(1)-)e-ak+xaa

=(1.67+5.93)e0.157k-5.93

=7.6e0.157k-5.93

求X(1)的模拟值

ˆ(1)={xˆ(1)(1),(1)(2),xˆ(1)(3),xˆ(1)(4),xˆ(1)(5)}Xx

=(1.67,2.962,4.474,6.202,8.311)还原出X(0)的模拟值，由

ˆ(0)(k+1)=xˆ(1)(k+1)-xˆ(1)(k)x

ˆ(0)={xˆ(0)(1),xˆ(0)(2),xˆ(0)(3),xˆ(0)(4),xˆ(0)(5)}得X

=（1.67，1.292，1.512，1.728，2.109）

ˆ的灰色关联度计算X与X

S=1(x(k)-x(1)+(x(5)-x(1))∑2k=2

=(-0.61-0.64+0.47+1⨯0.32=0.172

411ˆˆ(k)-xˆ(1)+(xˆ(5)-xˆ(1)=-0.378-0.1580.058+⨯0.S=∑(x22k=24

=0.2585

ˆ-S=S1ˆˆˆ(5)-xˆ(1))[][(x(5)-x(1))-(x(x(k)-x(1))-(x(k)-x(1))+∑k=242

=(-0.378+0.16+(-0.158+0.64)+(0.058-0.47+1(0.439-0.32)2

=-0.218+0.482-0.417+0.=0.0885

ε=ˆ1+S+S

ˆ+Sˆ-S1+S+S=1+0.17+0.25851.4285=1+0.17+0.2585+0.08851.517

=0.94＞0.90

精度为一级，关联度为一级，可以用

⎧xˆ(1)(k+1)=7.6e0.157k-5.93⎨ˆ(0)(1)(1)ˆˆx(k+1)=x(k+1)-x(k)⎩

进行预测。

ˆ(1)={xˆ(1)(6),xˆ(1)(7),xˆ(1)(8),xˆ(1)(9),xˆ(1)(10),X

=(10.732,13.565,16.879,20.756,25.293,ˆ(1)(11),xˆ10)(12),xˆ(1)(13),xˆ(1)(14),xˆ(1)(15)

x

30.601,36.811,44.076,52.577,62.523)}

,2.833,3.314,3.877,4.537,=(2.421

5.308,6.21,7.265,8.501,9.946)

ˆ(0)=(xˆ(0)(6),xˆ(0)(7),xˆ(0)(8),xˆ(0)(9),xˆ(0)(10),,xˆ(0)(15))X

(1,1)模型的适用范围

命题4.3.1当(n-1)∑⎡⎣z(k)⎤⎦→∑⎡⎣z(k)⎤⎦时，GM(1,1)模型无

k=2k=211n2n2

意义。因为此时a→∞,b→∞

命题4.3.2当GM(1,1)发展系数a≥2时，GM(1,1)模型无意义。即当a<2时，

GM(1,1)模型才有意义。但随着a的取值不同，预测效果也不同。要得到较高精

度，要求a∈⎛

2⎛-n2⎫也可以推出级比σ(k)=e+1,en+1⎪⎝⎭0∧∧-22⎫,⎪∈(-2,2)，⎝n+1n+1⎭

第五章灰色系统预测

预测就是借助于过去的探讨去推测，了解未来。灰色预测就是通过原始数据的处

理和灰色模型的建立，发现，掌握系统发展规律，对系统未来状态做出科学定理

预测。

定义5.1.1设原始数据序列

X0=(x0(1),x0(2),,x0(n))

相应的预测模型模拟序列

∧⎛∧x0=x0(1),x0(2),⎝∧⎫x0(n)⎪⎭∧

残差序列

ε0=(ε0(1),ε0(2),ε0(n))

∧∧⎛=x0(1)-x0(1),x0(2)-x0(2),⎝⎫x0(n)-x0(n)⎪⎭∧

相对误差序列⎛ε0(1)ε0(2)=,,x0(1)x0(2)⎝,ε0(n)⎫⎪={x0(n)⎪⎭k}n1则

1.对于k≤n，称⎛ε0(k)⎫=⎪kx0(k⎪⎝⎭，为k点的模拟相对误差，称1n

∆=∑∆k为平均相对误差。nk=1

2.称1-∆为平均相对精度，1-∆k为k点的模拟精度。

3.给定α，当∆<α且∆n<α成立时，称模型为残差合格模型。

定义5.1.2设X0为原始序列，X为相应的模拟序列，ε为X0与X的绝对关联

度，若对于给定的ε0>0，有ε>ε0，则称模型为关联度合格模型。

定义5.1.3设X0为原始序列，X为相应的模拟序列，ε0为X0与X的残差序列，

则

1n0x=∑x(k),nk=11n02S1=∑x(k)-xnk=10000()

22分别为X0的均值、方差；1n0=∑ε(k),nk=11n0S=∑ε(k)-nk=122()

分别为残差的均值、方差。

1.C=s2

s1称为均方差比值，对于给定的C0>0,当C

型为均方差比合格模型。

2.p=p(ε(k)-0

当p>p0，称模型为小误差概率合格模型。

精度检验等级参照表

5.2数列预测

数列预测是对系统变量的未来行为进行预测，GM（1，1）是较常用的数列预测

模型。根据实际情况，也可以考虑采用其他灰色模型，在定性分析的基础上，定

义适当的算子，对算子作用后的序列建立GM模型，通过精度检验后，即可用于

预测。

例5.2.1河南省长葛县乡镇企业产值（数据来源于长葛县统计局）。

解由统计资料查得产值序列为

X0=(x0(1),x0(2),x0(3),x0(4))=(1015512588，，23480，35388)引入二阶弱化算子

D2，令

X0D=(x0(1)d,x0(2)d,x0(3)d,x0(4)d)

其中x0(k)d=1(x0(k)+x0(k+1)+4-k+1+x0(4));以及

X0D2=(x0(1)d2,x0(2)d2,x0(3)d2,x0(4)d2)其中

x0(k)d2=1(x0(k)d+x0(k+1)d+4-k+1+x0(4)d);于是

X0D2=(27260,29547,32411,35388)X=(x(1),x(2),x(3),x(4))X的1-AGO序列为

X1=(27260,56807,89218,124606)设dx1+ax1=bdt

按最小二乘法求得参数a,b的估计值为

a=-0.089995,b=25790.28

得GM（1，1）模型白化方程

dx1-0.089995x1=25790.28dt其时间响应式为

0.089995k⎧-286574⎪x1(k+1)=313834e⎨⎪⎩x0(k+1)=x1(k+1)-x1(k)

得模拟序列

X=x(1),x(2),(,x(4)=(27260,29533,32337,35381))

残差序列ε0=(0,-6,74,7)

相对误差序列⎛ε0(1)ε0(2)=,,x0(1)x0(2)⎝,ε0(n)⎫⎪x0(n)⎪⎭

=(0,0.0002,0.00228,0.0002)(∆1,∆2,∆3,∆4)

平均相对误差

14

∆=∑∆k=0.00067=0.067%<0.014k=1

模拟误差∆4=0.0002<0.1，精度为一级。

计算X与x的灰色绝对关联度ε：

s=1⎡⎤xk-x1+∑⎣(()())⎦2⎡⎣x(4)-x(1)⎤⎦=11502

k=2

33

s=1⎡⎡⎤xk-x1+x(4)-x(1)⎤=11430.5()()∑⎣⎦⎣⎦2k=2()

s-s=1⎡⎡⎤xk-x1-xk-x1+x(4)-x(1)-x(4)-x(1)⎤=71.5()()()()∑⎣⎦2⎣⎦k=23()()从而

ε=1+|s|+|s|

1+|s|+|s|+|s-s|=1+11502+11430.5=0.997>0.901+11502+11430.5+71.5

关联度为一级

计算均方差比C14

x=∑x(k)=31151.5,4k=1

S12=1x(k)-x∑4k=14()2=37252465,S1=6103.48

14=∑ε4k=1

2S2=(k)=18.75,41∑ε4k=1((k)-ε)2=4154.75,S2=64.46

所以C=S264.46==0.01<0.35,均方差比值为一级。S16103.48

计算小误差概率：

0.6745S1=4116.80

ε(1)-=18.75,

ε(3)-=55.25,ε(2)-=24.75ε(4)-=11.75

所以p=p(ε(k)-0.95,小误差概率为一级，故可用

0.089995k⎧-286574⎪x1(k+1)=313834e⎨⎪⎩x0(k+1)=x1(k+1)-x1(k)

进行预测。这里给出5个预测值

X=x(5),x(6),(,x(9))

=(38714,42359,46349,50712,55488)