87号资源-源程序：《基于储能电站服务的冷热电多微网系统 双层优化配置》《考虑微网新能源经济消纳的共享储能优化》本人博客有解读

%========================================================================== % 分布式光伏消纳的微电网群共享储能配置策略研究 % 参考文献：谢雨龙，罗逸飏，李智威，等.考虑微网新能源经济消纳的共享储能优化配置 % 关键词：共享储能，微能源网，经济消纳，双层规划，KKT条件，Cplex % 场景2：微网考虑经济消纳，配置共享储能 %========================================================================== clc; clear; close all; %% 模型参数设定 W=4; %典型日类型数为4（春夏秋冬） Tw=91; %典型日天数为91 N=3; %微网数为3 M=1E8; %Big-M法中的M NT=24; %调度时段数为24； % 下层的三个微网的运行参数导入（这里将4个典型日数据合并，形成一组24*W=96列的数据，方便后续用矩阵形式列写约束） P_cool=zeros(24*W,3); %下层的3个微网的冷负荷功率 P_heat=zeros(24*W,3); %下层的3个微网的热负荷功率 P_load=zeros(24*W,3); %下层的3个微网的电负荷功率 P_PV_0=zeros(24*W,3); %下层的3个微网的光伏发电功率 %导入微网1的参数 P_cool(:,1)=[xlsread('微网数据1','I2:I25');xlsread('微网数据1','P2:P25');xlsread('微网数据1','W2:W25');xlsread('微网数据1','AD2:AD25')]; P_heat(:,1)=[xlsread('微网数据1','J2:J25');xlsread('微网数据1','Q2:Q25');xlsread('微网数据1','X2:X25');xlsread('微网数据1','AE2:AE25')]; P_load(:,1)=[xlsread('微网数据1','K2:K25');xlsread('微网数据1','R2:R25');xlsread('微网数据1','Y2:Y25');xlsread('微网数据1','AF2:AF25')]; P_PV_0(:,1)=1.85*[xlsread('微网数据1','L2:L25');xlsread('微网数据1','S2:S25');xlsread('微网数据1','Z2:Z25');xlsread('微网数据1','AG2:AG25')]; %导入微网2的参数(MG2中没有风电) P_cool(:,2)=[xlsread('微网数据2','J2:J25');xlsread('微网数据2','Q2:Q25');xlsread('微网数据2','X2:X25');xlsread('微网数据2','AE2:AE25')]; P_heat(:,2)=[xlsread('微网数据2','K2:K25');xlsread('微网数据2','R2:R25');xlsread('微网数据2','Y2:Y25');xlsread('微网数据2','AF2:AF25')]; P_load(:,2)=[xlsread('微网数据2','L2:L25');xlsread('微网数据2','S2:S25');xlsread('微网数据2','Z2:Z25');xlsread('微网数据2','AG2:AG25')]; P_PV_0(:,2)=1.85*[xlsread('微网数据2','M2:M25');xlsread('微网数据2','T2:T25');xlsread('微网数据2','AA2:AA25');xlsread('微网数据2','AH2:AH25')]; %导入微网3的参数 P_cool(:,3)=[xlsread('微网数据3','I2:I25');xlsread('微网数据3','P2:P25');xlsread('微网数据3','W2:W25');xlsread('微网数据3','AD2:AD25')]; P_heat(:,3)=[xlsread('微网数据3','J2:J25');xlsread('微网数据3','Q2:Q25');xlsread('微网数据3','X2:X25');xlsread('微网数据3','AE2:AE25')]; P_load(:,3)=[xlsread('微网数据3','K2:K25');xlsread('微网数据3','R2:R25');xlsread('微网数据3','Y2:Y25');xlsread('微网数据3','AF2:AF25')]; P_PV_0(:,3)=1.85*[xlsread('微网数据3','L2:L25');xlsread('微网数据3','S2:S25');xlsread('微网数据3','Z2:Z25');xlsread('微网数据3','AG2:AG25')]; %微网向储能电站售电单位电能电价矩阵 delt_s = [0.20,0.20,0.20,0.20,0.20,0.20,0.20,0.20,0.95,0.95,0.95,0.95,0.55,0.55,0.55,0.55,0.95,0.95,0.95,0.95,0.95,0.55,0.55,0.55]; delt_s = [delt_s,delt_s,delt_s,delt_s]; %微网从储能电站购电单位电能电价矩阵 delt_b = [0.40,0.40,0.40,0.40,0.40,0.40,0.40,0.40,1.15,1.15,1.15,1.15,0.75,0.75,0.75,0.75,1.15,1.15,1.15,1.15,1.15,0.75,0.75,0.75]; delt_b = [delt_b,delt_b,delt_b,delt_b]; %微网向储能电站缴纳单位功率服务费用矩阵 delt_f = 0.05*ones(1,24*W); %微网从电网的购电电价 delt_0 = [0.37,0.37,0.37,0.37,0.37,0.37,0.37,0.37,1.36,1.36,1.36,1.36,0.82,0.82,0.82,0.82,1.36,1.36,1.36,1.36,1.36,0.82,0.82,0.82]; delt_0 = [delt_0,delt_0,delt_0,delt_0]; %微网单位体积天然气费用矩阵，取江苏地区气价2.2元/m3 c_0 = 2.2*ones(1,24*W); %储能电站的容量成本参考某储能项目磷酸铁锂电池中标价格 r = 0.08;%资金年利率（取值为8%） gamma = 8;%装置寿命周期（年） delt_P = 1000;% 单位功率投资成本（元/kW） delt_E = 1897;% 单位容量投资成本（元/kWh） delt_M = 72;% 单位功率维护成本 belta = 2.74; % 储能电池能量倍率 yita_abs = 0.95;% 储能电站充电效率 yita_relea = 0.95;% 储能电站放电效率 k_min = 0.1;% 电站荷电状态下限 k_max = 0.9;% 电站荷电状态上限 P_ess_max = 5000;% 微网与共享储能电站最大交换功率 P_grid_max = 10000;% 从电网购电的最大功率 P_GT_min = 0;% 燃气轮机最小功率（kW） P_GT_max = 1500;% 燃气轮机最大功率（kW） P_EC_min = 0;% 电制冷机最小功率 P_EC_max = 4000;% 电制冷机最大功率 P_HX_min = 0;% 换热装置最小功率 P_HX_max = 4000;% 换热装置最大功率 Q_AC_min = 0;% 吸收式制冷机最小功率 Q_AC_max = 4000;% 吸收式制冷机最大功率 Q_GB_min = 0;% 燃气锅炉最小功率 Q_GB_max = 4000;% 燃气锅炉最大功率 Q_0 = 9.7;% 燃气热值（kWh/m3） yita_EC = 4;% 电制冷机能效比 yita_AC = 1.2;% 吸收式制冷机能效比 yita_HX = 0.9;% 换热装置效率 yita_WH = 0.8;% 余热锅炉效率 yita_GT = 0.3;% 燃气轮机发电效率 yita_GB = 0.9;% 燃气锅炉效率 gamma_GT = 1.47;% 燃气轮机热电比 %% 决策变量初始化 alpha = sdpvar(1); %经济消纳率 P_essmax=sdpvar(1); %储能电站最大充放电功率（即配置的储能功率） E_essmax=sdpvar(1); %储能电站的最大容量（即配置的储能容量） P_esss=sdpvar(3,24*W); %第i个微网向储能电站的售电功率 P_essb=sdpvar(3,24*W); %第i个微网向储能电站的购电功率 P_GT=sdpvar(3,24*W); %第i个微网的燃气轮机输出功率 P_grid=sdpvar(3,24*W); %第i个微网从电网的购电功率 P_EC=sdpvar(3,24*W); %第i个微网电制冷机消耗的电功率 P_PV = sdpvar(3,24*W); %微网光伏实际出力 Q_AC=sdpvar(3,24*W); %第i个微网制冷机的输出制冷功率 Q_GB=sdpvar(3,24*W); %第i个电网的燃气轮机输出热功率 P_HX=sdpvar(3,24*W); %第i个电网的换热装置输出制热功率 E_ess=sdpvar(1,24*W); %储能电站存储的能量 P_essabs=sdpvar(1,24*W); %储能电站充电功率 P_essrelea=sdpvar(1,24*W); %储能电站放电功率 U_abs=binvar(1,24*W); %储能电站充放电状态，0-1变量 U_relea=binvar(1,24*W); %储能电站充放电状态，0-1变量 %% 上层储能电站的目标函数 C_inv = ((r*(1+r)^gamma)/((1+r)^gamma-1))*(delt_P*P_essmax+delt_E*E_essmax)+delt_M*P_essmax; C_grid = Tw*sum(P_grid)*delt_0'; C_flue = Tw*sum(P_GT/(yita_GT*Q_0)+Q_GB/(yita_GB*Q_0))*c_0'; Obj = C_inv+C_grid+C_flue; %% 上层储能电站的约束条件 Con = []; %式（5）中储能电站容量和额定功率之间的倍率约束 Con=[Con E_essmax==belta*P_essmax E_ess([1,25,49,73])==0.2*E_essmax ]; %式（6）充放电约束，因为不能出现两个决策变量相乘，故这里用类似于大M法的思路代替 Con=[Con sum(P_essb)-sum(P_esss)==P_essrelea-P_essabs %式(6)储能电站充放电功率平衡约束 0<= P_essabs <= P_essmax*ones(1,24*W) 0<= P_essabs <= U_abs*M 0<= P_essrelea <= P_essmax*ones(1,24*W) 0<= P_essrelea <= U_relea*M ]; %式（7）储能电池荷电状态约束 Con=[Con E_ess(1,2:96)==E_ess(1,1:95)+yita_abs*P_essabs(1,2:96)-(1/yita_relea)*P_essrelea(1,2:96) % 式（7）储能电站与上一时段的能量变化约束 k_min*E_essmax<=E_ess<=k_max*E_essmax %式（7）储能电站容量上下限 U_abs+U_relea<=ones(1,24*W)%确定充放电状态不可能同时发生 ]; %% 进行下层的KKT条件处理 %下层等式约束导入 Con=[Con P_GT+P_PV+P_grid+P_essb-P_esss-P_EC-P_load'==0 %式(12)电功率平衡约束 P_HX/yita_HX+Q_AC/yita_AC-gamma_GT*yita_WH*P_GT==0 %式(15-1)余热锅炉余热平衡约束 Q_GB+P_HX-P_heat'==0 %式(15-2)热功率平衡约束 yita_EC*P_EC+Q_AC-P_cool'==0 %式(15-3)冷功率平衡约束 sum(sum(P_PV))==alpha*sum(sum(P_PV_0')); ]; %处理KKT条件 %定义等式约束的拉格朗日乘子 l1itw=sdpvar(3,24*W); l2itw=sdpvar(3,24*W); l3itw=sdpvar(3,24*W); l4itw=sdpvar(3,24*W); l5itw=sdpvar(1,24*W); l6itw=sdpvar(3,24*W); %定义不等式约束的拉格朗日乘子 u1itmin=sdpvar(3,24*W); u2itmin=sdpvar(3,24*W); u3itmin=sdpvar(3,24*W); u4itmin=sdpvar(3,24*W); u5itmin=sdpvar(3,24*W); u6itmin=sdpvar(3,24*W); u7itmin=sdpvar(3,24*W); u8itmin=sdpvar(3,24*W); u1itmax=sdpvar(3,24*W); u2itmax=sdpvar(3,24*W); u3itmax=sdpvar(3,24*W); u4itmax=sdpvar(3,24*W); u5itmax=sdpvar(3,24*W); u6itmax=sdpvar(3,24*W); u7itmax=sdpvar(3,24*W); u8itmax=sdpvar(3,24*W); u9itmax=sdpvar(3,24*W); u10itmin = sdpvar(3,24*W); u10itmax = sdpvar(3,24*W); U_sale=binvar(3,24*W); U_buy=binvar(3,24*W); %依次添加附录（5）-（14）KKT平衡条件约束 for i=1:3 for t=1:24*W Con=[Con Tw*c_0(t)/(yita_GT*Q_0)+l1