我国碳排放量预测模型

摘要

本文主要我国碳排放预测问题，同时根据预测结果提出合理性建议。以人口总量，城镇化，人均GDP,第三产业GDP比例,能源强度吨标准煤,煤炭消费比例 的数据，建立GM（1,1）预测模型、多元线性回归预测模型、BP神经网络预测模型，借助Matlab软件逐个对碳排放量和影响因素数据进行模拟与预测，然后采用绝对误差与相对误差两个参数对模型进行评价与对比，接着应用关联度分析法求得影响因素的重要性排序，最后结合重要性排序向相关部门提出建议。

对于GM（1,1）预测模型，通过对1986至2010年原始单变量数据进行生成处理，寻找系统的变化规律建立相应的微分方程预测模型，代入相关单变量数据用Matlab编程得到各单变量在2011至2015年的预测值。

对于多元线性回归预测模型，确定线性预测变量和因变量，即影响因素和测度指标，将数据代入Matlab统计软件，求得多元线性方程，将1986至2010年所有数据代入该方程，同时结合GM（1,1）预测模型对2011至2015年各单变量预测结果，用Matlab编程得到对应年份的碳排放量模拟值和预测值。

对于BP神经网络预测模型，首先根据碳排放量的排放趋势，确定输出层、中间隐层和输入层，然后把样本分为训练样本和测试样本两个部分，在以上基础,对样本数据进行归一化预处理，结合GM（1,1）预测模型对2011至2015年各单变量预测结果，采用Matlab软件中的神经网络计算功能,建立合理训练模型得到对应年份的旅游人数模拟值和预测值。

在模型求解过程中，将得到其对应的平均绝对误差值和相对误差值，通过比较知3个预测模型的精确度都合格。其中BP神经网络模型误差最小，预测效果最佳，三种模型2011-2015年预测数据如下表。

模型 2011 2012 2013 2014 2015 GM(1,1)模型 线性回归模型

Bp网络模型

77.8641 85.073 87.2029

83.4852 90.4646 95.4649

89.5121 96.1978 104.5097

95.9741 102.2945 114.4115

102.9026 108.7775 125.2514

对于影响因素重要性确定，本文应用关联度分析法建立因素排序模型，将数据代入关联系数公式得出影响因素数列对参考数列在每个年份的关联系数，关联度即各个关联系数之和的平均值，按关联度大小排序可得影响因素的重要性排序：人均GDP>人口>煤炭消费比例>城镇化>能源强度比例>第三产业GDP比例。 最后根据重要性排序，向有关部门提出一些减少碳排放量的建议。

关键词：碳排放量预测 GM（1,1）预测模型 BP神经网络预测模型 多元线性回归预测 关联度分析法 碳排放 Matlab软件

目录

1.问题重述…………………………………………………………………………4 2.问题分析…………………………………………………………………………4 3.问题假设…………………………………………………………………………4 4.符号说明…………………………………………………………………………5 5.模型建立与求解…………………………………………………………………5 5.1GM(1,1)预测模型……………………………………………………………6 5.1.1模型思路………………………………………………………………6 5.1.2模型建立………………………………………………………………6 5.1.3模型求解………………………………………………………………7 5.2多元线性回归预测模型……………………………………………………9 5.2.1模型思路………………………………………………………………9 5.2.2模型建立………………………………………………………………10 5.2.3模型求解………………………………………………………………10 5.3BP神经网络预测模型………………………………………………………11 5.3.1模型思路………………………………………………………………11 5.3.2模型建立………………………………………………………………11 5.3.3模型求解………………………………………………………………12 5.4因素排序模型………………………………………………………………13 5.4.1模型思路………………………………………………………………13 5.4.2模型建立………………………………………………………………13 5.4.3模型求解………………………………………………………………14 6.模型检验分析……………………………………………………………………14 7.建议………………………………………………………………………………15 8.模型评价与推广…………………………………………………………………15 8.1模型优点………………………………………………………………………15 8.2模型缺点………………………………………………………………………15 8.3模型推广………………………………………………………………………15 参考文献……………………………………………………………………………15 附录…………………………………………………………………………………16 附件一……………………………………………………………………………16 附件二……………………………………………………………………………16 附件三……………………………………………………………………………19 附件四……………………………………………………………………………22

1、问题重述

中国是世界上能源生产与消费大国。碳排放问题在我国已引起广泛的关注，“十二五”规划中明确提出要“节约能源，降低温室气体排放强度”。要实现这一目标，需要对碳排放的影响因素进行深入分析，构建科学的预测模型对未来碳排放进行预测，为制定有效的碳减排路径提供决策依据。

据此我们需要解决以下问题： 1、收集中国历年碳排放及其影响因素数据（收集至少近20年的相关数据）； 2、建立至少3种定量预测模型（其中GM（1,1）和BP神经网络模型必需，其它可考虑微分方程、多元回归分析等）对未来中国碳排放进行预测；

