前言：终于写一次本专业的博客了，前段时间有个同学问我作业，突然翻到以前写的重力模型的matlab demo,来此记录一下

重力模型主要有无约束，单约束，双约束三种，

doublecon.m文件：

%双约束重力模型标定及计算用
clear
clc
close all
%输入标定数据
qij=[17,7,4;7,38,6;4,5,17];
Oi=[28,51,26];
P=[38.6,91.9,36.0];
Dj=[28,50,27];
A=[39.3,90.3,36.9];
dij=[4,9,11;9,8,12;11,12,4];
Cij=[7,17,22;17,15,23;22,23,7];
Z=-0.5;
n=size(qij,1);
tic
while 2018
%gravity为自定义函数，用于双约束重力模型ki,kj的标定
[kim,kjm]=gravity(qij,Oi,Dj,Cij,Z);
%幂指数Z的标定
s1=0;
s2=0;
for i=1:n
        for j=1:n
            Qij(i,j)=kim(i)*kjm(j)*Oi(i)*Dj(j)*(Cij(i,j).^Z);
            s1=s1+ Qij(i,j)*Cij(i,j);
             s2=s2+ qij(i,j)*Cij(i,j);
        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(i,j)=kim(i)*kjm(j)*P(i)*A(j)*(dij(i,j).^Z);%生成初始分布
        end
end
cnt=2;
while 2018
    if mod(cnt,2)==0
rate=sum(Hij,2)./P';
Hij=Hij./repmat(rate,1,n);
    else
 rate=sum(Hij,1)./A;
Hij=Hij./repmat(rate,n,1);
    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

gravity.m文件：

function [ kim,kjm ] = gravity( qij,Oi,Dj,Cij,Z )
%用于双约束重力模型ki,kj的标定
%qij为当前已知的分布矩阵，Oi为qij按列求和而来，Dj为qij按行求和而来，Cij为阻抗矩阵，Z为阻抗函数的幂指数
n=size(qij,1);
kim=ones(1,n);
while 2018
    t=0;
    kjm=miao(qij,kim,Oi,Cij,Z);%miao为c语言编译而成为加速计算的函数
    kim1=miao(qij,kjm,Dj,Cij,Z);
   kjm1=miao(qij,kim1,Oi,Cij,Z);
%判断是否满足收敛条件
    for i=1:n
      r1(i)=kim1(i)/kim(i);
      r2(i)=kjm1(i)/kjm(i);
    end
  for i=1:n
      if abs(r1(i)-1)>0.01%收敛判别
          t=1;break;
      end
      if abs(r2(i)-1)>0.01%收敛判别
          t=1;break;
      end
  end
  if t==0
     break;
  end
   kim=miao(qij,kjm1,Dj,Cij,Z); %迭代
end
return;

 miao.mexw64文件是我当时为了加快运算速度，用c语言写的编译成matlab可执行得文件，点这里可以获取miao

对matlab与c语言混合编程感兴趣可以看我这一篇博客：matlab与c语言混合编程