目录
1 电力系统环境经济调度数学模型
2 知识回顾
3 Matlab代码实现
3.1 Matlab代码
主函数程序：
function MultiObj = GetProblemInfo(TestProblem) %1个多目标工程应用
switch TestProblem
case 1 % IEEE30Bus 环境经济负荷分配问题
numOfObj=3; %目标函数个数
nVar=6; %机组个数【dimension of problem】
VarMin=0.05*ones(1,nVar);
VarMax=1.5*ones(1,nVar);
name='IEEE 30-bus';
CostFunction = @Bus30;
end
%===结构体MultiObj参数=====
MultiObj.nVar=nVar;
MultiObj.var_min = VarMin;
MultiObj.var_max =VarMax;
MultiObj.fun=CostFunction;
MultiObj.numOfObj=numOfObj;
MultiObj.name=name;
end
%% 30节点数据及目标函数
function f=Bus30(x)
%=========发电机参数===============
%第1台机组参数
a(1,1) = 10; %ai
a(1,2) = 200; %bi
a(1,3) = 100; %ci
a(1,4) = 4.091; %aerfa
a(1,5)= -5.543; %beita
a(1,6)= 6.490; %gama
a(1,7) = 2.0 * 10.^(-4); %segama
a(1,8) = 2.857; %yita
%第2台机组参数
a(2,1) = 10;
a(2,2) = 150;
a(2,3) = 120;
a(2,4) = 2.543;
a(2,5)= -6.047;
a(2,6)= 5.638;
a(2,7) = 5.0 *10.^(-4);
a(2,8) = 3.333;
%第3台机组参数
a(3,1) = 20;
a(3,2) = 180;
a(3,3) = 40;
a(3,4) = 4.258;
a(3,5) = -5.094;
a(3,6) = 4.586;
a(3,7) = 1.0 * 10.^(-6);
a(3,8) = 8.000;
%第4台机组参数
a(4,1)= 10;
a(4,2) = 100;
a(4,3) = 60;
a(4,4)= 5.326;
a(4,5) = -3.550;
a(4,6) = 3.380;
a(4,7) = 2.0 * 10.^(-3);
a(4,8) = 2.000;
%第5台机组参数
a(5,1) = 20;
a(5,2) = 180;
a(5,3) = 40;
a(5,4) = 4.258;
a(5,5) = -5.094;
a(5,6) = 4.586;
a(5,7) = 1.0 * 10.^(-6);
a(5,8) = 8.000;
%第6台机组参数
a(6,1) = 10;
a(6,2) = 150;
a(6,3) = 100;
a(6,4) = 6.131;
a(6,5)= -5.555;
a(6,6) = 5.151;
a(6,7) = 1.0 * 10.^(-5);
a(6,8) = 6.667;
%第1台机组Bij参数
B(1,1) = 0.1382;
B(1,2) = -0.0299;
B(1,3) = 0.0044;
B(1,4) = -0.0022;
B(1,5) = -0.0010;
B(1,6) = -0.0008;
%第2台机组Bij参数
B(2,1) = -0.0299;
B(2,2) = 0.0487;
B(2,3) = -0.0025;
B(2,4) = 0.0004;
B(2,5) = 0.0016;
B(2,6) = 0.0041;
%第3台机组Bij参数
B(3,1) = 0.0044;
B(3,2)= -0.0025;
B(3,3) = 0.0182;
B(3,4) = -0.0070;
B(3,5) = -0.0066;
B(3,6) = -0.0066;
%第4台机组Bij参数
B(4,1) = -0.0022;
B(4,2) = 0.0004;
B(4,3) = -0.0070;
B(4,4) = 0.0137;
B(4,5) = 0.0050;
B(4,6) = 0.0033;
%第6台机组Bij参数
B(5,1) = -0.0010;
B(5,2) = 0.0016;
B(5,3) = -0.0066;
B(5,4) = 0.0050;
B(5,5) = 0.0109;
B(5,6) = 0.0005;
%第6台机组Bij参数
B(6,1) = -0.0008;
B(6,2) = 0.0041;
B(6,3) = -0.0066;
B(6,4) = 0.0033;
B(6,5) = 0.0005;
B(6,6) = 0.0244;
B0 = [-0.0107,0.0060,-0.0017,0.0009,0.0002,0.0030 ]; %线损参数B0
B00 = 9.8573 * 10.^(-4); %线损参数B00
f1 = 0.0;
f2 = 0.0;
f3 = 0.0;
%===============总燃料成本======================
for i=1:6
f1 =f1+ (a(i,1) + a(i,2) * x(i) + a(i,3) * x(i).^(2)); %sum(ai+bi*PGi+ci*PGi^2)【xi-PGi】
end
%===============总排放量========================
for i=1:6
f2 = f2 + 0.01 * (a(i,4) + a(i,5) *x(i) + a(i,6) * x(i).^(2)) + a(i,7) * 2.718281828.^(a(i,8) * x(i));
end
%===============总线损====================
for i=1:6
for j=1:6
f3 = f3 + x(i) * B(i,j) * x(j);
end
end
for i=1:6
f3 =f3 + B0(i) * x(i);
end
f3 =f3 + B00;
% %% 不等式约束
% for i=1:6
% if x(i)<PGmin(i)||x(i)>PGmax(i)
% x(i)=unifrnd(PGmin(i),PGmax(i))
% end
%% 等式约束
g = 0;
Pd = 2.834; %PD负荷量
for i=1:6
g =g+ x(i); %sum（PGi）
end
g = g - Pd - f3; %sum（PGi)-PD-PL
g=max(0,abs(g));
%% 目标函数+惩罚项
f(1)=f1+g;
f(2)=f2+0.1*g;
f(3)=f3+g;
end
3.2 结果
部分代码：
function MultiObj = GetProblemInfo(TestProblem) %1个多目标工程应用
switch TestProblem
case 1 % IEEE30Bus 环境经济负荷分配问题
numOfObj=3; %目标函数个数
nVar=6; %机组个数【dimension of problem】
VarMin=0.05*ones(1,nVar);
VarMax=1.5*ones(1,nVar);
name='IEEE 30-bus';
CostFunction = @Bus30;
end
%===结构体MultiObj参数=====
MultiObj.nVar=nVar;
MultiObj.var_min = VarMin;
MultiObj.var_max =VarMax;
MultiObj.fun=CostFunction;
MultiObj.numOfObj=numOfObj;
MultiObj.name=name;
end
%% 30节点数据及目标函数
function f=Bus30(x)
%=========发电机参数===============
%第1台机组参数
a(1,1) = 10; %ai
a(1,2) = 200; %bi
a(1,3) = 100; %ci
a(1,4) = 4.091; %aerfa
a(1,5)= -5.543; %beita
a(1,6)= 6.490; %gama
a(1,7) = 2.0 * 10.^(-4); %segama
a(1,8) = 2.857; %yita
%第2台机组参数
a(2,1) = 10;
a(2,2) = 150;
a(2,3) = 120;
a(2,4) = 2.543;
a(2,5)= -6.047;
a(2,6)= 5.638;
a(2,7) = 5.0 *10.^(-4);
a(2,8) = 3.333;
%第3台机组参数
a(3,1) = 20;
a(3,2) = 180;
a(3,3) = 40;
a(3,4) = 4.258;
a(3,5) = -5.094;
a(3,6) = 4.586;
a(3,7) = 1.0 * 10.^(-6);
a(3,8) = 8.000;
%第4台机组参数
a(4,1)= 10;
a(4,2) = 100;
a(4,3) = 60;
a(4,4)= 5.326;
a(4,5) = -3.550;
a(4,6) = 3.380;
a(4,7) = 2.0 * 10.^(-3);
a(4,8) = 2.000;
%第5台机组参数
a(5,1) = 20;
a(5,2) = 180;
a(5,3) = 40;
a(5,4) = 4.258;
a(5,5) = -5.094;
a(5,6) = 4.586;
a(5,7) = 1.0 * 10.^(-6);
a(5,8) = 8.000;
%第6台机组参数
a(6,1) = 10;
a(6,2) = 150;
a(6,3) = 100;
a(6,4) = 6.131;
a(6,5)= -5.555;
a(6,6) = 5.151;
a(6,7) = 1.0 * 10.^(-5);
a(6,8) = 6.667;%第1台机组Bij参数
B(1,1) = 0.1382;
B(1,2) = -0.0299;
B(1,3) = 0.0044;
B(1,4) = -0.0022;
B(1,5) = -0.0010;
B(1,6) = -0.0008;%第2台机组Bij参数
B(2,1) = -0.0299;
B(2,2) = 0.0487;
B(2,3) = -0.0025;
B(2,4) = 0.0004;
B(2,5) = 0.0016;
B(2,6) = 0.0041;%第3台机组Bij参数
B(3,1) = 0.0044;
B(3,2)= -0.0025;
B(3,3) = 0.0182;
B(3,4) = -0.0070;
B(3,5) = -0.0066;
B(3,6) = -0.0066;%第4台机组Bij参数
B(4,1) = -0.0022;
B(4,2) = 0.0004;
B(4,3) = -0.0070;
B(4,4) = 0.0137;
B(4,5) = 0.0050;
B(4,6) = 0.0033;%第6台机组Bij参数
B(5,1) = -0.0010;
B(5,2) = 0.0016;
B(5,3) = -0.0066;
B(5,4) = 0.0050;
B(5,5) = 0.0109;
B(5,6) = 0.0005;%第6台机组Bij参数
B(6,1) = -0.0008;
B(6,2) = 0.0041;
B(6,3) = -0.0066;
B(6,4) = 0.0033;
B(6,5) = 0.0005;
B(6,6) = 0.0244;
B0 = [-0.0107,0.0060,-0.0017,0.0009,0.0002,0.0030 ]; %线损参数B0
B00 = 9.8573 * 10.^(-4); %线损参数B00f1 = 0.0;
f2 = 0.0;
f3 = 0.0;%===============总燃料成本======================
for i=1:6
f1 =f1+ (a(i,1) + a(i,2) * x(i) + a(i,3) * x(i).^(2)); %sum(ai+bi*PGi+ci*PGi^2)【xi-PGi】
end
%===============总排放量========================
for i=1:6
f2 = f2 + 0.01 * (a(i,4) + a(i,5) *x(i) + a(i,6) * x(i).^(2)) + a(i,7) * 2.718281828.^(a(i,8) * x(i));
end
%===============总线损====================
for i=1:6
for j=1:6
f3 = f3 + x(i) * B(i,j) * x(j);
end
end
for i=1:6
f3 =f3 + B0(i) * x(i);
end
f3 =f3 + B00;% %% 不等式约束
% for i=1:6
% if x(i)<PGmin(i)||x(i)>PGmax(i)
% x(i)=unifrnd(PGmin(i),PGmax(i))
% end%% 等式约束
g = 0;
Pd = 2.834; %PD负荷量
for i=1:6
g =g+ x(i); %sum（PGi）
end
g = g - Pd - f3; %sum（PGi)-PD-PL
g=max(0,abs(g));
%% 目标函数+惩罚项
f(1)=f1+g;
f(2)=f2+0.1*g;
f(3)=f3+g;
end
function MultiObj = GetProblemInfo(TestProblem) %1个多目标工程应用
switch TestProblem
case 1 % IEEE30Bus 环境经济负荷分配问题
numOfObj=3; %目标函数个数
nVar=6; %机组个数【dimension of problem】
VarMin=0.05*ones(1,nVar);
VarMax=1.5*ones(1,nVar);
name='IEEE 30-bus';
CostFunction = @Bus30;
end
%===结构体MultiObj参数=====
MultiObj.nVar=nVar;
MultiObj.var_min = VarMin;
MultiObj.var_max =VarMax;
MultiObj.fun=CostFunction;
MultiObj.numOfObj=numOfObj;
MultiObj.name=name;
end
%% 30节点数据及目标函数
function f=Bus30(x)
%=========发电机参数===============
%第1台机组参数
a(1,1) = 10; %ai
a(1,2) = 200; %bi
a(1,3) = 100; %ci
a(1,4) = 4.091; %aerfa
a(1,5)= -5.543; %beita
a(1,6)= 6.490; %gama
a(1,7) = 2.0 * 10.^(-4); %segama
a(1,8) = 2.857; %yita
%第2台机组参数
a(2,1) = 10;
a(2,2) = 150;
a(2,3) = 120;
a(2,4) = 2.543;
a(2,5)= -6.047;
a(2,6)= 5.638;
a(2,7) = 5.0 *10.^(-4);
a(2,8) = 3.333;
%第3台机组参数
a(3,1) = 20;
a(3,2) = 180;
a(3,3) = 40;
a(3,4) = 4.258;
a(3,5) = -5.094;
a(3,6) = 4.586;
a(3,7) = 1.0 * 10.^(-6);
a(3,8) = 8.000;
%第4台机组参数
a(4,1)= 10;
a(4,2) = 100;
a(4,3) = 60;
a(4,4)= 5.326;
a(4,5) = -3.550;
a(4,6) = 3.380;
a(4,7) = 2.0 * 10.^(-3);
a(4,8) = 2.000;
%第5台机组参数
a(5,1) = 20;
a(5,2) = 180;
a(5,3) = 40;
a(5,4) = 4.258;
a(5,5) = -5.094;
a(5,6) = 4.586;
a(5,7) = 1.0 * 10.^(-6);
a(5,8) = 8.000;
%第6台机组参数
a(6,1) = 10;
a(6,2) = 150;
a(6,3) = 100;
a(6,4) = 6.131;
a(6,5)= -5.555;
a(6,6) = 5.151;
a(6,7) = 1.0 * 10.^(-5);
a(6,8) = 6.667;%第1台机组Bij参数
B(1,1) = 0.1382;
B(1,2) = -0.0299;
B(1,3) = 0.0044;
B(1,4) = -0.0022;
B(1,5) = -0.0010;
B(1,6) = -0.0008;%第2台机组Bij参数
B(2,1) = -0.0299;
B(2,2) = 0.0487;
B(2,3) = -0.0025;
B(2,4) = 0.0004;
B(2,5) = 0.0016;
B(2,6) = 0.0041;%第3台机组Bij参数
B(3,1) = 0.0044;
B(3,2)= -0.0025;
B(3,3) = 0.0182;
B(3,4) = -0.0070;
B(3,5) = -0.0066;
B(3,6) = -0.0066;%第4台机组Bij参数
B(4,1) = -0.0022;
B(4,2) = 0.0004;
B(4,3) = -0.0070;
B(4,4) = 0.0137;
B(4,5) = 0.0050;
B(4,6) = 0.0033;%第6台机组Bij参数
B(5,1) = -0.0010;
B(5,2) = 0.0016;
B(5,3) = -0.0066;
B(5,4) = 0.0050;
B(5,5) = 0.0109;
B(5,6) = 0.0005;%第6台机组Bij参数
B(6,1) = -0.0008;
B(6,2) = 0.0041;
B(6,3) = -0.0066;
B(6,4) = 0.0033;
B(6,5) = 0.0005;
B(6,6) = 0.0244;
B0 = [-0.0107,0.0060,-0.0017,0.0009,0.0002,0.0030 ]; %线损参数B0
B00 = 9.8573 * 10.^(-4); %线损参数B00f1 = 0.0;
f2 = 0.0;
f3 = 0.0;%===============总燃料成本======================
for i=1:6
f1 =f1+ (a(i,1) + a(i,2) * x(i) + a(i,3) * x(i).^(2)); %sum(ai+bi*PGi+ci*PGi^2)【xi-PGi】
end
%===============总排放量========================
for i=1:6
f2 = f2 + 0.01 * (a(i,4) + a(i,5) *x(i) + a(i,6) * x(i).^(2)) + a(i,7) * 2.718281828.^(a(i,8) * x(i));
end
%===============总线损====================
for i=1:6
for j=1:6
f3 = f3 + x(i) * B(i,j) * x(j);
end
end
for i=1:6
f3 =f3 + B0(i) * x(i);
end
f3 =f3 + B00;% %% 不等式约束
% for i=1:6
% if x(i)<PGmin(i)||x(i)>PGmax(i)
% x(i)=unifrnd(PGmin(i),PGmax(i))
% end%% 等式约束
g = 0;
Pd = 2.834; %PD负荷量
for i=1:6
g =g+ x(i); %sum（PGi）
end
g = g - Pd - f3; %sum（PGi)-PD-PL
g=max(0,abs(g));
%% 目标函数+惩罚项
f(1)=f1+g;
f(2)=f2+0.1*g;
f(3)=f3+g;
end
930
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