3、结合若干性能评价指标对模型进行分析比较；

4、指出影响碳排放的主要因素，向有关部门提出具体建议。

2、问题分析

本文整个过程主要解决问题是我国碳排放预测问题，通过分析确定测度指标为碳排放量y，相关影响因素为以下六个方面：人口总量x1，城镇化x2，人均GDPx3,第三产业GDP比例x4,能源强度吨标准煤x5,煤炭消费比例x6。 本文将建立建立GM（1,1）预测模型、多元线性回归预测模型、BP神经网络预测模型，借助Matlab软件逐个对碳排放量进行模拟与预测。3种定量预测模型各有各的优势与不足，故在几处采用了多个预测模型相结合的方法进行预测，使得模型进一步优化。

通过关联度分析得出影响因素的重要性排序，在此基础上结合我国碳排放量的发展趋势向有关部门提出合理建议，可提高模型建立的科学性。

3、模型假设

1、假设统计的数据真实科学，短期内稳定变化；

2、假设建立模型中，个别偏离太远的数据可据题适当调整；

3、假设碳排放量变化主要受人口总量，城镇化，人均GDP,第三产业GDP比例,能源强度吨标准煤,煤炭消费比例六方面因素的影响；

4、假设碳排放量y值大部分呈线性，满足多元线性回归条件；

5、假设影响碳排放量变化的各因素之间存在相关关系、各因素与碳排放量存在非线性关系，适用BP神经网络的性能。

4、符号说明

Xi：原始数据序列

a：发展恢数 b：内生控制恢数

xi：自变量i1,2,5

y：测度指标

i：参系数

xiT：样本实际输入值

xiO：网络模型输出值

i(k)：比较数列xi对参考数列x0在k时刻的关联系数

：分辨系数，且[0,1]

ri：数列xi对参考数列x0的关联度

5、模型建立与求解

本文探讨的是中国碳排放量的预测模型，根据中国近25年的碳排放量变化趋势进行解题。本文确定测度指标为碳排放量，相关影响因素为以下六方面：人口总量x1，城镇化x2，人均GDPx3,第三产业GDP比例x4,能源强度吨标准煤x5,煤炭消费比例x6。

根据2011年全国年鉴统计表得到有关真实数据如下：

表1 1986至2010年各自变量及测度指标有关数据

年份 1986 1987 1988 1989 1990

第三产业能源强度吨

人口/（万人均GDP/煤炭消费比碳排放量/

城镇化/% GDP比例标准煤/（万

人） （元/人） 例/% 亿吨

/% 元） 107507 109300 111026 112704 114333

24.52 25.32 25.81 26.21 26.41

963 1112 1366 1519 1644

29.14 29.64 30.51 32.06 31.54

9.8 9.4 9.1 9.1 8.9

75.8 76.2 76.1 76.1 76.2

19.70823 21.0278 22.40368 22.75338 22.69709

1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 115823 117171 118517 119850 121121 122389 123626 124761 125786 126743 127627 128453 129227 129988 130756 131448 132129 132802 133450 134091 26.94 27.46 27.99 28.51 29.04 30.48 31.91 33.35 34.78 36.22 37.66 39.09 40.53 41.76 42.99 44.34 45.89 46.99 48.34 49.95 1893 2311 2998 4044 5046 5846 6420 6796 7159 7858 8622 9398 10542 12336 14185 16500 20169 23708 25608 29992 33.69 34.76 33.72 33.57 32.86 32.77 34.17 36.23 37.77 39.02 40.46 41.47 41.23 40.38 40.51 40.94 41.89 41.82 43.43 43.14 8.6 7.9 7.4 6.9 6.6 6.2 5.7 5.3 5.1 4.8 4.6 4.5 4.7 5 4.9 4.8 4.5 4.3 4.2 4.2 76.1 75.7 74.7 75 74.6 73.5 71.4 70.9 70.6 69.2 68.3 68 69.8 69.5 70.8 71.1 71.1 70.3 70.4 68 23.69252 24.49162 26.26645 28.31547 28.61685 28.93377 30.81745 29.67256 28.85722 28.4975 29.69576 34.64843 40.69239 50.8978 55.12703 58.17144 62.56704 68.00468 74.63289 82.40958

◆5.1 GM（1,1）预测模型

▪5.1.1 模型思路

建立GM (1,1)灰色动态预测模型，通过对1986至2010年原始单变量数据进行生成处理，包括累加、转换、代入和还原等过程，以寻找系统的变化规律，在此基础上建立相应的微分方程预测模型。代入相关单变量数据用Matlab编程得到对应变量在2011至2015年的预测值。 ▪5.1.2 模型建立

第一步：假设X0为原始非负序列：

X0{x(0,1),x0,2,,x0,n}，

其中x0,k0，k1,2,,25。

利用累加生成序列可将序列X0{x(0,1),x0,2,,x0,n}生成序列X1：

X1{x(1,1),x1,2,,x1,n}，

