我国碳排放量预测模型

我国碳排放量预测模型
摘要
本文主要我国碳排放预测问题，同时根据预测结果提出合理性建议。

以人口总量，城镇化，人均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年预测数据如下表。

GM(1,1)模型77.8641 83.4852 89.5121 95.9741 102.9026 线性回归模
85.073 90.4646 96.1978 102.2945 108.7775 型
Bp网络模型87.2029 95.4649 104.5097 114.4115 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、问题分析
本文整个过程主要解决问题是我国碳排放预测问题，通过分析确定测度指标
x，城镇化2x，人均为碳排放量y，相关影响因素为以下六个方面：人口总量
1
x,第三产业GDP比例4x,能源强度吨标准煤5x,煤炭消费比例6x。

GDP
3
本文将建立建立GM（1,1）预测模型、多元线性回归预测模型、BP神经网络预测模型，借助Matlab软件逐个对碳排放量进行模拟与预测。

3种定量预测模型各有各的优势与不足，故在几处采用了多个预测模型相结合的方法进行预测，使得模型进一步优化。

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

3、模型假设
1、假设统计的数据真实科学，短期内稳定变化；
2、假设建立模型中，个别偏离太远的数据可据题适当调整；
3、假设碳排放量变化主要受人口总量，城镇化，人均GDP,第三产业GDP 比例,能源强度吨标准煤,煤炭消费比例六方面因素的影响；
4、假设碳排放量y值大部分呈线性，满足多元线性回归条件；
5、假设影响碳排放量变化的各因素之间存在相关关系、各因素与碳排放量存在非线性关系，适用BP神经网络的性能。

