[关闭]
@Lin-- 2019-09-04T16:10:54.000000Z 字数 5143 阅读 515

数学建模算法

数学建模

灰色系统预测

  1. #灰色预测
  3. import numpy as np
  4. import math
  5. import matplotlib.pyplot as plt
  6. #plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
  8. history_data = [724.57,746.62,778.27,800.8,827.75,871.1,912.37,954.28,995.01,1037.2]
  9. n = len(history_data)
  10. X0 = np.array(history_data)
  12. #可容覆盖
  13. Lamda=[]
  14. for i in range(n-1):
  15.     Lamda.append(history_data[i]/history_data[i+1])
  16. if min(Lamda)>math.exp(-2/(n+1)) and max(Lamda)<math.exp(2/(n+2)):
  17.     print("在可容覆盖范围内，无需处理数据")
  18. else:
  19.     pass
  20.     #print("不在可容覆盖范围内，需对数据做处理")
  22. #累加生成
  23. history_data_agg = [sum(history_data[0:i+1]) for i in range(n)]
  24. X1 = np.array(history_data_agg)
  26. #计算数据矩阵B和数据向量Y
  27. B = np.ones([n-1,2])
  28. Y = np.zeros([n-1,1])
  29. for i in range(0,n-1):
  30.     B[i][0] = -0.5*(X1[i] + X1[i+1])
  31.     Y[i][0] = X0[i+1]
  33. #计算GM(1,1)微分方程的参数a和u
  34. #A = np.zeros([2,1])
  35. A = np.linalg.inv(B.T.dot(B)).dot(B.T).dot(Y)
  36. a = A[0][0]
  37. u = A[1][0]
  39. #建立灰色预测模型
  40. XX0 = np.zeros(n)
  41. XX0[0] = X0[0]
  42. for i in range(1,n):
  43.     XX0[i] = (X0[0] - u/a)*(1-math.exp(a))*math.exp(-a*(i));
  45. #模型精度的后验差检验
  46. e = 0      #求残差平均值
  47. for i in range(0,n):
  48.     e += (X0[i] - XX0[i])
  49. e /= n
  51. '''
  52. #图像(曲线)
  53. plt.plot(np.arange(1,n+1),X0,color="#F08080")
  54. plt.plot(np.arange(1,n+1),XX0,color='b')
  55. plt.show()
  56. '''
  57. #点图
  58. plt.scatter(np.arange(1,n+1),X0,marker='*',color='r',label='历史数据')
  59. plt.scatter(np.arange(1,n+1),XX0,marker='+',color='b',label='预测数据')
  60. plt.legend(loc='upper right')
  61. plt.show()
  63. #求历史数据平均值
  64. aver = 0;
  65. for i in range(0,n):
  66.     aver += X0[i]
  67. aver /= n
  69. #求历史数据方差
  70. s12 = 0;
  71. for i in range(0,n):
  72.     s12 += (X0[i]-aver)**2;
  73. s12 /= n
  75. #求残差方差
  76. s22 = 0;
  77. for i in range(0,n):
  78.     s22 += ((X0[i] - XX0[i]) - e)**2;
  79. s22 /= n
  81. #求后验差比值
  82. C = s22 / s12
  84. #求小误差概率
  85. cout = 0
  86. for i in range(0,n):
  87.     if abs((X0[i] - XX0[i]) - e) < 0.6754*math.sqrt(s12):
  88.         cout = cout+1
  89.     else:
  90.         cout = cout
  91. P = cout / n
  93. if (C < 0.35 and P > 0.95):
  94.     #预测精度为一级
  95.     m = 10   #请输入需要预测的年数
  96.     #print('往后m各年负荷为：')
  97.     f = np.zeros(m)
  98.     for i in range(0,m):
  99.         f[i] = (X0[0] - u/a)*(1-math.exp(a))*math.exp(-a*(i+n))
  100.     print(f)
  101. else:
  102.     print('灰色预测法不适用')