其中x1,kxi1k0,i，k1,2,,25。

第二步：利用先生成的序列X1建立GM（1,1）模型的一般形式：

x0,kaz1,kb， （1）

用微分方程表示如下：

dx1ax1b， （2） dtZ1为X1的紧邻均值生成序列：

Z1{z(1,1),z1,2,,z1,n}，

其中

z1,k0.5*x1,k0.5*x1,k1，k1,2,,25， （3）

灰微分方程模型中a，b为待估参数，分别为发展恢数和内生控制恢数；灰微分方程的最小二乘估计参数列满足： a,bBTBT1BTY， （4）

Y为列向量，B为构造矩阵

0.5*x1,20.5*x1,11x0,20.5*x1,30.5*x1,21x0,3BY，， 0.5*0.5*1x1,nx1,n1x0,n第三步：构建灰色预测模型

结合上述（1）、（3）、（4）求解（2）式得：

bbx1,k1(x1,0)eak，k1,2,,n，

aa因为x1,0x0,1，所以建立的灰色模型如下

bbx1,k1(x0,1)eak，k1,2,,n，

aa第四步：求出原始数据的还原值即可得到GM（1,1）预测模型：

x0,k1x1,k1x1,k，k1,2,,n。

第五步：对模型进行后检验模型的精确性

▪5.1.3 模型求解

表一中各自变量构成原始非负序列，代入Matlab软件得GM（1,1）对y值模拟结果如下：

表5.1.3.1 1986至2010 年y值GM（1,1）模型模拟结果

年份 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 20005 2006 2007 2008 2009 2010

y原始数据 19.70823 21.0278 22.40368 22.75338 22.69709 23.69252 24.49162 26.26645 28.31547 28.61685 28.93377 30.81745 29.67256 28.85722 28.4975 29.69576 34.64843 40.69239 50.8978 55.12703 58.17144 62.56704 68.00468 74.63289 82.40958

y模拟值 19.7082 14.6152 15.6703 16.8015 18.0144 19.3149 20.7093 22.2043 23.8073 25.5259 27.3687 29.3445 31.4629 33.7342 36.1695 38.7807 41.5803 44.582 47.8004 51.2512 54.9511 58.9181 63.1714 67.7318 72.6215

绝对误差 0.00003 6.4126 6.73338 5.95188 4.68269 4.37762 3.78232 4.06215 4.50817 3.09095 1.56507 1.47295 -1.79034 -4.87698 -7.672 -9.08494 -6.93187 -3.88961 3.0974 3.87583 3.22034 3.64894 4.83328 6.90109 9.78808

相对误差 0.00% 30.50% 30.05% 26.16% 20.63% 18.48% 15.44% 15.47% 15.92% 10.80% 5.41% 4.78% -6.03% -16.90% -26.92% -30.59% -20.01% -9.56% 6.09% 7.03% 5.54% 5.83% 7.11% 9.25% 11.88%

模拟曲线如下：

110100908070605040302010 0510152025原始曲线GM(1,1)模拟曲线 30

图5.1.3.2 GM（1,1）模型y值模拟曲线

同时得到2011至2015年相关预测结果如下表：

表5.1.3.3 据GM（1,1）模型y值预测结果 年份 2011 2012 2013 2014 2015

◆5.2 多元线性回归预测模型

▪5.2.1 模型思路

由于1986至2010年各影响因素及测度指标有关数据呈现一定的递增形式，故此采用多元线性回归模型对其进行模拟与预测。确定线性预测变量和因变量，建立回归矩阵。预测变量即六个影响因素：x1人口总量，x2城镇化，x3人均GDP, x4第三产业GDP比例, x5能源强度吨标准煤, x6煤炭消费比例。 将实际数据用Matlab软件求解出符合该情况的多元线性方程，将1986至2010年所有数据代入该方程，同时结合GM（1,1）预测模型对2011至2015年

碳排放量/亿吨

77.8641 83.4852 89.5121 95.9741 102.9026

各单变量预测结果，用Matlab编程得到对应变量的y的模拟值和预测值。 ▪5.2.2 模型建立

设回归模型为：

y01x1据题可导：

5x56x6，

yi01xi15xi56xi6

▪5.2.3 模型求解

将实际数据用Matlab软件求解出多元线性方程为：

y0.0016x11.1399x20.0015x30.9702x45.8205x50.5353x6246.3984；

借助Matlab编程得到对应变量y的模拟值：

表5.2.3.1 1986至2010 年y值线性回归模拟结果

