基于Python实现粒子群PSO算法优化Prophet模型参数（源码+数据）.rar

# -*- coding: utf-8 -*- """ 目的：PSO优化 changepoint_range=0.9,changepoint_prior_scale=0.8。 流程：非等间隔优化预测， 根据pso-prophet_original_opt优化得出的参数，直接进行优化求解。并与未经优化的算法进行比较 # array([ 0.56699489, 0.87695148]) 6.4447443347379894 """ import numpy as np import matplotlib.pyplot as plt import pandas as pd from pylab import mpl mpl.rcParams['font.sans-serif'] = ['STZhongsong'] # 指定默认字体：解决plot不能显示中文问题 mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题 ## 读入原始数据并绘图 IO='D:\\spyder workspath\\LSTM\\04.边坡变形监测统计-累计变形.xlsx' # data=pd.read_excel(IO,sheet_name="P2",usecols=[0,1,2,3,4,5,6], nrows=np.arange(104), skiprows=1,header=0) data=pd.read_excel(IO,sheet_name="P2",usecols=[0,1,2,3,4,5,6], nrows=104, skiprows=1,header=0) data_use=data.iloc[:,[0,3]] # 读取有用的列 df=data_use.iloc[0:102,:] # 读取前104行，nrows有时不可用 df=df.dropna() # (先提取需要的，再去除空行) df.columns=['ds','y'] # 绘制三维累计变形曲线 fig=plt.figure('BX30变形时序',figsize=(4,3)) plt.plot(data.iloc[:,0],-data.iloc[:,2],'r-o',markersize=1.5,linewidth=1,label='E') plt.plot(data.iloc[:,0],data.iloc[:,3],'k-o',markersize=1.5,linewidth=1,label='N') plt.plot(data.iloc[:,0],-data.iloc[:,4],'b-o',markersize=1.5,linewidth=1,label='H') plt.xticks(rotation=15) plt.xlabel('时间',fontsize=10) plt.ylabel('累计变形量/mm',fontsize=10) plt.legend(loc='upper left')#必须加，否则不显示图例 ## 定义训练集、测试集、预测集 def mean_squared_error(t1,t2): dt=t1-t2 num=len(dt)-1 rmse=np.sqrt(sum(dt**2)/num) return rmse t1=72 t2=92 #[：t1] 训练集，[t1,t2] 测试集，[t2：] 预测集， dx=df[:t1] #训练集 dc=df[t1:t2] #测试集 s_dc=dc['y'].values #测试集上预测的时间 dp=df[t2:] #预测集上实测参数值，用以选取 future_dc = pd.DataFrame(dc['ds']) #训练集训练后，在测试集上预测的时间 fig=plt.figure('时间时序图',figsize=(4,3)) plt.plot(df['ds'],df['y'],'k-o',markersize=1.5) plt.axvline(df['ds'][t1],color='r',linestyle='dashdot') plt.axvline(df['ds'][t2],color='r',linestyle='dashdot') plt.text(df['ds'][50],35, "训练集", color = "k",size = 10, alpha = 0.9) plt.annotate(s='',xy=(df['ds'][t1],30),xytext=(-45,0),textcoords='offset points',fontsize=10, arrowprops=dict(arrowstyle='<-')) plt.text(df['ds'][77],35, "测试集", color = "k",size = 10, alpha = 0.9) plt.annotate(s='',xy=(df['ds'][t1],30),xytext=(60,0),textcoords='offset points',fontsize=10, arrowprops=dict(arrowstyle='<->')) plt.text(df['ds'][93],35, "预测集", color = "k",size = 10, alpha = 0.9) plt.annotate(s='',xy=(df['ds'][t2],30),xytext=(40,0),textcoords='offset points',fontsize=10, arrowprops=dict(arrowstyle='<-')) plt.xlabel('时间',fontsize=10) plt.ylabel('累计形变量/mm',fontsize=10) plt.xticks(rotation=20) #绘制图像完毕 # 绘制最佳粒子选择过程 IO2 = 'D:\\spyder workspath\\moder_aricmetic\\for_prophet\\modification-prophet\\xunlian1.xlsx' goal=pd.read_excel(IO2) fig=plt.figure('全局粒子搜索过程',figsize=(4,3)) #plt.axhline(y=min(goal.iloc[:,2]),color='r',linestyle='dashdot') plt.plot(sorted(goal.iloc[:,2],reverse=True)) plt.annotate(s='群体最优解 6.45',xy=(300,min(goal.iloc[:,2])),xytext=(-80,40),textcoords='offset points',fontsize=10, arrowprops=dict(arrowstyle='->',linestyle="dashed")) plt.scatter(300,min(goal.iloc[:,2]),c='g',marker='o') plt.xlabel('粒子数',fontsize=10) plt.ylabel('适应度',fontsize=10) # from fbprophet import Prophet from fbprophet.plot import add_changepoints_to_plot # 优化后的参数 dxc=df[:t2] #确定最佳参数后，将训预练集合和预测集一起作为训练数据 forecast_dp=pd.DataFrame(dp['ds']) #在测试集上 预测的时间 s_dp=dp['y'].values #在预测集上实测值，用于评价预测结果。 m = Prophet(changepoint_range=0.56699489,changepoint_prior_scale=0.87695148,weekly_seasonality=False,n_changepoints=6) m.fit(dxc) prediction = m.predict(pd.DataFrame(df['ds'])) #获取拟合值与预测值 def_pre=prediction['yhat'][t2:].values #预测集上的预测值 err_fina_pre1=mean_squared_error(def_pre,s_dp) ## 展示拟合与预测结果 #绘制提取的趋势项与周期项 # 突变点、原始点、与趋势项 fig=plt.figure('变形分解',figsize=(4,3)) plt.plot(dxc['ds'],dxc['y'],'ko',linewidth=1,markersize=1.5,label='原始时序') plt.plot(prediction['ds'][:t2],prediction['trend'][:t2],linewidth=1.5,label='趋势项') plt.plot(prediction['ds'][:t2],prediction['additive_terms'][:t2],linewidth=1.5,label='周期项') plt.xlabel('时间',fontsize=10) plt.ylabel('累计形变量/mm',fontsize=10) plt.xticks(rotation=20) plt.legend(loc='best',fontsize=8) # 突变点 v_df=pd.DataFrame(index=df['ds'],columns=['v']) for i in np.arange(1,len(df)): v_df['v'][i]=abs((df['y'][i]-df['y'][i-1])/(df['ds'][i]-df['ds'][i-1]).days) v_df=v_df.dropna() # (先提取需要的，再去除空行) #from fbprophet.plot import add_changepoints_to_plot #fig1 = m.plot(prediction) #plt.xlabel('time',fontsize=15) #plt.xticks(rotation=15) #plt.ylabel('形变量',fontsize=15) #for i in m.changepoints.values: # plt.axvline(i,color='r',linestyle='dashdot') #plt.axvline(m.changepoints.values[0],color='r',linestyle='dashdot',label='突变点') ## a = add_changepoints_to_plot(fig5.gca(), m, forecast,label='突变点') #plt.plot(v_df.index,v_df['v'],'k-',linewidth=2,label='速率') #plt.legend(loc='best')#必须加，否则不显示图例 fig=plt.figure("突变点示意",figsize=(4,3)) plt.plot(dxc['ds'],dxc['y'],'ko',markersize=1.5,linewidth=1,label='原始时序') plt.plot(prediction['ds'][:t2],prediction['trend'][:t2],linewidth=1,label='趋势项') #plt.plot(v_df.index,v_df['v']*5,'b-o',linewidth=1,markersize=1.5,label='变形速率') #plt.bar(v_df.index[:t2],v_df['v'][:t2]*50,width=1,color="blue",label='变形速率') for i in m.changepoints.values: plt.axvline(i,color='r',linestyle='dashdot') plt.axvline(m.changepoints.values[0],color='r',linestyle='dashdot',label='突变点时刻') plt.xlabel('时间',fontsize=10) plt.ylabel('累计形变量/mm',fontsize=10) plt.xticks(rotation=20) plt.legend(loc='best',fontsize=8) fig=plt.figure("突变点与变形速率示意",figsize=(4,3)) plt.plot(v_df.index,v_df['v'],'b-o',linewidth=1,markersize=1.5,label='变形速率') for i in m.changepoints.values: plt.scatter(i,v_df['v'][i],c='r',marker='o') plt.scatter(m.changepoints.values[-1],v_df['v'][m.changepoints.values].values[-1],c='r',marker='o',label='突变点时刻') plt.xlabel('时间',fontsize=10) plt.xticks(rotation=20) plt.ylabel('变形速率/mm',fontsize=10) plt.xticks(rotation=20) plt.legend(loc='best',fontsize=8) # 周期项示意图 fig=plt.figure("周期项示意图",figsize=(4,3)) plt.plot(prediction['ds'][:t2],prediction['additive_terms'][:t2],linewidth=1.5,label='周期项') plt.xlabel('时间',fontsize=10) plt.xticks(rotation=20) plt.ylabel('变形量/mm',fontsize=10) # 预测值 fig=plt.figure('组合预测',figsize=(4,3)) plt.plot(df['ds'],df['y'],'k-o',linewidth=1,markersize=1.5,label='实测值') plt.plot(prediction['ds'][t2:],prediction['yhat'][t2:],'r-^',linewidth=1,markersize=1.5,label='组合预测值') plt.plot(prediction['ds'][:t2],prediction['trend'][:t2],'k-o',linewidth=1,markersize=1.5,label='趋势项') plt.plot(prediction['ds'][t2:],prediction['trend'][t2:],'r-o',linewidth=1,markersize=1.5,label='趋势项预测值') plt.plot(prediction['ds'][:t2],prediction['additive_terms'][:t2],'k-^',linewidth=1,markersize=1.5,label='周期项') plt.plot(prediction['ds']