排队论算法的matlab实现

   朋友画了一天写的排队论算法，用来解决如下问题，不忍让它蒙尘，便发了出来让有缘人拜读拜读。

以下是某城区地理简图，直线表示可行路径——公路，每个方格表示一个区域（边长都为1公里），所有区域都进行了编号，共47个区域，其中46、47区为天然屏障，无法直接穿越。该城区仅有的送餐公司C总部现设在25区，有送餐员6名。现统计了某天该城区的顾客下单情况（见附件）。假设送餐员每次送餐结束后，总是先返回公司总部，（等待）领取下一个订单（先到先得，按下单时间顺序领取），然后再前往订单中指定的餐馆取餐，最后到下单顾客的地址送餐。送餐员在公路上的行驶速度为30公里每小时，在到达一个区的边界后停下，从打听顾客具体位置到最终完成送餐所需平均时间是3分钟、方差为1分钟。
1 2 3 4 5 6 7 8 9 10
11 46 46 14 15 16 17 18 19 20
21 22 46 24 25 26 27 47 47 30
31 32 33 34 35 36 37 38 39 40
41 42 43 44 45 12 13 23 28 29

1） 结合附件研究该城区顾客下单规律，估计顾客的平均等待时间。
（找到一个模型描述顾客等待时间与送餐的关系，然后得到一个平均关系）——（思路是通过模拟建模——蒙特卡洛法 即仿真法）

1. Step1 建立送餐地点和时间的概率模型——比较难弄这个。。。
2. Step2 ：建模 等待时间=F（x1,x2,x3,x4）
3. 模拟10000+次 ，从而确定平均时间

p=zeros(256,1);
c=6;
z=ones(c,1);
w=zeros(256,1);
f=zeros(256,1);
p(1:6,1)=0;
w(1:6,1)=0;
s=zeros(1000,1);
for i=1:c
    j=1;  
    M=10000;
    L=10000;  %运行次数
    x1=guwei(i,1);  
    y1=guwei(i,2);
    gw=zeros(M,2);   
    gw(:,1)=x1;
    gw(:,2)=y1;  %得到顾客的位置矩阵
    cw=gw+cgxd;  %得到该顾客位置下餐馆的随机位置矩阵
    while(j<=L)
        if cw(j,1)<1||cw(j,1)>10||cw(j,2)<1||cw(j,2)>5   %将餐馆随机位置矩阵中不符合实际条件的位置剔除
            cw=[cw(1:j-1,:);cw(j+1:L,:)];
            L=L-1;
        else
            j=j+1;
        end
    end
[m,~]=size(cw);   %得到处理后的餐馆随机位置矩阵的行数
zw=zeros(m,2);    
zw(:,1)=3;
zw(:,2)=5;       %得到总部位置矩阵
gw=zeros(m,2);   
gw(:,1)=x1;
gw(:,2)=y1;       %得到顾客位置矩阵
czxd=abs(cw-zw);  %得到餐馆与总部的相对距离矩阵（分x方向与y方向）
gcxd1=abs(gw-cw);   %得到顾客与餐馆的相对距离矩阵 
gzxd=abs(gw-zw);
juli=czxd(:,1)+czxd(:,2)+gcxd1(:,1)+gcxd1(:,2)+gzxd(:,1)+gzxd(:,2);   %总距离矩阵
f(i)=sum(juli)/m/0.5;   %得到该顾客位置下的一个平均距离，将其放到矩阵mj中
end
for j=1:c
    s(j)=f(j)+t(j);   %前c个人无需排队
end
for i=c+1:256;   %从第c+1个人开始模拟
    j=1;  
    M=10000;
    L=10000;  %运行次数
    x1=guwei(i,1);  
    y1=guwei(i,2);
    gw=zeros(M,2);   
    gw(:,1)=x1;
    gw(:,2)=y1;  %得到顾客的位置矩阵
    cw=gw+cgxd;  %得到该顾客位置下餐馆的随机位置矩阵
    while(j<=L)
        if cw(j,1)<1||cw(j,1)>10||cw(j,2)<1||cw(j,2)>5   %将餐馆随机位置矩阵中不符合实际条件的位置剔除
            cw=[cw(1:j-1,:);cw(j+1:L,:)];
            L=L-1;
        else
            j=j+1;
        end
    end
