粒子群算法求解三态马尔可夫过程参数

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)=0.5rand(num_particle,1)+0.1; % g的初始值为0.3 position(:,2)=0.5rand(num_particle,1)+0.1; % m的初始值为0.3 position(:,3)=0.5rand(num_particle,1)+0.1; % b的初始值为0.3 echo=0; pbestp=position;
for i=1:1:num_particle % 计算适应值 g=position(i,1); m=position(i,2); b=position(i,3); if (g+m+b)>1 fitness(i,1)=Inf; % 限制g+m+b的和小于等于1 else P=[g (1-g) 0; (1-m)/2 m (1-m)/2; 0 (1-b) b]; [~,d]=eig(P); fitness(i,1)=abs(d(1,1)); % d的最大特征值对应的适应值 end if abs(threshold-fitness(i,1))<abs(threshold-pbestf(i,1))
% 更新个体最优解 pbestf(i,1)=fitness(i,1); pbestp(i,:)=position(i,:); end end [~,id]=min(abs(threshold-pbestf)); % 找到个体最优解中适应值最接近阈值的位置 gbestf=pbestf(id); % 群体最优解的适应值 gbestp=pbestp(id,:); % 群体最优解的位置 for echo=1:1:itermax echo w=0.9-0.5echo/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,:)=position(i,:)+v(i,:); %位置更新 % 限制参数范围 if position(i,1)<0.1 position(i,1)=0.1; elseif position(i,1)>0.9 position(i,1)=0.9; end if position(i,2)<0.1 position(i,2)=0.1; elseif position(i,2)>0.9 position(i,2)=0.9; end if position(i,3)<0.1 position(i,3)=0.1; elseif position(i,3)>0.9 position(i,3)=0.9; end % 计算适应值 g=position(i,1); m=position(i,2); b=position(i,3); if (g+m+b)>1 fitness(i,1)=Inf; else P=[g (1-g) 0; (1-m)/2 m (1-m)/2; 0 (1-b) b]; [~,d]=eig(P); fitness(i,1)=abs(d(1,1)); end % 更新个体最优解和群体最优解 if abs(threshold-fitness(i,1))<abs(threshold-pbestf(i,1)) pbestf(i,1)=fitness(i,1); pbestp(i,:)=position(i,:); end if abs(threshold-fitness(i,1))<abs(threshold-gbestf) gbestf=fitness(i,1); gbestp=position(i,:); end end record(echo)=gbestf; end g=gbestp(1); m=gbestp(2); b=gbestp(3); arr_ser(1,j)=(1-b)/(2-g-b-m); rho=arr_rate(1,j)/arr_ser(1,j); end

function [threshold01]=Threshold_MMOO_01(p,q,a,Dmax,x,y) theta0 = [0.1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5]; f = @(theta) largedeviationMMOO_01(p,q,a,x,y,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; %IBP sp1=(q+p*(1-a)+paexp(theta1)+sqrt((q+p*(1-a)+pa.exp(theta1)).^2-4(q(p*(1-a)+pa.exp(theta1))-((1-q)(1-a)+(1-q)a.exp(theta1))(1-p))))./2; T1=[q (1-q)(1-a) (1-q)aexp(theta1);(1-p) p(1-a) paexp(theta1);(1-p) p*(1-a) paexp(theta1)]; [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); ha(1); ha(2); ha(3); 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=(y+x.exp(-theta1)+sqrt((y+x.exp(-theta1)).^2-4(x.yexp(-theta1)-(1-x).(1-y).exp(-theta1))))/2; T2=[xexp(-theta1) (1-x);(1-y)exp(-theta1) y]; [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-y)./(2-x-y)+hs(2)(1-x)./(2-x-y); Ks=-(log((y+x.exp(-theta1)+sqrt((y+x.exp(-theta1)).^2-4(x.yexp(-theta1)-(1-x).(1-y).exp(-theta1))))/2))/theta1; H=ha(3)hs(2); threshold01=(Eha0Ehs0exp(-theta1KsDmax))/H; end

function [fx]=largedeviationMMOO_01(p,q,a,x,y,theta) fx=((y+x.exp(-theta)+sqrt((y+x.exp(-theta)).^2-4(x.yexp(-theta)-(1-x).(1-y)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)aexp(theta))(1-p))))./2)-1;
end