79号资源-源程序：论文可在知网下载《考虑寿命损耗的微网电池储能容量优化配置》本人博客有解读

%计及电池储能寿命损耗的微电网经济调度+三类需求侧响应 clear all;warning off; ee_bat_int=300;%电池储能的容量 p_wt_int=200;%风机 p_pv_int=100;%光伏 p_g_int=200;% p_g_min=10;%燃气轮机出力下限 p_m_max=500;%联络线功率上限 p_bat_int=ee_bat_int*0.21;%假设储能的功率上限和容量上限有比值关系 ee0=0.55*ee_bat_int; %储能初始电量 eta=0.95;%储能充放电效率 mm=100000;%一个极大正实数 k_suo=1/1;%缩减系数 c_bat_int=3320/3*k_suo;%缩减储能单位成本为3320/3,乘上了缩减系数 c_wt_om=0.0296;c_pv_om=0.0096;c_g_om=0.059;c_bat_om=0.009;%运维成本系数 c_fuel=0.6;%燃料成本系数 %决策变量 p_ch=sdpvar(24,1);p_dis=sdpvar(24,1);ee=sdpvar(24,1); uu_bat=binvar(24,1);uu_m=binvar(24,1); p_buy=sdpvar(24,1);p_sell=sdpvar(24,1); p_g=sdpvar(24,1); %风光出力和电价（以春季典型日为例） p_l=xlsread('四个典型日数据.xlsx','0%','B3:B26')*900; p_wt=xlsread('四个典型日数据.xlsx','0%','H3:H26')*p_wt_int; p_pv=xlsread('四个典型日数据.xlsx','0%','N3:N26')*p_pv_int; price=xlsread('四个典型日数据.xlsx','电价','A2:A25');%分时电价 %储能损耗模型的参数和变量 kk1=0.44268;bb1=0.43071; kk2=0.59493;bb2=0.41454; kk3=0.65646;bb3=0.40992; kk4=0.69405;bb4=0.40761; kk5=0.72954;bb5=0.40551; ss_bat=binvar(24,1);%改设为0/1变量 g1=binvar(24,1);g2=binvar(24,1);g3=binvar(24,1);g4=binvar(24,1);g5=binvar(24,1); d1=sdpvar(24,1);d2=sdpvar(24,1);d3=sdpvar(24,1);d4=sdpvar(24,1);d5=sdpvar(24,1); c_bat=sdpvar(1,1); constraints=[]; %负荷响应模块 jz_pri=0.9.*ones(1,24);%基准电价 detapr=price'-jz_pri;%电价差 load=DR3(p_l',detapr,price');%负荷调整 %储能功率约束 for t=1:24 constraints=[constraints,0<=p_dis(t),p_dis(t)<=(1-uu_bat(t))*p_bat_int];%由于u_bat定义不同，所以公式和原文不同，1为充电，0为放电 constraints=[constraints,0<=p_ch(t),p_ch(t)<=uu_bat(t)*p_bat_int]; end %soc约束 for t=1:24 constraints=[constraints,ee(t)==ee0+sum(eta*p_ch(1:t)-1/eta*p_dis(1:t))]; constraints=[constraints,0.1*ee_bat_int<=ee(t),ee(t)<=0.9*ee_bat_int]; end %功率平衡约束 for t=1:24 constraints=[constraints, p_pv(t)+p_wt(t)+p_dis(t)-p_ch(t)+p_buy(t)-p_sell(t)+p_g(t)==load(t)]; end %购售电约束 for t=1:24 constraints=[constraints, 0<=p_buy(t),p_buy(t)<=uu_m(t)*p_m_max]; constraints=[constraints,0<=p_sell(t),p_sell(t)<=(1-uu_m(t))*p_m_max];%购电1，售电0 end %充放电量平衡约束 constraints=[constraints, ee0==ee(24)]; %微型燃气轮机出力约束 for t=1:24 constraints=[constraints,p_g_min<=p_g(t),p_g(t)<=p_g_int]; end %爬坡约束 for t=2:24 constraints=[constraints, -10<=(p_g(t)-p_g(t-1))<=10]; end %储能寿命损耗约束 for t=1:24 constraints=[constraints,c_bat>=kk1*sum(ss_bat(:))+bb1+(g1(t)-1)*mm];%这里的c_bat不是关于t的函数，而是一个最大值 constraints=[constraints,c_bat>=kk2*sum(ss_bat(:))+bb2+(g2(t)-1)*mm]; constraints=[constraints,c_bat>=kk3*sum(ss_bat(:))+bb3+(g3(t)-1)*mm]; constraints=[constraints,c_bat>=kk4*sum(ss_bat(:))+bb4+(g4(t)-1)*mm]; constraints=[constraints,c_bat>=kk5*sum(ss_bat(:))+bb5+(g5(t)-1)*mm]; end %放电深度约束 for t=1:24 constraints=[constraints,0<=d1(t),d1(t)<=0.35*g1(t)]; constraints=[constraints,0.35*g2(t)<=d2(t),d2(t)<=0.65*g2(t)]; constraints=[constraints,0.65*g3(t)<=d3(t),d3(t)<=0.8*g3(t)]; constraints=[constraints,0.8*g4(t)<=d4(t),d4(t)<=0.9*g4(t)]; constraints=[constraints,0.9*g5(t)<=d5(t),d5(t)<=1*g5(t)]; end for t=1:24 constraints=[constraints,1-ee(t)/ee_bat_int==d1(t)+d2(t)+d3(t)+d4(t)+d5(t)]; end %循环次数约束 constraints=[constraints,ss_bat(1)>=0;ss_bat(1)>=uu_bat(1)];%默认u_bat(0)=0 for t=2:24 constraints=[constraints,ss_bat(t)>=0;ss_bat(t)>=uu_bat(t)-uu_bat(t-1)]; end %春季典型日 obj_o=sum(price.*(p_buy(:,1)-p_sell(:,1)))+c_fuel*sum(p_g(:,1))+...%购售电成本和燃料成本 sum(c_wt_om*p_wt(:,1))+sum(c_pv_om*p_pv(:,1))+sum(c_g_om*p_g(:,1))+sum(c_bat_om*p_dis(:,1))+sum(c_bat_om*p_ch(:,1))+...%运维成本 c_bat(1,1)/3*ee_bat_int*k_suo;%储能寿命损耗成本 %Ⅱ类负荷的转移成本！！！！ ops=sdpsettings('solver','cplex'); reuslt=optimize(constraints,obj_o,ops); obj_o=double(obj_o); p_ch=double(p_ch);p_dis=double(p_dis);p_buy=double(p_buy);p_sell=double(p_sell);p_g=double(p_g); uu_bat=double(uu_bat);uu_m=double(uu_m); ee=double(ee); figure; plot(value(p_ch)-value(p_dis),'b-*') hold on plot(value(p_g),'r-o') plot(p_pv,'m->') plot(p_wt,'c-p') plot(value(p_buy-p_sell),'k--') legend('储能出力','发电机出力','光伏','风电','购售电') grid on figure; plot(load,'r-o') hold on plot(p_l,'b-*') xlabel('时间（h）');ylabel('功率（kw）'); yyaxis right plot(price) legend('响应后负荷','响应前负荷','价格'); xlabel('时间（h）');ylabel('价格（元/kwh）'); grid on