【优化求解】粒子群求解微电网多目标问题matlab源码

多目标粒子群(MOPSO)算法是由CarlosA. Coello Coello等在2004年提出来的，详细参考1。目的是将原来只能用在单目标上的粒子群算法(PSO)应用于多目标上。我们知道原来的单目标PSO流程很简单：

　　-->初始化粒子位置(一般都是随机生成均匀分布)

　　-->计算适应度值(一般是目标函数值-优化的对象)

　　-->初始化历史最优pbest为其本身和找出全局最优gbest

　　-->根据位置和速度公式进行位置和速度的更新

　　-->重新计算适应度

　　-->根据适应度更新历史最优pbest和全局最优gbest

　　-->收敛或者达到最大迭代次数则退出算法

　　速度的更新公式如下：

　　等式右边有三部分组成。第一部分是惯性量，是延续粒子上一次运动的矢量；第二部分是个体认知量，是向个体历史最优位置运动的量；第三部分是社会认知量，是粒子向全局最优位置运动的量。

有了速度，则位置更新自然出来了：

　　以上是对于多目标PSO算法的介绍。运用到多目标上去的话，出现的问题有以下几点：

1. 如何选择pbest。我们知道对于单目标优化来说选择pbest，只需要对比一下就可以选择出哪个较优。但是对于多目标来说两个粒子的对比，并不能对比出哪个好一些。如果粒子的每个目标都要好的话，则该粒子更优。若有些更好，有些更差的话，就无法严格的说哪个好些，哪个差一些。

2. 如何选择gbest。我们知道对于单目标在种群中只有一个最优的个体。而对于多目标来说，最优的个体有很多个。而对PSO来说，每个粒子只能选择一个作为最优的个体(领带者)。该如何选择呢？

   　　MOPSO对于第一个问题的做法是在不能严格对比出哪个好一些时随机选择一个其中一个作为历史最优。对于第二个问题，MOPSO则在最优集里面(存档中)根据拥挤程度选择一个领导者。尽量选择不那么密集位置的粒子(在这里用到了网格法)。

   　　MOPSO在选择领导者和对存档(也可以说是pareto临时最优断面)进行更新的时候应用了自适应网格法，详细参考2。

   　　如何选择领带者呢？

   　　MOPSO在存档中选择一个粒子跟随。如何选择呢？根据网格划分，假设每个网格中粒子数个，i代表第几个网格。该网格中的粒子被选择的概率为，即粒子越拥挤，则选择的概率越低。这是为了保证能够对未知的区域进行探索。

   　　如何进行存档呢？

   　　在种群更新完成之后，是如何进行存档的呢？MOPSO进行了三轮筛选。

   　　首先，根据支配关系进行第一轮筛选，将劣解去除，剩下的加入到存档中。

   　　其次，在存档中根据支配关系进行第二轮筛选，将劣解去除，并计算存档粒子在网格中的位置。

   　　最后，若存档数量超过了存档阀值，则根据自适应网格进行筛选，直到阀值限额为止。重新进行网格划分。

   　　refer:

3. Handling multiple objectives with particle swarm optimization

4. Approximating the non dominated front using the Pareto archivedevolution strateg