层次分析法

  1. #层次分析法
  3. import numpy as np
  4. import math
  5. import matplotlib.pyplot as plt
  7. #计算矩阵各行乘积
  8. def array_mul(a):
  9.     a1=[]
  10.     for i in range(a.shape[0]):
  11.         temp=1.0
  12.         for j in range(a.shape[1]):
  13.             temp*=a[i][j]
  14.         a1.append(temp)
  15.     return np.array(a1)
  17. #矩阵归一化
  18. def normalization(a):
  19.     a1=[]
  20.     sum=a.sum()
  21.     for i in a:
  22.         a1.append(i/sum)
  23.     return np.array(a1)
  25. #初始变量
  26. P1=[[1,0.3333,0.125]
  27.     ,[3,1,0.3333]
  28.     ,[8,3,1]]
  29. #计算每一行的乘积
  30. P1=np.array(P1)
  31. PP1=array_mul(P1)
  33. #三次方根
  34. PP1=np.power(PP1,0.3333)
  36. #归一化矩阵
  37. PP1=normalization(PP1)
  38. print(PP1)
  40. #求最大特征根
  41. AW=P1.dot(PP1.T)
  42. print(AW)
  43. Lamda=0
  44. for i in range(AW.size):
  45.     Lamda+=AW[i]/(AW.size*PP1[i])
  46. print(Lamda)
  48. #一致性检验
  49. CI=(Lamda-AW.size)/(AW.size-1)
  50. RI=[0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45,1.49,1.51]
  51. CR=CI/RI[AW.size-1]
  52. print(CR)

主成分分析法

  1. #主成分分析法
  3. import numpy as np
  4. from scipy import stats
  5. from numpy import linalg as LA
  7. #矩阵标准化函数
  8. def zscore(a):
  9.     a1=np.subtract(a,np.mean(a,axis=0))
  10.     return np.divide(a1,np.std(a,axis=0,ddof=1))
  12. #计算累计贡献率达到85%对应第几个特征值
  13. def contr(a):
  14.     sum=np.sum(a)
  15.     count,temp=1,0
  16.     for i in a:
  17.         temp+=i
  18.         if (temp/sum)>0.85:
  19.             break
  20.         else: count+=1
  21.     return count
  23. #初值
  24. P=[[10.2352,11.3220,12.3315,12.5644,12.6534,13.5541,13.6698,13.9964],
  25.    [10.1223,11.8110,11.9156,12.6654,9.8652,15.6634,10.0023,16.0034],
  26.    [15.2364,9.2563,10.1156,10.1156,11.9874,16.2345,14.2648,12.1238],
  27.    [7.5612,8.1234,9.6321,8.5678,9.6654,10.2322,11.2344,15.6652]]
  28. P=np.array(P)
  30. #标准化矩阵
  31. P1=zscore(P)
  33. #建立相关系数矩阵
  34. P2=np.corrcoef(P1,rowvar=False)
  36. #矩阵特征值以及主成分个数
  37. Lamda=np.linalg.eigvalsh(P2)
  38. Lamda=-np.sort(-Lamda)#对矩阵特征值从大到小排序
  39. m=contr(Lamda)#主成分个数
  41. #矩阵的特征向量，为特征值对应的那一列
  42. Lamda_vector=np.linalg.eig(P2)

