风力发电系统的随机调度研究（matlab代码实现）

👨‍🎓 个人主页： 研学社的博客

💥 💥 💞 💞 欢迎来到本博客 ❤️ ❤️ 💥 💥


🏆 博主优势： 🌞 🌞 🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。


⛳ 座右铭：行百里者，半于九十。

📋 📋 📋 本文目录如下： 🎁 🎁 🎁
目录
💥1 概述
📚2 运行结果
🎉3 参考文献
🌈4 Matlab代码实现

💥1 概述

文献来源：

时域调度模拟是预测风力发电系统运行成本最有效的工具，因为它可以表示限制热电厂、储能和需求侧措施平衡行动的跨时间约束。高风侵度要求及时调度决策能够反映风馈的不确定性，因此使用随机单元承诺(SUC)和滚动规划生成调度决策是可取的。然而，计算负担可能使这种方法在长时间的模拟中不切实际。我们提出了SUC问题的一个有效公式，设计用于单母线电力系统调度仿真。

与传统的SUC技术不同，该公式采用了基于分位数的情景树结构，避免了对外生运营储备的需求。在长达一年的大型系统模拟中，我们比较了各种树拓扑的性能。与固定分位数确定性方法相比，简单的基于分位数的树在统计上具有显著的成本改进，并且与基于蒙特卡罗生成场景的树相比更优。

任何国家的政府都致力于大幅降低其电力部门的碳强度。脱碳的电力系统可能具有以下特点:大的、相对不灵活的核基载;间歇性和部分不可预测的可再生能源的高渗透率，在许多国家将由风能主导;化石燃料车队的载客率比目前低得多。因此，大规模风力集成的一个核心挑战是利用灵活性降低的热机组吸收风力发电的能力。

衡量系统吸收风的能力的一个指标是风的不确定性成本，即由于短期风的不确定性而运行系统的额外成本。简单的统计分析无法告诉我们这一成本，因为平衡风力间歇性的可调度要素受制于限制其响应能力的跨时间约束。这些限制的例子是热发电机的启动时间，以及对可存储能量的限制。为了理解这些限制对系统吸收风的能力的影响，我们必须在时域内模拟整个系统的调度，在一定程度上考虑跨时间约束。如[1]所述，在许多研究中使用时域模拟来评估集成风的成本。

📚2 运行结果

部分代码：

%%

clearvars

clc

%% Inputs

load('Data\InputData.mat')

load('Data\SystemParameters.mat')

% The units are GW, k/GW and k to keep the range of variables closer

%% Define the optimisation model:

Load_balance = [];

Aux_constraints = [];

Constraints_Gen_limits = [];

Constraints_Ramp_limits = [];

Cost_node = [];

scenarios = size(Wind_scenarios,1);

% Create some matrices of DVs:

Wind_curtailed = sdpvar(scenarios,time_steps);

Load_curtailed = sdpvar(scenarios,time_steps);

% And some matrices of constraints:

Aux_constraints = [Aux_constraints,...

0 <= Wind_curtailed <= Wind_scenarios,...

0 <= Load_curtailed <= ones(scenarios,1)*Demand];

tic % Start time counter

for k=1:time_steps

%% Define DVs

x{k} = sdpvar(scenarios,num_Clusters);

% x{k} = [x1 x2 x3... (scenario 1)

% x1 x2 x3... (scenario 2)]; % Power generated by each Cluster in each scenario

y{k} = intvar(scenarios,num_Clusters);

% y{k} = [y1 y2 y3... (scenario 1)

% y1 y2 y3... (scenario 2)]; % Commitment decision for each Cluster in each scenario

% FORCE y{k} to be lower than number of CCGTs and number of OCGTs

Aux_constraints = [Aux_constraints,...

zeros(scenarios,1)*num_Gen <= y{k} <= ones(scenarios,1)*num_Gen];

% NO INITIAL COMMITMENT

% if k==1

% Aux_constraints = [Aux_constraints,...

% y{n,k}==Initial_commitment];

% end

% Fix commitment decision of inflexible_gen to be same as in the first

% scenario in this time-step "k"

Aux_constraints = [Aux_constraints,...

(ones(scenarios,1)*inflexible).*y{k} == ones(scenarios,1)*(inflexible.*y{k}(1,:))];

