论文复现：风电、光伏与储能（含电池和废弃矿井小型抽水蓄能）互补调度运行（MATLAB-Yalmip全代码）

论文复现：风电、光伏与储能（含电池和废弃矿井小型抽水蓄能）互补调度运行（MATLAB-Yalmip全代码）

针对风电、光伏与水电站互补运行的问题，已有大量通过启发式算法寻优的案例，但工程上更注重实用性和普适性。Yalmip工具箱则是一种基于MATLAB平台的优化软件工具箱,被广泛应用于工程界优化问题和控制理论问题求解。本文复现了论文中涉及的风-光-储互补调度运行代码，并提供模型说明（目标函数、约束条件）和能源套利调度策略。
参考论文：风-水电联合优化运行分析[J].太阳能学报,2008(01):80-84；面向光伏消纳的光伏-废弃矿井抽蓄-蓄电池联合发电系统优化调度策略[J].科技导报,2021,39(13):52-58.
订阅专栏可查看全代码。

1、风-光-电池-废弃矿井小型抽水蓄能经济运行模型

1.1 目标函数

考虑运行成本、弃电和缺电惩罚，利用分时电价的电价差，储能通过能源套利获利，实现互补调度系统收益最大化。目标函数数学模型参考第1篇博客和文献【面向光伏消纳的光伏-废弃矿井抽蓄-蓄电池联合发电系统优化调度策略】。

1.2 约束条件

除第1篇博客中的约束条件外，还需要满足电池SOC和充放电约束，详见文献【面向光伏消纳的光伏-废弃矿井抽蓄-蓄电池联合发电系统优化调度策略】。

2、实现步骤

2.1 环境准备：下载、安装Yalmip工具箱

参见本专栏文章：Yalmip工具箱下载及使用说明

2.2 主函数

clc
clear all
close all
%% 风-光-混合储能优化模型参数设置
% P_v=[334.9,262.5,219,168.9,146.8,148.9,166.8,199.6,255.3,275.4,276.6,297.8,310.6,273.7,278.7,338.4,254.2,205.4,125.2,64.6,73.0,127.8,149.8,142.7];
% P_pv1=[0,0,0,0,0,11.8,34.3,88.1,137.3,166.4,158.2,179.2,176.6,163.0,127.5,104.3,118.7,81.8,37.7,8.6,0,0,0,0];
% P_load=[542.3,512.9,508.0,490.4,496.5,550.0,633.3,668.7,697.4,729.1,755.7,753.2,690.3,709.8,721.1,718.1,7409,7531,7365,7402,7747,7307,6594,5892];
P_v=[196,218,243,237,243,187,191,172,193,243,295,321,341,286,293,273,281,303,229,240,213,216,221,202];
P_pv1=[0,0,0,0,0,0,0,0,44,93,143,161,174,191,193,176,125,51,36,0,0,0,0,0];
P_load=[498,460,461,429,435,442,467,480,494,543,569,583,572,582,590,592,580,628,596,576,618,619,578,543];

