基于粒子群算法的IBP-Markov数据包到达与服务模型阈值计算

clc; clear; p = 0.9; q = 0.8; a = 0.8; % DLee and PLRee are end-to-end delay,end-to-end packet loss probability. DLee=0; PLRee=0; Dmax=25; K=100; M=5; threshold=10^(-2);

itermax=100; %最大迭代次数100

num_particle=10;
c1=2;%学习因子,典型取值都取2，一般通过试凑调节这个值 c2=2; for j = 1:1:M arr_rate(1,j)=a*(1-q)/(2-p-q); v=zeros(num_particle,3);
position=zeros(num_particle,3); pbestf=zeros(num_particle,1); fitness=zeros(num_particle,1); record=zeros(itermax,1);
position(:,1)=rand(num_particle,1); %g position(:,2)=rand(num_particle,1)*0.5; %m position(:,3)=rand(num_particle,1); %b echo=0; pbestp=position; %所有粒子搜索到的最优位置（个体最优解） for i=1:1:num_particle pbestf(i,1)=Threshold_Markov_01(p,q,a,Dmax,position(i,1),position(i,2),position(i,3)); end gbestf(1,j)=pbestf(1,1); id=1; for i=2:1:num_particle if abs(threshold-pbestf(i,1))<abs(threshold-gbestf(1,j)) gbestf(1,j)=pbestf(i,1); id=i; end end

gbestp=position(id,:);         %群体搜索到的最优位置（群体最优解）
for echo=1:1:itermax
    echo;
    w=0.9-0.5*echo/itermax;%0.7298;   
    for i=1:1:num_particle         %速度更新公式
        v(i,:)=w.*v(i,:)+c1.*rand().*(gbestp-position(i,:))+c2.*rand().*(pbestp(i,:)-position(i,:));
        position(i,:)=position(i,:)+v(i,:);   %位置更新
        if position(i,1)<0
            position(i,1)=0;
        elseif  position(i,1)>1
            position(i,1)=1;
        end
        if position(i,2)<0
            position(i,2)=0;
        elseif  position(i,2)>0.5
            position(i,2)=0.5;
        end
        if position(i,3)<0
            position(i,3)=0;
        elseif  position(i,3)>1
            position(i,3)=1;
        end
        fitness(i,1)=Threshold_Markov_01(p,q,a,Dmax,position(i,1),position(i,2),position(i,3));
        if abs(threshold-fitness(i,1))<abs(threshold-pbestf(i,1))          %第i个粒子搜索到的最优位置的适应值（优化目标函数的值）
            pbestf(i,1)=fitness(i,1);
            pbestp(i,:)=position(i,:);
        end
        if abs(threshold-fitness(i,1))<abs(threshold-gbestf(1,j)) 
            gbestf(1,j)=fitness(i,1);          %群体搜索到的最优位置的适应值
            gbestp=position(i,:);
        end
    end
    record(echo)=gbestf(1,j);
end
x(1,j)=gbestp(1,1);
y(1,j)=gbestp(1,2);
z(1,j)=gbestp(1,3);
arr_ser(1,j)=(1-y(1,j))/(2-x(1,j)-y(1,j)-z(1,j));
rho=arr_rate(1,j)/arr_ser(1,j);%负载率

end g m b function [threshold01]=Threshold_Markov_01(p,q,a,Dmax,g,m,b) T=[1-p,p,0;1-q,0,q;0,(1-m)/2,(1-b)/2]; [z,y]=eig(T); [v,w] = size(y); max_ii=y(1,1); max=1; for i=2:v if (y(i,i)>max_ii) max_ii=y(i,i); max=i; end end h=z(:,max); Eh0=h(1)(1-p)/(2-p-q)+h(2)(1-q)/(2-p-q)+h(3)(1-m-b)/2; Ks=2/(2-p-q); threshold01=(Eh0exp(-Ks*Dmax))/(g+m+b); end

function [fx]=largedeviationMarkov_01(p,q,a,g,m,b,theta) fx=((1-g)(1-a)+g+(g+m)exp(theta)-(1-g-m-b)exp(-theta))/2((q+p(1-a)+paexp(theta)+sqrt((q+p(1-a)+paexp(theta)).^2-4*(q*(p*(1-a)+paexp(theta))-((1-q)*(1-a)+(1-q)*a.exp(theta))(1-p))))./2)-1; end