matlab-考虑新能源消纳的火电深度调峰—电网调度-gurobi-混合整数线性规划

matlab-考虑新能源消纳的火电深度调峰—电网调度-gurobi-混合整数线性规划

以下是一个基于 MATLAB + Gurobi 的简单混合整数线性规划（MILP）模型，用于模拟考虑新能源消纳的火电深度调峰与电网调度问题。这个模型的目标是：在满足负荷需求和系统约束的前提下，最小化运行成本，并尽可能多地消纳风电/光伏等新能源。

---

🧮 模型说明

✅ 决策变量：

• 火电机组出力 $ P_{t,i} $
• 新能源机组出力 $ W_t $
• 负荷不平衡量（弃风弃光）$ S_t $

✅ 目标函数：

• 最小化火电运行成本 + 弃风惩罚成本

✅ 主要约束：

1. 功率平衡约束
2. 火电机组上下限、爬坡约束
3. 新能源最大可发电量限制
4. 非负性约束

---

📦 MATLAB + Gurobi 代码示例

% 清除工作区
clear; clc;

% 参数设置
T = 24;            % 时间段数量（小时）
N = 3;             % 火电机组数量
C_fuel = [200; 250; 300];   % 火电单位燃料成本 (元/MWh)
P_min = [50; 60; 70];       % 火电最小出力 (MW)
P_max = [200; 220; 250];    % 火电最大出力 (MW)
RU = [50; 50; 50];          % 爬坡上限 (MW/h)
RD = [50; 50; 50];          % 爬坡下限 (MW/h)

% 新能源预测功率（如风电）
W_forecast = [0, 50, 80, 100, 120, 110, 90, 70, ...
              60, 50, 40, 30, 30, 40, 50, 70, ...
              90, 110, 120, 130, 120, 100, 80, 60]; % MW

% 负荷预测
Load = [200, 220, 230, 250, 260, 280, 270, 250, ...
        240, 230, 220, 210, 200, 210, 220, 240, ...
        260, 280, 300, 310, 300, 280, 260, 240]; % MW

% 创建Gurobi模型
model.modelname = 'Power Dispatch with Renewable Integration';
model.modelsense = 'min';

% 变量定义
% x(i,t): 火电机组i在t时刻的出力
% s(t): 弃风量
% 所有变量 >= 0

num_vars = N*T + T;
counter = 1;

for t = 1:T
    for i = 1:N
        varindex_x(i,t) = counter;
        counter = counter + 1;
    end
    varindex_s(t) = counter;
    counter = counter + 1;
end

% 初始化目标函数系数
obj = zeros(num_vars, 1);

% 设置目标函数：火电成本 + 弃风惩罚
penalty_wind = 1000; % 弃风惩罚成本
for t = 1:T
    for i = 1:N
        idx = varindex_x(i,t);
        obj(idx) = C_fuel(i);
    end
    idx_s = varindex_s(t);
    obj(idx_s) = penalty_wind;
end

model.obj = obj;

% 添加变量：x(i,t), s(t)
model.lb = zeros(num_vars, 1);
model.ub = Inf(num_vars, 1); % 上限后面通过约束处理

model.vtype = repmat('C', num_vars, 1); % 全部为连续变量

% 构造约束矩阵 A*x <= b
% 包括功率平衡、机组出力限制、爬坡限制等

% 使用 linprog 或者 sparse 构建 A 和 b
A = [];
b = [];

