飞道的博客

基于多目标遗传算法(NSGA-II)和多目标粒子群算法(MOPSO)的分布式仿真系统双目标负载平衡模型【Matlab代码实现】

514人阅读  评论(0)

 💥💥💥💞💞💞欢迎来到本博客❤️❤️❤️💥💥💥
 

📋📋📋本文目录如下:⛳️⛳️⛳️

目录

1 概述

2 数学模型

3 多目标遗传算法和多目标粒子群算法

4 运行结果

5 Matlab代码+数据+文档讲解 


1 概述

分布式系统具有高利用率、高可扩展性和并行性的优点,可以并行运行以处理计算密集型问题。此外,从分布式结构中衍生出几个难题,例如节点之间的通信问题、系统资源和计算节点的映射和分配问题。

高层体系结构(HLA)广泛用于分布式仿真开发。虽然HLA支持大规模分布式仿真,但仿真过程中的负载平衡机制尚未标准化。当仿真系统在非优化分配条件下运行时,由于计算量不足,会出现负载不平衡问题。

负载平衡的结果表现在以下两个方面:一是有效缩短任务的平均响应时间,即尽量缩短仿真过程的运行时间,使系统能够稳定、平稳地运行;其次,最大限度地利用整个系统的仿真资源,减少不必要的资源浪费。

高层体系结构(HLA)是分布式仿真系统的软件体系结构规范,不涉及负载平衡。因此,在HLA的仿真过程中,无法解决仿真时间长和仿真任务分布不均匀引起的失真问题。本文将研究基于组件级的仿真系统,并解决负载不平衡问题。目标是找到一组帕累托最优解,以最小化模型中具有负载限制约束的计算负载和总通信负载的不平衡。为了解决这个问题,提出了一个新的整数规划模型。采用带精英策略的非支配排序遗传算法(NSGA-II)和多目标粒子群优化算法(MOPSO)求解该问题。分析了MOPSO的不同全局最优选择方法(拥挤距离、自适应网格和综合排序)和扰动方法(快速下降和精英学习策略)。由于算法的参数对其性能有显著影响,因此使用具有新响应值的田口方法来调整所提出算法的参数。五个性能指标用于评估这些算法的结果。基于这些结果,将使用网络控制仿真平台来测试和验证这些方案。数值结果表明,提出的带ELS算子的多粒子群优化算法在求解该问题上优于其他提出的算法。此外,所提出的策略可以解决基于HLA的系统中的负载不平衡问题。

2 数学模型

本文考虑了一个具有计算和通信负载限制的双目标负载平衡问题。如前所述,该问题受到分布式仿真系统中负载分配问题的启发,即HLA。在基于HLA的分布式仿真系统中,仿真任务可以分为联邦成员级和组件级。

                   

等式(1)最小化了模拟节点计算负载的方差。等式(2)最小化了仿真节点之间的总通信负载。等式(3)规定每个模拟任务被精确地分配给一个模拟节点。等式(4)确保了计算负载限制。等式(5)保证了通信负载限制。等式(6)表示变量的非负性和完整性。更多细节参考第五部分,文件夹里面有详细讲解。

3 多目标遗传算法和多目标粒子群算法

见第五部分,详细文章讲解。

4 运行结果

