论文复现：风电、光伏与抽水蓄能电站互补调度运行（MATLAB-Yalmip全代码）

论文复现：风电、光伏与抽水蓄能电站互补调度运行（MATLAB-Yalmip全代码）

针对风电、光伏与抽水蓄能站互补运行的问题，已有大量通过启发式算法寻优的案例，但工程上更注重实用性和普适性。Yalmip工具箱则是一种基于MATLAB平台的优化软件工具箱,被广泛应用于工程界优化问题和控制理论问题求解，求解结果非常稳定。本文复现了论文中涉及的风-光-抽水蓄能经济运行代码，并提供模型说明（目标函数、约束条件）和能源套利调度策略。
参考论文：风-水电联合优化运行分析[J].太阳能学报,2008(01):80-84；风-光-水-火-抽蓄联合发电系统日前优化调度研究[J].太阳能学报,2020,41( 8) : 79-85．

订阅专栏可查看全代码。

1、风-光-抽水蓄能经济运行模型

1.1 目标函数

抽水蓄能电站作为能量型储能，在联合系统发电高峰期抽水，储存电能；在发电低谷期放水，发出电能。考虑弃电和缺电惩罚，利用分时电价的电价差，通过能源套利获利，实现互补调度系统收益最大化：

其中，C为分时电价，Pw/Ppv/Ph/Pp为风/光/抽水蓄能发电/抽水蓄能抽水功率，Fqd/Fqk为弃电惩罚和缺电惩罚。

1.2 约束条件

1.2.1 系统功率平衡

其中，PL为系统负荷。

1.2.2 发电功率限制约束

1.2.3 抽水功率限制约束

1.2.4 水库储能量约束

储能量不能超出范围

1.2.5 风/光出力限制约束

风光出力小于预测出力

1.2.6 工况限制约束

抽水-发电不可同时进行

2、实现步骤

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

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

2.2 主函数

clc
clear all
close all
%% 抽蓄-风-光优化模型参数设置
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;                  % 一天时长
[C,C_p]=price(n);        % 确定分时电价,调用price函数
C=C';                    % 上网分时电价
C_p=C_p';                % 抽蓄抽水电价
T=1:n;
t=1;                     % 尺度1h
P_hmax=200;P_hmin=0;     % 抽蓄最大、最小 发电 功率  MW
P_pmax=200;P_pmin=0;     % 抽蓄最大、最小 抽水 功率  MW
P_pps_r=200;             % 抽蓄额定功率
E_max=600;E_min=0;       % 水库储能
E_0=100;                 % 水库初始储能量 赋值
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=50; P_cmin=0;     % 电池 充放电  最大/小功率
P_dmax=50;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

%% 开始优化  抽蓄
% 申明变量
% 生成连续型变量空间
Chrom=sdpvar(1,8*n);
% 顺序 P_w   P_p  P_h    E   P_DL  P_pv   P_lack
% 顺序  风   抽水 发电 储能量 弃电    光     缺电
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_pv=Chrom(1,(5*n+1:6*n));
P_lack=Chrom(1,(6*n+1:7*n));
% 目标函数  经济效益最大  考虑分时电价、弃电惩罚、缺额惩罚成本
f=sum(C.*Chrom(1,1:n)+C.*Chrom(1,(2*n+1:3*n))+C.*Chrom(1,(5*n+1:6*n))-C_p.*Chrom(1,(n+1:2*n))-yibuxil_lack*Chrom(1,(7*n+1:8*n))-yibuxil_DL*Chrom(1,(4*n+1:5*n))-k_ps_h*Chrom(1,(2*n+1:3*n)));
% 约束条件
F=[];
for in=1:n
    F=[F P_w(1,in)+P_h(1,in)+P_pv(1,in)+P_lack(1,in)==P_load(in)];  % 抽蓄+风电+光伏+水电上网，功率平衡限制约束
    F=[F P_p(1,in)*P_h(1,in)==0];                           % 抽蓄抽水-发电不同时发生约束
    
    F=[F 0<=P_pv(1,in)<=P_pv1(1,in)];                       % 光伏上网限制约束
    F=[F 0<=P_w(1,in)<=P_v(1,in)];                          % 风电上网限制约束
    
    F=[F E_0<=E(1)<=E_0];                                   % 赋值 - 水库  初始储能E_0
    E(in+1)=E(in)+t*(eta_p*P_p(1,in)-P_h(1,in)/eta_h);      % 水库 储能变化
    F=[F 0<=E(in+1)<=E_max];                                % 水库 储能量限制约束
    
    F=[F P_pmin<=P_p(1,in)<=P_pmax];                        % 抽蓄 抽水约束
    F=[F P_hmin<=P_h(1,in)<=min(P_hmax,E(in)*eta_h/t)];     % 抽蓄 发电约束
      
    F=[F P_w(1,in)+P_pv(1,in)+P_p(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);
P_DL=double(P_DL);
P_pv=double(P_pv);
E=double(E);
P_lack=double(P_lack);
P_i=P_w+P_h+P_pv;
%% 结果可视化
figure(3)
grid on;
hold on;
plot(P_load,'k-+','linewidth',1);   
plot(P_w,'g-','linewidth',2);
plot(P_pv,'c-','linewidth',2);
plot(-P_p,'b-s','linewidth',2);
plot(P_h,'b-d','linewidth',2);
plot(P_i,'r-+','linewidth',1);
plot(-P_DL,'m-p','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]=price(n)
% CNY/MW
C=zeros(n,1);%先生成存储空间
C_p=zeros(n,1);
for i=1:7
    C(i)=324.9;
end

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

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

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

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

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

for i=24
    C(i)=324.9;
end
C_p=0.25*C; % 抽蓄 抽水 电费

2.4 输出结果