% 1. 功率平衡约束： sum(x(i,t)) + W_forecast(t) - s(t) == Load(t)
for t = 1:T
    row = zeros(num_vars, 1);
    for i = 1:N
        idx = varindex_x(i,t);
        row(idx) = 1;
    end
    row(varindex_s(t)) = -1;
    A = [A; row'];
    b = [b; Load(t) - W_forecast(t)];
end

% 2. 火电机组出力上下限
for t = 1:T
    for i = 1:N
        row = zeros(num_vars, 1);
        idx = varindex_x(i,t);
        row(idx) = 1;
        A = [A; row'];
        b = [b; P_max(i)];

        row = zeros(num_vars, 1);
        row(idx) = -1;
        A = [A; row'];
        b = [b; -P_min(i)];
    end
end

% 3. 爬坡约束
for t = 2:T
    for i = 1:N
        row = zeros(num_vars, 1);
        idx_t = varindex_x(i,t);
        idx_t1 = varindex_x(i,t-1);
        row(idx_t) = 1;
        row(idx_t1) = -1;
        A = [A; row'];
        b = [b; RU(i)];

        row = zeros(num_vars, 1);
        row(idx_t) = -1;
        row(idx_t1) = 1;
        A = [A; row'];
        b = [b; RD(i)];
    end
end

% 添加不等式约束到模型
model.A = sparse(A);
model.rhs = b;
model.sense = repmat('=', size(A,1),1); % 大部分是等式约束

% 求解
params.outputflag = 1;
result = gurobi(model, params);

% 输出结果
if strcmp(result.status, 'OPTIMAL')
    disp('找到最优解！');

    % 提取火电出力
    x = zeros(N,T);
    for t = 1:T
        for i = 1:N
            x(i,t) = result.x(varindex_x(i,t));
        end
    end

    % 提取弃风量
    s = zeros(1,T);
    for t = 1:T
        s(t) = result.x(varindex_s(t));
    end

    % 绘图展示结果
    figure;
    plot(1:T, x', 'LineWidth', 1.5);
    hold on;
    plot(1:T, Load, '--k', 'LineWidth', 1.5);
    plot(1:T, W_forecast, '-g', 'LineWidth', 1.5);
    plot(1:T, s, '-r', 'LineWidth', 1.5);
    legend('火电1', '火电2', '火电3', '负荷', '风电预测', '弃风');
    xlabel('时间 (小时)');
    ylabel('功率 (MW)');
    title('电网调度与新能源消纳优化');
    grid on;

else
    disp(['求解失败: ', result.status]);
end

---

📌 注意事项：

1. 数据准备：你需要根据实际系统提供 W_forecast（新能源预测）、Load（负荷预测）、火电机组参数。
2. Gurobi 安装配置：确保你已经正确安装了 Gurobi 并在 MATLAB 中设置了路径。
3. 扩展功能：你可以添加储能、抽水蓄能、电价响应、备用容量等高级特性。
4. 加速求解：可以将部分变量设为整数（如启停机状态），变成真正的混合整数线性规划（MILP）。

---

我们可以看到三条曲线分别代表了日负荷、一次调节后的净负荷以及二次调整后的净负荷。为了生成类似的图表，我们可以使用MATLAB编写代码来模拟这些数据并绘制相应的图形。以下是一个示例代码，

% 创建时间向量
t = 0:24;

% 模拟日负荷数据
daily_load = [900, 850, 750, 700, 720, 800, 900, 1000, 1100, 1200, 1300, 1400, 1500, 1400, 1300, 1200, 1100, 1000, 900, 1000, 1400, 1300, 1200, 1100, 1000];

% 模拟一次调节后净负荷数据
first_adjusted_load = daily_load - 100; % 假设每次调节减少100WM

% 模拟二次调整后净负荷数据
second_adjusted_load = first_adjusted_load - 50; % 再次假设每次调节减少50WM

% 绘制图形
figure;
plot(t, daily_load, '-r*', 'LineWidth', 1.5);
hold on;
plot(t, first_adjusted_load, '-g^', 'LineWidth', 1.5);
plot(t, second_adjusted_load, '-bd', 'LineWidth', 1.5);

% 添加图例
legend('日负荷', '一次调节后净负荷', '二次调整后净负荷');

% 设置坐标轴标签
xlabel('时间/t');
ylabel('高温储放热/WM');

% 设置标题
title('负荷调节示意图');

% 设置网格
grid on;

% 设置Y轴范围
ylim([600 1600]);

说明：

• daily_load：模拟的日负荷数据。
• first_adjusted_load：模拟的一次调节后净负荷数据，这里简单地假设每次调节减少100WM。
• second_adjusted_load：模拟的二次调整后净负荷数据，再次假设每次调节减少50WM。

---