4、符号说明
i X ：原始数据序列
a ：发展恢数
b ：内生控制恢数
i x ：自变量()5,2,1 =i
y ：测度指标
i β：参系数
T i x ：样本实际输入值
O i x ：网络模型输出值
)(k i ξ：比较数列i x 对参考数列0x 在k 时刻的关联系数
ρ：分辨系数，且[0,1]ρ∈
i r ：数列i x 对参考数列0x 的关联度
5、模型建立与求解
本文探讨的是中国碳排放量的预测模型，根据中国近25年的碳排放量变化趋势进行解题。

本文确定测度指标为碳排放量，相关影响因素为以下六方面：人口总量1x ，城镇化2x ，人均GDP 3x ,第三产业GDP 比例4x ,能源强度吨标准煤5x ,煤炭消费比例6x 。

根据2011年全国年鉴统计表得到有关真实数据如下：
表1 1986至2010年各自变量及测度指标有关数据
年份 人口/（万人） 城镇化/% 人均GDP/（元/人） 第三产业GDP 比例能源强度吨
标准煤/（万
煤炭消费比例/% 碳排放量/
亿吨
198610750724.5296329.149.875.819.708231987 109300 25.32 1112 29.64 9.4 76.2 21.0278 1988 111026 25.81 1366 30.51 9.1 76.1 22.40368 1989 112704 26.21 1519 32.06 9.1 76.1 22.75338 1990
114333
26.41
1644
31.54
8.9
76.2
22.69709
199111582326.94189333.698.676.123.692521992 117171 27.46 2311 34.76 7.9 75.7 24.49162 1993 118517 27.99 2998 33.72 7.4 74.7 26.26645 1994 119850 28.51 4044 33.57 6.9 75 28.31547 1995 121121 29.04 5046 32.86 6.6 74.6 28.61685 1996 122389 30.48 5846 32.77 6.2 73.5 28.93377 1997 123626 31.91 6420 34.17 5.7 71.4 30.81745 1998 124761 33.35 6796 36.23 5.3 70.9 29.67256 1999 125786 34.78 7159 37.77 5.1 70.6 28.85722 2000 126743 36.22 7858 39.02 4.8 69.2 28.4975 2001 127627 37.66 8622 40.46 4.6 68.3 29.69576 2002 128453 39.09 9398 41.47 4.5 68 34.64843 2003 129227 40.53 10542 41.23 4.7 69.8 40.69239 2004 129988 41.76 12336 40.38 5 69.5 50.8978 2005 130756 42.99 14185 40.51 4.9 70.8 55.12703 2006 131448 44.34 16500 40.94 4.8 71.1 58.17144 2007 132129 45.89 20169 41.89 4.5 71.1 62.56704 2008 132802 46.99 23708 41.82 4.3 70.3 68.00468 2009 133450 48.34 25608 43.43 4.2 70.4 74.63289 2010 134091 49.95 29992 43.14 4.2 68 82.40958
◆5.1 GM （1,1）预测模型
▪5.1.1 模型思路
建立GM (1,1)灰色动态预测模型，通过对1986至2010年原始单变量数据进行生成处理，包括累加、转换、代入和还原等过程，以寻找系统的变化规律，在此基础上建立相应的微分方程预测模型。

代入相关单变量数据用Matlab 编程得到对应变量在2011至2015年的预测值。

▪5.1.2 模型建立
第一步：假设X 0为原始非负序列：
()()},,,{,02,0)1,0(0n x x x X =，
其中()0,01,2,
,25k x k ≥=，。

利用累加生成序列可将序列()()},,,{,02,0)1,0(0n x x x X =生成序列1X ：
()()},,,{,12,1)1,1(1n x x x X =，
其中()(
)
1,0,1
1,2,,25k
k i i x x k ==
=∑，。

第二步：利用先生成的序列1X 建立GM （1,1）模型的一般形式：
()()b z a x k k =*+,1,0， （1）
用微分方程表示如下：
b ax dt
dx =+11
， （2） Z 1为X 1的紧邻均值生成序列：
()()},,,{,12,1)1,1(1n z z z Z =，
其中
()()()1,1,1,10.5*0.5*1,2,
,25k k k z x x k -=+=，， （3）
灰微分方程模型中a ，b 为待估参数，分别为发展恢数和内生控制恢数；灰微分方程的最小二乘估计参数列满足： ()()
Y B B
B b a T T T
***==-1
,σ， （4）
Y 为列向量，B 为构造矩阵
()()()⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎣⎡=n x x x Y ,03,02,0 ，()()()()()()1,21,11,31,21,1,10.5*0.5*10.5*0.5*10.5*0.5*1n n B x x x x x x ---⎡⎤
⎢⎥--⎢⎥=⎢⎥⎢⎥⎢⎥
--⎢⎥⎣⎦， 第三步：构建灰色预测模型
结合上述（1）、（3）、（4）求解（2）式得：
()()n k a
b
e a b x x ak k ,,2,1)(0,11,1 =+*-=-+，，
因为()()1,00,1x x =，所以建立的灰色模型如下
()()n k a
b
e a b x x ak k ,,2,1)(1,01,1 =+*-=-+，，
第四步：求出原始数据的还原值即可得到GM （1,1）预测模型：
()()()n k x x x k k k ,,2,1,11,11,0 =-=++，。

第五步：对模型进行后检验模型的精确性
▪5.1.3 模型求解
表一中各自变量构成原始非负序列，代入Matlab 软件得GM （1,1）对y 值模拟结果如下：
表5.1.3.1 1986至2010年y值GM（1,1）模型模拟结果
年份y原始数据y模拟值绝对误差相对误差
198619.7082319.70820.00003 0.00% 198721.027814.6152 6.412630.50% 198822.4036815.6703 6.7333830.05% 198922.7533816.8015 5.9518826.16% 199022.6970918.0144 4.6826920.63% 199123.6925219.3149 4.3776218.48% 199224.4916220.7093 3.7823215.44% 199326.2664522.2043 4.0621515.47% 199428.3154723.8073 4.5081715.92% 199528.6168525.5259 3.0909510.80% 199628.9337727.3687 1.56507 5.41% 199730.8174529.3445 1.47295 4.78% 199829.6725631.4629-1.79034-6.03% 199928.8572233.7342-4.87698-16.90% 200028.497536.1695-7.672-26.92% 200129.6957638.7807-9.08494-30.59% 200234.6484341.5803-6.93187-20.01% 200340.6923944.582-3.88961-9.56% 200450.897847.8004 3.0974 6.09% 2000555.1270351.2512 3.875837.03% 200658.1714454.9511 3.22034 5.54% 200762.5670458.9181 3.64894 5.83% 200868.0046863.1714 4.833287.11% 200974.6328967.7318 6.901099.25% 201082.4095872.62159.7880811.88%
模拟曲线如下：
102030405060708090100110
图5.1.3.2 GM （1,1）模型y 值模拟曲线
同时得到2011至2015年相关预测结果如下表：
表5.1.3.3 据GM （1,1）模型y 值预测结果 年份 碳排放量/亿吨
2011 77.8641 2012 83.4852 2013 89.5121 2014 95.9741 2015
102.9026
◆5.2 多元线性回归预测模型
▪5.2.1 模型思路
由于1986至2010年各影响因素及测度指标有关数据呈现一定的递增形式，故此采用多元线性回归模型对其进行模拟与预测。

确定线性预测变量和因变量，建立回归矩阵。

预测变量即六个影响因素：1x 人口总量，2x 城镇化，3x 人均GDP, 4x 第三产业GDP 比例, 5x 能源强度吨标准煤, 6x 煤炭消费比例。

将实际数据用Matlab 软件求解出符合该情况的多元线性方程，将1986至2010年所有数据代入该方程，同时结合GM （1,1）预测模型对2011至2015年
各单变量预测结果，用Matlab 编程得到对应变量的y 的模拟值和预测值。

▪5.2.2 模型建立
设回归模型为：
0115566y x x x ββββ=++
++，
据题可导：
0115566i i i i y x x x ββββ=++++
▪5.2.3 模型求解
将实际数据用Matlab 软件求解出多元线性方程为：
1234560.0016 1.13990.00150.9702 5.82050.5353246.3984y x x x x x x =++-++-；
借助Matlab 编程得到对应变量y 的模拟值：
表5.2.3.1 1986至2010 年y 值线性回归模拟结果
年份 y 原始数据 y 模拟值 绝对误差 相对误差 198619.7082324.3527-4.64447 -23.57%1987 21.0278 19.4675 1.56030 7.42% 1988 22.40368 20.7012 1.70248 7.60% 1989 22.75338 22.0132 0.74018 3.25% 1990 22.69709 23.4083 -0.71121 -3.13% 1991 23.69252 24.8918 -1.19928 -5.06% 1992 24.49162 26.4693 -1.97768 -8.07% 1993 26.26645 28.1469 -1.88045 -7.16% 1994 28.31547 29.9307 -1.61523 -5.70% 1995 28.61685 31.8276 -3.21075 -11.22% 1996 28.93377 33.8447 -4.91093 -16.97% 1997 30.81745 35.9896 -5.17215 -16.78% 1998 29.67256 38.2705 -8.59794 -28.98% 1999 28.85722 40.6959 -11.83868 -41.03% 2000 28.4975 43.275 -14.77750 -51.86% 2001 29.69576 46.0176 -16.32184 -54.96% 2002 34.64843 48.934 -14.28557 -41.23% 2003 40.69239 52.0325 -11.34011 -27.87% 2004 50.8978 55.333 -4.43520 -8.71% 20005 55.12703 58.8398 -3.71277 -6.73% 2006 58.17144 62.5688 -4.39736 -7.56% 2007 62.56704 66.5342 -3.96716 -6.34% 2008 68.00468 70.7508 -2.74612 -4.04% 2009 74.63289 75.2347 -0.60181 -0.81%
对2011至2015年区域旅游人数预测如下：
表5.2.3.2 据多元线性回归模型y值预测结果
年份碳排放量/亿吨
2011 85.073
2012 90.4646
2013 96.1978
2014 102.2945
2015 108.7775
◆5.3 BP神经网络预测模型
▪5.3.1 模型思路
由于神经网络能很好识别各相关参数之间的非线性关系，这符合影响中国碳排放量影响因素的特点，故下文将用3层BP神经网络模型来分析和预测我国碳排放量。

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

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

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

▪5.3.2 模型建立
*5.3.2. 1 确定输入层、输出层
x人口总量，2x城镇化，3x人均经分析确定输入层节点数为6,体现为：
1
x第三产业GDP比例, 5x能源强度吨标准煤, 6x煤炭消费比例。

GDP,
4
*5.3.2. 2数据预处理
将1986-2010年的数据用Matlab软件编程预测出2001-2010年的数据，把它与已知数据进行比较。

*5.3.2. 3 神经网络的训练
将所得资料数据从1986年至2000年的各因素数据作为训练样本,2001年至2010年的数据作为测试样本；
网络误差计算：网络模型输出值与已知训练、测试样本输出值之差平方和开根号，计算公式如下
∑=-
=
n
i
O
i T
i
x x
error
1
2
) (；
*5.3.2.4 Bp模型预测
对测试样本结果进行误差分析，然后结合GM（1,1）预测模型对2011至2015年六项因素预测结果，对2011～2015年碳排放量进行预测。

▪5.3.3 模型求解
表5.3.3.1 2001至2010年y值Bp网络训练结果
年份y原始数据y模拟值绝对误差相对误差
200129.6957629.67650.019260.06% 200234.6484338.6113-3.96287-11.44% 200340.6923942.2695-1.57711-3.88% 200450.897846.2737 4.62419.09% 200555.1270350.6586 4.468438.11% 200658.1714455.4582 2.71324 4.66% 200762.5670460.7126 1.85444 2.96% 200868.0046866.4648 1.53988 2.26% 200974.6328972.762 1.87089 2.51% 201082.4095879.6559 2.75368 3.34%
Bp模型训练模拟曲线如下：
图2 Bp 模型训y 值模拟曲线
对2011至2015年区域旅游人数预测如下：
表9 据Bp 网络模型y 值预测结果
年份
Y 碳排放量/亿吨 2011
87.2029 2012 95.4649 2013 104.5097 2014 114.4115 2015
125.2514
◆5.4 因素排序模型
▪5.4.1 模型思路
应用关联度分析法建立因素排序模型，对数据进行量化比较分析，将数据代入关联系数公式，得出影响因素数列对参考数列在每个年份的关联系数，从而得到关联度，按关联度大小排序即可得影响因素的重要性排序。

▪5.4.2 模型建立 选取参考数列：
{}y
x n y y y n k k y y ====0))(),2(),1((,,2,1)(，
其中k 表示时刻；
假设有m 个比较数列：
{()|1,2,...,}((1),(2),...,()),1,2,...,i i i i i x x k k n x x x n i m ==== 数列i x 对参考数列0x 在k 时刻的关联系数： 0000min min |()()|max max |()()|
()...........(1)|()()|max max |()()|
s s s
t
s
t
i i s s
t
x t x t x t x t k x k x k x t x t ρξρ-+-=
-+-
其中0min min |()()|s s
t
x t x t -与0max max |()()|s s
t
x t x t -分别为两级最小差及两级最大
差，分辨系数ρ越大，分辨率越大，ρ越小，分辨率越小。

根据上述可得数列i x 对参考数列0x 的关联度为：
1
1().............................................(2)n
i i k r k n ξ==∑
利用（2）公式，我们可以对各种问题进行因素分析，得到因素重要性排序。

▪5.4.3 模型求解
利用Matlab 编程实现上述语句，代入1986至2010年各自变量及测度指标有关数据，得影响因素的重要性排序：人均GDP>人口>煤炭消费比例>城镇化>能源强度比例>第三产业GDP 比例，具体结果如下：
表5.4.3.1 因素重要性排序
因素 x1人口 x2城镇化 x3人均GDP x4第三产业GDP 比例 x5能源强度吨标准煤 x6煤炭消费
比例 关联度 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 值预测结果
GM(1,1)模型 77.8641 83.4852 89.5121 95.9741 102.9026 线性回归模型
85.073
90.4646
96.1978
102.2945
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,12112 1,122389,123626,124761,125786,126743,127627,128453,129227,129988,130756,13 1448,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,785 8,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 24.8918 Columns 7 through 12
26.4693 28.1469 29.9307 31.8276 33.8447 35.9896 Columns 13 through 18
38.2705 40.6959 43.2750 46.0176 48.9340 52.0352 Columns 19 through 24
55.3330 58.8398 62.5688 66.5342 70.7508 75.2347 Columns 25 through 30
80.0028 85.0730 90.4646 96.1978 102.2945 108.7775
附件三：BP神经网络预测模型
（1）2001-2010年模拟预测程序
clc;clear all,close all
XX=[107507 24.52 963 29.14 9.8 75.8 19.70823
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 55.4582 60.7126
Columns 8 through 14
66.4648 72.7620 79.6559 87.2029 95.4649 104.5097 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
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。