基于混合规划求解的6机30节点的机组组合优化

一、问题重述

针对6机30节点的系统，求解计划时间内各支路的功率情况与机组的开停机情况，使得系统总成本达到最小。该问题的决策变量由两类，第一类是各时段机组的出力，为连续变量。第二类是各时段机组的启停状态，为整数变量，0表示关停，1表示启动。本问题属于混合整数规划（MIP）问题，即要在决策变量的可行解空间里找到一组最优解，使得目标函数尽可能取得极值。

二、模型建立

2.1 目标函数

目标函数即为最小化成本，包括发电带来的煤耗成本和机组启停产生的开停机成本。

2.2 约束条件

考虑了有功功率平衡约束、机组出力约束、机组爬坡约束、机组启停约束、潮流约束和启停费用约束。

（a）有功功率平衡约束

（b）热备用约束

（c）机组出力约束

（d）机组爬坡约束

（e）机组启停约束

（f）潮流约束

2.3 符号说明

三、部分代码展示

clear
clc
%%系统参数
%所有参数均用有名值表示
paragen=xlsread('excel2017','机组参数');
loadcurve=xlsread('excel2017','负荷曲线');
netpara=xlsread('excel2017','网络参数');
branch_num=size(netpara);%网络中的支路
branch_num=branch_num(1,1);
PL_max=netpara(:,6);%线路最大负荷
PL_min=netpara(:,7);%线路最小负荷
limit=paragen(:,3:4);%机组出力上下限//limit(:,1)表示上限，limit(:,2)表示下限
para=paragen(:,5:7);%成本系数//para(:,1)表示系数a,para(:,2)表示系数b,para(:,3)表示系数c。
price=100;
para=price*para;%价格换算
lasttime=paragen(:,9);%持续时间
Rud=paragen(:,8);%上下爬坡速率//因题中简化上坡下坡速度相同
H=paragen(:,10);%启动成本
J=paragen(:,11);%关停成本
u0=[1 1 1 1 1 1];%初始状态
%% 规模变量
%机组数
gennum=size(paragen);
gennum=gennum(1,1);
%节点数
numnodes=size(loadcurve);
numnodes=numnodes(1,1)-1;
%时间范围
T=size(loadcurve);
T=T(1,2)-1;
%线性化分段数(按需要更改)
m=4;
%各时刻节点总负荷
PL=loadcurve(numnodes+1,2:T+1);
%%
%决策变量
u=binvar(gennum,T,'full');%状态变量
p=sdpvar(gennum,T,'full');%即各机组实时功率p(i,t)
Ps=sdpvar(gennum,T,m,'full');%分段出力
costH=sdpvar(gennum,T,'full');%启动成本
costJ=sdpvar(gennum,T,'full');%关停成本
sum_PowerGSDF=sdpvar(T,branch_num,numnodes,'full');%发电机的输出功率转移总和
%% 目标函数
totalcost=0;%机组费用成本最小
%线性化的最优成本目标
for i=1:gennum
for t=1:T
for s=1:m
    totalcost=totalcost+K(i,s)*Ps(i,t,s);%线性化煤耗成本
end
    totalcost=totalcost+u(i,t)*(para(i,2)*limit(i,2)+para(i,1)*limit(i,2)^2+para(i,3));%加上表示机组开机并以最小出力 运行产生的煤耗
    totalcost=totalcost+costH(i,t)+costJ(i,t);%加上机组启停产生的开停机成本
end
end
%原二次函数式的最优成本目标
% for i=1:gennum
%     for t=1:T
%     totalcost=totalcost+para(i,1)*p(i,t).^2+para(i,2)*p(i,t)+para(i,3)*u(i,t);  %煤耗成本
%     totalcost=totalcost+costH(i,t);                                %启动成本
%     totalcost=totalcost+costJ(i,t);                                %关停成本
%     end
% end
%%
for t=1:T
st=st+[sum(p(:,t))==PL(1,t)];%负荷平衡约束;
end
%%
%% 机组爬坡约束
%按下式进行推导编程
% %启动最大升速率
% Su=(Pmax+Pmin)/2;
% %停机最大降速率
% Sd=(Pmax+Pmin)/2;
%Ru=Rud;Rd=Rud;
% %上爬坡约束
% for t=2:T
% st=st+[p(:,t)-p(:,t-1)<=u(:,t-1).*(Ru-Su)+Su];
% end
% %下爬坡约束
% for t=2:T
%st=st+[p(:,t-1)-p(:,t)<=u(:,t).*(Rd-Sd)+Sd];
% end
%展开表达式：
for t=2:T
    for i=1:gennum
    % st=st+[-Rud(i,1)*u(i,t)+(u(i,t)-u(i,t-1))*limit(i,2)-limit(i,1)*(1-u(i,t))<=p(i,t)-p(i,t-1)];
    % st=st+[p(i,t)-p(i,t-1)<=Rud(i,1)*u(i,t-1)+(u(i,t)-u(i,t-1))*limit(i,2)+limit(i,1)*(1-u(i,t))];
    %由于原式可能关机以后就无法再开动了，改用下式
    st=st+[p(i,t-1)-p(i,t)<=Rud(i,1)*u(i,t)+(1-u(i,t))*(limit(i,2)+limit(i,1))/2];%下坡
    st=st+[p(i,t)-p(i,t-1)<=Rud(i,1)*u(i,t-1)+(1-u(i,t-1))*(limit(i,2)+limit(i,1))/2];%上坡
    end
