一、问题描述
改进发电调度方式又是电力行业节能减排的主要环节。
改进发电调度方式需要在满足负荷需求和功率限制的条件下，使煤耗成本和购电成本尽可能降低，为解决此问题建立双目标优化模型。
二、双目标优化模型
三、模型求解
四、算例分析
以上两张表的来源是
李正哲,马燕峰,娄雅融,韩金铜.基于电力节能减排双目标调度优化模型及方法的研究[J].电力科学与工程,,28(06):44-50.
显然可以实现最初的双目标优化目的。
五、matlab实现
clcclearclose allF = zeros(2);H=zeros(120);%以购电成本最低为目标for i = 1:120if i>0 && i<25H(i,i) = 0.00212;endif i>24 && i<49H(i,i) = 0.00148;endif i>48 && i<73H(i,i) = 0.00289;endif i>72 && i<97H(i,i) = 0.00127;endif i>96 && i<121H(i,i) = 0.00261;endendf=zeros(120,1);for i = 1:120if i>0 && i<25f(i) = 1.8015;endif i>24 && i<49f(i) = 1.213;endif i>48 && i<73f(i) = 1.2643;endif i>72 && i<97f(i) = 1.1954;endif i>96 && i<121f(i) = 1.5354;endendAeq = zeros(24,120);for i = 1:24for j = 0:4Aeq(i,i+j*24) = 1;endendbeq = [1873.01;.01;1852.12;1880.80;1810.95;1786.95;2231.75;1898.98;.98;1849.90;2366.64;2055.99;2298.72;2272.76;2302.72;2195.77;2260.72;2303.72;2190.94;2723.50;2699.50;2693.79;2446.66;2576.55]; LB = zeros(1,120);for i = 1:120if i>0 && i<25LB(i) = 310;endif i>24 && i<49LB(i) = 300;endif i>48 && i<73LB(i) = 350;endif i>72 && i<97LB(i) = 325;endif i>96 && i<121LB(i) = 250;endendUB = zeros(1,120);for i = 1:120if i>0 && i<25UB(i) = 570;endif i>24 && i<49UB(i) = 610;endif i>48 && i<73UB(i) = 700;endif i>72 && i<97UB(i) = 660;endif i>96 && i<121UB(i) = 425;endend [x,fval,exitflag]=quadprog(H,f,[],[],Aeq,beq,LB,UB);plot(x(49:72))xlabel('时间/h')ylabel('负荷/MW')ylim([330 720])title('购电成本目标的机组3日有功负荷曲线')%f11为以购电成本为目标时的购电成本F(1,1) = fval;v=zeros(1,120);for i = 1:120if i>0 && i<25v(i) = 0.15;endif i>24 && i<49v(i) = 0.11;endif i>48 && i<73v(i) = 0.068;endif i>72 && i<97v(i) = 0.20;endif i>96 && i<121v(i) = 0.088;endend%f12为以购电成本为目标时的煤耗成本F(1,2) = v*x;%以煤耗成本为目标[x,fval,exitflag]=linprog(v,[],[],Aeq,beq,LB,UB);figureplot(x(49:72))xlabel('时间/h')ylabel('负荷/MW')ylim([330 720])title('煤耗成本目标的机组3日有功负荷曲线')%f22为以煤耗成本为目标时的煤耗成本F(2,2) = fval;%f21为以煤耗成本为目标时的购电成本F(2,1) = x'*H*x+f'*x;%选取购电成本为基本目标N = 30;epsilon = zeros(1,N+1);for i = 0:Nepsilon(i+1) = F(1,2)-(F(1,2)-F(2,2))/N*i;end%P为帕累托最优解集,PG记录不同解集下的电机出力情况P = zeros(2,31);PG = zeros(120,31);for i = 1:31[x,fval,exitflag]=quadprog(H,f,v,epsilon(i),Aeq,beq,LB,UB);PG(:,i) = x;P(1,i) = fval;P(2,i) = epsilon(i);endfigurescatter(P(1,:),P(2,:))title('帕累托最优解集')ylim([5900 6900])%R1为评价矩阵R1 = P;%R为规格化后的评价矩阵R = zeros(2,31);for i = 1:2for j = 1:31R(i,j) = (max(R1(i,:)-R1(i,j)))/(max(R1(i,:)-min(R1(i,:))));endend%计算熵权alphae = zeros(2,1);alpha = zeros(2,1);for i = 1:2for j = 1:Ne(i,1) = -(sum(R(i,j)/sum(R(i,:))*log(R(i,j)/sum(R(i,:)))))/(log(N));endendfor i = 1:2alpha(i,1) = (1-e(i,1))/(sum([1;1]-e));end%主观权值lambdalambda = [0.5;0.5];%修正权系数omegaomega = zeros(2,1);for i = 1:2omega(i) = (alpha(i)*lambda(i))/(sum(alpha.*lambda));end%建立加权的规格化评级矩阵RjRj = zeros(2,31);for i = 1:2for j = 1:NRj(i,j) = omega(i)*R(i,j);endend%确定双基点Fz = [max(Rj(1,:)),max(Rj(2,:))];Ff = [min(Rj(1,:)),min(Rj(2,:))];%计算帕累托解和理想点之间的距离Dz = zeros(1,31);Df = zeros(1,31);for i = 1:31Dz(i) = ((Rj(1,i)-Fz(1))*(Rj(1,i)-Fz(1))+(Rj(2,i)-Fz(2))*(Rj(2,i)-Fz(2)))^(0.5);Df(i) = ((Rj(1,i)-Ff(1))*(Rj(1,i)-Ff(1))+(Rj(2,i)-Ff(2))*(Rj(2,i)-Ff(2)))^(0.5);end%计算相对贴近度TjTj = zeros(1,31);for i = 1:31Tj(i) = Df(i)/(Dz(i)+Df(i));endpos = find(Tj == max(Tj));val = P(:,pos);%val存放最优解PGz = PG(:,pos);%PGz存放最优解时不同机组在不同时刻的出力情况figureplot(PGz(49:72))xlabel('时间/h')ylabel('负荷/MW')ylim([330 720])title('双目标优化的机组3日有功负荷曲线')插入代码片
后记
我不熟悉在csdn打数学公式，为了方便就直接都放图片了。
在论文基于电力节能减排双目标调度优化模型及方法的研究中所使用的模糊控制我也用代码实现了，虽然确实达到了优化的目的，但运行结束后lambda的值是0，我不清楚是哪里出了问题还是运行结果确实是这样。而且文中提到的对二次约束进行逐步线性化我也不明白是怎么实现的，所以直接用了matlab的fmincon函数进行求解。代码我也放在这里了，有了解的朋友可以私信或者发在评论区，谢谢。
论文中在以购电成本为单目标进行求解的时候得到的结果是115391，但我用matlab的函数quadprog进行求解后只有93828，差距较大，其他数据的差距都比较小，这里我也不清楚是不是有什么问题。
clcclearclose all%c01和c02是我直接用了论文中最后得到的数据c01 = 115391;c02 = 5899.5;%线性不等式约束A=ones(1,121);for i = 1:120if i>0 && i<25A(i) = 0.15;endif i>24 && i<49A(i) = 0.11;endif i>48 && i<73A(i) = 0.068;endif i>72 && i<97A(i) = 0.20;endif i>96 && i<121A(i) = 0.088;endendA(121) = 0.1*c02;b = 1.1*c02;%线性等式约束Aeq = zeros(24,121);for i = 1:24for j = 0:4Aeq(i,i+j*24) = 1;endendbeq = [1873.01;.01;1852.12;1880.80;1810.95;1786.95;2231.75;1898.98;.98;1849.90;2366.64;2055.99;2298.72;2272.76;2302.72;2195.77;2260.72;2303.72;2190.94;2723.50;2699.50;2693.79;2446.66;2576.55]; %上下界LB = zeros(1,121);for i = 1:120if i>0 && i<25LB(i) = 310;endif i>24 && i<49LB(i) = 300;endif i>48 && i<73LB(i) = 350;endif i>72 && i<97LB(i) = 325;endif i>96 && i<121LB(i) = 250;endendLB(121) = 0;UB = zeros(1,121);for i = 1:120if i>0 && i<25UB(i) = 570;endif i>24 && i<49UB(i) = 610;endif i>48 && i<73UB(i) = 700;endif i>72 && i<97UB(i) = 660;endif i>96 && i<121UB(i) = 425;endend UB(121) = 1;[x,fval]=fmincon('fun1',zeros(121,1),A,b,Aeq,beq,LB,UB,'fun2')H=zeros(121);for i = 1:120if i>0 && i<25H(i,i) = 0.00212;endif i>24 && i<49H(i,i) = 0.00148;endif i>48 && i<73H(i,i) = 0.00289;endif i>72 && i<97H(i,i) = 0.00127;endif i>96 && i<121H(i,i) = 0.00261;endendf=zeros(121,1);for i = 1:120if i>0 && i<25f(i) = 1.8015;endif i>24 && i<49f(i) = 1.213;endif i>48 && i<73f(i) = 1.2643;endif i>72 && i<97f(i) = 1.1954;endif i>96 && i<121f(i) = 1.5354;endendx'*H*x+f'*xv=zeros(1,121);for i = 1:120if i>0 && i<25v(i) = 0.15;endif i>24 && i<49v(i) = 0.11;endif i>48 && i<73v(i) = 0.068;endif i>72 && i<97v(i) = 0.20;endif i>96 && i<121v(i) = 0.088;endendv*xplot(x(49:72))
fun1函数
function f=fun1(x)f=-x(121);end
fun2函数
function [c,ceq] =fun2(x)c01 = 93828;%c01 = 115391;H=zeros(121);for i = 1:120if i>0 && i<25H(i,i) = 0.00212;endif i>24 && i<49H(i,i) = 0.00148;endif i>48 && i<73H(i,i) = 0.00289;endif i>72 && i<97H(i,i) = 0.00127;endif i>96 && i<121H(i,i) = 0.00261;endendf=zeros(121,1);for i = 1:120if i>0 && i<25f(i) = 1.8015;endif i>24 && i<49f(i) = 1.213;endif i>48 && i<73f(i) = 1.2643;endif i>72 && i<97f(i) = 1.1954;endif i>96 && i<121f(i) = 1.5354;endendc = x'*H*x+f'*x+11539.1*x(121)-1.1*c01;ceq = 0*x;end
