【最优潮流】基于分布式交变方向乘法器(ADMM)方法来求解带碳排放交易的直流动态最优潮流（Matlab代码实现）
💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。
⛳️座右铭：行百里者，半于九十。
📋📋📋本文目录如下：🎁🎁🎁
目录
💥1 概述
📚2 运行结果
🎉3 文献来源
🌈4 Matlab代码实现
💥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 文献来源
部分理论来源于网络，如有侵权请联系删除。
