第十五届电工杯5月26号就要开始啦,今天给大家回顾第十四届全国大学生电工数学建模竞赛A题,主要从赛题重述和问题分析与代码实战展开。第十五届全国大学生电工数学建模竞赛已经开始报名了哦,后续我也会分享对应的建模思路哦,大家记得多多关注。
一、赛题重述
“碳中和”目标驱动下未来电力系统必将是高比例可再生能源电力系统,可再生能源输出功率强随机波动性导致系统运行中功率实时平衡困难;储能被认为 是保障系统功率实时平衡的有效手段,由于储能成本相对昂贵,利用储能平衡系 统功率将增加系统运行成本;下面以高比例风电电力系统为例,探究“供给侧” 低碳化转型对电力系统运行经济性、可靠性影响。待研究系统包含火电、风电、储能和负荷,火电机组 3 台、装机容量 1050MW;某日风电、负荷归一化功率(1.0p.u.风电对应其装机容量,1.0p.u.负荷对应最大 负荷功率)数据见附件 1,风电渗透率(最大风电功率与最大负荷功率之比)递 增可能造成系统弃风、失负荷,影响系统功率平衡。定义:系统单位供电成本=系统发电总成本/系统总负荷电量,发电总成本= 火电成本+风电成本+储能成本+弃风损失+失负荷损失,其中:火电成本包括运行成本、碳捕集成本,其中火电运行成本由运行维护成 本和发电煤耗成本构成,,P 为机组出力/MW;运行维护成本按照 0.5 倍煤耗 成本考虑,碳捕集成本取决于碳排放量及碳捕集单价,火电机组相关参 数如附表 1 所示,电煤价格为 700 元/t。风电成本仅考虑运维成本,相关参数如附表 2。储能成本由投资成本、运维成本构成,相关参数如附表 3。注:在计算 每天成本时,需将投资成本平摊至每天,即平均每天投资成本=总投资 成本/运行年限/365 天。弃风损失按 0.3 元/kWh 计算,失负荷损失按 8 元/kWh 计算。基本题:假设系统日负荷功率最大值 900MW,单位碳捕集成本分别为 0 元/t、60 元 /t、80 元/t、100 元/t,摄动风电渗透率,分析计算以下问题:
1. 无风电接入,火电以最小成本运行,绘制机组日发电计划曲线,计算系 统单位供电成本,将结果填于表 1 相应栏(保留三位有效数字)。
2. 风电装机 300MW、替代机组 3 时,系统功率平衡发生什么变化?弃风电 量多少?在此场景下,为减少弃风又不失负荷,风电接入装机容量可以降低多少?
3. 风电装机 600MW、替代机组 2 时,系统功率平衡又发生什么变化?在此 场景下,为不失负荷,风电接入容量可增加多少?
4. 针对上述 2-3 风电替代场景,考虑上述 4 种碳捕捉成本,系统按最低发 电成本供电,计算系统单位供电成本,并将相关计算结果填入表 2、3。
5. 风电装机 900MW、替代机组 2、3 时,失负荷电量多少?为不失负荷,需 要配置的最小储能容量是多少(储能充放电效率 90%)?考虑储能成本、单位碳 捕捉成本(取 60 元/t),此时系统单位供电成本多大?
6. 当负荷功率不变,试分析风电替代容量递增給系统可靠供电带来哪些挑 战?为保障可靠供电,系统单位供电成本发生了什么变化?结合上述计算结果进 行定量分析。
二、问题分析
1、第一题的思路比较简单,系统没有风电接入,仅考虑火力发电,且当日负荷最大值小于三个火力机组最大技术出力之和,则系统成本没有储能成本、弃风成本和失负荷成本。所以以火电成本最小为目标函数,考虑系统功率平衡约束。对模型求解即可得到日发电计划曲线。求解方法可采用粒子群算法或采用CPLEX求解器求解,本文采用CPLEX求解器求解.
根据对问题1的分析,仅考虑火力发电,没有储能成本、弃风成本和失负荷成本,则发电总成本=火电成本=火电运行成本+碳捕集成本,其中火电运行成本=运行维护成+发电煤耗成本。建立以火电成本最小为优化目标,目标函数如下式所示:
考虑系统功率平衡约束和机组上、下限约束
采用Yalmip+Cplex对以上模型进行求解 ,代码如下
clc
clear
data = xlsread('附件1.xlsx');
T = 96;
N = 3; %火电机组的个数
%%提取负荷数据
pt = data(:,2);
PT_MAX = 900;
PT = PT_MAX*pt;
figure(1)
plot(PT)
%%火电机组最大最小出力
WT_MAX = [600 300 150];
WT_MIN = [180 90 45];
%% 发电耗能相关系数
a = [0.226 0 0;0 0.588 0;0 0 0.785];
b = [30.42 65.12 139.6];
c = [786.80 451.32 1049.50];
%%电煤单价
PRICE_WT = 700/1000;%%转换为1Kg电煤单价
%%碳捕集单价
d = 0;%依次修改为0、60、80、100
PRICE_C = d/1000;
%碳排放量
CARB = [0.72 0.75 0.79];
P = sdpvar(3,T);
Constraints = [];
%%功率平衡约束
for i = 1:T
Constraints = [Constraints; sum(P(:,i)) == PT(i)];
end
for i = 1:T
Constraints = [Constraints; WT_MIN(1) <= P(1,i) <= WT_MAX(1);
WT_MIN(2) <= P(2,i) <= WT_MAX(2);
WT_MIN(3) <= P(3,i) <= WT_MAX(3)];
end
%%目标函数
Objective = 0;
for i = 1:T
Objective = Objective +1.5*0.25*PRICE_WT*(P(:,i)'*a*P(:,i)+b*P(:,i)+c*ones(3,1))+1000*0.25*PRICE_C*CARB*P(:,i);
end
ops = sdpsettings('solver','cplex');%设定求解器
result = optimize(Constraints,Objective,ops);%求解
t = 1:96;
p = value(P);%火电机组出力
Objective = value(Objective);%目标函数最优解
figure(2)
plot(t,p(1,:),'r-','Linewidth',1);hold on;
plot(t,p(2,:),'k-','Linewidth',1);hold on;
plot(t,p(3,:),'b-','Linewidth',1);hold on;
plot(t,ones(length(t)),'k--');hold on;
legend('火电机组1','火电机组2','火电机组3');
P_c1=zeros(3,1);
for i = 1:3
P_c1(i,1) = d*0.25*CARB(i)*sum(p(i,:))/10000;
end
P_c=sum(P_c1) %碳捕捉成本,万元
P_om1=zeros(3,1);
for j = 1:3
for i = 1:T
P_om1(j) = P_om1(j) +1.5*0.25*PRICE_WT*(a(j,j)*p(j,i)^2+b(j)*p(j,i)+c(j)*ones(1,1));
end
end
P_om=sum(P_om1)/10000 %火电运行成本
P_gd=(P_c+P_om)/(sum(PT)/4)*10 %单位供电成本,元/KWh
P_s=P_om+P_c %总发电成本
2、对于问题2,风电装机300MW替代机组3,则要考虑风电成本、弃风成本和失负荷成本。当总发电功率超过负荷功率,形成弃风;总发电功率小于负荷功率时,形成失负荷。以最小总成本为目标函数,系统功率平衡为约束,使用MATLAB 求解模型。其中风功率数据由风功率归一化数据根据风机装机容量反归一化得到。而对于减少弃风又不失负荷, 建立以系统总成本最小,包括弃风成本最小为目标函数,采用差分进化算法优化风机接入容量,使用MATLAB 求解结果。
(1)当风电装机300MW替代机组3时,可能会造成系统弃风、失负荷、从而影响系统的功率平衡。此时系统总成本=火电机组1成本+火电机组2成本+风电成本+弃风损失+失负荷损失,建立以系统总成本最小为优化目标,目标函数如下:
考虑系统功率平衡约束,机组出力上、下限约束,弃风电量约束和切负荷约束
考虑不同单位碳捕集成本情况,利用Matlab采用Yalmip+Cplex对以上模型进行求解,代码如下:
clc
clear
close all
data = xlsread('附件1.xlsx');
T = 96;
N = 2; %火电机组的个数
%%提取负荷数据、风功率数据
pt = data(:,2);
pw = data(:,3);
PT_MAX = 900;
PW_MAX = 300;
PT = PT_MAX*pt;
PW = PW_MAX*pw;
figure(1)
plot(PT)
hold on
plot(PW)
%%火电机组最大最小出力
WT_MAX = [600 300];
WT_MIN = [180 90];
%% 发电耗能相关系数
a = [0.226 0;0 0.588];
b = [30.42 65.12];
c = [786.80 451.32];
%%电煤单价
PRICE_WT = 700/1000;%%转换为1Kg电煤单价
%%碳捕集单价
d = 0;%依次修改为0、60、80、100
PRICE_C = d/1000;
%碳排放量
CARB = [0.72 0.75];
P = sdpvar(2,T);
pw_loss = sdpvar(1,T);
Constraints = [];
%%功率平衡约束
for i = 1:T
Constraints = [Constraints; sum(P(:,i))+ PW(i) - pw_loss(i) == PT(i)];
end
for i = 1:T
Constraints = [Constraints; WT_MIN(1) <= P(1,i) <= WT_MAX(1);
WT_MIN(2) <= P(2,i) <= WT_MAX(2);
PW(i) - pw_loss(i) >= 0;
pw_loss(i) >= 0];
end
%目标函数
Objective = 0;
for i = 1:T
Objective = Objective +1.5*0.25*PRICE_WT*(P(:,i)'*a*P(:,i)+b*P(:,i)+c*ones(2,1))+1000*0.25*PRICE_C*CARB*P(:,i)+0.3*0.25*1000*pw_loss(i)+0.045*0.25*1000*(PW(i)-pw_loss(i));
end
ops = sdpsettings('solver','cplex');
result = optimize(Constraints,Objective,ops);
t = 1:96;
p = value(P);%火电机组1、2出力
p2 = value(pw_loss);%弃风
objective = value(Objective);%目标函数最优解
PW_REAL = PW - p2';%风机实际出力
figure(2)
plot(t,p(1,:),'r-','Linewidth',1);hold on;
plot(t,p(2,:),'k-','Linewidth',1);hold on;
plot(t,p2,'b-','Linewidth',1);hold on;
plot(t,PW_REAL,'g-','Linewidth',1);hold on;
legend('火电机组1','火电机组2','弃风量','风机实际出力');
P_c1=zeros(2,1);
for i = 1:2
P_c1(i,1) = d*0.25*CARB(i)*sum(p(i,:))/10000;
end
P_c=sum(P_c1) %碳捕捉成本,万元
P_om1=zeros(2,1);
for j = 1:2
for i = 1:T
P_om1(j) = P_om1(j) +1.5*0.25*PRICE_WT*(a(j,j)*p(j,i)^2+b(j)*p(j,i)+c(j)*ones(1,1));
end
end
P_om=sum(P_om1)/10000; %火电运行成本
P_w=sum(p2)*1000/4*0.3+sum(PW'-p2)*1000/4*0.045; %风电成本+弃风成本
P_total=P_c+P_om+P_w/10000 %总成本
P_gd=(P_c+P_om)/(sum(PT)/4) %单位供电成本,万元/MWh
p_qifeng=sum(p2)*1000/4 %弃风量
(2)在此场景下,为了保证减少弃风又不失负荷,本文考虑以系统总成本和弃风惩罚成本之和最小为目标建立相应的目标函数如下式所示:
考虑系统功率平衡约束,机组出力上、下限约束,弃风电量约束和切负荷约束
采用差分进化算法优化风电装机容量,差分进化算法如具体流程如下图,适应度函数为系统总成本和弃风惩罚成本之和。
clc
clear
close all
data = xlsread('附件1.xlsx');
pt = data(:,2);
pw = data(:,3);
PT_MAX = 900;
% PW_MAX = 300-40.7892000000111;
PT = PT_MAX*pt;
T = 96;
WT_MAX = [600 300];
WT_MIN = [180 90];
nPop=8; % 种群规模 Population Size
MaxIt=20; % 最大迭代次数Maximum Number of Iterations
nVar= 1; % 自变量维数,此例需要优化两个参数c和g Number of Decision Variables
VarSize= 1; % 决策变量矩阵大小 Decision Variables Matrix Size
beta_min=0.2; % 缩放因子下界 Lower Bound of Scaling Factor
beta_max=0.8; % 缩放因子上界 Upper Bound of Scaling Factor
pCR=0.3; % 交叉概率 Crossover Probability
lb= 0; % 参数取值下界
ub= 300; % 参数取值上界
%% 初始化 Initialization
empty_individual.Position=[]; % 种群初始化
empty_individual.Cost=[]; % 种群目标函数值初始化
BestSol.Cost=inf; % 最优值初始化
pop=repmat(empty_individual,nPop,1); % 将保存种群信息的结构体扩展为结构体矩阵,行数等于种群大小
for i=1:nPop % 遍历每个个体
pop(i).Position=init_individual(lb,ub,nVar,1); % 随机初始化个体
pop(i).Cost=fitness( pop(i).Position,pt,pw,PT_MAX); % 计算个体目标函数值
if pop(i).Cost<BestSol.Cost % 如果个体目标函数值优于当前最优值
BestSol=pop(i); % 更新最优值
end
end
BestCost=zeros(MaxIt,1); % 初始化迭代最优值
%% 主循环 DE Main Loop
for it=1:MaxIt
for i=1:nPop % 遍历每个个体
x=pop(i).Position; % 提取个体位置
% 随机选择三个个体以备变异使用
A=randperm(nPop); % 个体顺序重新随机排列
A(A==i)=[]; % 当前个体所排位置腾空(产生变异中间体时当前个体不参与)
a=A(1);
b=A(2);
c=A(3);
% 变异操作 Mutation
% beta=unifrnd(beta_min,beta_max,VarSize); % 随机产生缩放因子
beta = (beta_max-(beta_max-beta_min)*(i/MaxIt))*[1,1];
y=round(pop(a).Position+beta.*(pop(b).Position-pop(c).Position)); % 产生中间体
% 防止中间体越界
y=max(y,lb);
y=min(y,ub);
% 交叉操作 Crossover
z=zeros(size(x)); % 初始化一个新个体
j0=randi([1,numel(x)]); % 产生一个伪随机数,即选取待交换维度编号???
for j=1:numel(x) % 遍历每个维度
if j==j0 || rand<=pCR % 如果当前维度是待交换维度或者随机概率小于交叉概率
z(j)=y(j); % 新个体当前维度值等于中间体对应维度值
else
z(j)=x(j); % 新个体当前维度值等于当前个体对应维
end
end
NewSol.Position=z; % 交叉操作之后得到新个体
NewSol.Cost=fitness(NewSol.Position,pt,pw,PT_MAX); % 新个体目标函数值
if NewSol.Cost<pop(i).Cost % 如果新个体优于当前个体
pop(i)=NewSol; % 更新当前个体
if pop(i).Cost<BestSol.Cost % 如果当前个体(更新后的)优于最优个体
BestSol=pop(i); % 更新最优个体
end
end
end
% 保存当前迭代最优个体函数值 Update Best Cost
BestCost(it)=BestSol.Cost;
end
bestc=BestSol.Position(1);
% bestg=BestSol.Position(2);
%% 图示寻优过程
plot(BestCost);
xlabel('Iteration');
ylabel('Best Val');
PW_MAX = bestc;
PT = PT_MAX*pt;
PW = PW_MAX*pw;
T = 96;
%% 发电耗能相关系数
a = [0.226 0;0 0.588];
b = [30.42 65.12];
c = [786.80 451.32];
%%电煤单价
PRICE_WT = 700/1000;%%转换为1Kg电煤单价
%%碳捕集单价
d = 0;
PRICE_C = d/1000;
%碳排放量
CARB = [0.72 0.75];
P = sdpvar(2,T);
pw_loss = sdpvar(1,T);
p_loss = sdpvar(1,T);
Constraints = [];
%%功率平衡约束
for i = 1:T
Constraints = [Constraints; sum(P(:,i))+ PW(i) - pw_loss(i) == PT(i)];
end
for i = 1:T
Constraints = [Constraints; WT_MIN(1) <= P(1,i) <= WT_MAX(1);
WT_MIN(2) <= P(2,i) <= WT_MAX(2);
PW(i) - pw_loss(i) >= 0;
pw_loss(i) >= 0;
];
end
Objective = 0;
for i = 1:T
Objective = Objective +1.5*0.25*PRICE_WT*(P(:,i)'*a*P(:,i)+b*P(:,i)+c*ones(2,1))+1000*0.25*PRICE_C*CARB*P(:,i)+1000*0.3*0.25*pw_loss(i)+1000*0.045*0.25*PW(i);
end
ops = sdpsettings('solver','cplex');
result = optimize(Constraints,Objective,ops);
Target = value(Objective);
p = value(P);
p2 = value(pw_loss);
ss = sum(p2)
p3 = value(p_loss);
PW_REAL = PW - p2';
P_REAL = PT - p3';
figure(2)
t = 1:96;
plot(t,p(1,:),'r-','Linewidth',1);hold on;
plot(t,p(2,:),'k-','Linewidth',1);hold on;
plot(t,p2,'b-','Linewidth',1);hold on;
plot(t,PW_REAL,'g-','Linewidth',1);hold on;
plot(t,P_REAL,'y-','Linewidth',1);hold on;
legend('火电机组1','火电机组3','弃风量','风机实际出力','实际负荷');
差分进化算法迭代优化曲线如下图所示。
今天的分享就到这里了,后续想了解智能算法、机器学习、深度学习和信号处理相关理论的可以后台私信哦~希望大家多多转发点赞加收藏,你们的支持就是我源源不断的创作动力!