1 主要内容

该程序复现《计及调度经济性的光热电站储热容量配置方法》模型，综合考虑火电机组发电成本、光热发电并网消纳的环境效益和运行维护成本、系统旋转备用成本等调度经济性因素，求解光热电站最佳储热容量，并在此基础上探究储热装置配置成本与调度经济性的平衡点，确定光热电站储热容量配置。基于matlab+cplex求解器进行优化求解，通过IEEE-30节点系统验证所提方法的可行性与有效性。

目标函数

调度经济性

综合成本函数

光热电站能量传递过程

2 部分程序

clc
clear
tic
%% 参数定义
T=24;   %调度时刻
yita_c = 0.9; %热储能充热效率
yita_f = 0.9; %热储能放热效率
yita_d = 0.4; %储热系统热电转换效率
N_i=6; %火电机组台数
a=[0.2625 1.224 4.375 0.5838 1.75 1.75];%机组燃料成本系数
b=[140 122.5 70 227.5 210 210];
S=[150 141 121 133 111 141]; %启停成本
P_max=[100 80 50 35 30 40]; %火电机组出力上限
P_min=[53 41 26 18 16 21]; %火电机组出力下限
k_G=230; %环境效益系数230元/MW
alpha=157; %正旋转备用系数
beta=134;  %负旋转备用系数
gama=132;  %事故旋转备用系数
k_s=120; %集热装置供热发电的成本系数
k_Ts=100; %储热装置供热发电的成本系数
c=420000; %单位储热容量配置成本
idg=[1,2,22,27,23,13];
is=22;
mpc=case30;
l=sum(mpc.bus(:,3));
A=[0.95 0.965 0.925 0.945 0.965 0.976 0.953 0.942 0.934 0.943 0.934 0.98 0.96 0.93 0.94 0.967 0.935 0.95 0.98 0.99 0.97 0.94 0.95 0.97];
load = A*l;
loadi = repmat(A,30,1).*repmat(mpc.bus(:,3),1,24);
% load=[220 230 230 230 198 220 208 250 280 303 300 259 255 230 235 253 290 308 307 280 274 260 252 238];
%% 变量定义
P_cspr=[10.05 10.05 10.05 10.05 10.89 10.05 10.89 10.05 84.59 386.10 427.97 383.58 427.14 505.03 355.11 11.73 87.94 10.89 10.05 10.05 10.05 10.05 10.05 10.05]; %CSP集热装置吸收热功率
P_SFr=sdpvar(1,T,'full'); %集热装置直接发电的热功率
P_HT=sdpvar(1,T,'full');  %集热装置供给储热系统的热功率
P_SFd=sdpvar(1,T,'full');  %集热装置输出电功率
P_TH=sdpvar(1,T,'full');  %储热装置传递至传热流体的热功率
P_TSc=sdpvar(1,T,'full'); %储热装置的储热功率
P_TSf=sdpvar(1,T,'full'); %储热装置的放热功率
u_TSc=binvar(1,T,'full'); %储热装置的储热状态
u_TSf=binvar(1,T,'full'); %储热装置的放热状态
P_TSdf=sdpvar(1,T,'full'); %储热装置输出电功率
C_TSr=sdpvar(1,T,'full'); %储热量
P_G=sdpvar(1,T,'full'); %光热电站的输出功率
P=sdpvar(N_i,T,'full'); %火电机组功率
u=binvar(N_i,T,'full'); %火电机组运行状态
U=sdpvar(N_i,T,'full'); %火电机组正旋转备用容量
D=sdpvar(N_i,T,'full'); %火电机组负旋转备用容量
R=sdpvar(N_i,T,'full'); %火电机组事故旋转备用容量
x_theta=sdpvar(30,T,'full');
%% 光热电站运行模型
C=[];
%集热装置运行特性
C=[C,P_SFr+P_HT == P_cspr]; 
% C=[C,P_SFr>=0];
% C=[C,P_HT>=0];
%储热装置运行特性
C=[C,P_TSc==yita_c*P_HT, P_TSf==P_TH/yita_f];
%光热电站发电
C=[C,P_G==yita_d*(P_cspr-P_TSc/yita_c+yita_f*P_TSf)];
C=[C,P_SFd==yita_d*P_SFr];
C=[C,P_TSdf==(1-0.03)*yita_d*P_TSf];
C=[C,0<=P_SFd<=500,0<=P_TSdf<=500];
%% 约束条件
%功率平衡约束
% C=[C,sum(P)+P_G==load]; 
C = [C,getConsEQ2(P,P_G,loadi,case30, x_theta,idg,is)];
% 火电机组约束
for n=1:1:N_i
    for t=1:1:T
  C=[C,0<=P(n,t)<=P_max(n)]; %火电机组出力约束
    end
end