💥1 概述
文献来源：
摘要:本文提出了一种分布式交变方向乘法器(ADMM)方法来求解带碳排放交易的直流动态最优潮流问题(DC-DOPF-CET)。通常，基于adm的分布式方法公开了相邻子系统之间的边界总线和分支信息。与这些方法不同，本文提出的DC-ADMM-P方法采用了一种新颖的策略，即使用共识ADMM来解决DC-DOPF-CET的对偶问题，同时只公开相邻子系统之间的边界分支信息。此外，通过减少双乘法器数量和改进乘法器更新步骤，提高了DC-ADMM-P的收敛性能。DC-ADMM-P在6到1062总线的情况下进行了测试，并与其他分布式/分散方法进行了比较。仿真结果验证了DC-ADMM-P算法在求解具有复杂(非线性)因子的DC-DOPF问题时的有效性，这些复杂(非线性)因子可表示为凸可分函数。同时，通过减少双乘法器的数量和采用新的乘法器更新策略，提高了收敛性能。
📚2 运行结果
部分代码：
if isequal(includeCet,'yes') %包含碳排放约束 % -二次约束-start-% QCP.conQ = EconQ; QCP.conc = Econc; QCP.conb = Econb; %-二次约束-end-% end options = cplexoptimset; options.Display = 'off'; if isequal(includeCet,'yes') %包含碳排放约束 [x,fval,exitflag,output] = cplexqcp(QCP.Q,QCP.c,QCP.Aineq,QCP.bineq,QCP.Aeq,QCP.beq,QCP.conc,QCP.conQ,QCP.conb,QCP.lb,QCP.ub,[],options); disp(output.cplexstatusstring); else %不考虑CET [x,fval,exitflag,output] = cplexqp(QCP.Q,QCP.c,QCP.Aineq,QCP.bineq,QCP.Aeq,QCP.beq,QCP.lb,QCP.ub,[],options); disp(output.cplexstatusstring); end if isequal(isRTS,'yes') %RTS数据存在一个节点上有多台机组的情况 xx = sparse(PbusUnitsNumber(end,1)-PbusUnitsNumber(1,1),T);%P,行按照allNodes顺序排列 else xx = sparse(N,T);%P end st = sparse(N,T);%θ pf = 0;%整个系统的排放 pf_t = zeros(T,1);%排放量（按照时段划分） fd_t = zeros(T,1);%发电费用（按照时段划分） thpit = zeros(T,1);%机组出力（按照时段划分） eb = x(1);%买入碳排放deta_E_b es = x(2);%卖出碳排放deta_E_s if isequal(isRTS,'yes') %RTS数据存在一个节点上有多台机组的情况 dr = []; %弹性负荷变量dr hr = []; %辅助变量Hr units_number = 1; for i = 1:size(allNodes,1) bus_sequence_index = find(ismember(SCUC_data.busUnits.bus_sequence,allNodes(i,1))==1); %allNodes(i,1)在SCUC_data.busUnits.bus_sequence上的索引 P_start_index = (PbusUnitsNumber(i,1) - PbusUnitsNumber(1,1) + i - 1); %allNodes(i,1)对应变量P前面的所有变量P和θ的总数量 Seta_start_index = (PbusUnitsNumber(i+1,1) - PbusUnitsNumber(1,1) + i - 1); %allNodes(i,1)对应变量θ前面的所有变量P和θ的总数量 if ~isempty(bus_sequence_index) %节点上有发电机 for j = 1:size(SCUC_data.busUnits.unitIndex{bus_sequence_index,1},1) xx(units_number,:) = x(P_start_index*T + (j-1)*T + 2 + 1:P_start_index*T + (j-1)*T + 2 + T);%P units_number = units_number + 1; end else %节点上没有发电机 xx(units_number,:) = x(P_start_index*T + 2 + 1:P_start_index*T + 2 + T);%P units_number = units_number + 1; end st(i,:) = x(Seta_start_index*T + 2 + 1:Seta_start_index*T + 2 + T);%θ end elseif isequal(includeDR,'yes') dr = zeros(K,T); %弹性负荷变量dr hr = zeros(piecewiseNumber,T,K); %辅助变量Hr,第一个参数对应分段数，第二个参数对应时段，第三个参数对应节点编号 %按照片区顺序 for i = 1:n Pindex = 2*(PINumber{i}-1)*T + (EINumber{i}-1)*(piecewiseNumber+1)*T + 2; %+2为考虑碳排放的两个变量 Dindex = 2*(PINumber{i+1}-1)*T + (EINumber{i}-1)*(piecewiseNumber+1)*T + 2; %+2为考虑碳排放的两个变量 Hindex = Dindex + T; %取P和θ for j = 1:PINumber{i+1}-PINumber{i} %xx的行按照partitionData.allNodes(即allNodes)中节点编号的顺序 xx(PINumber{i}-1+j,:) = x(2*(j-1)*T+1+Pindex:2*(j-1)*T+T+Pindex,1);%P st(PINumber{i}-1+j,:) = x(2*(j-1)*T+T+1+Pindex:2*(j-1)*T+2*T+Pindex,1);%θ end %取dr和hr for k = 1:EINumber{i+1}-EINumber{i} %dr的行按照partitionData.allElasticityNodes中节点编号的顺序 dr(EINumber{i}-1+k,:) = x((k-1)*(piecewiseNumber+1)*T+1+Dindex:(k-1)*(piecewiseNumber+1)*T+T+Dindex,1); %dr for r = 1:piecewiseNumber hr(r,:,EINumber{i}-1+k) = x((k-1)*(piecewiseNumber+1)*T+(r-1)*T+1+Hindex:(k-1)*(piecewiseNumber+1)*T+(r-1)*T+T+Hindex,1); %hr end end end else dr = []; %弹性负荷变量dr hr = []; %辅助变量Hr for i = 1:N xx(i,:) = x((i-1)*2*T+2+1:(i-1)*2*T+2+T);%P st(i,:) = x((i-1)*2*T+2+T+1:(i-1)*2*T+2+2*T);%θ end end if isequal(isRTS,'yes') %RTS数据存在一个节点上有多台机组的情况 %发电费用（按照时段划分） for t = 1:T units_number = 1; for i = 1:size(allNodes,1) bus_sequence_index = find(ismember(SCUC_data.busUnits.bus_sequence,allNodes(i,1))==1); %allNodes(i,1)在SCUC_data.busUnits.bus_sequence上的索引 if ~isempty(bus_sequence_index) %节点上有发电机 for j = 1:size(SCUC_data.busUnits.unitIndex{bus_sequence_index,1},1) units_index = SCUC_data.busUnits.unitIndex{bus_sequence_index,1}(j,1); %SCUC_data.busUnits.unitIndex{bus_sequence_index,1}在SCUC_data.units.bus_G上的索引 fd_t(t,1) = fd_t(t,1) + xx(units_number,t).^2 * SCUC_data.units.gamma(units_index) ... + xx(units_number,t) * SCUC_data.units.beta(units_index) + SCUC_data.units.alpha(units_index); units_number = units_number + 1; end else units_number = units_number + 1; end end end else if isequal(includeCet,'yes') %包含碳排放约束 %整个系统的排放量（所有时段） for i = 1:unitN index = find(allNodes == SCUC_data.units.bus_G(i)); %发电机节点的在矩阵中的索引 for t = 1:T pf = pf + xx(index,t).^2 * SCUC_data.units.c(i) + xx(index,t) * SCUC_data.units.b(i) + SCUC_data.units.a(i); end end %排放量（按照时段划分） for t = 1:T for i = 1:unitN index = find(allNodes == SCUC_data.units.bus_G(i)); %发电机节点的在矩阵中的索引 pf_t(t,1) = pf_t(t,1) + xx(index,t).^2 * SCUC_data.units.c(i) + xx(index,t) * SCUC_data.units.b(i) + SCUC_data.units.a(i); end end end %发电费用（按照时段划分） for t = 1:T for i = 1:unitN index = find(allNodes == SCUC_data.units.bus_G(i)); %发电机节点的在矩阵中的索引 fd_t(t,1) = fd_t(t,1) + xx(index,t).^2 * SCUC_data.units.gamma(i) + xx(index,t) * SCUC_data.units.beta(i) + SCUC_data.units.alpha(i); end end end %所有时段的机组出力计划,比较发现，使用需求响应，确实能够削峰填谷 thpit = full(sum(xx)); disp('总费用='); disp((fval+QCP.b)*SCUC_data.baseparameters.standardP); %返回值赋值--start-% IEEE6:8.977148404647500e+04 ; IEEE30:2.136876446422936e+03 ;IEEE118:1.729592022969799e+06 fval = (fval+QCP.b)*SCUC_data.baseparameters.standardP; %返回值赋值--end-% end
🎉3 文献来源
部分理论来源于网络，如有侵权请联系删除。
