1 主要内容

这是一篇代码勘误，该程序复现文献《基于多目标粒子群算法的微电网优化调度》，代码实现了多目标优化调度模型——考虑微电网系统运行成本和环境保护成本的并网模式下微电网多目标优化调度模型，同时采用改进的粒子群算法对优化模型进行求解，具体改进的点包括：1）改进了惯性因子；2）在PSO算法中引入变异操作进行改进粒子群算法，在普通粒子群算法的基础上引入简单变异算子，在粒子每次更新之后，以一定概率重新初始化粒子。简单来说，该程序实现了风光、储能、柴油、燃气、电网交互等主体参与的微网优化调度，采用改进多目标粒子群进行优化，原代码是包括传统粒子群和改进粒子群两部分，但是多目标粒子群在实现过程中存在较大的问题，下面就来做一个说明。（下载网址包括原代码（3.9元）和【原代码+勘误】（99元）两种方式代码，大家可以按需下载！）

问题简述

运行主程序main_improve，得到如下的结果。

为了方便展示，将结果图采用subplot方式进行绘制。

可以注意到，最后一个子图是关于各个设备出力情况，折线图有点乱，改成柱状图+折线图来看一下效果。

上面这个图就能发现问题——功率严重不平衡，正常设备出力（在图中即柱状图部分）应该和负荷（折现）保持一致，这也是电力系统最重要的属性，但是结果和实际相去甚远，可见程序代码存在问题。

原因分析

多目标pareto解集定义如下：

即对于双目标模型而言，只要有一个目标值最小，另一个目标值再大也满足有效解的条件。回到程序上面，单目标算法对于约束处理的方式是采用罚函数形式，即在目标后面加一个较大的惩罚性的数值，当不满足约束条件时，罚函数起作用，在寻找最优解的过程中逐渐把不满足约束条件的解淘汰；多目标算法处理约束最好不采用罚函数的方式，如上面蓝字部分，不满足解集的也很有可能被作为有效解，虽然当前采用一个状态变量作为是否可行的标志，但这个方法仍有较大的偏差裕度和不准确性，程序的问题就出在这里。

%% 判断是否为可行解
if ele_sum>4500
    c=1;
else
    c=0;
end

该程序对于功率平衡可行的判断条件是小于4500，这个裕度太大了，导致该微网优化结果功率失衡。

result=C_grid+C_OM_F+f_DE*DE_sum+f_MT*MT_sum+f_SOC*SOC_sum_delt+f_ele*ele_sum;

其一目标包含有罚函数，该目标函数求解结果准确性低，导致成本数据计算不准确。对上述问题的勘误就需要修改约束条件编程方式，采用绝对约束方式进行，保证每一个解集均满足约束条件，这样才能保证最终结果的准确性。

2 部分代码

%程序初始化
clear;
clc;
close all;
%定义全局变量
global P_load; %电负荷
global WT;%风电
global PV;%光伏
global buy_price;
global sell_price;
%获取数据
data=xlsread('mopso_data');
P_load=data(:,1);
PV=data(:,2);
WT=data(:,3);
buy_price=data(:,4);
sell_price=data(:,5);

%蓄电池最大放电功率（正表示为电负荷供电，即放电）
BESSMax_dischar=30;
%蓄电池最大充电功率
BESSMax_char=-30;
%柴油机最大发电功率
DEMax=30;
%柴油机最小发电功率
DEMin=6;
%燃气轮机最大发电功率
MTMax=30;
%燃气轮机最小发电功率
MTMin=3;
%主网交互最大功率(正表示为电负荷供电)
GridMax=30;
%主网交互最小功率
GridMin=-30;
%% 调用mopso函数
mm=mopso; %调用mopso函数
nn=length(mm.swarm); %非支配解数目

%% 比较不同目标函数寻优对调度结果的影响
%% 第1种.将两个目标函数值归一化相加，取相加后最小的目标值的粒子，即寻找折衷解并画图
%将非支配解中的运行成本和环境保护成本分别赋值给yyy,xxx
% for i=1:nn
%    yyy(i)= mm.swarm(1,i).cost(1);
%    xxx(i)= mm.swarm(1,i).cost(2);
% end
% m1=max(yyy);
% m2=max(xxx);
% 
% for i=1:nn
%     object(i)= mm.swarm(1,i).cost(1)./m1+ mm.swarm(1,i).cost(2)./m2;
%     f1(i)=mm.swarm(1,i).cost(1)./m1;
%     f2(i)=mm.swarm(1,i).cost(2)./m2;
% end
% [m,p]=min(object); %得到有着最小目标值的微粒所在的行数P
% pg=mm.swarm(1,p).x; %pg为折衷解
% Title = sprintf('折衷解情况下');
%% 第2种寻找总成本最低时的解并画图
for i=1:nn
    object(i)= mm.swarm(1,i).cost(1)+ mm.swarm(1,i).cost(2);
end
[m,p]=min(object); %得到有着最小目标值的微粒所在的行数P
pg=mm.swarm(1,p).x; %pg为总成本最低时的解
Title = sprintf('总成本最低情况下');