[负荷预测]基于灰色GM(1,1)模型的中长期电力负荷预测

article/2025/10/31 0:32:34

目录

一、灰色模型GM（1,1）原理

二、模型构建前检验

三、预测精度的检验

3.1 残差检验

3.2 后验差检验

四、灰色模型GM(1,1)算法

五、Matlab编程实现

5.1 程序代码

5.2 输出结果

5.3 结果检验

5.4 预测泛化

六、Python编程实现

6.1 程序代码

6.2 程序输出

6.3 往后 10 年负荷为:（2021-2030年）

一、灰色模型GM（1,1）原理

灰色预测的主要特点是模型使用的不是原始数据序列，而是生成的数据序列。其核心体系是灰色模型（Grey Model，简称GM），即对原始数据作累加生成（或其它方法生成）得到近似的指数规律再进行建模的方法。灰色预测模型对于不同问题采用不同模型，GM（1，1）模型主要解决生成序列是有指数变化规律，只能描述单调的变化过程。

对于生成的数列，我们可以设想用一个指数曲线乃至一条直线来逼近这个生成的数列。于是我们构建一个一阶常微分方程来求解这个拟合曲线函数表达式。设：

给定观测数据列

经一次累加得

设满足一阶常微分方程

灰色预测模型求解的途径是一次累加序列通过最小二乘法来估计常数a与u.

最小二乘法的作用是计算未知数个数小于方程个数的情况下，未知数的使误差最小的取值，这里用于求u和a。从最后的结果可以看出，GM（1，1）只要针对指数变化且单调规律进行拟合，对于摆动序列就不行。

二、模型构建前检验

给定序列能否建立较高精度的GM(1,1)模型，一般用序列的级比lamda的大小与所属区间来判断。序列的级比定义为：

若满足：

例如n=4,

则认为序列x(0)是可作为GM（1,1）建模的。

三、预测精度的检验

模型对原始数据列进行预测，与原始结果比较，看看拟合度怎么样。残差指的是预测值与实际值的差值，即绝对误差。

3.1 残差检验

分为残差检验、后验差检验和关联度检验。

3.2 后验差检验

一个模型是否满足实际的需要，就要对模型进行检验。一个好的模型，要求C（后验差比值）越小越好；P（小误差概率）越大越好。按照P和C的大小，可将精度分为4个等级，如表1所示。

表4.1 预测精度等级

详细原理请阅读：https://heaven.blog.csdn.net/article/details/120639012

四、灰色模型GM(1,1)算法

五、Matlab编程实现

5.1 程序代码

%GM(1,1)灰色理论模型,2021.8.1
clc
clear
syms a u;
c=[a,u]';%构成矩阵
A_SG= [30.9250   32.5390   35.2270   34.6940   34.5060   36.0510   38.7450   42.3610   44.5360 45.783 47.52];%输入数据，可以修改
A = A_SG; 
Ago=cumsum(A);%原始数据一次累加,得到1-AGO序列xi(1)。
n=length(A);%原始数据个数
for k=1:(n-1)Z(k)=(Ago(k)+Ago(k+1))/2; %Z(i)为xi(1)的紧邻均值生成序列
end
Yn =A;%Yn为常数项向量
Yn(1)=[]; %从第二个数开始，即x(2),x(3)...
Yn=Yn';
E=[-Z;ones(1,n-1)]';%累加生成数据做均值
c=(E'*E)\(E'*Yn);%利用公式求出a，u
c= c';
a=c(1);%得到a的值
u=c(2);%得到u的值
F=[];
F(1)=A(1);
%n=n+1
for k=2:(n)F(k)=(A(1)-u/a)/exp(a*(k-1))+u/a;%求出GM(1,1)模型公式
end
G=[];
G(1)=A(1);
for k=2:(n)G(k)=F(k)-F(k-1);%两者做差还原原序列，得到预测数据
end
%t1=1:(n-1);
t1=1:n;
t2=1:n;
plot(t1,A,'bo--');
hold on;
plot(t2,G,'r*-'); 
title('predict result');
legend('real value','predict value')

5.2 输出结果

F(k)=(A(1)-u/a)/exp(a*(k-1))+u/a

G(k)=F(k)-F(k-1)

5.3 结果检验

%后验差检验
e=A-G;
q=e/A;%相对误差
s1=var(A);
s2=var(e);
c=s2/s1;%方差比
len=length(e);
p=0;  %小误差概率
for i=1:lenif(abs(e(i))<0.6745*s1)p=p+1;end
end
p=p/len

输出结果：

c = 0.045

p = 1#故预测模型可用

u= 29.6891

a= -0.0447