年份 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 20005 2006 2007 2008 2009 2010 y原始数据 19.70823 21.0278 22.40368 22.75338 22.69709 23.69252 24.49162 26.26645 28.31547 28.61685 28.93377 30.81745 29.67256 28.85722 28.4975 29.69576 34.64843 40.69239 50.8978 55.12703 58.17144 62.56704 68.00468 74.63289 82.40958 y模拟值 24.3527 19.4675 20.7012 22.0132 23.4083 24.8918 26.4693 28.1469 29.9307 31.8276 33.8447 35.9896 38.2705 40.6959 43.275 46.0176 48.934 52.0325 55.333 58.8398 62.5688 66.5342 70.7508 75.2347 80.0028 绝对误差 -4.64447 1.56030 1.70248 0.74018 -0.71121 -1.19928 -1.97768 -1.88045 -1.61523 -3.21075 -4.91093 -5.17215 -8.59794 -11.83868 -14.77750 -16.32184 -14.28557 -11.34011 -4.43520 -3.71277 -4.39736 -3.96716 -2.74612 -0.60181 2.40678 相对误差 -23.57% 7.42% 7.60% 3.25% -3.13% -5.06% -8.07% -7.16% -5.70% -11.22% -16.97% -16.78% -28.98% -41.03% -51.86% -54.96% -41.23% -27.87% -8.71% -6.73% -7.56% -6.34% -4.04% -0.81% 2.92% 对2011至2015年区域旅游人数预测如下：

表5.2.3.2 据多元线性回归模型y值预测结果 年份 2011 2012 2013 2014 2015

碳排放量/亿吨

85.073 90.4646 96.1978 102.2945 108.7775

◆5.3 BP神经网络预测模型 ▪5.3.1 模型思路

由于神经网络能很好识别各相关参数之间的非线性关系，这符合影响中国碳排放量影响因素的特点，故下文将用3层BP神经网络模型来分析和预测我国碳排放量。

首先根据我国碳排放量的需求趋势，确定输出层、中间隐层和输入层：输出层的输出节点为碳排放量y；用试凑法确定中间隐层节点数有数个，逐个做试验,通过BP网络训练结果比较,选取网络训练误差和网络训练次数组合最优所对应的隐层神经元数目；将碳排放量的影响因素作为确定对应输入层节点：x1人口总量，x2城镇化，x3人均GDP, x4第三产业GDP比例, x5能源强度吨标准煤,

x6煤炭消费比例。

然后把样本分为训练样本和测试样本两个部分，其中测试样本用于对所建的网络模型的检验和测试：将我国1986年至2000年的各因素数据作为训练样本,2001年至2010年的数据作为测试样本。

在以上基础,对样本数据进行归一化预处理，采用Matlab软件中的神经网络计算功能,结合GM（1,1）预测模型对2011至2015年各单变量预测结果，建立合理训练模型进行数据测试和预测,实现BP神经网络模型对碳排放的预测。

▪5.3.2 模型建立

*5.3.2. 1 确定输入层、输出层

经分析确定输入层节点数为6,体现为：x1人口总量，x2城镇化，x3人均GDP, x4第三产业GDP比例, x5能源强度吨标准煤, x6煤炭消费比例。 *5.3.2. 2数据预处理

将1986-2010年的数据用Matlab软件编程预测出2001-2010年的数据，把它与已知数据进行比较。 *5.3.2. 3 神经网络的训练

将所得资料数据从1986年至2000年的各因素数据作为训练样本,2001年至2010年的数据作为测试样本；

网络误差计算：网络模型输出值与已知训练、测试样本输出值之差平方和开根号，计算公式如下

errorTO2(xxii)； i1n*5.3.2.4 Bp模型预测

对测试样本结果进行误差分析，然后结合GM（1,1）预测模型对2011至2015年六项因素预测结果，对2011～2015年碳排放量进行预测。 ▪5.3.3 模型求解

表5.3.3.1 2001至2010 年y值Bp网络训练结果

年份 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010

y原始数据 29.69576 34.64843 40.69239 50.8978 55.12703 58.17144 62.56704 68.00468 74.63289 82.40958

y模拟值

29.6765 38.6113 42.2695 46.2737 50.6586 55.4582 60.7126 66.4648 72.762 79.6559

绝对误差

0.01926 -3.96287 -1.57711 4.6241 4.46843 2.71324 1.85444 1.53988 1.87089 2.75368

相对误差

0.06% -11.44% -3.88% 9.09% 8.11% 4.66% 2.96% 2.26% 2.51% 3.34%

Bp模型训练模拟曲线如下：

90检验期望数值BP神经网络预测检验值 80706050403020 2001200220032004200520062007200820092010图2 Bp模型训y值模拟曲线

对2011至2015年区域旅游人数预测如下：

表9 据Bp网络模型y值预测结果

Y 碳排放量/亿吨 年份

2011

2012 2013 2014 2015

◆5.4 因素排序模型

▪5.4.1 模型思路

应用关联度分析法建立因素排序模型，对数据进行量化比较分析，将数据代入关联系数公式，得出影响因素数列对参考数列在每个年份的关联系数，从而得到关联度，按关联度大小排序即可得影响因素的重要性排序。 ▪5.4.2 模型建立 选取参考数列：

87.2029 95.4649 104.5097 114.4115 125.2514

yy(k)k1,2,,n(y(1),y(2),y(n))，

