vip1.m

经典的排队论程序

function [Pv1,Pv2,Pv3,Pd1,Pd2,Pd3,PL]=vip1(av1,ad3,bv,bd,c1,c2,c3,d1,d2,d3,c,s)
%clear all;
%av1=0.07;%:0.005:0.11;
%ad3=0.18;
%avv=av1;
%avv1=av1;
%avv2=av1;
av2=0.25*av1;
av3=0.1*av1;
ad2=0.25*ad3;
ad1=0.1*ad3;
avv=av1+av2+av3;
avv1=av2+av3;
avv2=av3;
add=ad1+ad2+ad3;
add1=ad2+ad3;
add2=ad3;
%add=ad1;
%add1=add;
%add2=add;
for i=1:s
if  i>=1&i<=d1;  
    A1(i,1)=-add;
elseif i>=d1+1&i<=d2; 
    A1(i,1)=-add1;
elseif  i>=d2+1&i<=d3; 
    A1(i,1)=-add2;    
end 
end
A1=[A1;0];
for i=1:s
if  i>=1&i<=d1;  
    B1(i,1)=add+avv+min(c,i-1)*bd;
elseif i>=d1+1&i<=d2; 
    B1(i,1)=add1+avv+min(c,i-1)*bd;
elseif i>=d2+1&i<=d3; 
    B1(i,1)=add2+avv+min(c,i-1)*bd;    
end 
end
 B1(s+1,1)=avv+min(c,s)*bd;
 B1;
 for i=1:s
    C1(i,1)=-min(c,i)*bd;
 end
C1=[0;C1];
D1=spdiags([A1,B1,C1],[-1,0,1],s+1,s+1); 
for i=1:s+1
    E1(i,1)=-bv;
end
F1=spdiags([E1],[0],s+1,s+1);
F1;
H1=zeros(s+1,(c3-1)*(s+1));
H1;
J{1,1}=[D1,F1,H1];
J{1,1};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for x=1:c1-1
  J{x+1,1}=xunh1(x,av1,ad3,bv,bd,c1,c2,c3,d1,d2,d3,c,s);
end
%the numberof living channels  equal c1
A3=zeros(s+1,(c1-1)*(s+1));
 for  i=1:s+1
     B3(i,1)=-avv;
 end
 B3=spdiags([B3],[0],s+1,s+1);
   for i=1:s
if  i>=1&i<=d1;  
    C3(i,1)=-add;
elseif i>=d1+1&i<=d2; 
    C3(i,1)=-add1;
elseif  i>=d2+1&i<=d3; 
    C3(i,1)=-add2;    
end 
end
C3=[C3;0];
 for i=1:s
if  i>=1&i<=d1;  
    D3(i,1)=add+avv1+c1*bv+min(c-c1,i-1)*bd;
elseif i>=d1+1&i<=d2; 
    D3(i,1)=add1+avv1+c1*bv+min(c-c1,i-1)*bd;
elseif i>=d2+1&i<=d3; 
    D3(i,1)=add2+avv1+c1*bv+min(c-c1,i-1)*bd;    
end 
end
 D3(s+1,1)=avv1+c1*bv+min(c-c1,s)*bd;  
for i=1:s
    E3(i,1)=-min(c-c1,i)*bd;
end  
E3=[0;E3];
F3=spdiags([C3,D3,E3],[-1,0,1],s+1,s+1);
for i=1:s+1
    G3(i,1)=-(c1+1)*bv;
end
G3=spdiags([G3],[0],s+1,s+1);
H3=zeros(s+1,(c3-1-c1)*(s+1));
J{c1+1,1}=[A3,B3,F3,G3,H3];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for x=c1+1:c2-1
     J{x+1,1}=xunh2(x,av1,ad3,bv,bd,c1,c2,c3,d1,d2,d3,c,s);
end 
%the numberof living channels  equal c2
A5=zeros(s+1,(c2-1)*(s+1));
 for  i=1:s+1
     B5(i,1)=-avv1;
 end
 B5=spdiags([B5],[0],s+1,s+1);
  for i=1:s
if  i>=1&i<=d1;  
    C5(i,1)=-add;
elseif i>=d1+1&i<=d2; 
    C5(i,1)=-add1;
elseif  i>=d2+1&i<=d3; 
    C5(i,1)=-add2;    
