前言
在能源安全、环境污染和气候变化的大背景下，大力发展可再生能源是应对全球气候变化，实现“碳达峰、碳中和”和可持续发展的重大需求。截至2020年底，中国风电总装机容量为281GW，风力发电466.5TW·h，同比增长约15%；光伏发电总装机容量为253 GW，光伏发电260.5TW·h，同比增长约16%。风电、光伏等可再生能源出力具有波动性、不确定性以及不可调峰性等特点，大规模可再生能源接入对电力系统调峰能力提出了更高要求。现行的调峰成本分摊机制按照上网电量或电费的比例在发电侧内部平衡，无法反映负荷和可再生能源等主体对调峰成本引发的程度。明确调峰辅助服务的需求主体，制定公平合理的调峰费用分摊机制将有利于促进发电企业为消纳可再生能源提供辅助服务。
调峰成本分摊机制
调峰需求主体
调峰服务是指发电机组为了跟踪负荷的峰谷变化而有计划地、按照一定调节速度调整发电机组出力所提供的服务。区别于表现为静态容量特性的备用服务，调峰辅助服务更多地表现为应对负荷变化的动态调节特性。随着可再生能源上网电量比例的增加，系统净负荷的波动变化逐渐偏离系统负荷的波动变化。常规机组需要同时调节负荷的峰谷变化和可再生能源出力的波动，调峰服务的具体体现为对净负荷波动变化的调节。
风电、光伏等可再生能源受自然资源的影响，出力具有波动性和不可调峰性等特征，和负荷比较类似，可以看成是一种“负”负荷。
下图所示为含可再生能源消纳的电力系统调峰示意图。上调峰极限与净负荷之间的差值表示系统预留的正旋转备用，下调峰极限与净负荷之间的差值为系统的负旋转备用，系统的正、负旋转备用要求一般取负荷的一定比例。常规机组不仅要满足正负备用需求，还需要跟踪净负荷的峰谷变化，常规机组提供的最大下调峰深度随着净负荷峰谷差的加大而增加。当净负荷波动处于系统基本调峰范围内时，仅通过基本调峰就能满足净负荷的峰谷变化。当净负荷波动变化超过系统基本调峰临界值时，则需要根据经济性原则投入抽水蓄能、深度调峰、启停调峰等高成本调峰方式。随着净负荷峰谷差加大，投入的高成本调峰资源越多，调峰成本越大。当消纳可再生能源导致的调峰成本大于可再生能源的价值时，适当削减可再生能源上网电量有利于电力系统的经济运行。
调峰成本分析
在量化调峰需求主体导致的调峰成本时，借鉴了分析风电波动成本时用到的替代场景思路，但不同之处主要有两点：其一，研究对象不局限于风电，对每个调峰需求主体建立了无调峰需求的替代场景，将负荷和可再生能源出力曲线分别转换为无波动的均值线；其二，综合考虑了不同调峰需求主体对系统调峰成本的相互影响。单一主体有无调峰需求两种情况下系统调峰成本的差值仅能代表该主体导致的边际调峰成本。将调峰需求主体导致的加权平均边际调峰成本作为各自导致的调峰成本，这种方法既体现了不同主体之间对系统调峰成本的相互影响，也体现了公平性原则。对于负荷，其替代场景可以理解为负荷转移；对于可再生能源，其替代场景则可以理解为与储能或需求响应资源联合实现，如风电厂联合抽水蓄能电站实现电能的平缓输出。
调峰成本分摊方法
Shapley值是经济学上解决公共费用分摊与补偿问题最常用的方法之一，其思想为参与者应承担的成本等于参与者对其所参与联盟的平均边际贡献。为体现不同调峰需求主体对系统调峰成本的相互影响，采用Shapley值计算不同调峰需求主体导致的调峰成本，并以此为基础分摊有偿调峰成本。这种分摊方法体现了调峰需求主体对系统调峰成本的影响，符合“谁引发谁承担”的原则。由于目前可再生能源参与辅助服务的市场交易机制尚不完善，文中未考虑负荷和可再生能源主动参与调峰辅助服务的情况。另外，基于Shapley值的调峰成本量化方法需要2|N|次调度优化计算，面临一定的计算困难，因此这种方法适用于调峰需求主体有限的情况。为了降低计算难度，可采用分类或聚类等方法减少Shapley值计算的成员数量，然后对同一类别内的成员采用相同的规则分摊成本。
程序介绍
程序从净负荷波动的角度出发，建立了调峰成本的量化与分摊模型，构造了无调峰需求的替代场景，将负荷和可再生能源出力曲线分别转换为无波动的均值线。其次，建立了含深度调峰和抽水蓄能的调度优化模型，用于计算不同场景下的调峰成本，并将有无调峰需求两种场景下的系统调峰成本之差作为单一主体导致的边际调峰成本，然后采用Shapley值计算不同主体导致的调峰成本。最后，根据成本的引发程度分摊调峰成本。程序中算例丰富，注释清晰，干货满满，创新性和可扩展性很高，足以撑起一篇高水平论文！下面对程序做简要介绍！
程序适用平台：Matlab+Yalmip+Gurobi
参考文献：《高比例可再生能源电力系统的调峰成本量化与分摊模型》-电力系统自动化
程序结果
部分程序
%% 公式（6） 目标函数：调峰成本
Cost = sdpvar(1,1); %系统总调峰成本P300_Diot = sdpvar(18,4,24,'full');
P600_Diot = sdpvar(18,4,24,'full');
y300_it = binvar(18,24,'full'); % 机组i在时段t是否开启动作的0-1变量
z300_it = binvar(18,24,'full'); % 机组i在时段t是否停机动作的0-1变量
y600_it = binvar(18,24,'full'); z600_it = binvar(18,24,'full');
%% 公式（7） 电平衡约束
%需要注意的是，AA(i,j,t)，无法通过AA(1,1,1:24)的方式，调取出一个向量
%只能通过AA(1,1,1)的方式，调取出一个数值，范例如下：
% A = ones(18,4,24);
% AA = sum(sum(A(:,:,1)));
% AAA=AA(1,1,1:24);
P_PuG = sdpvar(4,24,'full'); %抽水蓄能发电
%% 公式（8-9）系统正负旋转备用约束
u200Git = binvar(20,24,'full'); %uit 表示机组运行状态的0-1变量,1开机，0关机
u300Git = binvar(18,24,'full'); u600Git = binvar(18,24,'full');
R_Ureq = (Eload*0.05+E_Wind*0.05+E_PV*0.05); %正旋转备用
R_Dreq = (Eload*0.05+E_Wind*0.05+E_PV*0.05); %负旋转备用
%% 公式（10-12）火电机组功率大小约束，这里换一种方法去写
uP300_Diot = binvar(18,4,24,'full');uP600_Diot = binvar(18,4,24,'full');
%公式（11-12）不用再写了，已经隐含在上式（10）里了。
%% 公式（13-14）
%约束的含义是，持续运行需要受到上下爬坡约束，但启停不受爬坡约束。火电机组爬坡率为2%/min,小时调度体现不出来
sum(P300_Diot(i,:,t))-sum(P300_Diot(i,:,t-1)) <= 0.3*300+(1-u300Git(i,t-1))*300,
P200_it(i,t)-P200_it(i,t-1) <= 0.3*200+(1-u200Git(i,t-1))*200,
P200_it(i,t-1)-P200_it(i,t) <= 0.3*200+(1-u200Git(i,t))*200,
%% 公式（15-16） 火电机组的启停状态标识位置
y200_it(:,2:24)-z200_it(:,2:24)==u200Git(:,2:24)-u200Git(:,1:23), %公式(15)
%% 公式（17-18） 最短运行时间，最短停机时间限值，取20h和10h
%需要注意的是，调度的维度只有24，因此，编写代码的时候，需要避免调用越界24
C = [C,sum(u200Git(i,t:(t+Ton-1))) >= Ton*y200_it(i,t)];C = [C,sum(u300Git(i,t:(t+Ton-1))) >= Ton*y300_it(i,t)];
%% 公式（19） 可再生能源出力约束
0<=P_PVcur,P_PVcur<=E_PV,P_PV+P_PVcur==E_PV,
0<=P_Wind,P_Wind<=E_Wind,
<=P_Windcur,P_Windcur<=E_Wind,
%% 公式（20-24） 抽水蓄能电站功率约束
U_PuG = binvar(4,24,'full');U_PuS = binvar(4,24,'full');
u_PuG = binvar(1,24,'full');u_PuS = binvar(1,24,'full');
%% 公式（25-27） 抽水蓄能电站容量约束
%抽水蓄能电站的综合效率为0.75，取抽水效率为0.9 ,发电效率为0.833 ，nG=1.2
%抽水蓄能可以直接用电量表示容量，也可以用蓄水量代表容量
%仔细看公式（25）可知，文中的抽水蓄能是四个抽水机，四个发电机，共用一个蓄水池部分内容源自网络，侵权联系删除！
欢迎感兴趣的小伙伴关注并私信获取完整版代码，小编会不定期更新高质量的学习资料、文章和程序代码，为您的科研加油助力！