x0y其中k表示时刻；

假设有m个比较数列：

xi{xi(k)|k1,2,...,n}(xi(1),xi(2),...,xi(n)),i1,2,...,m 数列xi对参考数列x0在k时刻的关联系数： i(k)minmin|x0(t)xs(t)|maxmax|x0(t)xs(t)|stst|x0(k)xi(k)|maxmax|x0(t)xs(t)|stst...........(1)

其中minmin|x0(t)xs(t)|与maxmax|x0(t)xs(t)|分别为两级最小差及两级最大

st差，分辨系数越大，分辨率越大，越小，分辨率越小。 根据上述可得数列xi对参考数列x0的关联度为：

1n rii(k).............................................(2)

nk1利用（2）公式，我们可以对各种问题进行因素分析，得到因素重要性排序。 ▪5.4.3 模型求解

利用Matlab编程实现上述语句，代入1986至2010年各自变量及测度指标有关数据，得影响因素的重要性排序：人均GDP>人口>煤炭消费比例>城镇化>能源强度比例>第三产业GDP比例，具体结果如下：

表5.4.3.1 因素重要性排序

因素 关联度 排序

x3人均x4第三产业x5能源强度x6煤炭消费

x1人口 x2城镇化 GDP GDP比例 吨标准煤 比例 0.9493 0.9662 0.7116 0.9564 0.918 0.943 2 4 1 6 5 3

6、模型检验分析

本文将利用绝对误差与相对误差两个参数对GM（1,1）预测模型、多元线性回归预测模型、BP神经网络预测模型3个模型进行检验，同时对其精确度进行对比评价。由表5.1.3.1，表5.2.3.1，表5.3.3.1三表对比其中的绝对误差和相对误差知BP神经网络预测模型误差最小故其预测合理性最高，预测效果最佳。

经模型求解可得相关数据如下：

表6.1 3个模型y值预测结果

2011 2012 2013 2014 2015 模型 GM(1,1)模型 线性回归模型

77.8641 85.073

83.4852 90.4646

89.5121 96.1978

95.9741 102.2945

102.9026 108.7775

Bp网络模型

87.2029 95.4649 104.5097 114.4115 125.2514

7、建议

根据影响因素的重要性排序：人均GDP>人口>煤炭消费比例>城镇化>能源强度比例>第三产业GDP比例，结合我国碳排放量变化趋势提出一些小建议： 1、国家适当调整人均GDP比率； 2、提倡计划生育，降低人口增长率；

3、合理开采煤炭，控制每年煤炭销售量，尽可能地使用氢能源、太阳能等无污染能源；

4、鼓励农村发展，先富带后富，城镇化比例需适中； 5、坚持走可持续发展道路，寻找开发新能源。

8、模型评价与推广

♦8.1模型优点

（1）3个预测模型总体来说，预测精度较高，所得结果预测结果真实性较高；

（2）BP神经网络预测模型,克服了一般回归分析的不够精确的局限性，使得结果更为精确，故其最适应碳排放量预测等复杂性较高的小样本型问题；

（3）用Matlab实现BP神经网络模型,使用方便且计算精确；

（4）使用Matlab统计软件完成多元线性回归的函数计算，简洁快速，值得推广；

（5）在几处采用了多个预测模型相结合的方法进行预测，使得模型进一步优化。 ♦8.2模型缺点

（1）GM (1, 1)模型适用于短期预测，对于近一两年的预测值很精确，而远期的数仅仅反应一种趋势；

（2）使用趋势移动平均预测模型存在很大的局限性，因为它只适用于时间序列呈直线趋势的情况。 ♦8.3模型推广

（1） 本文的结论是建立在数学模型基础上，有较高的可信度，故结论可以为

有关部门控制碳排放量提供科学依据.。

（2） 文中所用的3种预测模型精确度较高，预测效果较好，尤其是BP神经

网络预测模型，可推广应用到实际生活待预测的非线性问题上。

（3）建立的BP 神经网络预测模型，可以不断添加历史数据，进一步提高预测能力，提高网络的可信度。

参考文献

[1] 姜启源，谢金星著，数学建模案例选集.北京：高等教育出版社，2006. [2] 飞思科技产品研发中心，神经网络理论与Matlab7实现[M].北京：电子工业出版社，2005.

[3] 应玖茜，魏权龄，非线性规划及其理论.北京：中国人民大学出版社，1994.

附录

用Matlab编程如下：

附件一：GM(1,1)预测模型 clear;clc;

k=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25];

g=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30];

p=[19.70823,21.0278,22.40368,22.75338,22.69709,23.69252,24.49162,26.26645,28.31547,28.61658,28.93377,30.81745,29.67256,28.85722,28.4975,29.69567,34.64843,40.69239,50.8978,55.12703,58.17144,62.56704,68.00468,74.63289,82.40958]; sum=0; for j=1:25

q(j)=sum+p(j); sum=sum+p(j); end