end 
end
C5=[C5;0];
 for i=1:s
if  i>=1&i<=d1;  
    D5(i,1)=add+avv2+c2*bv+min(c-c2,i-1)*bd;
elseif i>=d1+1&i<=d2; 
    D5(i,1)=add1+avv2+c2*bv+min(c-c2,i-1)*bd;
elseif i>=d2+1&i<=d3; 
   D5(i,1)=add2+avv2+c2*bv+min(c-c2,i-1)*bd;    
end 
end
 D5(s+1,1)=avv2+c2*bv+min(c-c2,s)*bd;   
for i=1:s
    E5(i,1)=-min(c-c2,i)*bd;
end  
E5=[0;E5];
F5=spdiags([C5,D5,E5],[-1,0,1],s+1,s+1);
for i=1:s+1
    G5(i,1)=-(c2+1)*bv;
end
G5=spdiags([G5],[0],s+1,s+1);
H5=zeros(s+1,(c3-1-c2)*(s+1));
J{c2+1,1}=[A5,B5,F5,G5,H5];
%%%%%%%%%%%%%%%%%%%%%%%%%
for x=c2+1:c3-1
    J{x+1,1}=xunh3(x,av1,ad3,bv,bd,c1,c2,c3,d1,d2,d3,c,s);
end
%the numberof living channels  equal c3
 A7=zeros(s+1,(c3-1)*(s+1));
 for  i=1:s+1
     B7(i,1)=-avv2;
 end
 B7=spdiags([B7],[0],s+1,s+1);
  for i=1:s
if  i>=1&i<=d1;  
    C7(i,1)=-add;
elseif i>=d1+1&i<=d2; 
    C7(i,1)=-add1;
elseif  i>=d2+1&i<=d3; 
    C7(i,1)=-add2;    
end 
end
C7=[C7;0];
for i=1:s
if  i>=1&i<=d1;  
    D7(i,1)=add+c3*bv+min(c-c3,i-1)*bd;
elseif i>=d1+1&i<=d2; 
    D7(i,1)=add1+c3*bv+min(c-c3,i-1)*bd;
elseif i>=d2+1&i<=d3; 
   D7(i,1)=add2+c3*bv+min(c-c3,i-1)*bd;    
end 
end
 D7(s+1,1)=c3*bv+min(c-c3,s)*bd;   
for i=1:s
    E7(i,1)=-min(c-c3,i)*bd;
end  
E7=[0;E7];
F7=spdiags([C7,D7,E7],[-1,0,1],s+1,s+1);
H7=zeros(s+1,(c3-c3)*(s+1));
J{c3+1,1}=[A7,B7,F7,H7];
%%%%%%%%%%%%%%%%%%%%%%%%%555
 J{c3+2,1}=ones(1,(c3+1)*(s+1));
 J=cell2mat(J(:,1));
 J=full(J);
 L=zeros((c3+1)*(s+1)+1,1);
 L((c3+1)*(s+1)+1,1)=1;
 G=J\L;
 for i=1:c3+1
     for j=1:s+1
        p(i,j)=G((i-1)*s+j);
     end
 end
 %p=[];
 %p=reshape(M,c3+1,s+1);
 %vpa(p,10)
 Pv1=0;
 for i=c1+1:c3+1;
     for j=1:s+1
     Pv1=Pv1+p(i,j);
 end
end
Pv2=0;
for i=c2+1:c3+1
    for j=1:s+1
        Pv2=Pv2+p(i,j);
    end
end
Pv3=0;
for j=1:s+1
    Pv3=Pv3+p(c3+1,j);
end
 Pd1=0;
 for i=1:c3+1
     for j=d1+1:d3+1;
     Pd1=Pd1+p(i,j);
 end
end
Pd2=0;
for i=1:c3+1
    for j=d2+1:d3+1
        Pd2=Pd2+p(i,j);
    end
end
Pd3=0;
for i=1:c3+1
    Pd3=Pd3+p(i,d3+1);
end
Pd3;
%G=0.8*(Pvn+10*Pvh)+0.2*(Pbn+10*Pbh);
PL=0;
for i=1:c3+1
   for  j=1:s+1
  PL=PL+j*p(i,j);
end
end
%T=PL/(a3*(1-Pd));