%雙約束重力模型標定及計算用
clear
clc
close?all
%輸入標定數據
qij=[1774;7386;4517];
Oi=[285126];
P=[38.691.936.0];
Dj=[285027];
A=[39.390.336.9];
dij=[4911;9812;11124];
Cij=[71722;171523;22237];
Z=-0.5;
n=size（qij1）;
tic
while?2018
%gravity為自定義函數，用于雙約束重力模型kikj的標定
[kimkjm]=gravity（qijOiDjCijZ）;
%冪指數Z的標定
s1=0;
s2=0;
for?i=1:n
????????for?j=1:n
????????????Qij（ij）=kim（i）*kjm（j）*Oi（i）*Dj（j）*（Cij（ij）.^Z）;
????????????s1=s1+?Qij（ij）*Cij（ij）;
?????????????s2=s2+?qij（ij）*Cij（ij）;
????????end
end
R1=s1/sum（sum（Qij））;
R2=s2/sum（sum（qij））;
R=R1/R2;
if?abs（R-1）>0.01
????if?R>1
????????Z=1.5*Z;
????else
????????Z=0.5*Z;
????end
else
????break;
end
end
%利用底特律法進行OD均衡調整
for?i=1:n
????????for?j=1:n
????????????Hij（ij）=kim（i）*kjm（j）*P（i）*A（j）*（dij（ij）.^Z）;%生成初始分布
????????end
end
cnt=2;
while?2018
????if?mod（cnt2）==0
rate=sum（Hij2）./P‘;
Hij=Hij./repmat（rate1n）;
????else
?rate=sum（Hij1）./A;
Hij=Hij./repmat（raten1）;
????end???????
????t=0;
for?i=1:n
????if?abs（rate（i）-1）>0.01
????????t=1;
????????break;
????end
end
if?t==0
????break;
end
cnt=cnt+1;
end
toc
????fprintf（‘預測未來的分布為>>‘）;
????Hij

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件????????1255??2020-05-13?19:54??雙約束重力模型（1.0）\doublecon.m
?????文件?????????593??2020-05-13?19:59??雙約束重力模型（1.0）\gravity.m
?????文件???????20557??2020-05-13?20:03??雙約束重力模型（1.0）\miao.mexw64
?????文件??????????72??2020-05-13?19:57??雙約束重力模型（1.0）\使用說明.txt