5.4 预测泛化

（1） GM(1,1)模型

F(k)=(A(1)-u/a)/exp(a*(k-1))+u/a

F(12)=(30.925+29.6891/0.0447)*exp(0.0447*11)+29.6891/(-0.0447)= 472.3853

F(11)=(30.925+29.6891/0.0447)*exp(0.0447*10)+29.6891/(-0.0447)= 422.6993

（2）数据预测

G(k)=F(k)-F(k-1)

Y2021 =G(12)=F(12)-F(11)

=(30.925+29.6891/0.0447)/exp(-0.0447*11)-(30.925+29.6891/0.0477)/exp(-0.0477*10)

=(30.925+29.6891/0.0447)*(1/exp(-0.0447*11)-1/exp(-0.0447*10))

=49.6860

Y2021 = G(12) = F(12)-F(11)

= 472.3790 - 422.6942=49.6848

Y2022 = G(13) = F(13)-F(12)

= 524.3351-472.3790 = 51.9561

六、Python编程实现

6.1 程序代码

# -*- coding: utf-8 -*-
"""
Created on Sat Aug 21 20:51:43 2021
@author: zcq
"""
import matplotlib.pyplot as plt
import numpy as np
#载入math，提供对浮点数的数学运算函数
import math#年用电量
A_SG = [30.9250,32.5390,35.2270,34.6940,34.5060,36.0510,38.7450,42.3610,44.5360,45.783,47.52]
history_data = A_SG
#获取人口数据长度
n = len(history_data)
#将数据转换为数组
X0 = np.array(history_data)
#累加生成
history_data_agg = [sum(history_data[0:i+1]) for i in range(n)]
X1 = np.array(history_data_agg)
#计算数据矩阵B 和数据向量Y
B = np.zeros([n-1,2]) #初始化矩阵B
Y = np.zeros([n-1,1]) #初始化矩阵Y
for i in range(0,n-1):#用循环计算矩阵的每个元素B[i][0] = -0.5*(X1[i]+X1[i+1])B[i][1] = 1Y[i][0] = X0[i+1]
#计算GM(1,1)微分方程的参数a和u
A = np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)
#np.linalg.inv用于求逆矩阵,.T用于矩阵转置,使用最小二乘法估计微分方程的参数
'''
try:A = np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)#np.linalg.inv用于求逆矩阵,.T用于矩阵转置,使用最小二乘法估计微分方程的参数print(A)BB=np.linalg.inv(A)
except:print("矩阵不存在逆矩阵")
else:print(BB)
'''
a = A[0][0]
u = A[1][0]
#建立灰色预测模型
XX0 = np.zeros(n)#初始化预测结果矩阵XX0
XX0[0] = X0[0] #设置结果矩阵的第一项
for i in range(1,n): #循环预测结果矩阵的每一个值XX0[i] = (X0[0]-u/a)*(1-math.exp(a))*math.exp(-a*(i))
#模型精度的后验差检验
e = 0 #求残差平均值
for i in range(0,n):e += ( X0[i]- XX0[i])e /= n
print('残差均值为',np.round(e,2))
#求历史数据平均值
aver = 0;
for i in range( 0,n):aver += X0[i]aver /= n
print('历史数据平均值为',np.round(aver,2))
#求历史数据方差
s12 = 0;
for i in range( 0,n):s12 += (X0[i]-aver)**2s12 /= n
print('历史数据方差为',np.round(s12,2))
#求残差方差
s22 = 0;
for i in range(0,n):s22 += ( ( X0[i]-XX0[i]) - e)**2s22 /= n
print('残差方差为',np.round(s22,2))
#求后验差比值
C = s22 / s12
print('后验差比值为C=',np.round(C,2))
#求小误差概率
cout = 0
'''
服从标准正态分布，让误差概率落在25%到75%区域内。
由正态分布表可知，对应为25%和75%的值分别为－0．6745 和+0．6745。
'''
for i in range(0,n):    if abs((X0[i]-XX0[i])-e)<0.6754*math.sqrt(s12):cout = cout+1else:cout = cout
P = cout / n
print('小误差概率为P=',np.round(P,2))if ( C < 0.35 and P>0.95): #是否满足一级预测精度m = 10 #需要预测的年数print('往后',m,'年，各年负荷为:')f = np.zeros(m)for i in range(0,m):#预测结果的值f[i] = (X0[0]-u/a)*(1-math.exp(a))*math.exp(-a*(i+n))print(np.round(f[i],2)) #取两位小数
else:print('灰色预测法不适用')
#y = history_data
y = np.append(history_data,f)
x_num = m+n
plt.plot(y,"b.")
plt.grid()