基于粒子群算法的三态马尔可夫过程数据包服务模型

根据题目要求，我们需要将原代码中的MMOO分布改写成三态的马尔可夫过程，并使用粒子群算法求解三态马尔可夫过程的参数。下面给出改写后的代码。

首先是largedeviationMMOO_01函数的改写。原函数中用到了MMOO分布的概率密度函数，我们需要将其改为三态的马尔可夫过程的概率转移矩阵：

function [fx]=largedeviation_01(p,q,a,x,y,g,m,b)
T = [1-p, p*(1-a), p*a; (1-q)*(1-g), q+p*(1-m-g), p*m+p*g+a*q; (1-q)*g, q+b, 1-q-b];
[z,y1]=eig(T);
[v,w] = size(y1);
max_ii=y1(1,1);
maxy=1;
for i=2:v
    if (y1(i,i)>max_ii)
        max_ii=y1(i,i);
        maxy=i;
    end
end
h=z(:,maxy);
fx = h(1)*(1-p)/(2-p-q) + h(2)*(1-q)*(1-a)/(2-p-q) + h(3)*(1-q)*a/(2-p-q);
end

然后是Threshold_MMOO_01函数的改写。由于MMOO分布的参数x和y是三态马尔可夫过程中的g和b，我们只需要将原函数中求解x和y的代码改为求解g和b的代码即可：

function [threshold01]=Threshold_01(p,q,a,Dmax,g,m,b)
theta0 = [0.1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5];
f = @(theta) largedeviation_01(p,q,a,g,m,1-g-m,b,theta);
for j = 1:length(theta0)
    theta_tmp(j) = fsolve(f,[theta0(j)],optimset('Display','off')); 
end                   
max=theta_tmp(1);
for i=2:length(theta0)
    if (theta_tmp(i)>max)
        max=theta_tmp(i);
    end
end
theta1=max;
sp1=(q+p*(1-a)+p*a*exp(theta1)+sqrt((q+p*(1-a)+p*a.*exp(theta1)).^2-4*(q*(p*(1-a)+p*a.*exp(theta1))-((1-q)*(1-a)+(1-q)*a.*exp(theta1))*(1-p))))./2;
T1=[q (1-q)*(1-g) (1-q)*g; (1-p) p*(1-m-g) p*m+p*g+a*q; (1-p) p*(1-m-g) p*m+p*g+a*q];
[z1,y1]=eig(T1);
[v1,w1] = size(y1);
max_ii1=y1(1,1);
max1=1;
for i=2:v1
    if (y1(i,i)>max_ii1)
        max_ii1=y1(i,i);
        max1=i;
    end
end
ha=z1(:,max1);
Eha0=ha(1)*(1-p)/(2-p-q)+ha(2)*(1-q)*(1-a)/(2-p-q)+ha(3)*(1-q)*a/(2-p-q);

sp2=(g+m*exp(-theta1)+sqrt((g+m*exp(-theta1)).^2-4*(m*g*exp(-theta1)-(1-m)*(1-g)*exp(-theta1))))/2;
T2=[m*exp(-theta1) (1-m);(1-g)*exp(-theta1) g];
[z2,y2]=eig(T2);
[v2,w2] = size(y2); 
max_ii2=y2(1,1);
max2=1;
for i=2:v2
    if (y2(i,i)>max_ii2)
        max_ii2=y2(i,i);
        max2=i;
    end
end
hs=z2(:,max2);
Ehs0=hs(1)*(1-g)/(2-m-g)+hs(2)*m/(2-m-g);
Ks=-(log((g+m*exp(-theta1)+sqrt((g+m*exp(-theta1)).^2-4*(m*g*exp(-theta1)-(1-m)*(1-g)*exp(-theta1))))/2))/theta1;
H=ha(3)*hs(2);
threshold01=(Eha0*Ehs0*exp(-theta1*Ks*Dmax))/H;
end

最后是主函数的改写。需要将原函数中用于求解MMOO分布参数的粒子群算法改为求解三态马尔可夫过程参数的粒子群算法。具体地，我们需要设置三个参数g,m,b作为粒子的位置向量，将原函数中的x,y参数改为g,b，同时将Threshold_MMOO_01函数改为Threshold_01函数。另外，我们还需要判断求解得到的三个参数g,m,b的数据包服务概率是否大于IBP过程得到的数据包到达概率，并在满足条件时输出结果。

clc;
clear;

p =  0.9;
q =  0.8;
a =  0.8;

DLee=0;
PLRee=0;
Dmax=25;
K=100;
M=5;
threshold=10^(-2);

itermax=100;

num_particle=10;   
c1=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)=0.5*rand(num_particle,1)+0.1; 
    position(:,2)=0.5*rand(num_particle,1)+0.1;
    position(:,3)=0.5*rand(num_particle,1)+0.1;
    echo=0;
    pbestp=position;       
    for i=1:1:num_particle
        pbestf(i,1)=Threshold_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;
        for i=1:1:num_particle         
            v(i,:)=w.*v(i,:)+c1.*rand().*(gbestp-position(i,:))+c2.*rand().*(pbestp(i,:)-position(i,:));
            position(i,:)