y=zeros(24,1); y=(p(:,2:25))'; for i=2:25

z(1)=q(1);

z(i)=(0.5*q(i)+0.5*q(i-1)); end

b=zeros(24,2);

b(:,1)=(-1*z(:,2:25))'; b(:,2)=ones(24,1); c=zeros(2,1); c=inv(b'*b)*b'*y; for j=1:30

n(j)=(p(1)-c(2)/c(1))*exp(-c(1)*(j-1))+c(2)/c(1); end

for i=2:30

h(1)=z(1); h(i)=n(i)-n(i-1) end

plot(k,p,'r',g,h,'b')

legend('原始曲线','GM(1,1)模拟曲线')

运行结果： h =

Columns 1 through 6

19.7082 14.6152 15.6703 16.8015 18.0144 19.3149 Columns 7 through 12

20.7093 22.2043 23.8073 25.5259 27.3687 29.3445 Columns 13 through 18

31.4629 33.7342 36.1695 38.7807 41.5803 44.5820 Columns 19 through 24

47.8004 51.2512 54.9511 58.9181 63.1714 67.7318 Columns 25 through 30

72.6215 77.8641 83.4852 89.5121 95.9741 102.9026

附件二：多元线性回归预测模型 （1）多元线性方程求解代码：

>>x1=[107507,109300,111026,112704,114333,115823,117171,118517,119850,121121,122389,123626,124761,125786,126743,127627,128453,129227,129988,130756,131448,132129,132802,133450,134091];

x2=[24.52,25.32,25.81,26.21,26.41,26.94,27.46,27.99,28.51,29.04,30.48,31.91,33.35,34.78,36.22,37.66,39.09,40.53,41.76,42.99,44.34,45.89,46.99,48.34,49.95];

x3=[963,1112,1366,1519,1644,1893,2311,2998,4044,5046,5846,6420,6796,7159,7858,8622,9398,10542,12336,14185,16500,20169,23708,25608,29992];

x4=[29.14,29.64,30.51,32.06,31.54,33.69,34.76,33.72,33.57,32.86,32.77,34.17,36.23,37.77,39.02,40.46,41.47,41.23,40.38,40.51,40.94,41.89,41.82,43.43,43.14];

x5=[9.8,9.4,9.1,9.1,8.9,8.6,7.9,7.4,6.9,6.6,6.2,5.7,5.3,5.1,4.8,4.6,4.5,4.7,5,4.9,4.8,4.5,4.3,4.2,4.2];

x6=[75.8,76.2,76.1,76.1,76.2,76.1,75.7,74.7,75,74.6,73.5,71.4,70.9,70.6,69.2,68.3,68,69.8,69.5,70.8,71.1,71.1,70.3,70.4,68];

y=[19.70823,21.0278,22.40368,22.75338,22.69709,23.69252,24.49162,26.26645,28.31547,28.61685,28.93377,30.81745,29.67256,28.85722,28.4975,29.69576,34.64843,40.69239,50.8978,55.12703,58.17144,62.56704,68.00468,74.63289,82.40958]; [uuiijt]=regress(y',[ones(25,1) x1' x2' x3' x4' x5' x6'])

运行结果： uuiijt =

-246.3984 0.0016 1.1399 0.0015 -0.9702 5.8205 0.5353

（2）多元线性回归方程1986-2010年模拟预测代码

a=[107507 24.52 963 29.14 9.8 75.8 109300 25.32 1112 29.64 9.4 76.2 111026 25.81 1366 30.51 9.1 76.1 112704 26.21 1519 32.06 9.1 76.1 114333 26.41 1644 31.54 8.9 76.2

115823 26.94 1893 33.69 8.6 76.1 117171 27.46 2311 34.76 7.9 75.7 118517 27.99 2998 33.72 7.4 74.7 119850 28.51 4044 33.57 6.9 75 121121 29.04 5046 32.86 6.6 74.6 122389 30.48 5846 32.77 6.2 73.5 123626 31.91 6420 34.17 5.7 71.4 124761 33.35 6796 36.23 5.3 70.9 125786 34.78 7159 37.77 5.1 70.6 126743 36.22 7858 39.02 4.8 69.2 127627 37.66 8622 40.46 4.6 68.3 128453 39.09 9398 41.47 4.5 68 129227 40.53 10542 41.23 4.7 69.8 129988 41.76 12336 40.38 5 69.5 130756 42.99 14185 40.51 4.9 70.8 131448 44.34 16500 40.94 4.8 71.1 132129 45.89 20169 41.89 4.5 71.1 132802 46.99 23708 41.82 4.3 70.3 133450 48.34 25608 43.43 4.2 70.4 134091 49.95 29992 43.14 4.2 68]; for i=1:25

y(i)=0.0016*a(i,1)+1.1399*a(i,2)+0.0015*a(i,3)-0.9702*a(i,4)+5.8205*a(i,5)+0.5353*a(i,6)-246.3984 end

运行结果： y =

Columns 1 through 6

24.3527 25.7577 26.8151 28.6815 31.0974 30.5734 Columns 7 through 12

28.6234 29.9751 31.6655 34.5348 36.5753 35.6529 Columns 13 through 18

35.0799 36.0757 36.5885 37.7474 40.1405 47.0969 Columns 19 through 24

54.8178 60.2099 65.4898 71.1819 77.2967 80.6318 Column 25 89.06532

（3）多元线性回归方程1986-2015年数据预测

k=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25];