​部分代码: 


  
  1. function REP = MOPSOCD(params,MultiObj)
  2. % % 参数
  3. Np = params.Np;
  4. Nr = params.Nr;
  5. maxgen = params.maxgen;
  6. W = params.W;
  7. C1 = params.C1;
  8. C2 = params.C2;
  9. %ngrid = params.ngrid;
  10. maxvel = params.maxvel;
  11. %u_mut = params.u_mut;
  12. fun = MultiObj.fun;
  13. nVar = MultiObj.nVar;
  14. var_min = MultiObj.var_min(:);
  15. var_max = MultiObj.var_max(:);
  16. d =nVar; %搜索空间维数(未知数个数)
  17. xmin =var_min;
  18. xmax =var_max;
  19. xrange =xmax -xmin;
  20. n =Np; %初始化群体个体数目
  21. MaxDT =maxgen; %最大迭代次数
  22. rp =[]; %外部存档,档案里存的是每代保存下来的非支配解,也是算法求解的结果
  23. rpmax =Nr; %档案规模
  24. c1 =C1;c2 =C2;
  25. w =W; %惯性权重
  26. % pbx n *d 矩阵 ---个体最优解
  27. % gbx n *d 矩阵 ---全局最优解(每个粒子的全局最优解都不同)
  28. % pbf n *k 矩阵 ----个体最优值
  29. % gbf n *k 矩阵 ----全局最优值
  30. %产生随机 n 个粒子,初始化粒子的位置和速度
  31. %x =xmin +xrange. *rand(n,d);
  32. x =[];
  33. for i = 1:n
  34. x(i,:) =xmin '+xrange'. *rand( 1,d);
  35. end
  36. v =zeros(n,d);
  37. maxvel = (var_max -var_min). *maxvel. / 100;
  38. tmp = fun(x);
  39. % % 计算粒子的适应度
  40. k =size(tmp, 2); %目标函数个数
  41. x(:, 2 *d + 1: 2 *d +k) =tmp;
  42. x(:,d + 1: 2 *d) =v;
  43. % % 求个体极值和全局极值
  44. pbx =x(:, 1:d); %个体最优位置为初始粒子本身
  45. pbf =x(:, 2 *d + 1: 2 *d +k);
  46. gbx =x(:, 1:d); %全局最优位置为初始粒子本身
  47. gbf =x(:, 2 *d + 1: 2 *d +k);
  48. xk =x(:, 1:d);
  49. v =x(:,(d + 1): 2 *d);
  50. rp = x;
  51. % % 可视化
  52. if(k = = 2)
  53. h_fig = figure( 1);
  54. h_par = plot(x(:, 2 *d + 1),x(:, 2 *d + 2), 'or'); hold on;
  55. h_rep = plot(rp(:, 2 *d + 1),rp(:, 2 *d + 2), 'ok'); hold on;
  56. try
  57. set(gca, 'xtick',REP.hypercube_limits(:, 1) ','ytick ',REP.hypercube_limits(:,2)');
  58. axis([ min(REP.hypercube_limits(:, 1)) max(REP.hypercube_limits(:, 1)) ...
  59. min(REP.hypercube_limits(:, 2)) max(REP.hypercube_limits(:, 2))]);
  60. grid on; xlabel( 'f1'); ylabel( 'f2');
  61. end
  62. drawnow;
  63. end
  64. if(k = = 3)
  65. h_fig = figure( 1);
  66. h_par = plot3(POS_fit(:, 1),POS_fit(:, 2),POS_fit(:, 3), 'or'); hold on;
  67. h_rep = plot3(REP.pos_fit(:, 1),REP.pos_fit(:, 2),REP.pos_fit(:, 3), 'ok'); hold on;
  68. try
  69. set(gca, 'xtick',REP.hypercube_limits(:, 1) ','ytick ',REP.hypercube_limits(:,2)', 'ztick',REP.hypercube_limits(:, 3) ');
  70. axis([min(REP.hypercube_limits(:,1)) max(REP.hypercube_limits(:,1)) ...
  71. min(REP.hypercube_limits(:,2)) max(REP.hypercube_limits(:,2))]);
  72. end
  73. grid on; xlabel('f1 '); ylabel('f2 '); zlabel('f3 ');
  74. drawnow;
  75. axis square;
  76. end
  77. display(['Generation 0 - Repository size: ' num2str(size(rp,1))]);
  78. %开始循环
  79. for t=1:MaxDT
  80. %根据粒子群公式迭代,位置、速度更新
  81. gen = t;
  82. for i=1:n
  83. %v(i,:)=0.3.*randn(1,d)+c1.*rand.*(pbx(i,:)-xk(i,:))+c2.*rand.*(gbx(i,:)-xk(i,:));
  84. v(i,:)=w.*v(i,:)+c1.*rand.*(pbx(i,:)-xk(i,:))+c2.*rand.*(gbx(i,:)-xk(i,:));
  85. for j=1:d
  86. if v(i,j)>maxvel(j)
  87. v(i,j)=maxvel(j);
  88. elseif v(i,j)<-maxvel(j)
  89. v(i,j)=-maxvel(j);
  90. end
  91. end
  92. xk(i,:)=xk(i,:)+v(i,:);
  93. for j=1:d
  94. if xk(i,j)>xmax(j)
  95. xk(i,j)=xmax(j);
  96. elseif xk(i,j)<xmin(j)
  97. xk(i,j)=xmin(j);
  98. end
  99. end
  100. end
  101. for i=1:n
  102. mu=0.1; % Mutation Rate
  103. pm=(1-(t-1)/(MaxDT-1))^(1/mu);
  104. if rand<pm
  105. xk(i,:) = Mutate(xk(i,:),pm,xmin,xmax);
  106. end
  107. for j=1:d
  108. if xk(i,j)>xmax(j)
  109. xk(i,j)=xmax(j);
  110. elseif xk(i,j)<xmin(j)
  111. xk(i,j)=xmin(j);
  112. end
  113. end
  114. end
  115. x(:,1:d)=xk;
  116. x(:,(d+1):2*d)=v;
  117. %计算粒子的适应度
  118. x(:,2*d+1:2*d+k)=fun(xk);
  119. %求非支配解集 np(快速排序法)
  120. q=x;np=[];
  121. while isempty(q)==0
  122. xsb=q(1,:);
  123. q(1,:)=[];
  124. flag=1; %若 flag==1,则 xsb 是非支配解,就可以放到非支配解集中了
  125. [column,row]=size(q);
  126. for i=1:column
  127. dom_less=0;
  128. dom_equal=0;
  129. dom_more=0;
  130. for l=1:k
  131. if xsb(2*d+l)<q(i,2*d+l)
  132. dom_less=dom_less+1;
  133. elseif xsb(2*d+l)==q(i,2*d+l)
  134. dom_equal=dom_equal+1;
  135. else
  136. dom_more=dom_more+1;
  137. end
  138. end
  139. %若 xsb 支配 q(i),则将 q(i)所在的行标记为 0,用来到最后删除此行做标记
  140. if dom_more==0&dom_equal~=k
  141. q(i,:)=0;
  142. %若 q(i)支配 xsb,那么 xsb 不是非支配解,不能被放入非支配解集
  143. elseif dom_less==0&dom_equal~=k
  144. flag=0;break
  145. end
  146. end
  147. if flag==1 %若 xsb 是非支配解,则将其放入非支配解集 np 中
  148. np=[np;xsb];
  149. end
  150. %将 q 中标记为 0 的行删除,剩下不被 xsb 支配的粒子,即快速排序的简便之处
  151. q(~any(q,2),:)=[];
  152. end
  153. %更新个体极值(若当前位置支配其个体极值位置,则更新为当前位置)
  154. for i=1:n
  155. dom_less=0;
  156. dom_equal=0;
  157. dom_more=0;
  158. for l=1:k
  159. if x(i,2*d+l)<pbf(i,l)
  160. dom_less=dom_less+1;
  161. elseif x(i,2*d+l)==pbf(i,l)
  162. dom_equal=dom_equal+1;
  163. else
  164. dom_more=dom_more+1;
  165. end
  166. end
  167. if dom_more==0&dom_equal~=k
  168. pbf(i,:)=x(i,2*d+1:2*d+k);
  169. pbx(i,:)=x(i,1:d);
  170. end
  171. end
  172. %更新外部集 rp(将非支配集按支配关系插入外部集)
  173. [column,row]=size(rp);
  174. if column==0
  175. rp=np;
  176. else
  177. [column2,row2]=size(np);
  178. for i=1:column2
  179. flag=1; %若 flag==1,则 np(i,:)是就可以放到外部集中了
  180. [column1,row1]=size(rp);
  181. for j=1:column1
  182. dom_less=0;
  183. dom_equal=0;
  184. dom_more=0;
  185. for l=1:k
  186. if np(i,2*d+l)<rp(j,2*d+l)
  187. dom_less=dom_less+1;
  188. elseif np(i,2*d+l)==rp(j,2*d+l)
  189. dom_equal=dom_equal+1;
  190. else
  191. dom_more=dom_more+1;
  192. end
  193. end
  194. %若非支配集中的粒子 np(i,:)支配外部集中的粒子 rp(j,:),
  195. %则将 rp(j,:)所在的行标记为 0,用来到最后删除此行做标记
  196. if dom_more==0&dom_equal~=k
  197. rp(j,:)=0;
  198. %若 rp(j,:)支配 np(i,:),则表示 np(i,:)被 rp(j,:)支配,
  199. %那么 np(i,:)就一定不是非支配解,就不能被放入非支配解集中了
  200. elseif dom_less==0&dom_equal~=k
  201. flag=0;break
  202. end
  203. end
  204. if flag==1 %若 flag==1,则 np(i,:)是就可以放到外部集中了
  205. rp=[rp;np(i,:)];
  206. end
  207. rp(~any(rp,2),:)=[];
  208. end
  209. end
  210. np = rp;
  211. [column,row]=size(rp);
  212. for i=1:column
  213. for j=i:column
  214. if i~=j && sum(np(i,1:d)-np(j,1:d) == 0) == d
  215. rp(i,:)=0;
  216. break;
  217. end
  218. end
  219. end
  220. rp(~any(rp,2),:)=[];
  221. %基于拥挤距离控制档案规模
  222. %(计算外部集中粒子的拥挤距离,对拥挤距离进行降序排序,删除后面的多余个体)
  223. [column,row]=size(rp);
  224. indexrp=[1:column]';
  225. rp =[rp indexrp]; %索引
  226. rp(:, 2 *d +k + 2) =zeros( column, 1); %用来保存拥挤距离
  227. for i = 1:k
  228. rp =sortrows(rp, 2 *d +i); %将粒子进行对第 i 个目标函数值进行升序排列
  229. if column > 2
  230. for j = 2: column -1
  231. rp(j, 2 *d +k + 2) = rp(j, 2 *d +k + 2) +(rp(j + 1, 2 *d +i) -rp(j -1, 2 *d +i)) /(rp( column, 2 *d +i) -rp( 1, 2 *d +i));
  232. end
  233. end
  234. rp( 1, 2 *d +k + 2) =rp( 1, 2 *d +k + 2) +inf;
  235. rp( column, 2 *d +k + 2) =rp( column, 2 *d +k + 2) +inf;
  236. end
  237. rp =sortrows(rp, -( 2 *d +k + 2)); %外部集根据拥挤距离降序排序
  238. if column >rpmax
  239. rp =rp( 1:rpmax,:);
  240. end
  241. %更新全局极值(此时外部集已根据拥挤距离降序排序,从外部集的最前部分随机选择)
  242. for i = 1:n
  243. %产生一个随机数,从外部集前 10 %的粒子里面选
  244. randomnumber = ceil(rand * column / 10);
  245. gbx(i,:) =rp(randomnumber, 1:d); %全局最优位置
  246. gbf(i,:) =rp(randomnumber, 2 *d + 1: 2 *d +k);
  247. end
  248. rp(:, 2 *d +k + 1: 2 *d +k + 2) =[];
  249. % Plotting and verbose
  250. if(k = = 2)
  251. figure(h_fig); delete(h_par); delete(h_rep);
  252. h_par = plot(x(:, 2 *d + 1),x(:, 2 *d + 2), 'or'); hold on;
  253. h_rep = plot(rp(:, 2 *d + 1),rp(:, 2 *d + 2), 'ok'); hold on;
  254. try
  255. set(gca, 'xtick',REP.hypercube_limits(:, 1) ','ytick ',REP.hypercube_limits(:,2)');
  256. axis([ min(REP.hypercube_limits(:, 1)) max(REP.hypercube_limits(:, 1)) ...
  257. min(REP.hypercube_limits(:, 2)) max(REP.hypercube_limits(:, 2))]);
  258. end
  259. if(isfield(MultiObj, 'truePF'))
  260. try delete(h_pf); end
  261. h_pf = plot(MultiObj.truePF(:, 1),MultiObj.truePF(:, 2), '.', 'color', 0.8. *ones( 1, 3)); hold on;
  262. end
  263. grid on; xlabel( 'f1'); ylabel( 'f2');
  264. drawnow;
  265. axis square;
  266. end
  267. if(k = = 3) % 以后再改
  268. figure(h_fig); delete(h_par); delete(h_rep);
  269. h_par = plot3(POS_fit(:, 1),POS_fit(:, 2),POS_fit(:, 3), 'or'); hold on;
  270. h_rep = plot3(REP.pos_fit(:, 1),REP.pos_fit(:, 2),REP.pos_fit(:, 3), 'ok'); hold on;
  271. try
  272. set(gca, 'xtick',REP.hypercube_limits(:, 1) ','ytick ',REP.hypercube_limits(:,2)', 'ztick',REP.hypercube_limits(:, 3) ');
  273. axis([min(REP.hypercube_limits(:,1)) max(REP.hypercube_limits(:,1)) ...
  274. min(REP.hypercube_limits(:,2)) max(REP.hypercube_limits(:,2)) ...
  275. min(REP.hypercube_limits(:,3)) max(REP.hypercube_limits(:,3))]);
  276. end
  277. if(isfield(MultiObj,'truePF '))
  278. try delete(h_pf); end
  279. h_pf = plot3(MultiObj.truePF(:,1),MultiObj.truePF(:,2),MultiObj.truePF(:,3),'. ','color ',0.8.*ones(1,3)); hold on;
  280. end
  281. grid on; xlabel('f1 '); ylabel('f2 '); zlabel('f3 ');
  282. drawnow;
  283. axis square;
  284. end
  285. display(['Generation # ' num2str(gen) ' - Repository size: ' num2str(size(rp,1))]);
  286. end
  287. REP.pos = rp(:,1:d);
  288. REP.pos_fit = rp(:,2*d+1:2*d+k);
  289. hold off;
  290. % hold on;
  291. % rpcd=rp;
  292. % plot(rpcd(:,2*d+1),rpcd(:,2*d+2),' ^ ');
  293. end
  294. function xnew=Mutate(x,pm,VarMin,VarMax)
  295. nVar=numel(x);
  296. j=randi([1 nVar]);
  297. dx=pm*(VarMax(j)-VarMin(j));
  298. lb=x(j)-dx;
  299. ub=x(j)+dx;
  300. lb = max(lb,VarMin(j));
  301. lb = min(lb,VarMax(j));
  302. ub = max(ub,VarMin(j));
  303. ub = min(ub,VarMax(j));
  304. xnew=x;
  305. xnew(j)=unifrnd(lb,ub);
  306. end

5 Matlab代码+数据+文档讲解 

回复关键字


转载:https://blog.csdn.net/weixin_61181717/article/details/128094748
查看评论
* 以上用户言论只代表其个人观点,不代表本网站的观点或立场