[m,~]=size(cw);   %得到处理后的餐馆随机位置矩阵的行数
zw=zeros(m,2);    
zw(:,1)=3;
zw(:,2)=5;       %得到总部位置矩阵
gw=zeros(m,2);   
gw(:,1)=x1;
gw(:,2)=y1;       %得到顾客位置矩阵
czxd=abs(cw-zw);  %得到餐馆与总部的相对距离矩阵（分x方向与y方向）
gcxd1=abs(gw-cw);   %得到顾客与餐馆的相对距离矩阵
gzxd=abs(gw-zw);
juli=czxd(:,1)+czxd(:,2)+gcxd1(:,1)+gcxd1(:,2)+gzxd(:,1)+gzxd(:,2);   %总距离矩阵
f(i)=sum(juli)/m/0.5;   %得到该顾客位置下的一个平均距离，将其放到矩阵f中
     if p(i-1)>0    %当队长大于零时
        for j=1:c   %对服务台状态进行分析
            if s(j)<t(i)||s(j)==t(i)
                z(j)=0;
                s(j)=0;
            end         %挑出空闲服务台
        end
                k=p(i-1);  %队长赋值给k
                while k>0
                    j=1;
                    while j<=c
                        if z(j)==0
                           z(j)=1;
                           s(j)=t(i-k)+w(i-k)+f(i-k);
                           break;
                        end
                        j=j+1;
                    end
                    if j==c+1
                        break;
                    end
                    k=k-1;
                end       %处理第i个人之前等待人的状态
                if k>0
                    if p(i-1)-k>0
                    for j=1:k
                        s(c+j)=s(c+p(i-1)-k+j);
                    end
                    s((c+k+1):(c+p(i-1)))=0;
                    end
                    temp=s;
                    for l=1:k
                        temmin=temp(1);
                        tag=1;
                        for j=1:c
                        if temmin>temp(j)
                            temmin=temp(j);
                            tag=j;    
                        end
                        end
                        temp=[temp(1:(tag-1),1);temp((tag+1):(c+k-l+1),1)];
                    end
                    temmin=temp(1);
                    for j=1:c
                        if temmin>temp(j)
                            temmin=temp(j);
                        end
                    end
                    w(i)=temmin-t(i);
                    p(i)=k+1;
                    s(c+k+1)=t(i)+w(i)+f(i);      %得到等待时间，更新s,p,w
                else
                    s(7:c+p(i-1))=0;
                    for j=1:c
                    if s(j)<t(i)||s(j)==t(i)
                        z(j)=0;
                        s(j)=0;
                    end
                    end
                    j=1;
                    while j<=c
                        if z(j)==0
                            z(j)=1;
                            s(j)=t(i)+f(i);
                            w(i)=0;
                            p(i)=0;
                            break;
                        end
                        j=j+1;
                    end
                    if j==c+1
                        temmin=s(1);
                        for j=1:c
                            if temmin>s(j)
                                temmin=s(j);
                            end
                        end
                        w(i)=temmin-t(i);
                        p(i)=1;
                        s(c+1)=t(i)+w(i)+f(i);
                    end
                end
    else
         for j=1:c
            if s(j)<t(i)&&z(j)~=0
                z(j)=0;
                s(j)=0;
            end
         end
         j=1;
         while j<=c
             if z(j)==0
                 z(j)=1;
                 s(j)=t(i)+f(i);
                 w(i)=0;
                 break;
             end
               j=j+1;
         end
         if j==c+1
             p(i)=p(i-1)+1;
             temmin=s(1);
             for j=1:c
                 if temmin>s(j)
                     temmin=s(j);
                 end
             end
             w(i)=temmin-t(i);
             s(c+1)=t(i)+w(i)+f(i);
         else
             p(i)=0;
         end
    end
end
mu=sum(w)/256;
mutime=sum(f)/256+sum(w)/256;