n=24;                      % 一天24h
[C,C_p,C_c]=price(n);      % 确定分时电价,调用price函数，C_p为抽水费用、C_c为电池充电费用
C=C';                      % 上网分时电价
C_p=C_p';                  % 抽蓄抽水电价
C_c=C_c';                  % 电池运行成本 
T=1:n;
t=1;                       % 尺度1h
SOC_0=0.5;                 % SOC初始荷电状态 赋值
eta_p=0.87;                % 水泵抽水效率
eta_h=0.85;                % 水力发电效率
eta_c=0.9;                 % 电池 充电效率
eta_d=0.9;                 % 电池 放电效率
Emax_0=100;Emin=0;         % 电池 最大、最小容量MWh
P_cmax=100; P_cmin=0;      % 电池 充放电  最大/小功率
P_dmax=100;P_dmin=0;       
SOCmin=0.2;SOCmax=0.8;     % 电池 荷电状态 上下限
xgma=0.25/30/24;           % 电池每小时自放电率  20~30%/月
yibuxil_lack=100;          % 缺电惩罚系数  元/MWh
yibuxil_DL=30;             % 弃电惩罚系数  元/MWh
M_co2=0.877;               % 火电厂发单位电量产生的CO2量  tco2/MWh
k_ps_h=46;                 % 抽蓄发电运行成本    元/MWh
k_ba_d=28.7;               % 电池发电运行成本    元/MWh
P_hmax=150;P_hmin=0;       % 抽蓄最大、最小 发电 功率  MW
P_pmax=150;P_pmin=0;       % 抽蓄最大、最小 抽水 功率  MW
P_pps_r=150;               % 抽蓄额定功率
E_max_0=300;E_min=0;       % 水库储能
E_0=150;                   % 水库初始储能量 赋值
%% 开始优化  抽蓄+电池
% 申明变量
Chrom=sdpvar(1,10*n); % 整形 intvar      实形 sdpvar  
%P_w  P_p  P_h  E  P_DL  P_c  P_d  SOC  P_lack   P_pv 
% 1    2    3   4   5     6    7    8     9       10   
Chrom_E=sdpvar(1,2); % 抽蓄和电池容量
I_q=binvar(1,n);    %抽蓄启动抽水  0/1
I_h=binvar(1,n);    %抽蓄启动发电  0/1
% 目标函数
f=sum(C.*Chrom(1,1:n)+C.*Chrom(1,(2*n+1:3*n))+C.*Chrom(1,(6*n+1:7*n))+C.*Chrom(1,(9*n+1:10*n))-C_p.*Chrom(1,(n+1:2*n))-C_c.*Chrom(1,(5*n+1:6*n))-yibuxil_lack*Chrom(1,(8*n+1:9*n))-yibuxil_DL*Chrom(1,(4*n+1:5*n))-k_ps_h*Chrom(1,(2*n+1:3*n))-k_ba_d*Chrom(1,(6*n+1:7*n)));

% 约束条件
F=[];
P_w = Chrom(1,1:n);
P_p = Chrom(1,(n+1:2*n));
P_h = Chrom(1,(2*n+1:3*n));
E =   Chrom(1,(3*n+1:4*n));
P_DL =Chrom(1,(4*n+1:5*n));
P_c = Chrom(1,(5*n+1:6*n));
P_d = Chrom(1,(6*n+1:7*n));
SOC = Chrom(1,(7*n+1:8*n));
P_lack = Chrom(1,(8*n+1:9*n));
P_pv =Chrom(1,(9*n+1:10*n));
E_max=Chrom_E(1,1);           % 申明 水库最大储能量 便于-混合储能-计算
Emax= Chrom_E(1,2);           % 申明 电池最大容量   便于-混合储能-计算
% 得出 混合储能量 的上下界
P_net_max = P_load-P_v;
P_pc=zeros(1,n);P_hd=zeros(1,n);% 生成存储空间
for im=1:n
   if P_net_max(im) >= 0
       P_pc(im)=P_net_max(im); % 最大充电 上限
   end
   if P_net_max(im) < 0
       P_hd(im)=P_net_max(im); % 最大发电 上限
   end