end
%% 热备用约束
hp=0.05;%热备用系数
for t=1:T
st=st+[sum(u(:,t).*limit(:,1)-p(:,t))>=hp*PL(1,t)];
end
%% 启停时间约束
%启动约束
for t=2:T
    for i=1:gennum
        indicator=u(i,t)-u(i,t-1);%启停时间约束的简化表达式（自己推导的）,indicator为1表示启动，为0表示停止
        range=t:min(T,t+lasttime(i)-1);
        st=st+[u(i,range)>=indicator];
    end
end
%停机约束
for t=2:T
    for i=1:gennum
        indicator=u(i,t-1)-u(i,t);%启停时间约束
        range=t:min(T,t+lasttime(i)-1);%特别限制时间上限
        st=st+[u(i,range)<=1-indicator];
    end
end
%% 启停成本约束
for t=1:T   %启停成本零限约束
    for i=1:gennum
        st=st+[costH(i,t)>=0]; 
        st=st+[costJ(i,t)>=0];
    end
end
for i=1:gennum  %启停成本条件约束
   for t=2:T
         st=st+[costH(i,t)>=H(i,1)*(u(i,t)-u(i,t-1))];
         st=st+[costJ(i,t)>=J(i,1)*(u(i,t-1)-u(i,t))];
   end
    st=st+[costH(i,1)>=H(i,1)*(u(i,1)-u0(1,i))];%初始状态下的启停成本
    st=st+[costJ(i,1)>=J(i,1)*(u0(1,i)-u(i,1))];
end
%% 直流潮流约束
%% 直流潮流下的导纳矩阵节点参数初始化
netpara(:,4)=1./netpara(:,4);%电抗求倒数成电纳
slack_bus=26;%按不同的平衡节点号更改
Y=zeros(numnodes,numnodes);
%% 直流潮流的导纳矩阵计算
for k=1:branch_num 
    i=netpara(k,2);%首节点
    j=netpara(k,3);%尾节点
    Y(i,j)=-netpara(k,4);%导纳矩阵中非对角元素
    Y(j,i)= Y(i,j);
end
for k=1:numnodes
       Y(k,k)=-sum(Y(k,:)); %导纳矩阵中的对角元素 
end
%再删除掉平衡节点所在的行与列
Y(slack_bus,:)=[];
Y(:,slack_bus)=[];
%% 求解
ops=sdpsettings('solver', 'cplex');
result=solvesdp(st,totalcost);
double(totalcost) 
figure(1)
bar(value(p)','stack')%阶梯图
legend('Unit 1','Unit 2','Unit 3','Unit 4','Unit 5','Unit 6');  %在坐标轴上添加图例
figure(2)
stairs(value(p)')
legend('Unit 1','Unit 2','Unit 3','Unit 4','Unit 5','Unit 6');  %在坐标轴上添加图例
xlswrite('机组组合问题求解结果',double(u),'机组各时段启停计划');
P=(sum(sum_PowerGSDF(:,:,:),3)-sum_nodeGSDF(:,:))';%各段支路的实时潮流
P_sp=zeros(numnodes,T);%各个节点的直流潮流功率
for i=1:numnodes
for k=1:branch_num
    m=netpara(k,2);%首端节点
    n=netpara(k,3);%末端节点
    if m==i
        P_sp(i,:)=P_sp(i,:)+P(k,:);
    end
    if n==i
        P_sp(i,:)=P_sp(i,:)-P(k,:);
    end
end
end
dot_theta=zeros(numnodes,T);
dot_theta=X*P_sp;       
xlswrite('机组组合问题求解结果',double(P),'支路各时段的直流潮流');
xlswrite('机组组合问题求解结果',double(P_sp),'节点各时段的潮流功率');
xlswrite('机组组合问题求解结果',double(dot_theta),'节点各时段的潮流相角');


set(0,'ShowHiddenHandles','On')
set(gcf,'menubar','figure')

四、结果展示