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

clearclc%%系统参数%所有参数均用有名值表示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:gennumfor t=1:Tfor 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);%加上机组启停产生的开停机成本endend%原二次函数式的最优成本目标% 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:Tst=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];%上坡    endend%% 热备用约束hp=0.05;%热备用系数for t=1:Tst=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];    endend%停机约束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];    endend%% 启停成本约束for t=1:T   %启停成本零限约束    for i=1:gennum        st=st+[costH(i,t)>=0];         st=st+[costJ(i,t)>=0];    endendfor 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);endfor 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:numnodesfor 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,:);    endendenddot_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')