@Lin--
2019-09-04T16:10:54.000000Z
字数 5143
阅读 392
数学建模
#灰色预测
import numpy as np
import math
import matplotlib.pyplot as plt
#plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
history_data = [724.57,746.62,778.27,800.8,827.75,871.1,912.37,954.28,995.01,1037.2]
n = len(history_data)
X0 = np.array(history_data)
#可容覆盖
Lamda=[]
for i in range(n-1):
Lamda.append(history_data[i]/history_data[i+1])
if min(Lamda)>math.exp(-2/(n+1)) and max(Lamda)<math.exp(2/(n+2)):
print("在可容覆盖范围内,无需处理数据")
else:
pass
#print("不在可容覆盖范围内,需对数据做处理")
#累加生成
history_data_agg = [sum(history_data[0:i+1]) for i in range(n)]
X1 = np.array(history_data_agg)
#计算数据矩阵B和数据向量Y
B = np.ones([n-1,2])
Y = np.zeros([n-1,1])
for i in range(0,n-1):
B[i][0] = -0.5*(X1[i] + X1[i+1])
Y[i][0] = X0[i+1]
#计算GM(1,1)微分方程的参数a和u
#A = np.zeros([2,1])
A = np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)
a = A[0][0]
u = A[1][0]
#建立灰色预测模型
XX0 = np.zeros(n)
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
'''
#图像(曲线)
plt.plot(np.arange(1,n+1),X0,color="#F08080")
plt.plot(np.arange(1,n+1),XX0,color='b')
plt.show()
'''
#点图
plt.scatter(np.arange(1,n+1),X0,marker='*',color='r',label='历史数据')
plt.scatter(np.arange(1,n+1),XX0,marker='+',color='b',label='预测数据')
plt.legend(loc='upper right')
plt.show()
#求历史数据平均值
aver = 0;
for i in range(0,n):
aver += X0[i]
aver /= n
#求历史数据方差
s12 = 0;
for i in range(0,n):
s12 += (X0[i]-aver)**2;
s12 /= n
#求残差方差
s22 = 0;
for i in range(0,n):
s22 += ((X0[i] - XX0[i]) - e)**2;
s22 /= n
#求后验差比值
C = s22 / s12
#求小误差概率
cout = 0
for i in range(0,n):
if abs((X0[i] - XX0[i]) - e) < 0.6754*math.sqrt(s12):
cout = cout+1
else:
cout = cout
P = cout / n
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(f)
else:
print('灰色预测法不适用')
#层次分析法
import numpy as np
import math
import matplotlib.pyplot as plt
#计算矩阵各行乘积
def array_mul(a):
a1=[]
for i in range(a.shape[0]):
temp=1.0
for j in range(a.shape[1]):
temp*=a[i][j]
a1.append(temp)
return np.array(a1)
#矩阵归一化
def normalization(a):
a1=[]
sum=a.sum()
for i in a:
a1.append(i/sum)
return np.array(a1)
#初始变量
P1=[[1,0.3333,0.125]
,[3,1,0.3333]
,[8,3,1]]
#计算每一行的乘积
P1=np.array(P1)
PP1=array_mul(P1)
#三次方根
PP1=np.power(PP1,0.3333)
#归一化矩阵
PP1=normalization(PP1)
print(PP1)
#求最大特征根
AW=P1.dot(PP1.T)
print(AW)
Lamda=0
for i in range(AW.size):
Lamda+=AW[i]/(AW.size*PP1[i])
print(Lamda)
#一致性检验
CI=(Lamda-AW.size)/(AW.size-1)
RI=[0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45,1.49,1.51]
CR=CI/RI[AW.size-1]
print(CR)
#主成分分析法
import numpy as np
from scipy import stats
from numpy import linalg as LA
#矩阵标准化函数
def zscore(a):
a1=np.subtract(a,np.mean(a,axis=0))
return np.divide(a1,np.std(a,axis=0,ddof=1))
#计算累计贡献率达到85%对应第几个特征值
def contr(a):
sum=np.sum(a)
count,temp=1,0
for i in a:
temp+=i
if (temp/sum)>0.85:
break
else: count+=1
return count
#初值
P=[[10.2352,11.3220,12.3315,12.5644,12.6534,13.5541,13.6698,13.9964],
[10.1223,11.8110,11.9156,12.6654,9.8652,15.6634,10.0023,16.0034],
[15.2364,9.2563,10.1156,10.1156,11.9874,16.2345,14.2648,12.1238],
[7.5612,8.1234,9.6321,8.5678,9.6654,10.2322,11.2344,15.6652]]
P=np.array(P)
#标准化矩阵
P1=zscore(P)
#建立相关系数矩阵
P2=np.corrcoef(P1,rowvar=False)
#矩阵特征值以及主成分个数
Lamda=np.linalg.eigvalsh(P2)
Lamda=-np.sort(-Lamda)#对矩阵特征值从大到小排序
m=contr(Lamda)#主成分个数
#矩阵的特征向量,为特征值对应的那一列
Lamda_vector=np.linalg.eig(P2)
#遗传算法
#求y = 10 * sin(5x) + 7 * cos(4x)的最大值
import numpy as np
import matplotlib.pyplot as plt
import math
import random
pop_size=500#种群规模
max_value=10#基因的搜索空间的上界
chrom_length=10#染色体长度,即二进制位数
pc=0.6#交配概率
pm=0.01#基因突变概率
results=[[]]#每一代的最优解,N个二元组
fit_value=[]#个体适应度
fit_mean=[]#平均适应度
#产生随机生成序列
def GeneEncoding(pop_size, chrom_length):
pop=[[]]
for i in range(pop_size):
temp=[]
for j in range(chrom_length):
temp.append(random.randint(0,1))
pop.append(temp)
return pop[1:]
pop=GeneEncoding(pop_size,chrom_length)#生成初始数据
#解码计算每一个个体的表现型
def DecodeChrom(pop,chrom_length):#解码过程
temp=[]
for i in range(len(pop)):
t=0
for j in range(chrom_length):
t+=pop[i][j]*(math.pow(2,j))
temp.append(t)
return temp
def CalObjValue(pop,chrom_length,max_value):#计算表现型
temp1=[]
obj_value=[]
temp1=DecodeChrom(pop,chrom_length)
for i in temp1:
x=i*max_value/(math.pow(2,chrom_length)-1)
obj_value.append(10*math.sin(5*x)+7*math.cos(4*x))
return obj_value
#淘汰负值
def CalFitValue(obj_value):
fit_value=[]
for i in obj_value:
if i>0:
fit_value.append(i)
else:
fit_value.append(0.0)
return fit_value
#选择
#适应度总和
def SumFitValue(fit_value):
total=0
for i in fit_value:
total+=i
return total
#累计概率
def cumsum(newfit_value):
for i in range(len(fit_value-2,-1,-1)):
t,j=0,0
while j<i:
t+=newfit_value[j]
j+=1
newfit_value[i]=t
newfit_value[-1]=1
return newfit_value
#选择函数
def selection(pop,fit_value):
newfit_value=[]
total_fit=SumFitValue(fit_value)#总适应度
for i in fit_value:
newfit_value.append(i/total_fit)#归一化
newfit_value=cumsum(newfit_value)#累计概率
#转轮盘算法选择
ms=[]
for i in range(len(pop)):
ms.append(random.random())
ms.sort()
fitin,newin=0,0
newpop=pop
while newin<len(pop):
if ms[newin]<newfit_value[fitin]:
newpop=pop[fitin]
newin+=1
else:
fitin+=1
pop=newpop
#交配/交换染色体
def crossover(pop,pc):
for i in range(len(pop)-1):
if random.random()<pc:
cpoint=random.randint(0,len(pop[0]))
temp1,temp2=[],[]
temp1.extend(pop[i][0:cpoint])
temp1.extend(pop[i+1][cpoint:-1])
temp2.extend(pop[i+1][0:cpoint])
temp2.extend(pop[i][cpoint:-1])
pop[i]=temp1
pop[i+1]=temp2
#变异/基因突变
def mutation(pop,pm):
px=len(pop)
py=len(pop[0])
for i in range(px):
if random.random()<pm:
mpoint=random.randint(0,py-1)
if pop[i][mpoint]==0:
pop[i][mpoint]=1
else:
pop[i][mpoint]=0