import pandas as pd
import numpy as np

od=np.array([17,7,4,28,38.6,7,38,6,51,91.9,4,5,17,26,36,28,50,27,105,0,39.3,90.3,36.9,0,156.5])
od=od.reshape(5,5)
od=pd.DataFrame(od)
od

初始OD矩阵

image.png

time_tom=[4,9,11,9,8,12,11,12,4]
time_tom=np.array(time_tom)
time_tom=time_tom.reshape(3,3)
time_tom=pd.DataFrame(time_tom)
time_tom

将来时间

image.png

time_now=[7,17,22,17,15,23,22,23,7]
time_now=np.array(time_now).reshape([3,3])
time_now=pd.DataFrame(time_now)
time_now

现在时间

image.png

temp=pd.DataFrame(np.zeros([len(time_now)**2,9]))
temp.columns=['样本点','Qij','Oi','Dj','Oi+Dj','Cij','y','X1','X2']
temp['样本点']=['i='+str(i)+',j= '+str(j) for i in range(1,4) for j in range(1,4)]
temp['Qij']=od.iloc[0:3,0:3].values.flatten()
temp['Oi']=np.array([x for x in od[3].values[0:len(od)-2]]).repeat(3)
temp['Dj']=[x for x in od.iloc[3,:].values[0:len(od)-2]]*3
temp['Oi+Dj']=temp['Oi']*temp['Dj']
temp['Cij']=time_now.values.flatten()
temp['y']=np.log(temp.Qij)
temp['X1']=np.log(temp.Oi*temp.Dj)
temp['X2']=np.log(temp.Cij)
temp

重力法矩阵

image.png

# #最小二乘法
# from sklearn.linear_model import LinearRegression 
# clf=LinearRegression()
# linear=clf.fit(temp.iloc[:,-2:].values,temp.y.values)
# A1,A2=linear.coef_
# A0=linear.intercept_
# print("三参数A0,A1,A2为: ",A0,A1,A2)
# a=np.exp(A0)
# b=A1
# r=-A2
# print("三参数a,b,r为: ",a,b,r)

三参数A0,A1,A2为:  -2.08379785868 1.17268924579 -1.45531274106
三参数a,b,r为:  0.124456644748 1.17268924579 1.45531274106

##转使用题目中给定的参数
a=0.183
b=1.152
r=1.536

迭代用函数

#计算第一次得到的OD表
def getod(df):
    #OD=df.values.reshape(-3,3)
    OD=pd.DataFrame(df.reshape(-1,3))
    OD['合计']=pd.DataFrame(OD).apply(sum,axis=1)
    OD.loc['合计'] = OD.apply(lambda x: x.sum())
    return OD
#计算Foi和Fdj
#福莱特法
#求发生增长系数和吸引增长系数
def F(i):
    Foi.append(od.iloc[:-2,-1].values/OD[i].iloc[:-1,-1].values)
    Fdj.append(od.iloc[-1,:-2].values/OD[i].loc['合计'][:-1].values)
    Loi=od.iloc[:-2,3].values/np.sum(Fdj[i]*od.iloc[:3,:-2].values,axis=1)
    Ldj=od.iloc[3,:-2].values/np.sum(Foi[i]*od.iloc[:3,:-2].T.values,axis=1)
    return Loi,Ldj

def jianyan(i,c=0.03):
    count=np.sum([abs(Foi[i]-1)>c,abs(Fdj[i]-1)>c])
    if count!=0:
        print(count,'个不满足模型分布的约束条件，将使用福莱特法进行修正。')
        return count
    if count==0:
        print("迭代结束")
        return count
        
#求qij 求Q转换为od
def find_Q(i):
    q=OD[i].iloc[:-1,:-1].values.flatten()
    foi=Foi[i].repeat(3)
    fdj=Fdj[i].tolist()*3
    li=Loi.repeat(3)
    lj=Ldj.tolist()*3
    Q=q*foi*fdj*(li+lj)/2
    return Q

迭代过程

OD={}
Foi=[]
Fdj=[]
i=0
OD[0]=getod(df.result)
while 1:
    Loi,Ldj=F(i)
    print("Loi,Ldj分别为：",Loi,Ldj)
    if jianyan(i,c=0.01)==0:
        print("第 " ,i+1 ," 次收敛")
        print("最后一次For,Fdj分别为为：",np.hstack([Foi[i],Fdj[i]]).tolist())
        break
    p=find_Q(i)
    i+=1
    OD[i]=getod(p)
    

结果

image.png