end
for in=1:n
    F=[F P_w(1,in)+P_pv(1,in)+P_h(1,in)+P_d(1,in)+P_lack(1,in)==P_load(in)]; % 抽蓄+电池+风电+光伏上网，功率平衡限制约束
    
    F=[F 0<=P_pv(1,in)<=P_pv1(1,in)];                            % 光伏上网限制约束
    F=[F 0<=P_w(1,in)<=P_v(1,in)];                               % 风电上网限制约束
    
    F=[F P_p(1,in)*P_h(1,in)==0];                                 % 抽蓄抽水-发电   不同时发生约束
    F=[F P_c(1,in)*P_d(1,in)==0];                                 % 电池充-放电     不同时发生约束
    
    F=[F (E_max_0<=E_max<=E_max_0)];                              % 赋值 - 水库 最大储能量
    F=[F (Emax_0<=Emax<=Emax_0)];                                 % 赋值 - 电池 最大储能量
    F=[F E_0<=E(1)<=E_0];                                         % 赋值 - 水库     初始储能,最大值现为10MWh
    
    E(in+1)=E(in)+t*(eta_p*P_p(1,in)-P_h(1,in)/eta_h);            % 水库 储能变化
    F=[F E_min<=E(in+1)<=E_max];                                  % 水库 储能量限制约束
    
    F=[F SOC_0<=SOC(1)<=SOC_0];                                    % 赋值 - 电池荷电状态  初始值 E_0=0.5
    SOC(in+1)=SOC(in)*(1-xgma)+t*(eta_c*P_c(1,in)/Emax)-t*P_d(1,in)/Emax/eta_d;   % 电池荷电状态 约束
    F=[F SOCmin<=SOC(1,in+1)<=SOCmax];                                            % 电池 储能量限制约束
    F=[F P_cmin<=P_c(1,in)<=min(P_cmax,Emax*(SOCmax-SOC(in)*(1-xgma))/t/eta_c)];  % 电池 充电功率限制 #--- m ---#
    F=[F P_dmin<=P_d(1,in)<=min(P_dmax,(SOC(in)*(1-xgma)-SOCmin)*eta_d*Emax/t)];  % 电池 放电功率限制 #--- p ---#
    
    
    F=[F P_pmin<=P_p(1,in)<=min(P_pmax,(E_max-E(in))/t/eta_p)];                   % 抽蓄 抽水约束     #--- n ---#
    F=[F P_hmin<=P_h(1,in)<=min(P_hmax,E(in)*eta_h/t)];                           % 抽蓄 发电约束     #--- q ---#
    
    F=[F P_w(1,in)+P_pv(1,in)+P_p(1,in)+P_c(1,in)+P_DL(1,in)==P_v(1,in)+P_pv1(1,in)];  % 能量守恒约束
    
    F=[F 0<=P_DL(1,in)];                                                          % 功率舍弃量 约束
    F=[F 0<=P_lack(1,in)];                                                        % 缺口量 约束
end
% 运算求解
output=solvesdp(F,-f);
z=double(f);
P_w=double(P_w);
P_p=double(P_p);P_h=double(P_h);E=double(E);       % 抽蓄 参数
P_c=double(P_c);P_d=double(P_d);SOC=double(SOC);   % 电池 参数
P_pv=double(P_pv);
P_DL=double(P_DL);                                 % 弃电值
P_lack=double(P_lack);                             % 缺口值
E_max=double(E_max);Emax=double(Emax);             % 抽蓄、电池最大容量
P_i=P_w+P_pv+P_h+P_d;%联合系统上网功率
%% 输出结果 
% 绘图
figure(4)
grid on;
hold on;
plot(P_load,'k-+','linewidth',1);  
plot(P_w,'g-','linewidth',2);
plot(P_pv,'k-+','linewidth',2);
plot(-P_p,'b-','linewidth',2);plot(P_h,'r-','linewidth',2);
plot(-P_c,'g--*','linewidth',2);plot(P_d,'g--+','linewidth',2);
plot(P_i,'r-','linewidth',1);
plot(P_DL,'y-','linewidth',1)
plot(P_lack,'m-h','linewidth',1);
legend('负荷','风电','光伏','抽蓄抽水','抽蓄发电','电池充电','电池放电','联合上网','弃电值','缺电值');
set(gca,'xtick',[1 6 12 18 24],'xticklabel',{'01:00' ,'06:00' ,'12:00' ,'18:00', '24:00'  });
xlabel('时刻')
ylabel('功率/MW')
axis([1,24,-inf,inf]);

2.3 电价函数

分时电价源于山东省分时电价：

function [C,C_p,C_c]=price(n)
% 输入参数
C=zeros(n,1);%先生成存储空间
C_p=zeros(n,1);
C_c=zeros(n,1);
for i=1:7
    C(i)=312.1;
end

for i=8:9
   C(i)=617.7;
end

for i=10:15
   C(i)=905.7;
end

for i=16:20
   C(i)=617.7;
end

for i=21:22
   C(i)=1087.9;
end

for i=23
   C(i)=617.7;
end

for i=24
    C(i)=312.1;
end
C_p=0.25*C; % 抽蓄 抽水 电费
C_c=0.1*C;  % 电池 抽水 电费

2.4 输出结果