5. ``` %%注意，目标函数的编写有可能会对粒子存档数产生很大影响，可能会使粒子存档数不会正常增加，产生波动，而达不到设置的你Rep

   clc;

   clear ;

   close all;

   tic

   %% Problem Definition

   %% MOPSO Settings

   % number_obj=3;%目标函数为3个

   nPop=100; % Population Size种群大小100

   nRep=100; % Repository Size存档大小为100

   MaxIt=150; % Maximum Number of Iterations原始最大迭代次数100

   phi1=2.05;%phi1=phi2=2.05

   phi2=2.05;

   phi=phi1+phi2;%4.1

   chi=2/(phi-2+sqrt(phi^2-4*phi));%0.7298

   w=chi; % Inertia Weight惯性权重

   wdamp=1; % Inertia Weight Damping Ratio惯性权重阻尼比1

   c1=chi*phi1; % Personal Learning Coefficient个体学习因子1.4962

   c2=chi*phi2; % Global Learning Coefficient群体学习因子1.4962

   nVar=55;%可以理解为粒子的维数,24+24+7

   % VarMin=[0 0 0 0.0033];%粒子能取到的最小值

   % VarMax=[0.00195 0.0009 0.00365 0.0037]; %粒子能取到的最大值

   load data_load

   load data_sun

   % VarMax=repmat(2.2,1,24);VarMax=[VarMax,repmat(0.1,1,24)];VarMax=[VarMax,dataload(4,1:12)];VarMax=[VarMax,dataload(7,1:12)];

   % VarMin=repmat(0,1,24);VarMin=[VarMin,repmat(-0.1,1,24)];VarMin=[VarMin,repmat(0,1,12)];VarMin=[VarMin,repmat(0,1,12)];%可平移负荷设置为节点7

   VarMax=repmat(0.1,1,24);VarMax=[VarMax,repmat(0.05,1,24)];VarMax=[VarMax,dataload(7,11:14),dataload(7,19:21)];

   VarMin=repmat(-0.1,1,24);VarMin=[VarMin,repmat(-0.05,1,24)];VarMin=[VarMin,repmat(0,1,7)];%可平移负荷设置为节点7

   end

   %% Results

   figure

   plot(record);

   xlabel('迭代次数');ylabel('存档最优目标函数值');

   title('存档最优目标函数值迭代过程图');

   hold on

   %测试rep中粒子位置是否符合储能限制条件

   for i=1:100

   repflag(i)=batteryflag(rep(i).Position);

   end

   costs=GetCosts(particle);

   rep_costs=GetCosts(rep);

   [final final_num]=min(integrate);

   bestposition=rep(finalnum).Best.Position;

   % plot3(bestposition(1),bestposition(2),best_position(3),'rP','MarkerSize',12);

   % hold on

   figure;

   subplot(2,1,1)

   plot(best_position(1,1:24),'-*');

   title('节点10储能运行计划图')

   xlabel('时间/h');ylabel('充放电功率/Mw');

   grid on

   xlim([1 24]);

   subplot(2,1,2)

   plot(best_position(1,25:48),'-*');

   grid on

   xlim([1 24]);

   title('节点11储能运行计划图')

   xlabel('时间/h');ylabel('充放电功率/Mw');

   % figure

   % temp=dataload(7,1:12)-bestposition(1,49:60);

   % load7=[temp(1,5:12)+dataload(7,1:8),bestposition(1,49:60),temp(1,1:4)+data_load(7,21:24)];

   % plot(data_load(7,1:24));hold on

   % plot(load7,'-*');hold on;

   % xlim([1 24]);

   % legend('平移前','平移后')

   % title('负荷7平移前与平移后日负荷曲线')

   % xlabel('时间/h');ylabel('功率/Mw');

   % grid on

   figure

   temp1=dataload(7,11:14)-bestposition(1,49:52);

   temp2=dataload(7,19:21)-bestposition(1,53:55);

   % temp=[temp1,temp2];

   load7=[temp1+dataload(7,1:4),temp2+dataload(7,5:7),dataload(7,8:11),bestposition(1,49:52),dataload(7,16:18),bestposition(1,53:55),data_load(7,22:24)];

   plot(data_load(7,1:24));hold on

   plot(load7,'-*');hold on;

   xlim([1 24]);

   legend('平移前','平移后')

   title('负荷7平移前与平移后日负荷曲线')

   xlabel('时间/h');ylabel('功率/Mw');

   grid on

   check=zeros(100,1);

   for number=1:100

   check(number)=particle(number).Dominated;

   end

   num_position=find(check==0);

   num=numel(num_position);

   clear i;

   figure;

   bestcost=rep(finalnum).Best.Cost;

   plot3(bestcost(1),bestcost(2),best_cost(3),'rP','MarkerSize',16)

   hold on

   plot3(costs(1,:),costs(2,:),costs(3,:),'b.');

   hold on;

   plot3(repcosts(1,:),repcosts(2,:),rep_costs(3,:),'gx');

   % legend('Best Repository','Main Population','Repository');

   legend('存档最优值','主要种群值','存档值');

   xlabel('用于负荷平移的用户电费/元');ylabel('电网侧成本/元');zlabel('电压偏移');

   grid on;

   title('种群和存档最优解对应目标函数的值')

   if it>MaxIt

   disp(['迭代次数达到限制最大值100,','跳出迭代']);

   end

   disp(['最优解为:',num2str(best_position)]);

   disp(['最优值为:',num2str(best_cost)]);

   toc ```