91号资源-源程序：论文可在知网下载《计及调度经济性的光热电站储热容量配置方法》本人博客有解读

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 for t=1:1:T-1 C=[C,-20<=P(:,t+1)-P(:,t)<=30]; %火电爬坡约束 end for n=1:1:N_i for t=2:1:T C=[C,implies(u(n,t-1)==0&u(n,t)==1,P(n,t)==P_min(n))]; end for t=1:1:T-1 %火电机组启停时刻出力约束 C=[C,implies(u(n,t+1)==0&u(n,t)==1,P(n,t)==P_min(n))]; end end for n=1:1:N_i for t=1:1:T C=[C,U(n,t)==min(P_max(n)-P(n,t),30)]; %火电机组正旋转备用容量约束 C=[C,0.05*load(t)<=U(n,t)]; C=[C,D(n,t)==min(P(n,t)-0,20)]; %火电机组负旋转备用容量约束 C=[C,0.05*load(t)<=D(n,t)]; C=[C,0.05*max(load)<=R(n,t)<=0.1*max(load)] %%%事故备用容量约束 end end %%%%%%光热电站出力约束%%%%%%% C=[C,0<=P_G<=100]; %光热电站出力约束 for t=2:1:T C=[C,-20<=P_G(t)-P_G(t-1)<=30]; %爬坡约束 end %%%%%光热电站储热约束%%%%%% for t=1:1:T C=[C,10<=C_TSr(t)<=1000]; C=[C,0<=P_TSc(t)<=u_TSc(t)*500]; C=[C,0<=P_TSf(t)<=u_TSf(t)*500]; C=[C,u_TSc(t)+u_TSf(t)<=1]; end for t=1:1:T-1 C=[C,C_TSr(t+1)==C_TSr(t)+P_TSc(t)-P_TSf(t)]; %储热装置储热量约束 end C=[C,C_TSr(1)==C_TSr(24)]; %% 目标函数 %%%%火电机组发电成本%%%%% e1=0; e2=0; for n=1:1:N_i for t=1:1:T e1=e1+a(n)*P(n,t)^2+b(n)*P(n,t); end for t=2:1:T e2=e2+u(n,t)*(1-u(n,t-1))*S(n); end end E1=e1+e2; %%%%%光热发电环境效益%%%% E2=0; E2=k_G*sum(P_G); %%%%%旋转备用成本%%%%%% E3=0; for n=1:1:N_i for t=1:1:T E3=E3+alpha*U(n,t)+beta*D(n,t)+gama*R(n,t); end end %%%%%运行维护费用%%%%%% E4=0; E4=sum(k_s*P_SFd+k_Ts*P_TSdf); %%%%储能容量配置成本%%%%% F=0; C_set=max(C_TSr) ; %储能配置容量 F=c*C_set; %%%%%总费用 E=(E1-E2+E3+E4)*365; f=F+E; %% 求解 ops = sdpsettings('solver','cplex', 'verbose', 1);%参数指定程序用cplex求解器 ops.cplex.mip.tolerances.mipgap=0.1; % 0-1 default 0.0001 sol=optimize(C,f,ops); toc %% if sol.problem == 0 disp('succcessful solved'); else disp('error'); yalmiperror(sol.problem) end % [model,recoveryalmip,diagnostic,internalmodel]=export(C,f,ops);%转为cplex模型 % milpt=Cplex('milp for htc'); % milpt.Model.sense='minimize'; % milpt.Model.obj=model.f; % milpt.Model.lb=model.lb; % milpt.Model.ub=model.ub; % milpt.Model.A=[model.Aineq;model.Aeq]; % milpt.Model.lhs=[-inf*ones(size(model.bineq,1),1);model.beq]; % milpt.Model.rhs=[model.bineq;model.beq]; % milpt.Model.ctype=model.ctype; % milpt.writeModel('ab.lp');%输出cplex模型（注意大小写） % milpt.solve();%模型求解 P_cspr=value(P_cspr); P_SFr=value(P_SFr); P_HT=value(P_HT); P_TSc=value(P_TSc); P_TSf=value(P_TSf); C_TSr=value(C_TSr); P_G=value(P_G); P=value(P); u=value(u); U=value(U); D=value(D); R=value(R); C_set=value(C_set); %% 结果输出%%% disp([' ******************************** ']); disp([' ']); disp([' 最 优 热 储 能 容 量 为 ',num2str(C_set),' M W h ']); disp([' ']); disp([' ******************************** ']); figure(1) plot(P(1,:),'-x','linewidth',1.8); hold on plot(P(2,:),'-','linewidth',1.8); hold on plot(P(3,:),'-*','linewidth',1.8); hold on plot(P(4,:),'-^','linewidth',1.8); hold on plot(P(5,:),'-.','linewidth',1.8); hold on plot(P(6,:),'-+','linewidth',1.8); hold on axis([1,24,0,100]); xlabel('时刻/h'); ylabel('有功功率/kW'); grid on legend('机组1出力','机组2出力','机组3出力','机组4出力','机组5出力','机组6出力'); figure(2) t=1:1:24; b=bar(t,C_TSr,'linewidth',1); b.FaceColor=[0.3 0.2 0.5]; ax = gca; ax.YColor = 'k'; axis([0.5,24.5,0,1000]); xlabel('时刻/h'); ylabel('热储能储热量/MWh'); % 给左y轴添加轴标签 yyaxis right; % 激活右边的轴 hold on plot(P_TSc-P_TSf,'-^g','linewidth',2); ax = gca; ax.YColor = 'k'; axis([1,24,-300,300]); ylabel('热功率/MW'); % 给右y轴添加轴标签 legend('储热量','储热功率'); grid on figure(3) F=zeros(2,24); F(1,:)=P_SFr; F(2,:)=P_HT; AA=bar(F(:,:)','stack'); set(AA(1),'FaceColor',[0.1 0.5 0.9]); set(AA(2),'FaceColor',[0.9 0.1 0.5]); hold on plot(t,P_cspr,'-*b','linewidth',2); axis([0.5,24.5,0,600]); xlabel('时刻/h'); ylabel('热功率/MW'); legend('集热装置直接发电的热功率','集热装置供给储热系统的热功率','集热装置吸收热功率'); grid on figure(4) plot(t,load,'-*r','linewidth',1.8); hold on plot(t,load-P_G,'--xb','linewidth',1.8); hold on plot(t,P_G,':.k','linewidth',1.8); axis([0.5,24.5,0,250]); xlabel('时刻/h'); ylabel('电功率/MW'); legend('原始负荷','光热电站接入后的等效负荷','光热电站出力'); grid on figure(5) subplot(2,1,1) plot(U(1,:),'-x','linewidth',1.8); hold on plot(U(2,:),'-','linewidth',1.8); hold on plot(U(3,:),'-*','linewidth',1.8); hold on plot(U(4,:),'-^','linewidth',1.8); hold on plot(U(5,:),'-.','linewidth',1.8); hold on plot(U(6,:),'-+','linewidth',1.8); axis([1,24,5,35]); xlabel('时刻/h'); ylabel('向上旋转备用容量/MW'); grid on legend('机组1','机组2','机组3','机组4','机组5','机组6'); subplot(2,1,2) plot(D(1,:),'-x','linewidth',1.8); hold on plot(D(2,:),'-','linewidth',1.8); hold on plot(D(3,:),'-*','linewidth',1.8); hold on plot(D(4,:),'-^','linewidth',1.8); hold on plot(D(5,:),'-.','linewidth',1.8); hold on plot(D(6,:),'-+','linewidth',1.8); axis([1,24,5,25]); xlabel('时刻/h'); ylabel('向下旋转备用容量/MW'); grid on