基本遗传算法

  1. #遗传算法
  2. #求y = 10 * sin(5x) + 7 * cos(4x)的最大值
  3. import numpy as np
  4. import matplotlib.pyplot as plt
  5. import math
  6. import random
  7. pop_size=500#种群规模
  8. max_value=10#基因的搜索空间的上界
  9. chrom_length=10#染色体长度，即二进制位数
  10. pc=0.6#交配概率
  11. pm=0.01#基因突变概率
  12. results=[[]]#每一代的最优解，N个二元组
  13. fit_value=[]#个体适应度
  14. fit_mean=[]#平均适应度
  16. #产生随机生成序列
  17. def GeneEncoding(pop_size, chrom_length):
  18.     pop=[[]]
  19.     for i in range(pop_size):
  20.         temp=[]
  21.         for j in range(chrom_length):
  22.             temp.append(random.randint(0,1))
  23.         pop.append(temp)
  24.     return pop[1:]
  26. pop=GeneEncoding(pop_size,chrom_length)#生成初始数据
  28. #解码计算每一个个体的表现型
  29. def DecodeChrom(pop,chrom_length):#解码过程
  30.     temp=[]
  31.     for i in range(len(pop)):
  32.         t=0
  33.         for j in range(chrom_length):
  34.             t+=pop[i][j]*(math.pow(2,j))
  35.         temp.append(t)
  36.     return temp
  38. def CalObjValue(pop,chrom_length,max_value):#计算表现型
  39.     temp1=[]
  40.     obj_value=[]
  41.     temp1=DecodeChrom(pop,chrom_length)
  42.     for i in temp1:
  43.         x=i*max_value/(math.pow(2,chrom_length)-1)
  44.         obj_value.append(10*math.sin(5*x)+7*math.cos(4*x))
  45.     return obj_value
  47. #淘汰负值
  48. def CalFitValue(obj_value):
  49.     fit_value=[]
  50.     for i in obj_value:
  51.         if i>0:
  52.             fit_value.append(i)
  53.         else:
  54.             fit_value.append(0.0)
  55.     return fit_value
  57. #选择
  58. #适应度总和
  59. def SumFitValue(fit_value):
  60.     total=0
  61.     for i in fit_value:
  62.         total+=i
  63.     return total
  64. #累计概率
  65. def cumsum(newfit_value):
  66.     for i in range(len(fit_value-2,-1,-1)):
  67.         t,j=0,0
  68.         while j<i:
  69.             t+=newfit_value[j]
  70.             j+=1
  71.         newfit_value[i]=t
  72.         newfit_value[-1]=1
  73.     return newfit_value
  74. #选择函数
  75. def selection(pop,fit_value):
  76.     newfit_value=[]
  77.     total_fit=SumFitValue(fit_value)#总适应度
  78.     for i in fit_value:
  79.         newfit_value.append(i/total_fit)#归一化
  80.     newfit_value=cumsum(newfit_value)#累计概率
  81.     #转轮盘算法选择
  82.     ms=[]
  83.     for i in range(len(pop)):
  84.         ms.append(random.random())
  85.     ms.sort()
  86.     fitin,newin=0,0
  87.     newpop=pop
  88.     while newin<len(pop):
  89.         if ms[newin]<newfit_value[fitin]:
  90.             newpop=pop[fitin]
  91.             newin+=1
  92.         else:
  93.             fitin+=1
  94.         pop=newpop
  96. #交配/交换染色体
  97. def crossover(pop,pc):
  98.     for i in range(len(pop)-1):
  99.         if random.random()<pc:
  100.             cpoint=random.randint(0,len(pop[0]))
  101.             temp1,temp2=[],[]
  102.             temp1.extend(pop[i][0:cpoint])
  103.             temp1.extend(pop[i+1][cpoint:-1])
  104.             temp2.extend(pop[i+1][0:cpoint])
  105.             temp2.extend(pop[i][cpoint:-1])
  106.             pop[i]=temp1
  107.             pop[i+1]=temp2
  109. #变异/基因突变
  110. def mutation(pop,pm):
  111.     px=len(pop)
  112.     py=len(pop[0])
  113.     for i in range(px):
  114.         if random.random()<pm:
  115.             mpoint=random.randint(0,py-1)
  116.             if pop[i][mpoint]==0:
  117.                 pop[i][mpoint]=1
  118.             else:
  119.                 pop[i][mpoint]=0
添加新批注
在作者公开此批注前，只有你和作者可见。
回复批注