g=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30];

p=[24.3527,25.7577,26.8151,28.6815,31.0974,30.5734,28.6234,29.9751,31.6655,34.5348,36.5753,35.6529,35.0799,36.0757,36.5885,37.7474,40.1405,47.0969,54.8178,6

0.2099,65.4898,71.1819,77.2967,80.6318,89.06532]; sum=0; for j=1:25

q(j)=sum+p(j); sum=sum+p(j); end

y=zeros(24,1); y=(p(:,2:25))'; for i=2:25

z(1)=q(1);

z(i)=(0.5*q(i)+0.5*q(i-1)); end

b=zeros(24,2);

b(:,1)=(-1*z(:,2:25))'; b(:,2)=ones(24,1); c=zeros(2,1); c=inv(b'*b)*b'*y; for j=1:30

n(j)=(p(1)-c(2)/c(1))*exp(-c(1)*(j-1))+c(2)/c(1); end

for i=2:30

h(1)=z(1); h(i)=n(i)-n(i-1) end

plot(k,p,'r',g,h,'b')

legend('原始曲线','GM(1,1)模拟曲线')

运行结果： h =

Columns 1 through 6

24.3527 19.4675 20.7012 22.0132 23.4083 Columns 7 through 12

26.4693 28.1469 29.9307 31.8276 33.8447 Columns 13 through 18

38.2705 40.6959 43.2750 46.0176 48.9340 Columns 19 through 24

55.3330 58.8398 62.5688 66.5342 70.7508 Columns 25 through 30

80.0028 85.0730 90.4646 96.1978 102.2945

附件三：BP神经网络预测模型 （1）2001-2010年模拟预测程序 clc;clear all,close all

XX=[107507 24.52 963 29.14 9.8 75.8 19.70823

24.8918 35.9896 52.0352 75.2347 108.7775 109300 25.32 1112 29.64 9.4 76.2 21.02780 111026 25.81 1366 30.51 9.1 76.1 22.40368 112704 26.21 1519 32.06 9.1 76.1 22.75338 114333 26.41 1644 31.54 8.9 76.2 22.69709 115823 26.94 1893 33.69 8.6 76.1 23.69252 117171 27.46 2311 34.76 7.9 75.7 24.49162 118517 27.99 2998 33.72 7.4 74.7 26.26645 119850 28.51 4044 33.57 6.9 75.0 28.31547 121121 29.04 5046 32.86 6.6 74.6 28.61685 122389 30.48 5846 32.77 6.2 73.5 28.93377 123626 31.91 6420 34.17 5.7 71.4 30.81745 124761 33.35 6796 36.23 5.3 70.9 29.67256 125786 34.78 7159 37.77 5.1 70.6 28.85722 126743 36.22 7858 39.02 4.8 69.2 28.49750 127627 37.66 8622 40.46 4.6 68.3 29.69576 128453 39.09 9398 41.47 4.5 68.0 34.64843 129227 40.53 10542 41.23 4.7 69.8 40.69239 129988 41.76 12336 40.38 5.0 69.5 50.89780 130756 42.99 14185 40.51 4.9 70.8 55.12703 131448 44.34 16500 40.94 4.8 71.1 58.17144 132129 45.8 20169 41.89 4.5 71.1 62.56704 132802 46.99 23708 41.82 4.3 70.3 68.00468 133450 48.34 25608 43.4 4.2 70.4 74.63289 134091 49.95 29992 43.14 4.2 68.0 82.40958]; [u,v]=size(XX); for i=1:u for j=1:v

X0(i,j)=(XX(i,j)-min(XX(:,j)))/(max(XX(:,j))-min(XX(:,j))); end end X0

P_train=X0(1:u-9,1:v-1)'; T_train=X0(1:u-9,v)';

net=newff(minmax(P_train),[20,1],{'tansig' 'logsig'}, 'traingdx'); net.trainParam.epochs=2000; net.trainParam.goal=0.00001; net=train(net,P_train,T_train); y=sim(net,P_train)

y=y*(max(XX(:,v))-min(XX(:,v)))+min(XX(:,v))

P_test=X0(u-9:u,1:v-1)'; O_test=sim(net,P_test);

O_test=O_test*(max(XX(:,v))-min(XX(:,v)))+min(XX(:,v)) x1=2001:2010;

f=[29.69576 34.64843 40.69239 50.8978 55.12703 58.17144 62.56704 68.00468 74.63289 82.40958]; g=O_test;

