请优化以下代码,在原有代码上不改变原有逻辑上进行修改:clc
clear
close all
%% 参数初始化
iter=0; % 迭代次数初始化
maxIter=100; % 设定最大迭代次数
tolerant=0.001; % 设定收敛阈值
X=zeros(maxIter+1,48); % 历史交易价格记录矩阵
Power_buy=zeros(maxIter,24); % 下层微网的历史交易功率保存
Power_sell=zeros(maxIter,24);
% 共享储能从主网购电价格
u_Db=[0.4,0.4,0.4,0.4,0.4,0.4,0.79,0.79,0.79,1.2,1.2,1.2,1.2,1.2,0.79,0.79,0.79,1.2,1.2,1.2,0.79,0.79,0.4,0.4];
% 共享储能向主网售电价格
u_Ds=[0.35,0.35,0.35,0.35,0.35,0.35,0.68,0.68,0.68,1.12,1.12,1.12,1.12,1.12,0.68,0.68,0.68,1.12,1.12,1.12,0.79,0.79,0.35,0.35];
% 共享储能与综合能源系统的交易价格上下限
u_Pbmax=1.5*u_Db; % 综合能源系统的购电价上限
u_Pbmin=u_Db; % 综合能源系统的购电价下限
u_Psmax=u_Ds; % 综合能源系统的售电价上限
u_Psmin=u_Psmax-0.35*ones(1,24); % 综合能源系统的售电价下限
% 价格上限
ub=[u_Pbmax,u_Psmax];
% 价格下限
lb=[u_Pbmin,u_Psmin];
% 初始价格赋值
X(1,:)=(ub+lb)/1.2;
% 二分法添加约束的历史上、下限保存
X_up=zeros(maxIter+1,48);
X_low=zeros(maxIter+1,48);
for k=1:maxIter+1
X_up(k,:)=ub;
X_low(k,:)=lb;
end
U=0; % 是否增加约束的标志,如果U=1,则向Fun_DSO增加约束
%% 迭代
while 1
if iter==maxIter % 限制迭代次数
disp('迭代不收敛,参数有误');
break;
end
disp(['迭代还未收敛,当前迭代第 ', num2str(iter),' 次']);
% 直接完成前两次迭代,省略步骤5
if iter==0
% 第一次迭代
[Power_buy(1,1:24),Power_sell(1,1:24),Obj_IES(1,:),p_tran,obj_user]=Fun_LA(X(1,1:48));
[X(2,1:48),Obj_ESO(1),X_up,X_low]=Fun_ESO(U,Power_buy(1,1:24),Power_sell(1,1:24),zeros(1,48),zeros(1,48),X_up,X_low,iter);
% 第二次迭代
[Power_buy(2,1:24),Power_sell(2,1:24),Obj_IES(2,:),p_tran,obj_user]=Fun_LA(X(2,1:48));
[X(3,1:48),Obj_ESO(2),X_up,X_low]=Fun_ESO(U,Power_buy(2,1:24),Power_sell(2,1:24),zeros(1,48),zeros(1,48),X_up,X_low,iter);
iter=2; % 完成两次迭代+2
end
% 步骤4,判断收敛条件
if sum(abs(X(iter+1,:)-X(iter,:)))/sum(X(iter,:))<=tolerant
disp(['迭代收敛,在第 ', num2str(iter),' 次收敛']);
break;
else
iter=iter+1; % 迭代次数自增
end
% 步骤6
if isequal(X(iter,1:48),X(iter-2,1:48))
iter=iter+1; % 迭代次数自增,相比于原文这里自增提前了,所以后面判断部分的iter-1
X(iter,1:48)=(X(iter-1,1:48)+X(iter-2,1:48))/2; % 二分法赋值
Obj_ESO(iter-1)=Obj_ESO(iter-2);
Obj_IES(iter-1,:)=Obj_IES(iter-2,:);
disp(['迭代还未收敛,当前迭代第 ', num2str(iter-1),' 次']);
if iter-1==maxIter % 限制迭代次数
disp('迭代不收敛,参数有误');
break;
end
U=0; % 重置U,避免添加约束
else
U=1; % 是否增加约束的标志,如果U=1,则向Fun_ESO增加约束
end
% 此处开始真正的一般迭代步骤,先执行步骤2和步骤3,然后转到上方执行步骤4和步骤6
% 步骤2,优化购能策略
[Power_buy(iter,1:24),Power_sell(iter,1:24),Obj_IES(iter,:),p_tran,obj_user]=Fun_LA(X(iter,1:48));
% 步骤3,优化交易价格
[X(iter+1,1:48),Obj_ESO(iter),X_up,X_low]=Fun_ESO(U,Power_buy(iter,1:24),Power_sell(iter,1:24),X(iter,1:48),X(iter-1,1:48),X_up,X_low,iter);
end
%% 画图输出运行结果
% 分布式优化迭代收敛过程
figure
yyaxis left
plot(Obj_IES(1:iter-1),'r->','LineWidth',1.5);
ylabel('综合能源系统收益/元');
xlabel('迭代次数');
yyaxis right
plot(Obj_ESO(1:iter-1),'c-d','LineWidth',1.5);
ylabel('共享储能收益/元');
title('分布式优化迭代收敛过程');
hold on
legend('综合能源系统','共享储能');
legend('boxoff');
grid on
box off
% 共享储能制定的交易电价
figure
plot(X(iter,1:24),'b-*','LineWidth',1.5);
hold on
plot(X(iter,25:48),'r-s','LineWidth',1.5);
xlabel('时段');
ylabel('交易电价/元');
legend('购电价','售电价');
title('共享储能给综合能源系统制定的电价');
box off
function [Power_buy,Power_sell,Obj_IES,p_tran,obj_user]=Fun_LA(lambda)
% 综合能源系统和用户聚合商分布式优化模型
%% 导入基础数据
% 综合能源系统向共享储能的购售电价(上层给下层的购售电价)
Pri_buy=lambda(1:24); % 综合能源系统从共享储能的购电价格
Pri_sell=lambda(25:48); % 综合能源系统向共享储能的售电价格
% 电网购/售电价格
pri_buy=[0.40*ones(1,7),0.75*ones(1,4),1.20*ones(1,3),0.75*ones(1,4),1.20*ones(1,4),0.40*ones(1,2)];
pri_sell=[0.20*ones(1,7),0.40*ones(1,4),0.60*ones(1,3),0.40*ones(1,4),0.60*ones(1,4),0.20*ones(1,2)];
% 综合能源系统可再生功率
Pre_e_rew=[875.530000000000,846.018000000000,862.711000000000,870.271000000000,836.205000000000,841.486000000000,782.689000000000,698.204000000000,701.427000000000,621.098000000000,581.150000000000,596.790000000000,580.380000000000,594.440000000000,566.580000000000,574.680000000000,585.510000000000,552.640000000000,670.590000000000,783.630000000000,804.215000000000,850.555000000000,884.575000000000,897.079000000000];
% 机会约束置信水平0.8,处理新能源不确定
margin =norminv(0.8);
% 3个用户的屋顶光伏功率、电负荷功率、热负荷功率
p_PV_1=[0,0,0,0,0,0,129.249800400000,214.269080100000,406.062737200000,187.890089100000,348.089547200000,489.493149800000,715.739240300000,535.311515700000,419.338819500000,106.062817800000,110.394000800000,0,0,0,0,0,0,0];
p_PV_2=[0,0,0,0,0,0,120.055543400000,220.397293500000,415.706359500000,418.698799700000,491.053826400000,481.272437700000,625.935805400000,626.515006000000,242.663224400000,139.605285200000,100.423479600000,0,0,0,0,0,0,0];
p_PV_3=[0,0,0,0,0,0,133.403152700000,202.166004400000,330.858423400000,322.621504700000,563.029096300000,730.455581800000,500.117797900000,517.608432000000,531.138746500000,201.970629000000,99.6730838300000,0,0,0,0,0,0,0];
p_L_1=[438.916916300000,449.056859100000,384.870009500000,383.558344100000,555.779696000000,829.075341000000,695.912178100000,634.751870600000,824.460547100000,534.335901900000,811.547316800000,886.409315700000,730.531820100000,887.845537000000,541.619819400000,696.932235100000,769.473943000000,699.578910600000,572.864157800000,540.616560300000,435.745547100000,396.302703800000,474.343413500000,437.774595000000];
p_L_2=[311.348637500000,496.957953500000,422.784317900000,490.746788800000,501.394977100000,743.940698300000,606.486646600000,663.093908500000,700.896210900000,848.333544800000,800.266967400000,776.039508700000,651.840857300000,983.690799900000,671.581982400000,618.503165900000,512.865450300000,509.703972100000,756.912181800000,687.646096900000,554.832135700000,461.364115700000,344.505043300000,451.308587500000];
p_L_3=[392.652248000000,394.794708200000,398.316433800000,350.809513300000,375.836497300000,690.772657300000,627.585185200000,806.969387600000,545.598216500000,717.080418000000,874.862341300000,782.362704200000,809.257738000000,819.566430800000,549.259711900000,802.076021500000,831.482225700000,612.444478600000,597.203128300000,564.869066100000,544.120910100000,579.738373700000,481.816023100000,412.592805000000];
h_L_1=[281.524850700000,210.022691300000,303.033573100000,218.448903100000,263.472930800000,350.300699900000,375.959274400000,330.863083900000,337.319821000000,444.792387000000,436.016969100000,503.658902500000,312.272982400000,366.629675500000,462.303817700000,500.065901100000,438.912277000000,331.586449400000,454.634115400000,356.664546200000,321.477668300000,363.056798100000,246.404780100000,320.909875400000];
h_L_2=[292.692635300000,279.790642700000,249.415523200000,335.280209500000,324.713923800000,311.205789900000,440.457761000000,505.704909100000,346.729019700000,328.935902000000,360.900341100000,296.288544300000,434.498292300000,444.006202200000,359.390592300000,373.673058500000,554.704060400000,503.057494500000,296.993898500000,336.521125600000,328.445557900000,356.387830000000,328.938629800000,293.093364400000];
h_L_3=[260.079123200000,247.818540100000,262.922674500000,282.560914900000,402.309088300000,310.230575500000,394.913254600000,441.181108700000,295.731606800000,348.935949700000,413.055395700000,447.172973900000,439.955735100000,355.544164000000,496.099300600000,402.265911800000,402.001222900000,557.987919600000,482.691546700000,437.399880500000,307.046253500000,243.881573900000,303.084737900000,306.684971700000];
% Big-M法中的大M
M=1E8;
%% 决策变量
% 综合能源系统决策变量
P_e_rew=sdpvar(1,24); % 综合能源系统中的可再生的实际出力
P_e_cur=sdpvar(1,24); % 综合能源系统中的弃可再生的量
P_e_CHP=sdpvar(1,24); % CHP的产电功率
P_e_P2G=sdpvar(1,24); % CHP的供给P2G的功率
P_e_CCS=sdpvar(1,24); % CHP的供给CCS的功率
P_eJ_CHP=sdpvar(1,24); % CHP的净出力
P_h_CHP=sdpvar(1,24); % CHP的产热功率
CO2_CHP=sdpvar(1,24); % CHP的碳排放量
Gas_CHP=sdpvar(1,24); % CHP的耗气量
CO2_CCS=sdpvar(1,24); % CCS的碳捕集量/P2G所用的二氧化碳量
Gas_P2G=sdpvar(1,24); % P2G的产气量
P_h_GB=sdpvar(1,24); % GB的产热功率
Gas_GB=sdpvar(1,24); % GB的耗气量
P_buy=sdpvar(1,24); % 综合能源系统向共享储能的购电量
P_sell=sdpvar(1,24); % 综合能源系统向共享储能的售电量
miu_ies=binvar(1,24); % 综合能源系统与共享储能交互的状态变量
E_CO2=sdpvar(1,1); % 综合能源系统的实际碳排放量
E_c=sdpvar(1,1); % 综合能源系统的碳排放配额
E_1=sdpvar(1,1); % 阶梯碳交易引入的辅助区间变量
E_2=sdpvar(1,1);
E_3=sdpvar(1,1);
pri_e_1=sdpvar(1,24); % 综合能源系统给用户1制定的电价
pri_h_1=sdpvar(1,24); % 综合能源系统给用户1制定的热价
pri_e_2=sdpvar(1,24); % 综合能源系统给用户2制定的电价
pri_h_2=sdpvar(1,24); % 综合能源系统给用户2制定的热价
pri_e_3=sdpvar(1,24); % 综合能源系统给用户3制定的电价
pri_h_3=sdpvar(1,24); % 综合能源系统给用户3制定的热价
% 用户决策变量
p_cut_1=sdpvar(1,24); % 用户1的可削减电负荷
h_cut_1=sdpvar(1,24); % 用户1的可削减热负荷
p_load_1=sdpvar(1,24); % 用户1的购电量
h_load_1=sdpvar(1,24); % 用户1的购热量
p_cut_2=sdpvar(1,24); % 用户2的可削减电负荷
h_cut_2=sdpvar(1,24); % 用户2的可削减热负荷
p_load_2=sdpvar(1,24); % 用户2的购电量
h_load_2=sdpvar(1,24); % 用户2的购热量
p_cut_3=sdpvar(1,24); % 用户3的可削减电负荷
h_cut_3=sdpvar(1,24); % 用户3的可削减热负荷
p_load_3=sdpvar(1,24); % 用户3的购电量
h_load_3=sdpvar(1,24); % 用户3的购热量
p_tran_12=zeros(1,24); % 用户1与用户2的交互电量
p_tran_13=zeros(1,24); % 用户1与用户3的交互电量
p_tran_23=zeros(1,24); % 用户2与用户3的交互电量
% KKT条件部分引入的变量
lambda_e_1=sdpvar(1,24);
lambda_h_1=sdpvar(1,24);
miu_11_lb=sdpvar(1,24);
miu_11_ub=sdpvar(1,24);
miu_12_lb=sdpvar(1,24);
miu_12_ub=sdpvar(1,24);
miu_13_lb=sdpvar(1,24);
miu_13_ub=sdpvar(1,24);
miu_14_lb=sdpvar(1,24);
miu_14_ub=sdpvar(1,24);
lambda_e_2=sdpvar(1,24);
lambda_h_2=sdpvar(1,24);
miu_21_lb=sdpvar(1,24);
miu_21_ub=sdpvar(1,24);
miu_22_lb=sdpvar(1,24);
miu_22_ub=sdpvar(1,24);
miu_23_lb=sdpvar(1,24);
miu_23_ub=sdpvar(1,24);
miu_24_lb=sdpvar(1,24);
miu_24_ub=sdpvar(1,24);
lambda_e_3=sdpvar(1,24);
lambda_h_3=sdpvar(1,24);
miu_31_lb=sdpvar(1,24);
miu_31_ub=sdpvar(1,24);
miu_32_lb=sdpvar(1,24);
miu_32_ub=sdpvar(1,24);
miu_33_lb=sdpvar(1,24);
miu_33_ub=sdpvar(1,24);
miu_34_lb=sdpvar(1,24);
miu_34_ub=sdpvar(1,24);
miu_1_lb=sdpvar(1,24);
miu_1_ub=sdpvar(1,24);
miu_2_lb=sdpvar(1,24);
miu_2_ub=sdpvar(1,24);
miu_3_lb=sdpvar(1,24);
miu_3_ub=sdpvar(1,24);
b_11_lb=binvar(1,24);
b_11_ub=binvar(1,24);
b_12_lb=binvar(1,24);
b_12_ub=binvar(1,24);
b_13_lb=binvar(1,24);
b_13_ub=binvar(1,24);
b_14_lb=binvar(1,24);
b_14_ub=binvar(1,24);
b_21_lb=binvar(1,24);
b_21_ub=binvar(1,24);
b_22_lb=binvar(1,24);
b_22_ub=binvar(1,24);
b_23_lb=binvar(1,24);
b_23_ub=binvar(1,24);
b_24_lb=binvar(1,24);
b_24_ub=binvar(1,24);
b_31_lb=binvar(1,24);
b_31_ub=binvar(1,24);
b_32_lb=binvar(1,24);
b_32_ub=binvar(1,24);
b_33_lb=binvar(1,24);
b_33_ub=binvar(1,24);
b_34_lb=binvar(1,24);
b_34_ub=binvar(1,24);
b_1_lb=binvar(1,24);
b_1_ub=binvar(1,24);
b_2_lb=binvar(1,24);
b_2_ub=binvar(1,24);
b_3_lb=binvar(1,24);
b_3_ub=binvar(1,24);
% 针对双线性项的McCormick线性化引入的变量
chi_e_1=sdpvar(1,24); % 用来等效购电成本的辅助变量
chi_h_1=sdpvar(1,24); % 用来等效购热成本的辅助变量
chi_e_2=sdpvar(1,24); % 用来等效购电成本的辅助变量
chi_h_2=sdpvar(1,24); % 用来等效购热成本的辅助变量
chi_e_3=sdpvar(1,24); % 用来等效购电成本的辅助变量
chi_h_3=sdpvar(1,24); % 用来等效购热成本的辅助变量
%% 导入约束条件
C=[];
% 综合能源系统约束部分
C=[C,
% 综合能源系统给用户的购售电/热价定价约束部分
pri_sell<=pri_e_1<=pri_buy, % 综合能源系统给用户1的售电价格上下限约束
sum(pri_e_1)/24<=0.6, % 综合能源系统给用户1的售电价格平均值约束
0.15<=pri_h_1<=0.6, % 综合能源系统给用户1的售热价格上下限约束
sum(pri_h_1)/24<=0.5, % 综合能源系统给用户1的售热价格平均值约束
pri_sell<=pri_e_2<=pri_buy, % 综合能源系统给用户2的售电价格上下限约束
sum(pri_e_2)/24<=0.6, % 综合能源系统给用户2的售电价格平均值约束
0.15<=pri_h_2<=0.6, % 综合能源系统给用户2的售热价格上下限约束
sum(pri_h_2)/24<=0.5, % 综合能源系统给用户2的售热价格平均值约束
pri_sell<=pri_e_3<=pri_buy, % 综合能源系统给用户3的售电价格上下限约束
sum(pri_e_3)/24<=0.6, % 综合能源系统给用户3的售电价格平均值约束
0.15<=pri_h_3<=0.6, % 综合能源系统给用户3的售热价格上下限约束
sum(pri_h_3)/24<=0.5, % 综合能源系统给用户3的售热价格平均值约束
% 综合能源系统的机组运行约束部分
P_e_cur+P_e_rew==Pre_e_rew, % 可再生出力平衡约束
0<=P_e_cur<=Pre_e_rew, % 可再生实际出力上下限约束
0<=P_e_rew<=Pre_e_rew, % 弃可再生量上下限约束
P_e_CHP==0.3*9.7*Gas_CHP, % CHP的发电和耗气功率的关系约束
P_h_CHP==0.45*9.7*Gas_CHP, % CHP的发热和耗气功率的关系约束
0<=P_e_CHP<=2000, % CHP的发电功率上下限约束
P_h_GB==0.9*9.7*Gas_GB, % GB的发热和耗气功率的关系约束
0<=P_h_GB<=2000, % GB的发热功率上下限约束
% 碳捕集和电转气运行约束部分
CO2_CHP==0.390*9.7*Gas_CHP, % CHP的碳排放强度计算式约束
0<=CO2_CCS<=0.9*CO2_CHP, % CHP中的碳捕集设备的捕碳量上下限约束
P_e_CCS==0.109*CO2_CCS, % CHP中的碳捕集设备的耗电量与捕碳量的关系约束
0<=P_e_CCS<=500, % 碳捕集设备的运行上下限约束
CO2_CCS==0.8*P_e_P2G, % 生产甲烷所需要的二氧化碳的质量与P2G消耗的电功率的关系约束
Gas_P2G==0.95/3.9*P_e_P2G, % P2G生成CH4的体积与消耗的电功率的等式约束
0<=P_e_P2G<=500, % 电转气设备的运行上下限约束
P_eJ_CHP==P_e_CHP-P_e_CCS-P_e_P2G, % CHP机组的净出力关系式约束
P_eJ_CHP>=0, % CHP机组的净出力非负性约束
% 综合能源系统与共享储能交互约束部分
0<=P_buy<=1000*miu_ies, % 综合能源系统向共享储能的购电量上下限约束
0<=P_sell<=1000*(1-miu_ies), % 综合能源系统向共享储能的售电量上下限约束
% 综合能源系统功率平衡约束部分
(P_buy-P_sell)+P_eJ_CHP+P_e_rew-p_load_1-p_load_2-p_load_3>=margin*sqrt((0.2*Pre_e_rew).^2), % 电功率平衡约束
P_h_CHP+P_h_GB==h_load_1+h_load_2+h_load_3, % 热功率平衡约束
% 阶梯碳交易约束部分
E_CO2==sum(CO2_CHP-CO2_CCS)+0.968*sum(Gas_GB), % 综合能源系统的实时碳排放量
E_c==0.047*sum(1.67*P_e_CHP+P_h_CHP+P_h_GB), % 综合能源系统的碳排放配额
E_CO2-E_c==E_1+E_2+E_3, % 阶梯碳交易等效区间赋值
0<=E_1<=5000, % 分段区间上下限约束
0<=E_2<=5000,
0<=E_3,
];
% 用户约束部分
% 用户的原始约束部分
C=[C,
% 用户1的约束部分
0<=p_cut_1<=0.15*p_L_1, % 可削减电负荷的上下限约束,miu_11_lb、miu_11_ub
p_load_1==p_L_1-p_PV_1-p_cut_1+p_tran_12+p_tran_13, % 负荷侧购电量功率平衡约束,lambda_e_1
0<=h_cut_1<=0.15*h_L_1, % 可削减热负荷的上下限约束,miu_12_lb、miu_12_ub
h_load_1==h_L_1-h_cut_1, % 负荷侧购热量功率平衡约束,lambda_h_1
0<=p_load_1<=1000, % 购电量上下限约束,miu_13_lb、miu_13_ub
0<=h_load_1<=1000, % 购热量上下限约束,miu_14_lb、miu_14_ub
% 用户2的约束部分
0<=p_cut_2<=0.15*p_L_2, % 可削减电负荷的上下限约束,miu_21_lb、miu_21_ub
p_load_2==p_L_2-p_PV_2-p_cut_2-p_tran_12+p_tran_23, % 负荷侧购电量功率平衡约束,lambda_e_2
0<=h_cut_2<=0.15*h_L_2, % 可削减热负荷的上下限约束,miu_22_lb、miu_22_ub
h_load_2==h_L_2-h_cut_2, % 负荷侧购热量功率平衡约束,lambda_h_2
0<=p_load_2<=1000, % 购电量上下限约束,miu_23_lb、miu_23_ub
0<=h_load_2<=1000, % 购热量上下限约束,miu_24_lb、miu_24_ub
% 用户3的约束部分
0<=p_cut_3<=0.15*p_L_3, % 可削减电负荷的上下限约束,miu_31_lb、miu_31_ub
p_load_3==p_L_3-p_PV_3-p_cut_3-p_tran_13-p_tran_23, % 负荷侧购电量功率平衡约束,lambda_e_3
0<=h_cut_3<=0.15*h_L_3, % 可削减热负荷的上下限约束,miu_32_lb、miu_32_ub
h_load_3==h_L_1-h_cut_3, % 负荷侧购热量功率平衡约束,lambda_h_3
0<=p_load_3<=1000, % 购电量上下限约束,miu_33_lb、miu_33_ub
0<=h_load_3<=1000, % 购热量上下限约束,miu_34_lb、miu_34_ub
% 各个用户之间的电能交易约束部分
-200<=p_tran_12<=200, % 用户1与用户2之间的电能交互上下限约束,miu_1_lb、miu_1_ub
-200<=p_tran_13<=200, % 用户1与用户3之间的电能交互上下限约束,miu_2_lb、miu_2_ub
-200<=p_tran_23<=200, % 用户2与用户3之间的电能交互上下限约束,miu_3_lb、miu_3_ub
];
% 用户的KKT等效约束部分
C=[C,
% KKT平衡条件约束部分
0.5-lambda_e_1-miu_11_lb+miu_11_ub==0, % 拉格朗日函数关于p_cut_1的偏导部分
0.5-lambda_e_2-miu_21_lb+miu_21_ub==0, % 拉格朗日函数关于p_cut_2的偏导部分
0.5-lambda_e_3-miu_31_lb+miu_31_ub==0, % 拉格朗日函数关于p_cut_3的偏导部分
% lambda_e_1-lambda_e_2-miu_1_lb+miu_1_ub==0, % 拉格朗日函数关于p_tran_12的偏导部分
% lambda_e_1-lambda_e_3-miu_2_lb+miu_2_ub==0, % 拉格朗日函数关于p_tran_13的偏导部分
% lambda_e_2-lambda_e_3-miu_3_lb+miu_3_ub==0, % 拉格朗日函数关于p_tran_23的偏导部分
0.4-lambda_h_1-miu_12_lb+miu_12_ub==0, % 拉格朗日函数关于h_cut_1的偏导部分
0.4-lambda_h_2-miu_22_lb+miu_22_ub==0, % 拉格朗日函数关于h_cut_2的偏导部分
0.4-lambda_h_3-miu_32_lb+miu_32_ub==0, % 拉格朗日函数关于h_cut_3的偏导部分
pri_e_1-lambda_e_1-miu_13_lb+miu_13_ub==0, % 拉格朗日函数关于p_load_1的偏导部分
pri_e_2-lambda_e_2-miu_23_lb+miu_23_ub==0, % 拉格朗日函数关于p_load_2的偏导部分
pri_e_3-lambda_e_3-miu_33_lb+miu_33_ub==0, % 拉格朗日函数关于p_load_3的偏导部分
pri_h_1-lambda_h_1-miu_14_lb+miu_14_ub==0, % 拉格朗日函数关于h_load_1的偏导部分
pri_h_2-lambda_h_2-miu_24_lb+miu_24_ub==0, % 拉格朗日函数关于h_load_2的偏导部分
pri_h_3-lambda_h_3-miu_34_lb+miu_34_ub==0, % 拉格朗日函数关于h_load_3的偏导部分
% 经过Big-M法线性化后的KKT方向和互补条件约束部分
0<=p_cut_1<=M*b_11_lb,
0<=miu_11_lb<=M*(1-b_11_lb),
0<=0.15*p_L_1-p_cut_1<=M*b_11_ub,
0<=miu_11_ub<=M*(1-b_11_ub),
0<=p_cut_2<=M*b_21_lb,
0<=miu_21_lb<=M*(1-b_21_lb),
0<=0.15*p_L_2-p_cut_2<=M*b_21_ub,
0<=miu_21_ub<=M*(1-b_21_ub),
0<=p_cut_3<=M*b_31_lb,
0<=miu_31_lb<=M*(1-b_31_lb),
0<=0.15*p_L_3-p_cut_3<=M*b_31_ub,
0<=miu_31_ub<=M*(1-b_31_ub),
% 0<=p_tran_12+200<=M*b_1_lb,
% 0<=miu_1_lb<=M*(1-b_1_lb),
% 0<=200-p_tran_12<=M*b_1_ub,
% 0<=miu_1_ub<=M*(1-b_1_ub),
% 0<=p_tran_13+200<=M*b_2_lb,
% 0<=miu_2_lb<=M*(1-b_2_lb),
% 0<=200-p_tran_13<=M*b_2_ub,
% 0<=miu_2_ub<=M*(1-b_2_ub),
% 0<=p_tran_23+200<=M*b_3_lb,
% 0<=miu_3_lb<=M*(1-b_3_lb),
% 0<=200-p_tran_23<=M*b_3_ub,
% 0<=miu_3_ub<=M*(1-b_3_ub),
0<=h_cut_1<=M*b_12_lb,
0<=miu_12_lb<=M*(1-b_12_lb),
0<=0.15*h_L_1-h_cut_1<=M*b_12_ub,
0<=miu_12_ub<=M*(1-b_12_ub),
0<=h_cut_2<=M*b_22_lb,
0<=miu_22_lb<=M*(1-b_22_lb),
0<=0.15*h_L_2-h_cut_2<=M*b_22_ub,
0<=miu_22_ub<=M*(1-b_22_ub),
0<=h_cut_3<=M*b_32_lb,
0<=miu_32_lb<=M*(1-b_32_lb),
0<=0.15*h_L_3-h_cut_3<=M*b_32_ub,
0<=miu_32_ub<=M*(1-b_32_ub),
0<=p_load_1<=M*b_13_lb,
0<=miu_13_lb<=M*(1-b_13_lb),
0<=1000-p_load_1<=M*b_13_ub,
0<=miu_13_ub<=M*(1-b_13_ub),
0<=p_load_2<=M*b_23_lb,
0<=miu_23_lb<=M*(1-b_23_lb),
0<=1000-p_load_2<=M*b_23_ub,
0<=miu_23_ub<=M*(1-b_23_ub),
0<=p_load_3<=M*b_33_lb,
0<=miu_33_lb<=M*(1-b_33_lb),
0<=1000-p_load_3<=M*b_33_ub,
0<=miu_33_ub<=M*(1-b_33_ub),
0<=h_load_1<=M*b_14_lb,
0<=miu_14_lb<=M*(1-b_14_lb),
0<=1000-h_load_1<=M*b_14_ub,
0<=miu_14_ub<=M*(1-b_14_ub),
0<=h_load_2<=M*b_24_lb,
0<=miu_24_lb<=M*(1-b_24_lb),
0<=1000-h_load_2<=M*b_24_ub,
0<=miu_24_ub<=M*(1-b_24_ub),
0<=h_load_3<=M*b_34_lb,
0<=miu_34_lb<=M*(1-b_34_lb),
0<=1000-h_load_3<=M*b_34_ub,
0<=miu_34_ub<=M*(1-b_34_ub),
];
% 针对双线性项的McCormick线性化约束部分
for t=1:24
C=[C,
% 用户1的双线性项线性化
% 电部分
chi_e_1(t)>=pri_sell(t)*p_load_1(t)+pri_e_1(t)*0-pri_sell(t)*0,
chi_e_1(t)>=pri_buy(t)*p_load_1(t)+pri_e_1(t)*1000-pri_buy(t)*1000,
chi_e_1(t)<=pri_buy(t)*p_load_1(t)+pri_e_1(t)*0-pri_buy(t)*0,
chi_e_1(t)<=pri_e_1(t)*1000+pri_sell(t)*p_load_1(t)-pri_sell(t)*1000,
% 热部分
chi_h_1(t)>=0.15*h_load_1(t)+pri_h_1(t)*0-0.15*0,
chi_h_1(t)>=0.6*h_load_1(t)+pri_h_1(t)*1000-0.6*1000,
chi_h_1(t)<=0.6*h_load_1(t)+pri_h_1(t)*0-0.6*0,
chi_h_1(t)<=pri_h_1(t)*1000+0.15*h_load_1(t)-0.15*1000,
% 用户2的双线性项线性化
% 电部分
chi_e_2(t)>=pri_sell(t)*p_load_2(t)+pri_e_2(t)*0-pri_sell(t)*0,
chi_e_2(t)>=pri_buy(t)*p_load_2(t)+pri_e_2(t)*1000-pri_buy(t)*1000,
chi_e_2(t)<=pri_buy(t)*p_load_2(t)+pri_e_2(t)*0-pri_buy(t)*0,
chi_e_2(t)<=pri_e_2(t)*1000+pri_sell(t)*p_load_2(t)-pri_sell(t)*1000,
% 热部分
chi_h_2(t)>=0.15*h_load_2(t)+pri_h_2(t)*0-0.15*0,
chi_h_2(t)>=0.6*h_load_2(t)+pri_h_2(t)*1000-0.6*1000,
chi_h_2(t)<=0.6*h_load_2(t)+pri_h_2(t)*0-0.6*0,
chi_h_2(t)<=pri_h_2(t)*1000+0.15*h_load_2(t)-0.15*1000,
% 用户3的双线性项线性化
% 电部分
chi_e_3(t)>=pri_sell(t)*p_load_3(t)+pri_e_3(t)*0-pri_sell(t)*0,
chi_e_3(t)>=pri_buy(t)*p_load_3(t)+pri_e_3(t)*1000-pri_buy(t)*1000,
chi_e_3(t)<=pri_buy(t)*p_load_3(t)+pri_e_3(t)*0-pri_buy(t)*0,
chi_e_3(t)<=pri_e_3(t)*1000+pri_sell(t)*p_load_3(t)-pri_sell(t)*1000,
% 热部分
chi_h_3(t)>=0.15*h_load_3(t)+pri_h_3(t)*0-0.15*0,
chi_h_3(t)>=0.6*h_load_3(t)+pri_h_3(t)*1000-0.6*1000,
chi_h_3(t)<=0.6*h_load_3(t)+pri_h_3(t)*0-0.6*0,
chi_h_3(t)<=pri_h_3(t)*1000+0.15*h_load_3(t)-0.15*1000,
];
end
%% 目标函数
Obj=-(3.2*sum(Gas_CHP+Gas_GB-Gas_P2G)+sum(Pri_buy.*P_buy-Pri_sell.*P_sell)+0.5*sum(P_e_cur)+0.2*E_1+0.3*E_2+0.45*E_3...
-sum(chi_e_1+chi_h_1)-sum(chi_e_2+chi_h_2)-sum(chi_e_3+chi_h_3));
%% 求解器配置与求解
ops=sdpsettings('solver','cplex','verbose',0,'usex0',0);
ops.cplex.mip.tolerances.mipgap=1e-6;
result=solvesdp(C,-Obj,ops); % 目标函数为综合能源系统的收益函数,所以有'负号'
%% 导出数据
% 综合能源系统决策变量
P_e_rew=double(P_e_rew); % 综合能源系统中的可再生的实际出力
P_e_cur=double(P_e_cur); % 综合能源系统中的弃可再生的量
P_e_CHP=double(P_e_CHP); % CHP的产电功率
P_e_P2G=double(P_e_P2G); % CHP的供给P2G的功率
P_e_CCS=double(P_e_CCS); % CHP的供给CCS的功率
P_eJ_CHP=double(P_eJ_CHP); % CHP的净出力
P_h_CHP=double(P_h_CHP); % CHP的产热功率
CO2_CHP=double(CO2_CHP); % CHP的碳排放量
Gas_CHP=double(Gas_CHP); % CHP的耗气量
Gas_P2G=double(Gas_P2G); % P2G的产气量
CO2_CCS=double(CO2_CCS); % CCS的碳捕集量/P2G所用的二氧化碳量
P_h_GB=double(P_h_GB); % GB的产热功率
Gas_GB=double(Gas_GB); % GB的耗气量
P_buy=double(P_buy); % 综合能源系统向共享储能的购电量
P_sell=double(P_sell); % 综合能源系统向共享储能的售电量
miu_ies=double(miu_ies); % 综合能源系统与共享储能交互的状态变量
E_CO2=double(E_CO2); % 综合能源系统的实际碳排放量
E_c=double(E_c); % 综合能源系统的碳排放配额
E_1=double(E_1); % 阶梯碳交易引入的辅助区间变量
E_2=double(E_2);
E_3=double(E_3);
pri_e_1=double(pri_e_1); % 综合能源系统给用户1制定的电价
pri_h_1=double(pri_h_1); % 综合能源系统给用户1制定的热价
pri_e_2=double(pri_e_2); % 综合能源系统给用户2制定的电价
pri_h_2=double(pri_h_2); % 综合能源系统给用户2制定的热价
pri_e_3=double(pri_e_3); % 综合能源系统给用户3制定的电价
pri_h_3=double(pri_h_3); % 综合能源系统给用户3制定的热价
% 用户决策变量
p_cut_1=double(p_cut_1); % 用户1的可削减电负荷
h_cut_1=double(h_cut_1); % 用户1的可削减热负荷
p_load_1=double(p_load_1); % 用户1的购电量
h_load_1=double(h_load_1); % 用户1的购热量
p_cut_2=double(p_cut_2); % 用户2的可削减电负荷
h_cut_2=double(h_cut_2); % 用户2的可削减热负荷
p_load_2=double(p_load_2); % 用户2的购电量
h_load_2=double(h_load_2); % 用户2的购热量
p_cut_3=double(p_cut_3); % 用户3的可削减电负荷
h_cut_3=double(h_cut_3); % 用户3的可削减热负荷
p_load_3=double(p_load_3); % 用户3的购电量
h_load_3=double(h_load_3); % 用户3的购热量
p_tran_12=double(p_tran_12); % 用户1与用户2的交互电量
p_tran_13=double(p_tran_13); % 用户1与用户3的交互电量
p_tran_23=double(p_tran_23); % 用户2与用户3的交互电量
% 针对双线性项的McCormick线性化引入的变量
chi_e_1=double(chi_e_1); % 用来等效购电成本的辅助变量
chi_h_1=double(chi_h_1); % 用来等效购热成本的辅助变量
chi_e_2=double(chi_e_2); % 用来等效购电成本的辅助变量
chi_h_2=double(chi_h_2); % 用来等效购热成本的辅助变量
chi_e_3=double(chi_e_3); % 用来等效购电成本的辅助变量
chi_h_3=double(chi_h_3); % 用来等效购热成本的辅助变量
% 用户目标函数
obj_user_1=0.5*sum(p_cut_1)+0.4*sum(h_cut_1)+sum(chi_e_1+chi_h_1); % 用户1成本
obj_user_2=0.5*sum(p_cut_2)+0.4*sum(h_cut_2)+sum(chi_e_2+chi_h_2); % 用户2成本
obj_user_3=0.5*sum(p_cut_3)+0.4*sum(h_cut_3)+sum(chi_e_3+chi_h_3); % 用户3成本
% 传出数据
Power_buy=double(P_buy);
Power_sell=double(P_sell);
Obj_IES=double(Obj);
p_tran=[p_tran_12;p_tran_13;p_tran_23];
obj_user=[obj_user_1;obj_user_2;obj_user_3];
end
function [lambda,Obj,X_up,X_low]=Fun_ESO(U,Power_buy,Power_sell,X_v1,X_v2,X_up,X_low,iter)
% 共享储能分布式优化模型
%% 导入基础数据
% 综合能源系统向共享储能的购售电量(下层给上层的购售电量)
P_buy=Power_buy; % 综合能源系统向共享储能的购电量
P_sell=Power_sell; % 综合能源系统向共享储能的售电量
% 共享储能从主网购电价格
u_Db=[0.4,0.4,0.4,0.4,0.4,0.4,0.79,0.79,0.79,1.2,1.2,1.2,1.2,1.2,0.79,0.79,0.79,1.2,1.2,1.2,0.79,0.79,0.4,0.4];
% 共享储能向主网售电价格
u_Ds=[0.35,0.35,0.35,0.35,0.35,0.35,0.68,0.68,0.68,1.12,1.12,1.12,1.12,1.12,0.68,0.68,0.68,1.12,1.12,1.12,0.79,0.79,0.35,0.35];
% 共享储能与综合能源系统的交易价格上下限
u_Pbmax=1.5*u_Db; % 综合能源系统的购电价上限
u_Pbmin=u_Db; % 综合能源系统的购电价下限
u_Psmax=u_Ds; % 综合能源系统的售电价上限
u_Psmin=u_Psmax-0.35*ones(1,24); % 综合能源系统的售电价下限
u_Ps_ave=0.85; % 最大平均售电价
u_Pb_ave=1.20; % 最大平均购电价
%% 决策变量
u_Ps=sdpvar(1,24); % 综合能源系统向共享储能的售电价格(上层给下层的购售电价)
u_Pb=sdpvar(1,24); % 综合能源系统从共享储能的购电价格
P_ch=sdpvar(1,24); % 共享储能的充电功率
P_dis=sdpvar(1,24); % 共享储能的放电功率
miu_eso=binvar(1,24); % 共享储能的充放电状态变量
S_e=sdpvar(1,24); % 共享储能的储电量
P_grid_buy=sdpvar(1,24); % 共享储能向主网的购电量
P_grid_sell=sdpvar(1,24); % 共享储能向主网的售电量
%% 约束条件
C=[];
C=[C,
% 共享储能给综合能源系统的购售电价定价约束部分
u_Pbmin<=u_Pb<=u_Pbmax, % 共享储能给综合能源系统的购电价格上下限约束
u_Psmin<=u_Ps<=u_Psmax, % 共享储能给综合能源系统的售电价格上下限约束
sum(u_Pb)/24<=u_Pb_ave, % 共享储能给综合能源系统的购电价格平均值约束
sum(u_Ps)/24<=u_Ps_ave, % 共享储能给综合能源系统的售电价格平均值约束
% 共享储能运行约束部分
S_e(1)==0.2*5000+P_ch(1)*0.98-P_dis(1)/0.98, % 共享储能0-1时段的储电量变化约束
S_e(2:24)==S_e(1:23)+P_ch(2:24)*0.98-P_dis(2:24)/0.98, % 共享储能1-24时段的储电量变化约束
0.1*5000<=S_e<=0.9*5000, % 共享储能的储电量上下限约束
0<=P_ch<=3000*miu_eso, % 共享储能的充电功率上下限约束
0<=P_dis<=3000*(1-miu_eso), % 共享储能的放电功率上下限约束
0<=P_grid_buy<=3000, % 共享储能向主网的购电功率上下限约束
0<=P_grid_sell<=3000, % 共享储能向主网的售电功率上下限约束
% 共享储能的功率平衡约束部分
P_ch-P_dis==P_grid_buy-P_grid_sell+P_sell-P_buy, % 电功率平衡约束
];
lambda=[u_Pb,u_Ps]; % 给lambda赋值,方便写约束
% 执行步骤6,判断是否增加约束
if U==1
% 确实需要增加约束
% 将X_v1和X_v2的低值导入X_low,高值导入X_up
for t=1:48
if X_v1(t)<=X_v2(t)
X_low(iter,t)=X_v1(t);
X_up(iter,t)=X_v2(t);
else
X_low(iter,t)=X_v2(t);
X_up(iter,t)=X_v1(t);
end
end
% 添加约束
for k=1:iter
C=[C,X_low(k,:)<=lambda<=X_up(k,:),];
end
end
%% 目标函数
Obj=u_Ds*P_grid_sell'-u_Db*P_grid_buy'-0.01*sum(P_ch+P_dis)+sum(u_Pb.*P_buy-u_Ps.*P_sell);
%% 求解器配置与求解
ops=sdpsettings('solver','cplex','verbose',0,'usex0',0);
ops.cplex.mip.tolerances.mipgap=1e-6;
result=solvesdp(C,-Obj,ops); % 目标函数为共享储能的收益函数,所以有'负号'
%% 数据输出
lambda=double(lambda);
Obj=double(Obj);
最新发布