%% 请先确保YALMIP工具箱和CPLEX求解器正确安装
%% 初始化
clc;
clear all;
close all;
%% 导入微网
Pw=sdpvar(1,24,'full'); %风机出力
Ppv=sdpvar(1,24,'full');%光伏出力
Pbat=sdpvar(1,24,'full');%蓄电池出力
Pnet=sdpvar(1,24,'full');%交换功率
Pbuy=sdpvar(1,24,'full');%从电网购电电量
Psell=sdpvar(1,24,'full');%向电网售电电量
Temp_net=binvar(1,24,'full'); % 购|售电标志
Temp_cha=binvar(1,24,'full'); %充电标志
Temp_dis=binvar(1,24,'full'); %放电标志
Temp_static=binvar(1,24,'full'); %电池静置标志
Pcha=sdpvar(1,24);
Pdis=sdpvar(1,24);
Constraints = [];
load load.mat %微网基础数据读入
load fengji_predict.mat %风机预测出力
load pv_predict.mat %光伏预测出力
load price.mat %分时电价
load chuneng.mat %储能电池参数定义
figure(1);
bar(Load);
hold on
plot(Pw,'-gd','Color',[1,0,0],'LineWidth',1,'MarkerEdgeColor','R','MarkerFaceColor',[1 0 0],'MarkerSize',5);
hold on
plot(Ppv,'-mo','Color',[0,1,0],'LineWidth',1,'MarkerEdgeColor','R','MarkerFaceColor',[1 0 0],'MarkerSize',5);
hold off
axis([0 24 0 240]) ;
xlabel('时间/h');
set(gca,'xTick',(0:2:24),'yTick',(0:40:240));
ylabel('功率/kW');
legend('电负荷','风机预测出力','光伏预测出力');
legend('Location', 'northwest');
title('微网基本情况');
%% 约束
for k = 1:24
Constraints = [Constraints, -POWER<=Pnet(1,k)<=POWER,0<=Pbuy(1,k)<=POWER, -POWER<=Psell(1,k)<=0]; %主网功率交换约束，不大于160kW
Constraints = [Constraints,Pnet(1,k)+Pw(1,k)+Ppv(1,k)==Load(1,k)+Pbat(1,k)]; %功率平衡约束，电网+风电+光伏 = 负载 + 电池充电
Constraints = [Constraints, implies(Temp_net(1,k),[Pnet(1,k)>=0,Pbuy(1,k)==Pnet(1,k),Psell(1,k)==0])]; %购电情况约束，Pnet>0是购电，Pnet<0是售电
Constraints = [Constraints, implies(1-Temp_net(1,k),[Pnet(1,k)<=0,Psell(1,k)==Pnet(1,k),Pbuy(1,k)==0])]; %售电情况约束
%蓄电池约束
Constraints = [Constraints, -Pcs<=Pbat(1,k)<=Pcs,0<=Pcha(1,k)<=Pcs,-Pcs<=Pdis(1,k)<=0]; %电池充放电约束,PCS功率是40kW
Constraints = [Constraints, implies(Temp_cha(1,k),[Pbat(1,k)>=0,Pcha(1,k)==Pbat(1,k),Pdis(1,k)==0])]; %充电情况约束
Constraints = [Constraints, implies(Temp_dis(1,k),[Pbat(1,k)<=0,Pdis(1,k)==Pbat(1,k),Pcha(1,k)==0])]; %放电情况约束
Constraints = [Constraints, implies(Temp_static(1,k),[Pbat(1,k)==0,Pdis(1,k)==0,Pcha(1,k)==0])]; %静置情况约束
Constraints = [Constraints,Temp_cha(1,k)+Temp_dis(1,k)+Temp_static(1,k)==1];
Constraints=[Constraints,Ebattery*(socmin-soc0)<=sum(Pdis(1,1:k)+Pcha(1,1:k))<=Ebattery*(socmax-soc0)]; %SOC约束，电池容量300kwh，初始S0C为0.4，0.3<=SOC<=0.95
end
Constraints=[Constraints,sum(Pdis+Pcha)==0]; %ST=S0，始末SOC相等约束
%% 目标函数
F=0;
%总费用计算
for k = 1:24
F = F + 0.52*Pw(:,k)+0.72*Ppv(:,k)+C_buy(:,k)*Pbuy(:,k)+C_sell(:,k)*Psell(1,k)+0.2*Pdis(1,k);
end
ops=sdpsettings('solver', 'cplex');%参数指定程序用cplex求解器
optimize(Constraints,F,ops);
disp(['总费用=']);
value(F)
%% 画图
x=1:24;
PP=[Pbuy;-Pdis;Pw;Ppv];
PP_neg=[Psell;-Pcha];
figure(2)
bar(PP','stack');
hold on
h=legend('电网交换功率','蓄电池充电功率','风机出力','光伏出力','Location','northwest');
set(h,'Orientation','horizon');
bar(PP_neg','stack');
hold on
plot(x,value(Load),'r','linewidth',2);
hold off
xlabel('时间/h');
set(gca,'xTick',(1:2:24),'yTick',(0:50:300));
ylabel('功率/kW');
figure(3)
plot(x,value(Pbuy+Psell),'-ro','Color',[1,0,0],'LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor',[0 0 1],'MarkerSize',5);
hold on
plot(x,value(Pdis+Pcha),'-gd','Color',[1,0,1],'LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor',[0 0 1],'MarkerSize',5);
hold off
xlabel('时间/h');
set(gca,'xTick',(0:2:24),'yTick',(-120:20:60));
ylabel('功率/kW');
h=legend('电网购售电功率 >0是购电 <0是售电','蓄电池充放电功率 >0是充电 <0是放电');
soc = zeros(1,25);
s = zeros(1,25);
soc(1)=soc0;
for k=1:24
s(k)=value(sum(Pdis(1,1:k)+Pcha(1,1:k)))/Ebattery+soc0;
soc(k+1)=s(k);
end
soc(1)=soc0;
xx=0:24;
figure(4)
plot(xx,soc,'-ro','Color',[1,0,0],'LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor',[0 0 1],'MarkerSize',5);
xlabel('时间/h');
set(gca,'xTick',(0:2:24),'yTick',(0.2:0.1:1));
ylabel('SOC值');
title('蓄电池SOC状态');