plot(x1,f,'r+',x1,g,'b*')

legend('检验期望数值 ','BP神经网络预测检验值')

运行结果： y =

Columns 1 through 6

0.0101 0.0205 0.0383 0.0488 0.0505 0.0604 Columns 7 through 12

0.0791 0.1037 0.1362 0.1441 0.1459 0.1771 Columns 13 through 16

0.1598 0.1447 0.1411 0.1590 y =

Columns 1 through 6

20.3386 20.9964 22.1097 22.7712 22.8718 23.4944 Columns 7 through 12

24.6689 26.2084 28.2502 28.7430 28.8572 30.8107 Columns 13 through 16

29.7297 28.7795 28.5561 29.6765 O_test =

Columns 1 through 6

29.6765 32.5109 44.5326 42.6808 51.0667 57.3003 Columns 7 through 10

66.5810 70.4552 72.1191 74.0750

（2）2001-2015年模拟预测程序 clear;clc;

k=[1,2,3,4,5,6,7,8,9,10];

g=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15];

p=[29.6765,32.5109,44.5326,42.6808,51.0667,57.3003,66.5810,70.4552,72.1191,74.0750]; sum=0; for j=1:10

q(j)=sum+p(j); sum=sum+p(j); end

y=zeros(9,1); y=(p(:,2:10))'; for i=2:10

z(1)=q(1);

z(i)=(0.5*q(i)+0.5*q(i-1)); end

b=zeros(9,2);

b(:,1)=(-1*z(:,2:10))'; b(:,2)=ones(9,1); c=zeros(2,1); c=inv(b'*b)*b'*y; for j=1:15

n(j)=(p(1)-c(2)/c(1))*exp(-c(1)*(j-1))+c(2)/c(1); end

for i=2:15

h(1)=z(1); h(i)=n(i)-n(i-1) end

plot(k,p,'r',g,h,'b')

legend('原始曲线','模拟曲线') 运行结果： h =

Columns 1 through 7

29.6765 38.6113 42.2695 46.2743 50.6586 60.7126

Columns 8 through 14

66.4648 72.7620 79.6559 87.2029 95.4649 114.4115 Column 15 125.2514

附件四：因素排序模型 clc,clear

x=[19.70823 107507 24.52 963 29.14 9.8 75.8 21.0278 109300 25.32 1112 29.64 9.4 76.2

22.40368 111026 25.81 1366 30.51 9.1 76.1 22.75338 112704 26.21 1519 32.06 9.1 76.1 22.69709 114333 26.41 1644 31.54 8.9 76.2 23.69252 115823 26.94 1893 33.69 8.6 76.1 24.49162 117171 27.46 2311 34.76 7.9 75.7 26.26645 118517 27.99 2998 33.72 7.4 74.7 28.31547 119850 28.51 4044 33.57 6.9 75 28.61685 121121 29.04 5046 32.86 6.6 74.6 28.93377 122389 30.48 5846 32.77 6.2 73.5 30.81745 123626 31.91 6420 34.17 5.7 71.4 29.67256 124761 33.35 6796 36.23 5.3 70.9 28.85722 125786 34.78 7159 37.77 5.1 70.6 28.4975 126743 36.22 7858 39.02 4.8 69.2

29.69576 127627 37.66 8622 40.46 4.6 68.3 34.64843 128453 39.09 9398 41.47 4.5 68 40.69239 129227 40.53 10542 41.23 4.7 69.8

55.4582 104.5097 50.8978 129988 41.76 12336 40.38 5 69.5 55.12703 130756 42.99 14185 40.51 4.9 70.8 58.17144 131448 44.34 16500 40.94 4.8 71.1 62.56704 132129 45.89 20169 41.89 4.5 71.1 68.00468 132802 46.99 23708 41.82 4.3 70.3 74.63289 133450 48.34 25608 43.43 4.2 70.4 82.40958 134091 49.95 29992 43.14 4.2 68]'; for i=1:6

x(i,:)=x(i,:)/x(i,1); %标准化数据 end for i=7

x(i,:)=x(i,1)./x(i,:); %标准化数据 end data=x;

n=size(data,2); %求矩阵的列数，即观测时刻的个数 ck=data(1,:); %提出参考数列 bj=data(2:end,:); %提出比较数列

m2=size(bj,1); %求出比较数列的个数 for j=1:m2 t(j,:)=bj(j,:)-ck; end

mn=min(min(abs(t'))); %求最小差 mx=max(max(abs(t'))); %求最大差 rho=0.5; %分辨系数设置

ksi=(mn+rho*mx)./(abs(t)+rho*mx); %求关联系数 r=sum(ksi')/n %求关联度

[rs,rind]=sort(r,'descend') %对关联度进行排序

运行结果： r =

0.9493 0.9662 0.7116 0.9564 0.9180 0.9430 rs =

0.9662 0.9564 0.9493 0.9430 0.9180 0.7116 rind =

2 4 1 6 5 3