% Define DVs for generators started up

startup_DV{k} = sdpvar(scenarios,num_Clusters);

if k==1

Aux_constraints = [Aux_constraints,...

startup_DV{k} == y{k}];

else

% These 2 constraints model the startup, it's easy to verify

% that they indeed do:

Aux_constraints = [Aux_constraints,...

startup_DV{k} >= y{k} - y{k-1},...

startup_DV{k} >= zeros(scenarios,num_Clusters)];

end

%% Constraints:

Load_balance = [Load_balance,...

sum(x{k},2)+Wind_scenarios(:,k)-Wind_curtailed(:,k) == ones(scenarios,1)*Demand(k)-Load_curtailed(:,k)];

% Generation limits:

for i=1:num_Clusters

Constraints_Gen_limits = [Constraints_Gen_limits,...

y{k}.*(ones(scenarios,1)*Gen_limits(:,1)') <= x{k} <= y{k}.*(ones(scenarios,1)*Gen_limits(:,2)')];

end

% % Ramp limits:

% if k>1

%

% Constraints_Ramp_limits = [Constraints_Ramp_limits,...

% -tau*(ones(scenarios,1)*Ramp_limits(:,1)') <= x{k}-x{k-1} <= tau*(ones(scenarios,1)*Ramp_limits(:,2)')];

%

% end

%% Finally, define the nodal costs:

Cost_node = horzcat(Cost_node,...

sum((ones(scenarios,1)*stc').*startup_DV{k},2)... % Startup costs

+ tau*(sum((ones(scenarios,1)*NLHR').*y{k},2)... % fixed generation costs

+ sum((ones(scenarios,1)*HRS').*x{k},2)... % variable generation costs

+ VOLL*Load_curtailed(:,k))); % load shed cost

end

time.BuildingModel = toc;

tic

Objective = sum(sum(prob_nodes.*Cost_node));

Constraints = [Aux_constraints,...

Load_balance,...

Constraints_Gen_limits

Constraints_Ramp_limits];

options = sdpsettings('solver','gurobi','gurobi.MIPGap',0.1e-2);

solution = optimize(Constraints,Objective,options)

time.Optimisation = toc;

%% Analyse results

tic

sol.Wind_curtailed = value(Wind_curtailed);

sol.Load_curtailed = value(Load_curtailed);

sol.Cost = [];

sol.Cost_StartUp = [];

sol.Cost_NLHR = [];

sol.Cost_HRS = [];

sol.Cost_VOLL = [];

sol.Check_Load_balance(scenarios,time_steps) = 0;

for k=1:time_steps

sol.x{k} = value(x{k});

sol.y{k} = value(y{k});

sol.startup_DV{k} = value(startup_DV{k});

sol.Cost = horzcat(sol.Cost,...

sum((ones(scenarios,1)*stc').*sol.startup_DV{k},2)... % Startup costs

+ tau*(sum((ones(scenarios,1)*NLHR').*sol.y{k},2)... % fixed generation costs

+ sum((ones(scenarios,1)*HRS').*sol.x{k},2)... % variable generation costs

+ VOLL*sol.Load_curtailed(:,k)));

sol.Cost_StartUp = horzcat(sol.Cost_StartUp,...

sum((ones(scenarios,1)*stc').*sol.startup_DV{k},2));

sol.Cost_NLHR = horzcat(sol.Cost_NLHR,...

tau*(sum((ones(scenarios,1)*NLHR').*sol.y{k},2)));

sol.Cost_HRS = horzcat(sol.Cost_HRS,...

tau*(sum((ones(scenarios,1)*HRS').*sol.x{k},2)));

sol.Cost_VOLL = horzcat(sol.Cost_VOLL,...

tau*(VOLL*sol.Load_curtailed(:,k)));

sol.Check_Load_balance(:,k) = sum(sol.x{k},2)+Wind_scenarios(:,k)-sol.Wind_curtailed(:,k)-(ones(scenarios,1)*Demand(k)-sol.Load_curtailed(:,k));

end

time.GetSolution = toc;

clearvars -except Wind_scenarios Demand options sol solution time

save('Solution.mat')

🎉3 参考文献

部分理论来源于网络，如有侵权请联系删除。

🌈4 Matlab代码实现