七 dis_max--best fields策略

本文探讨了在搜索引擎中如何通过布尔模型和最佳字段策略计算文档与查询的相关性得分,详细解析了不同策略下文档排序的变化,以及如何优化查询以获取更相关的结果。

1. 假设数据

doc1:

title: black cat.

content: The cat like eat fish.

doc2:

title: red dog.

content: The dog like eat bone.

doc3:

title: yellow cat.

content: yellow cat look good.

doc4:

title: white dog.

content: People like white dog.

doc5:

title: animal.

content: cat is animal, and dog is also animal.

 

2. 查询需求

查询title或content中包含cat或is的帖子。

GET /forum/article/_search
{
    "query": {
        "bool": {
            "should": [
                { "match": { "title": "cat is" }},
                { "match": { "content":  "cat is" }}
            ]
        }
    }
}

 

3.结果分析

查询结果排序可能是:doc1>doc3>doc5

 

4. 原因分析

计算每个document的relevance score:每个query的分数相加,乘以matched query数量,除以总query数量。

doc1---->  title查询分数:1.1,content查询分数:1.2,两者相加得2.3;matched query数量为2(title和content各1次),总查询数量为2(title和content各被查询了1次)。relevance score:2.3 * 2 / 2 = 2.3

doc3与doc1类似。

doc5---->  title查询分数:没有match到,为0,content查询分数:2.3(cat和is都匹配到了,分数较高),两者相加得2.3;matched query数量为1(content 1次),总查询数量为2(title和content各被查询了1次)。relevance score:2.3 * 1 / 2 = 1.15

 

5. best fields策略

如果我们想要到查询结果是:一个filed中能够匹配到尽可能多的关键字,则排在前面,即让doc5排在前面,则可以使用下面的查询方式:

GET /forum/article/_search
{
    "query": {
        "dis_max": {
            "queries": [
                { "match": { "title": "cat is" }},
                { "match": { "content":  "cat is" }}
            ]
        }
    }
}

这样就会直接以最大的filed分数作为relevance score,那么doc1的relevance score分数为1.2,doc5的relevance score分数为2.3

请修改以下代码:%% 基于电-氢混合储能的电力系统优化调度模型(日内滚动优化) clc clear %% 参数设置 T = 96; % 时间周期数 (15分钟间隔) N_Gens = 6; % 常规机组数量 W_Size = 16; % 滚动窗口大小 (4小时) R_Step = 4; % 滚动步长 (1小时) N_Windows = floor((T - W_Size)/R_Step) + 1; % 滚动次数 %% 数据加载和预处理 load('P_load_24.mat', 'P_load_e'); % 基础负荷数据 load('P_wt_24', 'P_wt_e'); % 风电数据 mpc = loadcase('case30'); baseMVA = mpc.baseMVA; % 将24小时数据扩展到96个15分钟时段 P_wt_96 = repelem(P_wt_e, 1, 4); P_load_96 = repelem(P_load_e, 1, 4); % 常规机组参数 [机组ID, P_max, P_min, a, b, c, ramp, 碳排放系数] Gen_Params = [1 250 50 0.0375 20 372.5 72 0.85; 2 80 20 0.175 17.5 352.3 48 0.92; 3 50 15 0.625 10 316.5 30 1.05; 4 35 10 0.0834 32.5 329.2 21 0.78; 5 30 10 0.25 30 276.4 18 0.95; 6 40 12 0.25 30 232.2 24 0.88]; % 提取发电机参数 a2 = Gen_Params(:,4); % 二次项系数 a1 = Gen_Params(:,5); % 一次项系数 a0 = Gen_Params(:,6); % 常数项 Ramp_e = Gen_Params(:,7); % 爬坡率 (MW/15min) P_GT_Min = Gen_Params(:,3); % 最小出力 P_GT_Max = Gen_Params(:,2); % 最大出力 % 启动成本参数 (单位: 元) C_Start_qt = [2768; 2874; 2136; 2677; 2262; 2869]; % 电储能参数 P_Cha_Max = 45; % 电储能最大充电功率 P_Dis_Max = 90; % 电储能最大放电功率 K_ess_Eff = 0.92; % 电储能转换效率 C_ess_Max = 300; % 储能容量限制 SOC_ESS_Init = 10; % 初始电SOC (MWh) % 氢系统参数 P_FC_Max = 400; % 燃料电池最大功率 K_H2P_xs = 9.1; % 氢转电系数 (kg/MWh) P_EL_Max = 400; % 电解槽最大功率 K_P2H_xs = 0.22; % 电转氢系数 (kg/MWh) H2_cha_Max = 100; % 最大充氢速率 (kg/15min) H2_dis_Max = 100; % 最大放氢速率 (kg/15min) C_HSS_Max = 400; % 最大储氢容量 (kg) C_HSS_Min = 0; % 最小储氢容量 (kg) SOC_HSS_Init = 20; % 初始氢SOC (kg) % 调整氢系统功率上限 P_FC_Max_adj = min(P_FC_Max, H2_dis_Max/(0.25*K_H2P_xs)); P_EL_Max_adj = min(P_EL_Max, H2_cha_Max/(0.25*K_P2H_xs)); % 节点位置 Bus_Gen = [1 2 5 8 11 13]; % 火电机组节点 Bus_Wt = 2; % 风电节点 Bus_Ess = 11; % 电储能节点 Bus_Hss = 8; % 氢储能节点 % 环境成本系数 C_carbon = 150; % 碳排放成本 (元/吨CO2) Carbon_Factor = Gen_Params(:,8); % 碳排放系数 (tCO2/MWh) %% ==================== 网络拓扑初始化 ==================== [Ybus, ~, ~] = makeYbus(mpc); B = imag(Ybus); % 电纳矩阵虚部 ref_node = 1; % 参考节点 non_ref = 2:30; % 非参考节点 Bred = B(non_ref, non_ref); % 降阶B矩阵 % 创建设备-节点映射矩阵 N = size(mpc.bus, 1); G_map = sparse(Bus_Gen, 1:N_Gens, 1, N, N_Gens); % 常规机组 W_map = sparse(Bus_Wt, 1, 1, N, 1); % 风电机组 E_map = sparse(Bus_Ess, 1, 1, N, 1); % 电储能 H_map = sparse(Bus_Hss, 1, 1, N, 1); % 氢系统 % 初始化t=0时刻状态 B_ug_z0 = ones(N_Gens,1); % 初始机组状态 (全运行) P_GT_z0 = P_GT_Min; % 初始机组出力 (最小出力) %% ==================== 主循环滚动优化 ==================== % 初始化结果存储结构 Results = struct(); Results.P_GT = repmat(P_GT_Min, 1, T); % 初始化为最小出力 Results.P_wt = P_wt_96; % 初始化为预测值 Results.P_cha = zeros(1, T); Results.P_dis = zeros(1, T); Results.SOC_ESS = SOC_ESS_Init * ones(1, T+1); Results.P_FC = zeros(1, T); Results.P_EL = zeros(1, T); Results.SOC_HSS = SOC_HSS_Init * ones(1, T+1); B_ug_96_full = ones(N_Gens, T); % 初始状态全为运行 % 成本跟踪 total_fuel_cost = 0; total_start_cost = 0; total_carbon_cost = 0; % YALMIP选项设置 Opts = sdpsettings('solver', 'cplex', 'verbose', 0, 'showprogress', 0, ... 'warning', 1, 'cplex.timelimit', 300); % 滚动优化主循环 for win_idx = 1:N_Windows k_start_t = (win_idx-1)*R_Step + 1; k_end_t = min(k_start_t + W_Size - 1, T); % 防止越界 k_win_len = k_end_t - k_start_t + 1; % 实际窗口长度 fprintf('\n====== 处理窗口 %d/%d: 时段 %d-%d (长度=%d) ======\n', win_idx, N_Windows, k_start_t, k_end_t, k_win_len); % 获取当前窗口初始状态 if k_start_t == 1 SOC_ESS_init = SOC_ESS_Init; SOC_HSS_init = SOC_HSS_Init; prev_status = B_ug_z0; prev_power = P_GT_z0; else SOC_ESS_init = Results.SOC_ESS(k_start_t); SOC_HSS_init = Results.SOC_HSS(k_start_t); prev_status = B_ug_96_full(:, k_start_t-1); prev_power = Results.P_GT(:, k_start_t-1); end fprintf('初始状态: SOC_ESS=%.2f MWh, SOC_HSS=%.2f kg\n', SOC_ESS_init, SOC_HSS_init); % 定义当前窗口变量 % 机组变量 P_GT_e = sdpvar(N_Gens, k_win_len, 'full'); % 火电机组出力 B_U_GT = binvar(N_Gens, k_win_len, 'full'); % 机组运行状态 B_Start_GT = binvar(N_Gens, k_win_len, 'full'); % 启动状态 % 可再生能源 P_wt_e = sdpvar(1, k_win_len); % 风电出力 % 电储能 P_cha_e = sdpvar(1, k_win_len); % 充电功率 P_dis_e = sdpvar(1, k_win_len); % 放电功率 SOC_ESS = sdpvar(1, k_win_len+1); % SOC状态 % 氢系统 P_FC_e = sdpvar(1, k_win_len); % 燃料电池出力 P_e_EL = sdpvar(1, k_win_len); % 电解槽用电 SOC_HSS = sdpvar(1, k_win_len+1); % 氢SOC状态 % 网络变量 Bus_Theta = sdpvar(30, k_win_len, 'full'); % 节点相角 % 约束初始化 Cons = []; % ==================== 初始状态约束 ==================== Cons = [Cons, SOC_ESS(1) == SOC_ESS_init]; Cons = [Cons, SOC_HSS(1) == SOC_HSS_init]; % ==================== 机组约束 ==================== for t_wd = 1:k_win_len t_global = k_start_t + t_wd - 1; % 机组出力上下限 Cons = [Cons, B_U_GT(:, t_wd).*P_GT_Min <= P_GT_e(:, t_wd) <= B_U_GT(:, t_wd).*P_GT_Max]; % 启动状态逻辑 if t_global == 1 Cons = [Cons, B_Start_GT(:, t_wd) == B_U_GT(:, t_wd) - B_ug_z0]; else if t_wd == 1 prev_win_status = prev_status; else prev_win_status = B_U_GT(:, t_wd-1); end Cons = [Cons, B_Start_GT(:, t_wd) >= B_U_GT(:, t_wd) - prev_win_status]; Cons = [Cons, B_Start_GT(:, t_wd) <= B_U_GT(:, t_wd)]; Cons = [Cons, B_Start_GT(:, t_wd) <= 1 - prev_win_status]; end % 爬坡约束 if t_wd == 1 % 第一个时间步 Cons = [Cons, P_GT_e(:, t_wd) - prev_power <= Ramp_e .* prev_status + (B_U_GT(:, t_wd) - prev_status) .* P_GT_Min]; Cons = [Cons, prev_power - P_GT_e(:, t_wd) <= Ramp_e .* prev_status + (prev_status - B_U_GT(:, t_wd)) .* P_GT_Min]; else % 后续时间步 Cons = [Cons, P_GT_e(:, t_wd) - P_GT_e(:, t_wd-1) <= Ramp_e .* B_U_GT(:, t_wd) + (B_U_GT(:, t_wd) - B_U_GT(:, t_wd-1)) .* P_GT_Min]; Cons = [Cons, P_GT_e(:, t_wd-1) - P_GT_e(:, t_wd) <= Ramp_e .* B_U_GT(:, t_wd-1) + (B_U_GT(:, t_wd-1) - B_U_GT(:, t_wd)) .* P_GT_Min]; end end % ==================== 风电机组约束 ==================== Cons = [Cons, 0 <= P_wt_e <= P_wt_96(k_start_t:k_end_t)]; % ==================== 电储能约束 ==================== bin_cha = binvar(1, k_win_len); % 充电状态指示器 bin_dis = binvar(1, k_win_len); % 放电状态指示器 for t_wd = 1:k_win_len % SOC动态更新 Cons = [Cons, SOC_ESS(t_wd+1) == SOC_ESS(t_wd) + K_ess_Eff * P_cha_e(t_wd) - P_dis_e(t_wd)/K_ess_Eff]; % SOC上下限 Cons = [Cons, 0 <= SOC_ESS(t_wd+1) <= C_ess_Max]; % 充放电功率限制 Cons = [Cons, 0 <= P_cha_e(t_wd) <= bin_cha(t_wd)*P_Cha_Max]; Cons = [Cons, 0 <= P_dis_e(t_wd) <= bin_dis(t_wd)*P_Dis_Max]; Cons = [Cons, bin_cha(t_wd) + bin_dis(t_wd) <= 1]; % 互斥约束 end % ==================== 氢系统约束 ==================== for t_wd = 1:k_win_len % 功率约束 Cons = [Cons, 0 <= P_FC_e(t_wd) <= P_FC_Max_adj]; Cons = [Cons, 0 <= P_e_EL(t_wd) <= P_EL_Max_adj]; % 氢流量计算 (kg/15min) H2_FC = P_FC_e(t_wd) * 0.25 * K_H2P_xs; % 燃料电池耗氢量 H2_EL = P_e_EL(t_wd) * 0.25 * K_P2H_xs; % 电解槽产氢量 % 氢平衡约束 Cons = [Cons, SOC_HSS(t_wd+1) == SOC_HSS(t_wd) - H2_FC + H2_EL]; % SOC上下限 Cons = [Cons, C_HSS_Min <= SOC_HSS(t_wd+1) <= C_HSS_Max]; % 充放氢速率限制 Cons = [Cons, 0 <= H2_FC <= H2_dis_Max]; Cons = [Cons, 0 <= H2_EL <= H2_cha_Max]; end % ==================== 网络约束 ==================== for t_wd = 1:k_win_len t_global = k_start_t + t_wd - 1; % 节点功率平衡 total_injection = G_map * P_GT_e(:, t_wd) + ... W_map * P_wt_e(t_wd) + ... E_map * (P_dis_e(t_wd) - P_cha_e(t_wd)) + ... H_map * (P_FC_e(t_wd) - P_e_EL(t_wd)); % 基础负荷 load_vector = P_load_96(:, t_global); % 节点净注入功率 P_net = total_injection - load_vector; % 直流潮流约束 Cons = [Cons, Bred * Bus_Theta(non_ref, t_wd) == P_net(non_ref) / baseMVA]; end % 参考节点相角约束 Cons = [Cons, Bus_Theta(ref_node, :) == 0]; % ==================== 目标函数 ==================== Obj = 0; for t_wd = 1:k_win_len % 1. 燃料成本 (乘以0.25转换为15分钟成本) for g = 1:N_Gens fuel_cost = a2(g)*(P_GT_e(g,t_wd)^2) + a1(g)*P_GT_e(g,t_wd) + a0(g)*B_U_GT(g,t_wd); Obj = Obj + 0.25 * fuel_cost; end % 2. 启动成本 for g = 1:N_Gens Obj = Obj + C_Start_qt(g) * B_Start_GT(g, t_wd); end % 3. 碳排放成本 for g = 1:N_Gens carbon_cost = C_carbon * Carbon_Factor(g) * P_GT_e(g,t_wd) * 0.001 * 0.25; Obj = Obj + carbon_cost; end end % 求解优化问题 fprintf('>> 求解窗口 %d 优化问题...\n', win_idx); Sol = optimize(Cons, Obj, Opts); % 保存结果 if Sol.problem == 0 % 提取当前窗口结果 Results.P_GT(:, k_start_t:k_end_t) = value(P_GT_e); Results.P_wt(k_start_t:k_end_t) = value(P_wt_e); Results.P_cha(k_start_t:k_end_t) = value(P_cha_e); Results.P_dis(k_start_t:k_end_t) = value(P_dis_e); Results.SOC_ESS(k_start_t+1:k_end_t+1) = value(SOC_ESS(2:end)); Results.P_FC(k_start_t:k_end_t) = value(P_FC_e); Results.P_EL(k_start_t:k_end_t) = value(P_e_EL); Results.SOC_HSS(k_start_t+1:k_end_t+1) = value(SOC_HSS(2:end)); % 保存机组状态到全局数组 B_ug_96_full(:, k_start_t:k_end_t) = value(B_U_GT); fprintf('>>> 优化成功 | 成本: %.2f 元\n', value(Obj)); % 计算并记录成本 win_obj = value(Obj); total_fuel_cost = total_fuel_cost + win_obj * 0.7; % 示例比例 total_start_cost = total_start_cost + win_obj * 0.2; total_carbon_cost = total_carbon_cost + win_obj * 0.1; else fprintf('>>> 优化失败: %s\n', Sol.info); yalmiperror(Sol.problem); % 使用保守策略:最小出力运行 Results.P_GT(:, k_start_t:k_end_t) = repmat(P_GT_Min, 1, k_win_len); Results.P_wt(k_start_t:k_end_t) = P_wt_96(k_start_t:k_end_t); fprintf('>>> 使用默认出力值\n'); end % 添加系统平衡检查 for t = k_start_t:k_end_t total_gen = sum(Results.P_GT(:,t)) + Results.P_wt(t) + ... (Results.P_dis(t) - Results.P_cha(t)) + ... (Results.P_FC(t) - Results.P_EL(t)); total_load = sum(P_load_96(:,t)); imbalance = total_gen - total_load; if abs(imbalance) > 1e-3 fprintf('警告: 时段 %d 功率不平衡: %.4f MW\n', t, imbalance); end end end
最新发布
11-27
请检查以下代码,并整理出一个修改后的完整代码:%% 基于电-氢混合储能的电力系统优化调度模型(最终修正版) clc clear close all %% ==================== 参数设置 ==================== T = 96; % 时间周期数 (15分钟间隔) N_Gens = 6; % 常规机组数量 W_Size = 16; % 滚动窗口大小 (4小时) R_Step = 4; % 滚动步长 (1小时) N_Windows = floor((T - W_Size)/R_Step) + 1; % 滚动次数 %% ==================== 数据加载和预处理 ==================== % 安全加载风电数据 try load(‘P_wt_24.mat’, ‘P_wt_e’); fprintf(‘成功加载风电数据\n’); if isstruct(P_wt_e) fields = fieldnames(P_wt_e); P_wt_e = P_wt_e.(fields{1}); end P_wt_e = double(P_wt_e(😃’); % 确保行向量 catch warning(‘使用模拟风电数据’); % 生成模拟风电数据 (24点) P_wt_e = 50 + 20 * sin(linspace(0, 2*pi, 24)); end % 扩展为96点 if numel(P_wt_e) == 24 P_wt_96 = repelem(P_wt_e, 1, 4); else P_wt_96 = P_wt_e(1:min(96, numel(P_wt_e))); if numel(P_wt_96) < 96 P_wt_96 = [P_wt_96, zeros(1, 96 - numel(P_wt_96))]; end end P_wt_96 = reshape(P_wt_96, 1, []); % 确保1×96 % 安全加载负荷数据 try load(‘P_load_24.mat’, ‘P_load_e’); fprintf(‘成功加载负荷数据\n’); if isstruct(P_load_e) fields = fieldnames(P_load_e); P_load_e = P_load_e.(fields{1}); end P_load_e = double(P_load_e); catch warning(‘使用模拟负荷数据’); % 生成模拟负荷数据 (30节点×24小时) P_load_e = 100 + 50 * rand(30, 24); end % 标准化负荷数据格式并扩展为96点 if size(P_load_e, 2) == 24 P_load_96 = repelem(P_load_e, 1, 4); elseif size(P_load_e, 2) < 96 P_load_96 = [P_load_e, repmat(P_load_e(:, end), 1, 96 - size(P_load_e, 2))]; else P_load_96 = P_load_e(:, 1:96); end % 确保30×96维度 if size(P_load_96, 1) ~= 30 if size(P_load_96, 1) == 1 P_load_96 = repmat(P_load_96, 30, 1); else P_load_96 = repmat(P_load_96(1, 😃, 30, 1); end end % 加载电网拓扑 try mpc = loadcase(‘case30’); fprintf(‘成功加载电网拓扑数据\n’); catch warning(‘使用默认30节点系统’); % 创建默认30节点系统 mpc = case30; end baseMVA = mpc.baseMVA; % 常规机组参数 [机组ID, P_max, P_min, a, b, c, ramp, 碳排放系数] Gen_Params = [1 250 50 0.0375 20 372.5 72 0.85; 2 80 20 0.175 17.5 352.3 48 0.92; 3 50 15 0.625 10 316.5 30 1.05; 4 35 10 0.0834 32.5 329.2 21 0.78; 5 30 10 0.25 30 276.4 18 0.95; 6 40 12 0.25 30 232.2 24 0.88]; % 提取发电机参数 a2 = Gen_Params(:,4); % 二次项系数 a1 = Gen_Params(:,5); % 一次项系数 a0 = Gen_Params(:,6); % 常数项 Ramp_e = Gen_Params(:,7); % 爬坡率 (MW/15min) P_GT_Min = Gen_Params(:,3); % 最小出力 P_GT_Max = Gen_Params(:,2); % 最大出力 % 启动成本参数 (单位: 元) C_Start_qt = [2768; 2874; 2136; 2677; 2262; 2869]; % 电储能参数 P_Cha_Max = 45; % 电储能最大充电功率 P_Dis_Max = 90; % 电储能最大放电功率 K_ess_Eff = 0.92; % 电储能转换效率 C_ess_Max = 300; % 储能容量限制 SOC_ESS_Init = 10; % 初始电SOC (MWh) % 氢系统参数 P_FC_Max = 400; % 燃料电池最大功率 K_H2P_xs = 9.1; % 氢转电系数 (kg/MWh) P_EL_Max = 400; % 电解槽最大功率 K_P2H_xs = 0.22; % 电转氢系数 (kg/MWh) H2_cha_Max = 100; % 最大充氢速率 (kg/15min) H2_dis_Max = 100; % 最大放氢速率 (kg/15min) C_HSS_Max = 400; % 最大储氢容量 (kg) C_HSS_Min = 0; % 最小储氢容量 (kg) SOC_HSS_Init = 20; % 初始氢SOC (kg) % 调整氢系统功率上限 P_FC_Max_adj = min(P_FC_Max, H2_dis_Max/(0.25K_H2P_xs)); P_EL_Max_adj = min(P_EL_Max, H2_cha_Max/(0.25K_P2H_xs)); % 节点位置 Bus_Gen = [1 2 5 8 11 13]; % 火电机组节点 Bus_Wt = 2; % 风电节点 Bus_Ess_E = 11; % 电储能节点 Bus_Hss_H = 8; % 氢储能节点 % 环境成本系数 C_carbon = 150; % 碳排放成本 (元/吨CO2) Carbon_Factor = Gen_Params(:,8); % 碳排放系数 (tCO2/MWh) %% ==================== 网络拓扑初始化 ==================== [Ybus, ~, ~] = makeYbus(mpc); B = imag(Ybus); % 电纳矩阵虚部 ref_node = 1; % 参考节点 non_ref = 2:30; % 非参考节点 Bred = B(non_ref, non_ref); % 降阶B矩阵 % 创建设备-节点映射矩阵 N = size(mpc.bus, 1); G_map = sparse(Bus_Gen, 1:N_Gens, 1, N, N_Gens); % 常规机组 W_map = sparse(Bus_Wt, 1, 1, N, 1); % 风电机组 E_map = sparse(Bus_Ess_E, 1, 1, N, 1); % 电储能 H_map = sparse(Bus_Hss_H, 1, 1, N, 1); % 氢系统 % 初始化t=0时刻状态 B_ug_z0 = ones(N_Gens,1); % 初始机组状态 (全运行) P_GT_z0 = P_GT_Min; % 初始机组出力 (最小出力) %% ==================== 主循环滚动优化 ==================== % 初始化结果存储结构 Results = struct(); Results.P_GT = repmat(P_GT_Min, 1, T); % 初始化为最小出力 Results.P_wt = P_wt_96; % 初始化为预测值 Results.P_cha = zeros(1, T); Results.P_dis = zeros(1, T); Results.SOC_ESS = SOC_ESS_Init * ones(1, T+1); Results.P_FC = zeros(1, T); Results.P_EL = zeros(1, T); Results.SOC_HSS = SOC_HSS_Init * ones(1, T+1); B_ug_96_full = ones(N_Gens, T); % 初始状态全为运行 % 成本跟踪 window_costs = zeros(1, N_Windows); fuel_costs = zeros(1, N_Windows); start_costs = zeros(1, N_Windows); carbon_costs = zeros(1, N_Windows); % YALMIP选项设置 Opts = sdpsettings(‘solver’, ‘cplex’, ‘verbose’, 1, ‘showprogress’, 1, … ‘warning’, 1, ‘debug’, 1, ‘cplex.timelimit’, 300); % 滚动优化主循环 for win_idx = 1:N_Windows k_start_t = (win_idx-1)*R_Step + 1; k_end_t = min(k_start_t + W_Size - 1, T); % 防止越界 k_win_len = k_end_t - k_start_t + 1; % 实际窗口长度 fprintf('\n====== 处理窗口 %d/%d: 时段 %d-%d (长度=%d) ======\n', ... win_idx, N_Windows, k_start_t, k_end_t, k_win_len); %% 获取当前窗口初始状态 if k_start_t == 1 SOC_ESS_init = SOC_ESS_Init; SOC_HSS_init = SOC_HSS_Init; prev_status = B_ug_z0; prev_power = P_GT_z0; else SOC_ESS_init = Results.SOC_ESS(k_start_t); SOC_HSS_init = Results.SOC_HSS(k_start_t); prev_status = B_ug_96_full(:, k_start_t-1); prev_power = Results.P_GT(:, k_start_t-1); end fprintf('初始状态: SOC_ESS=%.2f MWh, SOC_HSS=%.2f kg\n', ... SOC_ESS_init, SOC_HSS_init); %% 定义当前窗口变量 % 机组变量 P_GT_e = sdpvar(N_Gens, k_win_len, 'full'); % 火电机组出力 B_U_GT = binvar(N_Gens, k_win_len, 'full'); % 机组运行状态 B_Start_GT = binvar(N_Gens, k_win_len, 'full'); % 启动状态 % 可再生能源 P_wt_cur = sdpvar(1, k_win_len); % 风电出力 % 电储能 P_cha_e = sdpvar(1, k_win_len); % 充电功率 P_dis_e = sdpvar(1, k_win_len); % 放电功率 SOC_ESS = sdpvar(1, k_win_len+1); % SOC状态 % 氢系统 P_FC_e = sdpvar(1, k_win_len); % 燃料电池出力 P_e_EL = sdpvar(1, k_win_len); % 电解槽用电 SOC_HSS = sdpvar(1, k_win_len+1); % 氢SOC状态 % 网络变量 Bus_Theta = sdpvar(30, k_win_len, 'full'); % 节点相角 %% 约束初始化 Cons = []; % ==================== 初始状态约束 ==================== Cons = [Cons, SOC_ESS(1) == SOC_ESS_init]; Cons = [Cons, SOC_HSS(1) == SOC_HSS_init]; % ==================== 机组约束 ==================== for t_wd = 1:k_win_len t_global = k_start_t + t_wd - 1; % 机组出力上下限 Cons = [Cons, B_U_GT(:, t_wd).*P_GT_Min <= P_GT_e(:, t_wd) <= B_U_GT(:, t_wd).*P_GT_Max]; % 启动状态逻辑 if t_global == 1 Cons = [Cons, B_Start_GT(:, t_wd) == B_U_GT(:, t_wd) - B_ug_z0]; else if t_wd == 1 prev_win_status = prev_status; else prev_win_status = B_U_GT(:, t_wd-1); end Cons = [Cons, B_Start_GT(:, t_wd) >= B_U_GT(:, t_wd) - prev_win_status]; Cons = [Cons, B_Start_GT(:, t_wd) <= B_U_GT(:, t_wd)]; Cons = [Cons, B_Start_GT(:, t_wd) <= 1 - prev_win_status]; end % 爬坡约束 if t_wd == 1 % 第一个时间步 Cons = [Cons, P_GT_e(:, t_wd) - prev_power <= Ramp_e .* prev_status + ... (B_U_GT(:, t_wd) - prev_status) .* P_GT_Min]; Cons = [Cons, prev_power - P_GT_e(:, t_wd) <= Ramp_e .* prev_status + ... (prev_status - B_U_GT(:, t_wd)) .* P_GT_Min]; else % 后续时间步 Cons = [Cons, P_GT_e(:, t_wd) - P_GT_e(:, t_wd-1) <= Ramp_e .* B_U_GT(:, t_wd) + ... (B_U_GT(:, t_wd) - B_U_GT(:, t_wd-1)) .* P_GT_Min]; Cons = [Cons, P_GT_e(:, t_wd-1) - P_GT_e(:, t_wd) <= Ramp_e .* B_U_GT(:, t_wd-1) + ... (B_U_GT(:, t_wd-1) - B_U_GT(:, t_wd)) .* P_GT_Min]; end end % ==================== 风电机组约束 ==================== Cons = [Cons, 0 <= P_wt_cur <= P_wt_96(k_start_t:k_end_t)]; % ==================== 电储能约束 ==================== bin_cha = binvar(1, k_win_len); % 充电状态指示器 bin_dis = binvar(1, k_win_len); % 放电状态指示器 for t_wd = 1:k_win_len % SOC动态更新 Cons = [Cons, SOC_ESS(t_wd+1) == SOC_ESS(t_wd) + ... K_ess_Eff * P_cha_e(t_wd) - P_dis_e(t_wd)/K_ess_Eff]; % SOC上下限 Cons = [Cons, 0 <= SOC_ESS(t_wd+1) <= C_ess_Max]; % 充放电功率限制 Cons = [Cons, 0 <= P_cha_e(t_wd) <= bin_cha(t_wd)*P_Cha_Max]; Cons = [Cons, 0 <= P_dis_e(t_wd) <= bin_dis(t_wd)*P_Dis_Max]; Cons = [Cons, bin_cha(t_wd) + bin_dis(t_wd) <= 1]; % 互斥约束 end % ==================== 氢系统约束 ==================== for t_wd = 1:k_win_len % 功率约束 Cons = [Cons, 0 <= P_FC_e(t_wd) <= P_FC_Max_adj]; Cons = [Cons, 0 <= P_e_EL(t_wd) <= P_EL_Max_adj]; % 氢流量计算 (kg/15min) H2_FC = P_FC_e(t_wd) * 0.25 * K_H2P_xs; % 燃料电池耗氢量 H2_EL = P_e_EL(t_wd) * 0.25 * K_P2H_xs; % 电解槽产氢量 % 氢平衡约束 Cons = [Cons, SOC_HSS(t_wd+1) == SOC_HSS(t_wd) - H2_FC + H2_EL]; % SOC上下限 Cons = [Cons, C_HSS_Min <= SOC_HSS(t_wd+1) <= C_HSS_Max]; % 充放氢速率限制 Cons = [Cons, 0 <= H2_FC <= H2_dis_Max]; Cons = [Cons, 0 <= H2_EL <= H2_cha_Max]; end % ==================== 网络约束 ==================== for t_wd = 1:k_win_len t_global = k_start_t + t_wd - 1; % 节点功率平衡 total_injection = G_map * P_GT_e(:, t_wd) + ... W_map * P_wt_cur(t_wd) + ... E_map * (P_dis_e(t_wd) - P_cha_e(t_wd)) + ... H_map * (P_FC_e(t_wd) - P_e_EL(t_wd)); % 基础负荷 load_vector = P_load_96(:, t_global); % 节点净注入功率 P_net = total_injection - load_vector; % 直流潮流约束 Cons = [Cons, Bred * Bus_Theta(non_ref, t_wd) == P_net(non_ref) / baseMVA]; end % 参考节点相角约束 Cons = [Cons, Bus_Theta(ref_node, :) == 0]; % ==================== 目标函数 ==================== Obj = 0; for t_wd = 1:k_win_len % 1. 燃料成本 (乘以0.25转换为15分钟成本) for g = 1:N_Gens fuel_cost = a2(g)*(P_GT_e(g,t_wd)^2) + a1(g)*P_GT_e(g,t_wd) + a0(g)*B_U_GT(g,t_wd); Obj = Obj + 0.25 * fuel_cost; end % 2. 启动成本 for g = 1:N_Gens Obj = Obj + C_Start_qt(g) * B_Start_GT(g, t_wd); end % 3. 碳排放成本 (吨到千克转换: 元/吨CO2 * 吨CO2/MWh * MW * 0.25小时 * 0.001吨/千克) for g = 1:N_Gens carbon_cost = C_carbon * Carbon_Factor(g) * P_GT_e(g,t_wd) * 0.001 * 0.25; Obj = Obj + carbon_cost; end end %% 求解优化问题 fprintf('>> 求解窗口 %d 优化问题...\n', win_idx); Sol = optimize(Cons, Obj, Opts); %% 保存结果 if Sol.problem == 0 % 提取当前窗口结果 Results.P_GT(:, k_start_t:k_end_t) = value(P_GT_e); Results.P_wt(k_start_t:k_end_t) = value(P_wt_cur); Results.P_cha(k_start_t:k_end_t) = value(P_cha_e); Results.P_dis(k_start_t:k_end_t) = value(P_dis_e); Results.SOC_ESS(k_start_t+1:k_end_t+1) = value(SOC_ESS(2:end)); Results.P_FC(k_start_t:k_end_t) = value(P_FC_e); Results.P_EL(k_start_t:k_end_t) = value(P_e_EL); Results.SOC_HSS(k_start_t+1:k_end_t+1) = value(SOC_HSS(2:end)); % 保存机组状态到全局数组 B_ug_96_full(:, k_start_t:k_end_t) = value(B_U_GT); % 计算并保存窗口成本 window_costs(win_idx) = value(Obj); fprintf('>>> 优化成功 | 成本: %.2f 元\n', value(Obj)); % 分解成本并输出 fuel_cost_win = 0; start_cost_win = 0; carbon_cost_win = 0; for t_wd = 1:k_win_len for g = 1:N_Gens % 燃料成本 fuel_cost = a2(g)*(value(P_GT_e(g,t_wd))^2) + ... a1(g)*value(P_GT_e(g,t_wd)) + ... a0(g)*value(B_U_GT(g,t_wd)); fuel_cost_win = fuel_cost_win + 0.25 * fuel_cost; % 启动成本 if value(B_Start_GT(g, t_wd)) == 1 start_cost_win = start_cost_win + C_Start_qt(g); end % 碳排放成本 carbon_cost = C_carbon * Carbon_Factor(g) * ... value(P_GT_e(g,t_wd)) * 0.001 * 0.25; carbon_cost_win = carbon_cost_win + carbon_cost; end end fuel_costs(win_idx) = fuel_cost_win; start_costs(win_idx) = start_cost_win; carbon_costs(win_idx) = carbon_cost_win; fprintf(' - 燃料成本: %.2f 元\n', fuel_cost_win); fprintf(' - 启动成本: %.2f 元\n', start_cost_win); fprintf(' - 碳排成本: %.2f 元\n', carbon_cost_win); else fprintf('>>> 优化失败: %s\n', Sol.info); yalmiperror(Sol.problem); % 使用保守策略:最小出力运行 Results.P_GT(:, k_start_t:k_end_t) = repmat(P_GT_Min, 1, k_win_len); Results.P_wt(k_start_t:k_end_t) = P_wt_96(k_start_t:k_end_t); fprintf('>>> 使用默认出力值\n'); end % 添加系统平衡检查 for t = k_start_t:k_end_t total_gen = sum(Results.P_GT(:,t)) + Results.P_wt(t) + ... (Results.P_dis(t) - Results.P_cha(t)) + ... (Results.P_FC(t) - Results.P_EL(t)); total_load = sum(P_load_96(:,t)); imbalance = total_gen - total_load; if abs(imbalance) > 1e-3 fprintf('警告: 时段 %d 功率不平衡: %.4f MW\n', t, imbalance); end end end %% ==================== 结果后处理与可视化 ==================== fprintf(‘\n====== 优化完成,开始结果分析 ======\n’); % 时间向量 (15分钟间隔) time_vector = linspace(0, 24, T); % 1. 机组出力曲线 figure(‘Name’, ‘机组出力曲线’, ‘Position’, [100, 100, 1200, 600]); for i = 1:N_Gens subplot(2, 3, i); plot(time_vector, Results.P_GT(i, 😃, ‘LineWidth’, 1.5); hold on; plot([time_vector(1), time_vector(end)], [P_GT_Min(i), P_GT_Min(i)], ‘r–’); plot([time_vector(1), time_vector(end)], [P_GT_Max(i), P_GT_Max(i)], ‘r–’); title(sprintf(‘机组 %d 出力 (%.0f-%.0f MW)’, i, P_GT_Min(i), P_GT_Max(i))); xlabel(‘时间 (小时)’); ylabel(‘出力 (MW)’); grid on; xlim([0, 24]); end % 2. 新能源与储能出力 figure(‘Name’, ‘新能源与储能出力’, ‘Position’, [100, 100, 1200, 800]); subplot(3, 1, 1); plot(time_vector, Results.P_wt, ‘b’, ‘LineWidth’, 1.5); hold on; plot(time_vector, P_wt_96, ‘r–’); title(‘风电实际出力 vs 预测’); xlabel(‘时间 (小时)’); ylabel(‘出力 (MW)’); legend(‘实际’, ‘预测’); grid on; xlim([0, 24]); subplot(3, 1, 2); net_ess = Results.P_dis - Results.P_cha; plot(time_vector, net_ess, ‘r’, ‘LineWidth’, 1.5); title(‘电储能净出力’); xlabel(‘时间 (小时)’); ylabel(‘出力 (MW)’); grid on; xlim([0, 24]); subplot(3, 1, 3); net_h2 = Results.P_FC - Results.P_EL; plot(time_vector, net_h2, ‘g’, ‘LineWidth’, 1.5); title(‘氢系统净出力’); xlabel(‘时间 (小时)’); ylabel(‘出力 (MW)’); grid on; xlim([0, 24]); % 3. 储能状态变化 figure(‘Name’, ‘储能状态变化’, ‘Position’, [100, 100, 1200, 500]); subplot(1, 2, 1); plot(time_vector, Results.SOC_ESS(1:end-1), ‘b’, ‘LineWidth’, 1.5); hold on; plot([time_vector(1), time_vector(end)], [0, 0], ‘k–’); plot([time_vector(1), time_vector(end)], [C_ess_Max, C_ess_Max], ‘k–’); title(‘电储能SOC’); xlabel(‘时间 (小时)’); ylabel(‘SOC (MWh)’); grid on; xlim([0, 24]); subplot(1, 2, 2); plot(time_vector, Results.SOC_HSS(1:end-1), ‘r’, ‘LineWidth’, 1.5); hold on; plot([time_vector(1), time_vector(end)], [C_HSS_Min, C_HSS_Min], ‘k–’); plot([time_vector(1), time_vector(end)], [C_HSS_Max, C_HSS_Max], ‘k–’); title(‘氢储能SOC’); xlabel(‘时间 (小时)’); ylabel(‘SOC (kg)’); grid on; xlim([0, 24]); % 4. 系统总成本分析 total_fuel_cost = 0; total_start_cost = 0; total_carbon_cost = 0; for t = 1:T for g = 1:N_Gens % 仅当机组运行时计算燃料成本 if B_ug_96_full(g,t) == 1 % 燃料成本 P_gt_val = Results.P_GT(g,t); fuel_cost = a2(g)*(P_gt_val^2) + a1(g)*P_gt_val + a0(g); total_fuel_cost = total_fuel_cost + 0.25 * fuel_cost; % 碳排放成本 carbon_cost = C_carbon * Carbon_Factor(g) * P_gt_val * 0.001 * 0.25; total_carbon_cost = total_carbon_cost + carbon_cost; end % 启动成本 (仅在状态变化时计算) if t > 1 if B_ug_96_full(g,t) == 1 && B_ug_96_full(g,t-1) == 0 total_start_cost = total_start_cost + C_Start_qt(g); end end end end % 添加窗口成本汇总 figure(‘Name’, ‘窗口成本分析’, ‘Position’, [100, 100, 1200, 600]); subplot(2,2,1); bar(window_costs); title(‘各窗口总成本’); xlabel(‘窗口索引’); ylabel(‘成本(元)’); grid on; subplot(2,2,2); bar([fuel_costs; start_costs; carbon_costs]', ‘stacked’); title(‘成本分解’); xlabel(‘窗口索引’); ylabel(‘成本(元)’); legend(‘燃料’, ‘启动’, ‘碳排’); grid on; % 成本饼图 total_cost = total_fuel_cost + total_start_cost + total_carbon_cost; if abs(total_cost) < 1e-3 error(‘总成本为零,请检查成本计算逻辑’); end figure(‘Name’, ‘系统总成本分析’); cost_data = [total_fuel_cost, total_start_cost, total_carbon_cost]; pie(cost_data); legend({‘燃料成本’, ‘启动成本’, ‘碳排放成本’}, ‘Location’, ‘bestoutside’); title(sprintf(‘总成本: %.2f 元’, total_cost)); % 机组运行状态统计 figure(‘Name’, ‘机组运行状态’); for g = 1:N_Gens subplot(2, 3, g); plot(time_vector, B_ug_96_full(g, 😃, ‘LineWidth’, 2); title(sprintf(‘机组 %d 状态’, g)); ylim([-0.1, 1.1]); xlabel(‘时间(小时)’); ylabel(‘状态(0/1)’); grid on; end fprintf(‘====== 结果分析完成 ======\n’); fprintf(’ - 总燃料成本: %.2f 元\n’, total_fuel_cost); fprintf(’ - 总启动成本: %.2f 元\n’, total_start_cost); fprintf(’ - 总碳排放成本: %.2f 元\n’, total_carbon_cost); fprintf(’ - 系统总成本: %.2f 元\n’, total_cost); % 保存优化结果 save(‘optimization_results.mat’, ‘Results’, ‘B_ug_96_full’, ‘total_cost’, … ‘total_fuel_cost’, ‘total_start_cost’, ‘total_carbon_cost’); %% 基于电-氢混合储能的电力系统优化调度模型(最终稳定版) clc clear close all %% ==================== 参数设置 ==================== T = 96; % 时间周期数 (15分钟间隔) N_Gens = 6; % 常规机组数量 W_Size = 16; % 滚动窗口大小 (4小时) R_Step = 4; % 滚动步长 (1小时) N_Windows = floor((T - W_Size)/R_Step) + 1; % 滚动次数 %% ==================== 数据加载和预处理 ==================== % 安全加载风电数据 try % 尝试直接加载 load(‘P_wt_24.mat’, ‘P_wt_e’); catch try % 尝试UTF-8编码 load(‘P_wt_24.mat’, ‘-mat’, ‘-ascii’, ‘-nocompression’, ‘P_wt_e’); catch try % 尝试Latin1编码 load(‘P_wt_24.mat’, ‘-mat’, ‘-ascii’, ‘-latin1’, ‘P_wt_e’); catch ME warning(‘使用模拟风电数据: %s’, ME.message); % 生成模拟风电数据 P_wt_e = 50 + 20 * sin(linspace(0, 2*pi, 24)); end end end % 标准化风电数据格式 if isstruct(P_wt_e) fields = fieldnames(P_wt_e); P_wt_e = P_wt_e.(fields{1}); end if iscolumn(P_wt_e) P_wt_e = P_wt_e’; % 转置为行向量 end % 扩展为96点 if numel(P_wt_e) == 24 P_wt_96 = repelem(P_wt_e, 1, 4); else P_wt_96 = P_wt_e(1:min(96, numel(P_wt_e))); if numel(P_wt_96) < 96 P_wt_96 = [P_wt_96, zeros(1, 96 - numel(P_wt_96))]; end end P_wt_96 = reshape(P_wt_96, 1, []); % 确保1×96 % 安全加载负荷数据 try load(‘P_load_24.mat’, ‘P_load_e’); catch warning(‘使用模拟负荷数据’); % 生成模拟负荷数据 (30节点×24小时) P_load_e = 100 + 50 * rand(30, 24); end % 标准化负荷数据格式 if isvector(P_load_e) if size(P_load_e, 1) == 1 % 行向量 P_load_96 = repmat(P_load_e’, 1, 4); else % 列向量 P_load_96 = repmat(P_load_e, 1, 4); end else P_load_96 = repelem(P_load_e, 1, 4); end % 确保30×96维度 if size(P_load_96, 1) ~= 30 if size(P_load_96, 1) == 1 P_load_96 = repmat(P_load_96, 30, 1); else P_load_96 = repmat(P_load_96(1, 😃, 30, 1); end end if size(P_load_96, 2) < 96 P_load_96 = [P_load_96, repmat(P_load_96(:, end), 1, 96 - size(P_load_96, 2))]; end % 加载电网拓扑 try mpc = loadcase(‘case30’); catch % 创建默认30节点系统 mpc = case30; end baseMVA = mpc.baseMVA; % 常规机组参数 [机组ID, P_max, P_min, a, b, c, ramp, 碳排放系数] Gen_Params = [1 250 50 0.0375 20 372.5 72 0.85; 2 80 20 0.175 17.5 352.3 48 0.92; 3 50 15 0.625 10 316.5 30 1.05; 4 35 10 0.0834 32.5 329.2 21 0.78; 5 30 10 0.25 30 276.4 18 0.95; 6 40 12 0.25 30 232.2 24 0.88]; % 提取发电机参数 a2 = Gen_Params(:,4); % 二次项系数 a1 = Gen_Params(:,5); % 一次项系数 a0 = Gen_Params(:,6); % 常数项 Ramp_e = Gen_Params(:,7); % 爬坡率 (MW/15min) P_GT_Min = Gen_Params(:,3); % 最小出力 P_GT_Max = Gen_Params(:,2); % 最大出力 % 启动成本参数 (单位: 元) C_Start_qt = [2768; 2874; 2136; 2677; 2262; 2869]; % 电储能参数 P_Cha_Max = 45; % 电储能最大充电功率 P_Dis_Max = 90; % 电储能最大放电功率 K_ess_Eff = 0.92; % 电储能转换效率 C_ess_Max = 300; % 储能容量限制 SOC_ESS_Init = 10; % 初始电SOC (MWh) % 氢系统参数 P_FC_Max = 400; % 燃料电池最大功率 K_H2P_xs = 9.1; % 氢转电系数 (kg/MWh) P_EL_Max = 400; % 电解槽最大功率 K_P2H_xs = 0.22; % 电转氢系数 (kg/MWh) H2_cha_Max = 100; % 最大充氢速率 (kg/15min) H2_dis_Max = 100; % 最大放氢速率 (kg/15min) C_HSS_Max = 400; % 最大储氢容量 (kg) C_HSS_Min = 0; % 最小储氢容量 (kg) SOC_HSS_Init = 20; % 初始氢SOC (kg) % 调整氢系统功率上限 P_FC_Max_adj = min(P_FC_Max, H2_dis_Max/(0.25K_H2P_xs)); P_EL_Max_adj = min(P_EL_Max, H2_cha_Max/(0.25K_P2H_xs)); % 节点位置 Bus_Gen = [1 2 5 8 11 13]; % 火电机组节点 Bus_Wt = 2; % 风电节点 Bus_Ess_E = 11; % 电储能节点 Bus_Hss_H = 8; % 氢储能节点 % 环境成本系数 C_carbon = 150; % 碳排放成本 (元/吨CO2) Carbon_Factor = Gen_Params(:,8); % 碳排放系数 (tCO2/MWh) %% ==================== 网络拓扑初始化 ==================== [Ybus, ~, ~] = makeYbus(mpc); B = imag(Ybus); % 电纳矩阵虚部 ref_node = 1; % 参考节点 non_ref = 2:30; % 非参考节点 Bred = B(non_ref, non_ref); % 降阶B矩阵 % 创建设备-节点映射矩阵 N = size(mpc.bus, 1); G_map = sparse(Bus_Gen, 1:N_Gens, 1, N, N_Gens); % 常规机组 W_map = sparse(Bus_Wt, 1, 1, N, 1); % 风电机组 E_map = sparse(Bus_Ess_E, 1, 1, N, 1); % 电储能 H_map = sparse(Bus_Hss_H, 1, 1, N, 1); % 氢系统 % 初始化t=0时刻状态 B_ug_z0 = ones(N_Gens,1); % 初始机组状态 (全运行) P_GT_z0 = P_GT_Min; % 初始机组出力 (最小出力) %% ==================== 主循环滚动优化 ==================== % 初始化结果存储结构 Results = struct(); Results.P_GT = zeros(N_Gens, T); Results.P_wt = zeros(1, T); Results.P_cha = zeros(1, T); Results.P_dis = zeros(1, T); Results.SOC_ESS = SOC_ESS_Init * ones(1, T+1); Results.P_FC = zeros(1, T); Results.P_EL = zeros(1, T); Results.SOC_HSS = SOC_HSS_Init * ones(1, T+1); B_ug_96_full = zeros(N_Gens, T); % 全时段机组状态 B_ug_96_full(:,1) = B_ug_z0; % 设置初始机组状态 % 调试信息 fprintf(‘====== 系统初始化完成 ======\n’); fprintf(‘数据维度验证:\n’); fprintf(’ - 风电数据: %d×%d\n’, size(P_wt_96)); fprintf(’ - 负荷数据: %d×%d\n’, size(P_load_96)); fprintf(’ - 机组参数: %d×%d\n’, size(Gen_Params)); fprintf(‘开始滚动优化 (%d个窗口)…\n’, N_Windows); % YALMIP选项设置 Opts = sdpsettings(‘solver’, ‘cplex’, ‘verbose’, 1, ‘showprogress’, 1, … ‘warning’, 1, ‘debug’, 1); % 滚动优化主循环 for win_idx = 1:N_Windows k_start_t = (win_idx-1)*R_Step + 1; k_end_t = min(k_start_t + W_Size - 1, T); % 防止越界 k_win_len = k_end_t - k_start_t + 1; % 实际窗口长度 fprintf('\n====== 处理窗口 %d/%d: 时段 %d-%d (长度=%d) ======\n', ... win_idx, N_Windows, k_start_t, k_end_t, k_win_len); %% 获取当前窗口初始状态 if k_start_t == 1 SOC_ESS_init = SOC_ESS_Init; SOC_HSS_init = SOC_HSS_Init; else SOC_ESS_init = Results.SOC_ESS(k_start_t); SOC_HSS_init = Results.SOC_HSS(k_start_t); end %% 定义当前窗口变量 % 机组变量 P_GT_e = sdpvar(N_Gens, k_win_len, 'full'); % 火电机组出力 B_U_GT = binvar(N_Gens, k_win_len, 'full'); % 机组运行状态 B_Start_GT = binvar(N_Gens, k_win_len, 'full'); % 启动状态 % 可再生能源 P_wt_cur = sdpvar(1, k_win_len); % 风电出力 % 电储能 P_cha_e = sdpvar(1, k_win_len); % 充电功率 P_dis_e = sdpvar(1, k_win_len); % 放电功率 SOC_ESS = sdpvar(1, k_win_len+1); % SOC状态 % 氢系统 P_FC_e = sdpvar(1, k_win_len); % 燃料电池出力 P_e_EL = sdpvar(1, k_win_len); % 电解槽用电 SOC_HSS = sdpvar(1, k_win_len+1); % 氢SOC状态 % 网络变量 Bus_Theta = sdpvar(30, k_win_len, 'full'); % 节点相角 %% 约束初始化 Cons = []; % ==================== 初始状态约束 ==================== Cons = [Cons, SOC_ESS(1) == SOC_ESS_init]; Cons = [Cons, SOC_HSS(1) == SOC_HSS_init]; % ==================== 机组约束 ==================== for t_wd = 1:k_win_len t_global = k_start_t + t_wd - 1; % 获取前一时刻状态 if t_global == 1 prev_status = B_ug_z0; prev_power = P_GT_z0; else prev_status = B_ug_96_full(:, t_global-1); prev_power = Results.P_GT(:, t_global-1); end % 机组出力上下限 Cons = [Cons, B_U_GT(:, t_wd).*P_GT_Min <= P_GT_e(:, t_wd) <= B_U_GT(:, t_wd).*P_GT_Max]; % 启动状态逻辑 if t_global == 1 Cons = [Cons, B_Start_GT(:, t_wd) == B_U_GT(:, t_wd)]; else Cons = [Cons, B_Start_GT(:, t_wd) >= B_U_GT(:, t_wd) - prev_status]; Cons = [Cons, B_Start_GT(:, t_wd) <= B_U_GT(:, t_wd)]; Cons = [Cons, B_Start_GT(:, t_wd) <= 1 - prev_status]; end % 爬坡约束 if t_global >= 1 % 上爬坡约束 Cons = [Cons, P_GT_e(:, t_wd) - prev_power <= Ramp_e .* prev_status + ... (B_U_GT(:, t_wd) - prev_status) .* P_GT_Min]; % 下爬坡约束 Cons = [Cons, prev_power - P_GT_e(:, t_wd) <= Ramp_e .* B_U_GT(:, t_wd) + ... (prev_status - B_U_GT(:, t_wd)) .* P_GT_Min]; end end % ==================== 风电机组约束 ==================== Cons = [Cons, 0 <= P_wt_cur <= P_wt_96(k_start_t:k_end_t)]; % ==================== 电储能约束 ==================== bin_cha = binvar(1, k_win_len); % 充电状态指示器 bin_dis = binvar(1, k_win_len); % 放电状态指示器 for t_wd = 1:k_win_len % SOC动态更新 Cons = [Cons, SOC_ESS(t_wd+1) == SOC_ESS(t_wd) + ... K_ess_Eff * P_cha_e(t_wd) - P_dis_e(t_wd)/K_ess_Eff]; % SOC上下限 Cons = [Cons, 0 <= SOC_ESS(t_wd+1) <= C_ess_Max]; % 充放电功率限制 Cons = [Cons, 0 <= P_cha_e(t_wd) <= bin_cha(t_wd)*P_Cha_Max]; Cons = [Cons, 0 <= P_dis_e(t_wd) <= bin_dis(t_wd)*P_Dis_Max]; Cons = [Cons, bin_cha(t_wd) + bin_dis(t_wd) <= 1]; % 互斥约束 end % ==================== 氢系统约束 ==================== for t_wd = 1:k_win_len % 功率约束 Cons = [Cons, 0 <= P_FC_e(t_wd) <= P_FC_Max_adj]; Cons = [Cons, 0 <= P_e_EL(t_wd) <= P_EL_Max_adj]; % 氢流量计算 (kg/15min) H2_FC = P_FC_e(t_wd) * 0.25 * K_H2P_xs; % 燃料电池耗氢量 H2_EL = P_e_EL(t_wd) * 0.25 * K_P2H_xs; % 电解槽产氢量 % 氢平衡约束 Cons = [Cons, SOC_HSS(t_wd+1) == SOC_HSS(t_wd) - H2_FC + H2_EL]; % SOC上下限 Cons = [Cons, C_HSS_Min <= SOC_HSS(t_wd+1) <= C_HSS_Max]; % 充放氢速率限制 Cons = [Cons, 0 <= H2_FC <= H2_dis_Max]; Cons = [Cons, 0 <= H2_EL <= H2_cha_Max]; end % ==================== 网络约束 ==================== for t_wd = 1:k_win_len t_global = k_start_t + t_wd - 1; % 节点功率平衡 total_injection = G_map * P_GT_e(:, t_wd) + ... W_map * P_wt_cur(t_wd) + ... E_map * (P_dis_e(t_wd) - P_cha_e(t_wd)) + ... H_map * (P_FC_e(t_wd) - P_e_EL(t_wd)); % 基础负荷 load_vector = P_load_96(:, t_global); % 节点净注入功率 P_net = total_injection - load_vector; % 直流潮流约束 Cons = [Cons, Bred * Bus_Theta(non_ref, t_wd) == P_net(non_ref) / baseMVA]; end % 参考节点相角约束 Cons = [Cons, Bus_Theta(ref_node, :) == 0]; % ==================== 目标函数 ==================== Obj = 0; for t_wd = 1:k_win_len % 1. 燃料成本 (乘以0.25转换为15分钟成本) for g = 1:N_Gens fuel_cost = a2(g)*(P_GT_e(g,t_wd)^2) + a1(g)*P_GT_e(g,t_wd) + a0(g)*B_U_GT(g,t_wd); Obj = Obj + 0.25 * fuel_cost; end % 2. 启动成本 for g = 1:N_Gens Obj = Obj + C_Start_qt(g) * B_Start_GT(g, t_wd); end % 3. 碳排放成本 for g = 1:N_Gens carbon_cost = C_carbon * Carbon_Factor(g) * P_GT_e(g,t_wd); Obj = Obj + 0.25 * carbon_cost; end end %% 求解优化问题 fprintf('>> 求解窗口 %d 优化问题...\n', win_idx); Sol = optimize(Cons, Obj, Opts); %% 保存结果 if Sol.problem == 0 % 提取当前窗口结果 Results.P_GT(:, k_start_t:k_end_t) = value(P_GT_e); Results.P_wt(k_start_t:k_end_t) = value(P_wt_cur); Results.P_cha(k_start_t:k_end_t) = value(P_cha_e); Results.P_dis(k_start_t:k_end_t) = value(P_dis_e); Results.SOC_ESS(k_start_t+1:k_end_t+1) = value(SOC_ESS(2:end)); Results.P_FC(k_start_t:k_end_t) = value(P_FC_e); Results.P_EL(k_start_t:k_end_t) = value(P_e_EL); Results.SOC_HSS(k_start_t+1:k_end_t+1) = value(SOC_HSS(2:end)); % 保存机组状态到全局数组 B_ug_96_full(:, k_start_t:k_end_t) = value(B_U_GT); fprintf('>>> 优化成功 | 成本: %.2f 元\n', value(Obj)); else fprintf('>>> 优化失败: %s\n', Sol.info); yalmiperror(Sol.problem); % 尝试部分保存结果 try Results.P_GT(:, k_start_t:k_end_t) = value(P_GT_e); catch warning('无法保存机组出力结果'); end end end %% ==================== 结果后处理与可视化 ==================== fprintf(‘\n====== 优化完成,开始结果分析 ======\n’); % 时间向量 (15分钟间隔) time_vector = linspace(0, 24, T); % 1. 机组出力曲线 figure(‘Name’, ‘机组出力曲线’, ‘Position’, [100, 100, 1200, 600]); for i = 1:N_Gens subplot(2, 3, i); plot(time_vector, Results.P_GT(i, 😃, ‘LineWidth’, 1.5); title(sprintf(‘机组 %d 出力’, i)); xlabel(‘时间 (小时)’); ylabel(‘出力 (MW)’); grid on; xlim([0, 24]); end % 2. 新能源与储能出力 figure(‘Name’, ‘新能源与储能出力’, ‘Position’, [100, 100, 1200, 800]); subplot(3, 1, 1); plot(time_vector, Results.P_wt, ‘b’, ‘LineWidth’, 1.5); title(‘风电出力’); xlabel(‘时间 (小时)’); ylabel(‘出力 (MW)’); grid on; xlim([0, 24]); subplot(3, 1, 2); plot(time_vector, Results.P_dis - Results.P_cha, ‘r’, ‘LineWidth’, 1.5); title(‘电储能净出力’); xlabel(‘时间 (小时)’); ylabel(‘出力 (MW)’); grid on; xlim([0, 24]); subplot(3, 1, 3); plot(time_vector, Results.P_FC - Results.P_EL, ‘g’, ‘LineWidth’, 1.5); title(‘氢系统净出力’); xlabel(‘时间 (小时)’); ylabel(‘出力 (MW)’); grid on; xlim([0, 24]); % 3. 储能状态变化 figure(‘Name’, ‘储能状态变化’, ‘Position’, [100, 100, 1200, 500]); subplot(1, 2, 1); plot(time_vector, Results.SOC_ESS(1:end-1), ‘b’, ‘LineWidth’, 1.5); title(‘电储能SOC’); xlabel(‘时间 (小时)’); ylabel(‘SOC (MWh)’); grid on; xlim([0, 24]); subplot(1, 2, 2); plot(time_vector, Results.SOC_HSS(1:end-1), ‘r’, ‘LineWidth’, 1.5); title(‘氢储能SOC’); xlabel(‘时间 (小时)’); ylabel(‘SOC (kg)’); grid on; xlim([0, 24]); % 4. 系统总成本分析 total_fuel_cost = 0; total_start_cost = 0; total_carbon_cost = 0; for t = 1:T for g = 1:N_Gens % 燃料成本 fuel_cost = a2(g)*(Results.P_GT(g,t)^2) + a1(g)*Results.P_GT(g,t) + a0(g); total_fuel_cost = total_fuel_cost + 0.25 * fuel_cost; % 启动成本 (仅在状态变化时计算) if t > 1 if B_ug_96_full(g,t) == 1 && B_ug_96_full(g,t-1) == 0 total_start_cost = total_start_cost + C_Start_qt(g); end end % 碳排放成本 carbon_cost = C_carbon * Carbon_Factor(g) * Results.P_GT(g,t); total_carbon_cost = total_carbon_cost + 0.25 * carbon_cost; end end total_cost = total_fuel_cost + total_start_cost + total_carbon_cost; % 成本饼图 figure(‘Name’, ‘成本分析’); cost_data = [total_fuel_cost, total_start_cost, total_carbon_cost]; pie(cost_data); legend({‘燃料成本’, ‘启动成本’, ‘碳排放成本’}, ‘Location’, ‘bestoutside’); title(sprintf(‘总成本: %.2f 元’, total_cost)); fprintf(‘====== 结果分析完成 ======\n’); fprintf(’ - 总燃料成本: %.2f 元\n’, total_fuel_cost); fprintf(’ - 总启动成本: %.2f 元\n’, total_start_cost); fprintf(’ - 总碳排放成本: %.2f 元\n’, total_carbon_cost); fprintf(’ - 系统总成本: %.2f 元\n’, total_cost); %% 基于电-氢混合储能的电力系统优化调度模型(日内滚动优化) - 修正版 clc clear %% 参数设置 T = 96; % 时间周期数 (15分钟间隔) N_Gens = 6; % 常规机组数量 W_Size = 16; % 滚动窗口大小 (4小时) R_Step = 4; % 滚动步长 (1小时) N_Windows = floor((T - W_Size)/R_Step) + 1; % 滚动次数 %% 数据加载和预处理 load(‘P_load_24.mat’, ‘P_load_e’); % 基础负荷数据 load(‘P_wt_24’, ‘P_wt_e’); % 风电数据 mpc = loadcase(‘case30’); baseMVA = mpc.baseMVA; % 将24小时数据扩展到96个15分钟时段 P_wt_96 = repelem(P_wt_e, 1, 4); P_load_96 = repelem(P_load_e, 1, 4); % 确保负荷数据为列向量 (30×96) if size(P_load_96, 1) == 1 P_load_96 = P_load_96’; % 转置为列向量 end % 常规机组参数 [机组ID, P_max, P_min, a, b, c, ramp, 碳排放系数] Gen_Params = [1 250 50 0.0375 20 372.5 72 0.85; 2 80 20 0.175 17.5 352.3 48 0.92; 3 50 15 0.625 10 316.5 30 1.05; 4 35 10 0.0834 32.5 329.2 21 0.78; 5 30 10 0.25 30 276.4 18 0.95; 6 40 12 0.25 30 232.2 24 0.88]; % 提取发电机参数 a2 = Gen_Params(:,4); % 二次项系数 a1 = Gen_Params(:,5); % 一次项系数 a0 = Gen_Params(:,6); % 常数项 Ramp_e = Gen_Params(:,7); % 爬坡率 (MW/15min) P_GT_Min = Gen_Params(:,3); % 最小出力 P_GT_Max = Gen_Params(:,2); % 最大出力 % 启动成本参数 (单位: 元) C_Start_qt = [2768; 2874; 2136; 2677; 2262; 2869]; % 电储能参数 P_Cha_Max = 45; % 电储能最大充电功率 P_Dis_Max = 90; % 电储能最大放电功率 K_ess_Eff = 0.92; % 电储能转换效率 C_ess_Max = 300; % 储能容量限制 SOC_ESS_Init = 10; % 初始电SOC (MWh) % 氢系统参数 P_FC_Max = 400; % 燃料电池最大功率 K_H2P_xs = 9.1; % 氢转电系数 (kg/MWh) P_EL_Max = 400; % 电解槽最大功率 K_P2H_xs = 0.22; % 电转氢系数 (kg/MWh) H2_cha_Max = 100; % 最大充氢速率 (kg/15min) H2_dis_Max = 100; % 最大放氢速率 (kg/15min) C_HSS_Max = 400; % 最大储氢容量 (kg) C_HSS_Min = 0; % 最小储氢容量 (kg) SOC_HSS_Init = 20; % 初始氢SOC (kg) % 调整氢系统功率上限 (考虑氢流量约束) P_FC_Max_adj = min(P_FC_Max, H2_dis_Max/(0.25K_H2P_xs)); P_EL_Max_adj = min(P_EL_Max, H2_cha_Max/(0.25K_P2H_xs)); % 节点位置 Bus_Gen = [1 2 5 8 11 13]; % 火电机组节点 Bus_Wt = 2; % 风电节点 Bus_Ess_E = 11; % 电储能节点 Bus_Hss_H = 8; % 氢储能节点 % 环境成本系数 C_carbon = 150; % 碳排放成本 (元/吨CO2) Carbon_Factor = Gen_Params(:,8); % 碳排放系数 (tCO2/MWh) %% 网络拓扑初始化 [Ybus, ~, ~] = makeYbus(mpc); B = imag(Ybus); % 电纳矩阵虚部 ref_node = 1; % 参考节点 non_ref = 2:30; % 非参考节点 Bred = B(non_ref, non_ref); % 降阶B矩阵 % 创建设备-节点映射矩阵 N = size(mpc.bus, 1); G_map = sparse(Bus_Gen, 1:N_Gens, 1, N, N_Gens); % 常规机组 W_map = sparse(Bus_Wt, 1, 1, N, 1); % 风电机组 E_map = sparse(Bus_Ess_E, 1, 1, N, 1); % 电储能 H_map = sparse(Bus_Hss_H, 1, 1, N, 1); % 氢系统 % 初始化t=0时刻状态 (确保列向量格式) B_ug_z0 = ones(N_Gens,1); % 初始机组状态 (全运行) P_GT_z0 = P_GT_Min; % 初始机组出力 (最小出力) %% 主循环滚动优化 Results = struct(); Results.P_GT = zeros(N_Gens, T); Results.P_wt = zeros(1, T); Results.P_cha = zeros(1, T); Results.P_dis = zeros(1, T); Results.SOC_ESS = SOC_ESS_Init * ones(1, T+1); Results.P_FC = zeros(1, T); Results.P_EL = zeros(1, T); Results.SOC_HSS = SOC_HSS_Init * ones(1, T+1); B_ug_96_full = zeros(N_Gens, T); % 全时段机组状态 (列向量存储) B_ug_96_full(:,1) = B_ug_z0; % 设置初始机组状态 disp(‘====== 开始滚动优化 ======’); for win_idx = 1:N_Windows k_start_t = (win_idx-1)*R_Step + 1; k_end_t = min(k_start_t + W_Size - 1, T); % 防止越界 k_win_len = k_end_t - k_start_t + 1; % 实际窗口长度 fprintf('>> 处理窗口 %d: 时段 %d-%d (长度=%d)\n', ... win_idx, k_start_t, k_end_t, k_win_len); %% 获取当前窗口初始状态 % 电储能初始SOC if k_start_t == 1 SOC_ESS_init = SOC_ESS_Init; else SOC_ESS_init = Results.SOC_ESS(k_start_t); end % 氢储能初始SOC if k_start_t == 1 SOC_HSS_init = SOC_HSS_Init; else SOC_HSS_init = Results.SOC_HSS(k_start_t); end %% 定义当前窗口变量 % 机组变量 P_GT_e = sdpvar(N_Gens, k_win_len, 'full'); % 火电机组出力 B_U_GT = binvar(N_Gens, k_win_len, 'full'); % 机组运行状态 B_Start_GT = binvar(N_Gens, k_win_len, 'full'); % 启动状态 % 可再生能源 P_wt_cur = sdpvar(1, k_win_len); % 风电出力 % 电储能 P_cha_e = sdpvar(1, k_win_len); % 充电功率 P_dis_e = sdpvar(1, k_win_len); % 放电功率 SOC_ESS = sdpvar(1, k_win_len+1); % SOC状态 % 氢系统 P_FC_e = sdpvar(1, k_win_len); % 燃料电池出力 P_e_EL = sdpvar(1, k_win_len); % 电解槽用电 SOC_HSS = sdpvar(1, k_win_len+1); % 氢SOC状态 % 网络变量 Bus_Theta = sdpvar(30, k_win_len, 'full'); % 节点相角 %% 约束初始化 Cons = []; % ==================== 初始状态约束 ==================== Cons = [Cons, SOC_ESS(1) == SOC_ESS_init]; % 电储能初始状态 Cons = [Cons, SOC_HSS(1) == SOC_HSS_init]; % 氢储能初始状态 % ==================== 机组约束 ==================== for t_wd = 1:k_win_len t_global = k_start_t + t_wd - 1; % 获取前一时刻状态 if t_global == 1 prev_status = B_ug_z0; % N_Gens×1 列向量 prev_power = P_GT_z0; % N_Gens×1 列向量 else prev_status = B_ug_96_full(:, t_global-1); % N_Gens×1 列向量 prev_power = Results.P_GT(:, t_global-1); % N_Gens×1 列向量 end % 机组出力上下限 Cons = [Cons, B_U_GT(:, t_wd).*P_GT_Min <= P_GT_e(:, t_wd) <= B_U_GT(:, t_wd).*P_GT_Max]; % 启动状态逻辑 (统一使用列向量) if t_global == 1 Cons = [Cons, B_Start_GT(:, t_wd) == B_U_GT(:, t_wd)]; else Cons = [Cons, B_Start_GT(:, t_wd) >= B_U_GT(:, t_wd) - prev_status]; Cons = [Cons, B_Start_GT(:, t_wd) <= B_U_GT(:, t_wd)]; Cons = [Cons, B_Start_GT(:, t_wd) <= 1 - prev_status]; end % 爬坡约束 (统一使用列向量计算) if t_global >= 1 % 上爬坡约束 Cons = [Cons, P_GT_e(:, t_wd) - prev_power <= Ramp_e .* prev_status + ... (B_U_GT(:, t_wd) - prev_status) .* P_GT_Min]; % 下爬坡约束 Cons = [Cons, prev_power - P_GT_e(:, t_wd) <= Ramp_e .* B_U_GT(:, t_wd) + ... (prev_status - B_U_GT(:, t_wd)) .* P_GT_Min]; end end % ==================== 风电机组约束 ==================== Cons = [Cons, 0 <= P_wt_cur <= P_wt_96(k_start_t:k_end_t)]; % ==================== 电储能约束 ==================== % 添加二进制变量实现充放电互斥 bin_cha = binvar(1, k_win_len); % 充电状态指示器 bin_dis = binvar(1, k_win_len); % 放电状态指示器 for t_wd = 1:k_win_len % SOC动态更新 Cons = [Cons, SOC_ESS(t_wd+1) == SOC_ESS(t_wd) + ... K_ess_Eff * P_cha_e(t_wd) - P_dis_e(t_wd)/K_ess_Eff]; % SOC上下限 Cons = [Cons, 0 <= SOC_ESS(t_wd+1) <= C_ess_Max]; % 充放电功率限制 (使用二进制变量实现互斥) Cons = [Cons, 0 <= P_cha_e(t_wd) <= bin_cha(t_wd)*P_Cha_Max]; Cons = [Cons, 0 <= P_dis_e(t_wd) <= bin_dis(t_wd)*P_Dis_Max]; Cons = [Cons, bin_cha(t_wd) + bin_dis(t_wd) <= 1]; % 互斥约束 end % ==================== 氢系统约束 ==================== for t_wd = 1:k_win_len % 功率约束 (使用调整后的上限) Cons = [Cons, 0 <= P_FC_e(t_wd) <= P_FC_Max_adj]; Cons = [Cons, 0 <= P_e_EL(t_wd) <= P_EL_Max_adj]; % 氢流量计算 (kg/15min) H2_FC = P_FC_e(t_wd) * 0.25 * K_H2P_xs; % 燃料电池耗氢量 H2_EL = P_e_EL(t_wd) * 0.25 * K_P2H_xs; % 电解槽产氢量 % 氢平衡约束 Cons = [Cons, SOC_HSS(t_wd+1) == SOC_HSS(t_wd) - H2_FC + H2_EL]; % SOC上下限 Cons = [Cons, C_HSS_Min <= SOC_HSS(t_wd+1) <= C_HSS_Max]; % 充放氢速率限制 Cons = [Cons, 0 <= H2_FC <= H2_dis_Max]; Cons = [Cons, 0 <= H2_EL <= H2_cha_Max]; end % ==================== 网络约束 ==================== for t_wd = 1:k_win_len t_global = k_start_t + t_wd - 1; % 节点功率平衡 total_injection = G_map * P_GT_e(:, t_wd) + ... W_map * P_wt_cur(t_wd) + ... E_map * (P_dis_e(t_wd) - P_cha_e(t_wd)) + ... H_map * (P_FC_e(t_wd) - P_e_EL(t_wd)); % 基础负荷 (确保维度匹配) load_vector = P_load_96(:, t_global); % 节点净注入功率 P_net = total_injection - load_vector; % 直流潮流约束 Cons = [Cons, Bred * Bus_Theta(non_ref, t_wd) == P_net(non_ref) / baseMVA]; end % 参考节点相角约束 Cons = [Cons, Bus_Theta(ref_node, :) == 0]; % ==================== 目标函数 ==================== Obj = 0; for t_wd = 1:k_win_len % 1. 燃料成本 (乘以0.25转换为15分钟成本) for g = 1:N_Gens fuel_cost = a2(g)*(P_GT_e(g,t_wd)^2) + a1(g)*P_GT_e(g,t_wd) + a0(g)*B_U_GT(g,t_wd); Obj = Obj + 0.25 * fuel_cost; end % 2. 启动成本 (单次启动成本,无需时间转换) for g = 1:N_Gens Obj = Obj + C_Start_qt(g) * B_Start_GT(g, t_wd); end % 3. 碳排放成本 (乘以0.25转换为15分钟成本) for g = 1:N_Gens carbon_cost = C_carbon * Carbon_Factor(g) * P_GT_e(g,t_wd); Obj = Obj + 0.25 * carbon_cost; end end %% CPLEX求解优化问题 Opts = sdpsettings('solver', 'cplex', 'verbose', 1, 'showprogress', 1); Sol = optimize(Cons, Obj, Opts); %% 保存结果 if Sol.problem == 0 % 提取当前窗口结果 Results.P_GT(:, k_start_t:k_end_t) = value(P_GT_e); Results.P_wt(k_start_t:k_end_t) = value(P_wt_cur); Results.P_cha(k_start_t:k_end_t) = value(P_cha_e); Results.P_dis(k_start_t:k_end_t) = value(P_dis_e); Results.SOC_ESS(k_start_t+1:k_end_t+1) = value(SOC_ESS(2:end)); Results.P_FC(k_start_t:k_end_t) = value(P_FC_e); Results.P_EL(k_start_t:k_end_t) = value(P_e_EL); Results.SOC_HSS(k_start_t+1:k_end_t+1) = value(SOC_HSS(2:end)); % 保存机组状态到全局数组 B_ug_96_full(:, k_start_t:k_end_t) = value(B_U_GT); fprintf('>>>窗口 %d/%d 优化完成 | 成本: %.2f 元\n', win_idx, N_Windows, value(Obj)); % 调试输出:检查状态传递 fprintf(' 电储能初始SOC: %.2f -> 结束SOC: %.2f\n', ... SOC_ESS_init, value(SOC_ESS(end))); fprintf(' 氢储能初始SOC: %.2f -> 结束SOC: %.2f\n', ... SOC_HSS_init, value(SOC_HSS(end))); else fprintf('>>>窗口 %d 优化失败: %s\n', win_idx, Sol.info); yalmiperror(Sol.problem); error('优化失败,请检查约束条件'); end end %% 结果后处理与可视化 Total_Load = sum(P_load_96, 1); % 系统总负荷 (1×96) % 创建时间标签 (15分钟间隔) time_labels = cell(1, T); for h = 0:23 for q = 1:4 idx = h*4 + q; time_labels{idx} = sprintf(‘%02d:%02d’, h, (q-1)*15); end end %% 功率平衡图 figure(‘Name’, ‘日内调度功率平衡’, ‘Position’, [100, 100, 1200, 600]); subplot(2,1,1); % 计算各组件贡献 total_generation = sum(Results.P_GT, 1) + Results.P_wt + Results.P_dis + Results.P_FC; net_storage = Results.P_dis - Results.P_cha + Results.P_FC - Results.P_EL; % 创建堆叠区域图 components = [ sum(Results.P_GT, 1); Results.P_wt; Results.P_dis; Results.P_FC; -Results.P_cha; -Results.P_EL ]'; area(1:T, components(:,1:5), ‘FaceAlpha’, 0.7); hold on; plot(1:T, Total_Load, ‘k-’, ‘LineWidth’, 2); plot(1:T, total_generation - Results.P_cha - Results.P_EL, ‘r–’, ‘LineWidth’, 1.5); % 设置图形属性 legend(‘常规机组’, ‘风电’, ‘储能放电’, ‘燃料电池’, ‘储能充电’, ‘电解槽’, ‘总负荷’, ‘净发电量’, … ‘Location’, ‘northoutside’, ‘Orientation’, ‘horizontal’); title(‘日内调度功率平衡’); xlabel(‘时间 (15分钟间隔)’); ylabel(‘功率 (MW)’); xlim([1 T]); xticks(1:4:T); xticklabels(time_labels(1:4:end)); grid on; %% 储能状态图 subplot(2,1,2); yyaxis left; plot(1:T, Results.SOC_ESS(2:end), ‘b-o’, ‘LineWidth’, 1.5); ylabel(‘电储能SOC (MWh)’); hold on; plot(1:T, Results.SOC_HSS(2:end), ‘r-s’, ‘LineWidth’, 1.5); ylabel(‘氢储能SOC (kg)’); yyaxis right; bar(1:T, Results.P_cha, ‘FaceColor’, [0 0.5 0], ‘FaceAlpha’, 0.5); hold on; bar(1:T, -Results.P_dis, ‘FaceColor’, [0.5 0 0], ‘FaceAlpha’, 0.5); ylabel(‘储能充放电功率 (MW)’); % 设置图形属性 legend(‘电储能SOC’, ‘氢储能SOC’, ‘充电功率’, ‘放电功率’, … ‘Location’, ‘northoutside’, ‘Orientation’, ‘horizontal’); title(‘储能系统状态’); xlabel(‘时间 (15分钟间隔)’); xlim([1 T]); xticks(1:4:T); xticklabels(time_labels(1:4:end)); grid on; %% 机组出力详情图 figure(‘Name’, ‘机组出力详情’, ‘Position’, [100, 100, 1200, 800]); % 机组出力堆叠图 subplot(2,1,1); plot_data = Results.P_GT’; area(1:T, plot_data, ‘FaceAlpha’, 0.7); hold on; plot(1:T, sum(Results.P_GT, 1), ‘k-’, ‘LineWidth’, 1.5); % 设置图形属性 legend(arrayfun(@(x) sprintf(‘机组%d’, x), 1:N_Gens, ‘UniformOutput’, false), … ‘总出力’, ‘Location’, ‘best’); title(‘机组出力分布’); xlabel(‘时间 (15分钟间隔)’); ylabel(‘出力 (MW)’); xlim([1 T]); xticks(1:4:T); xticklabels(time_labels(1:4:end)); grid on; % 机组启停状态热图 subplot(2,1,2); commitment = zeros(N_Gens, T); for g = 1:N_Gens commitment(g, 😃 = (Results.P_GT(g, 😃 > 0); end imagesc(commitment); colormap([0.9 0.9 0.9; 0.2 0.6 0.2]); % 灰色-停机, 绿色-运行 % 设置坐标轴 colorbar(‘Ticks’, [0.25, 0.75], ‘TickLabels’, {‘停机’, ‘运行’}); set(gca, ‘YTick’, 1:N_Gens, … ‘YTickLabel’, arrayfun(@(x) sprintf(‘机组%d’, x), 1:N_Gens, ‘UniformOutput’, false), … ‘XTick’, 1:4:T, … ‘XTickLabel’, time_labels(1:4:end)); title(‘机组启停状态’); xlabel(‘时间 (15分钟间隔)’); ylabel(‘机组编号’); %% 氢系统运行状态 figure(‘Name’, ‘氢系统运行状态’, ‘Position’, [100, 100, 1200, 600]); % 氢流量图 subplot(1,2,1); H2_FC_flow = Results.P_FC * 0.25 * K_H2P_xs; % kg/15min H2_EL_flow = Results.P_EL * 0.25 * K_P2H_xs; % kg/15min bar_data = [H2_EL_flow; -H2_FC_flow]'; bar(1:T, bar_data, ‘stacked’); legend(‘电解槽产氢’, ‘燃料电池耗氢’, ‘Location’, ‘best’); title(‘氢流量平衡’); xlabel(‘时间 (15分钟间隔)’); ylabel(‘氢流量 (kg/15min)’); xlim([1 T]); xticks(1:4:T); xticklabels(time_labels(1:4:end)); grid on; % 储氢状态图 subplot(1,2,2); yyaxis left; plot(1:T, Results.SOC_HSS(2:end), ‘b-o’, ‘LineWidth’, 1.5); ylabel(‘储氢量 (kg)’); yyaxis right; bar(1:T, Results.P_FC, ‘FaceColor’, [0.5 0 0.5], ‘FaceAlpha’, 0.5); hold on; bar(1:T, -Results.P_EL, ‘FaceColor’, [0 0.5 0.5], ‘FaceAlpha’, 0.5); ylabel(‘功率 (MW)’); legend(‘储氢量’, ‘燃料电池出力’, ‘电解槽用电’, ‘Location’, ‘best’); title(‘氢系统功率与储量’); xlabel(‘时间 (15分钟间隔)’); xlim([1 T]); xticks(1:4:T); xticklabels(time_labels(1:4:end)); grid on; disp(‘====== 日内滚动优化完成 ======’); %% 基于电-氢混合储能的电力系统优化调度模型(日内滚动优化) clc clear %% 参数设置 T = 96; % 时间周期数 (15分钟间隔) N_Gens = 6; % 常规机组数量 W_Size = 16; % 滚动窗口大小 (4小时) R_Step = 4; % 滚动步长 (1小时) N_Windows = floor((T - W_Size)/R_Step) + 1; % 滚动次数 %% 数据加载和预处理 load(‘P_load_24.mat’, ‘P_load_e’); % 基础负荷数据 load(‘P_wt_24’, ‘P_wt_e’); % 风电数据 mpc = loadcase(‘case30’); baseMVA = mpc.baseMVA; % 将24小时数据扩展到96个15分钟时段 P_wt_96 = repelem(P_wt_e, 1, 4); P_load_96 = repelem(P_load_e, 1, 4); % 常规机组参数 [机组ID, P_max, P_min, a, b, c, ramp, 碳排放系数] Gen_Params = [1 250 50 0.0375 20 372.5 72 0.85; 2 80 20 0.175 17.5 352.3 48 0.92; 3 50 15 0.625 10 316.5 30 1.05; 4 35 10 0.0834 32.5 329.2 21 0.78; 5 30 10 0.25 30 276.4 18 0.95; 6 40 12 0.25 30 232.2 24 0.88]; % 提取发电机参数 a2 = Gen_Params(:,4); % 二次项系数 a1 = Gen_Params(:,5); % 一次项系数 a0 = Gen_Params(:,6); % 常数项 Ramp_e = Gen_Params(:,7); % 爬坡率 (MW/15min) P_GT_Min = Gen_Params(:,3); % 最小出力 P_GT_Max = Gen_Params(:,2); % 最大出力 % 启动成本参数 (单位: 元) C_Start_qt = [2768; 2874; 2136; 2677; 2262; 2869]; % 电储能参数 P_Cha_Max = 45; % 电储能最大充电功率 P_Dis_Max = 90; % 电储能最大放电功率 K_ess_Eff = 0.92; % 电储能转换效率 C_ess_Max = 300; % 储能容量限制 SOC_ESS_Init = 10; % 初始电SOC (MWh) % 氢系统参数 P_FC_Max = 400; % 燃料电池最大功率 K_H2P_xs = 9.1; % 氢转电系数 (kg/MWh) P_EL_Max = 400; % 电解槽最大功率 K_P2H_xs = 0.22; % 电转氢系数 (kg/MWh) H2_cha_Max = 100; % 最大充氢速率 (kg/15min) H2_dis_Max = 100; % 最大放氢速率 (kg/15min) C_HSS_Max = 400; % 最大储氢容量 (kg) C_HSS_Min = 0; % 最小储氢容量 (kg) SOC_HSS_Init = 20; % 初始氢SOC (kg) % 调整氢系统功率上限 (考虑氢流量约束) P_FC_Max_adj = min(P_FC_Max, H2_dis_Max/(0.25K_H2P_xs)); P_EL_Max_adj = min(P_EL_Max, H2_cha_Max/(0.25K_P2H_xs)); % 节点位置 Bus_Gen = [1 2 5 8 11 13]; % 火电机组节点 Bus_Wt = 2; % 风电节点 Bus_Ess_E = 11; % 电储能节点 Bus_Hss_H = 8; % 氢储能节点 % 环境成本系数 C_carbon = 150; % 碳排放成本 (元/吨CO2) Carbon_Factor = Gen_Params(:,8); % 碳排放系数 (tCO2/MWh) %% 网络拓扑初始化 [Ybus, ~, ~] = makeYbus(mpc); B = imag(Ybus); % 电纳矩阵虚部 ref_node = 1; % 参考节点 non_ref = 2:30; % 非参考节点 Bred = B(non_ref, non_ref); % 降阶B矩阵 % 创建设备-节点映射矩阵 N = size(mpc.bus, 1); G_map = sparse(Bus_Gen, 1:N_Gens, 1, N, N_Gens); % 常规机组 W_map = sparse(Bus_Wt, 1, 1, N, 1); % 风电机组 E_map = sparse(Bus_Ess_E, 1, 1, N, 1); % 电储能 H_map = sparse(Bus_Hss_H, 1, 1, N, 1); % 氢系统 % 初始化t=0时刻状态 (确保列向量格式) B_ug_z0 = ones(N_Gens,1); % 初始机组状态 (全运行) P_GT_z0 = P_GT_Min; % 初始机组出力 (最小出力) %% 主循环滚动优化 Results = struct(); Results.P_GT = zeros(N_Gens, T); Results.P_wt = zeros(1, T); Results.P_cha = zeros(1, T); Results.P_dis = zeros(1, T); Results.SOC_ESS = SOC_ESS_Init * ones(1, T+1); Results.P_FC = zeros(1, T); Results.P_EL = zeros(1, T); Results.SOC_HSS = SOC_HSS_Init * ones(1, T+1); B_ug_96_full = zeros(N_Gens, T); % 全时段机组状态 (列向量存储) disp(‘====== 开始滚动优化 ======’); for win_idx = 1:N_Windows k_start_t = (win_idx-1)*R_Step + 1; k_end_t = min(k_start_t + W_Size - 1, T); k_win_len = k_end_t - k_start_t + 1; %% 定义当前窗口变量 % 机组变量 P_GT_e = sdpvar(N_Gens, k_win_len, 'full'); % 火电机组出力 B_U_GT = binvar(N_Gens, k_win_len, 'full'); % 机组运行状态 B_Start_GT = binvar(N_Gens, k_win_len, 'full'); % 启动状态 % 可再生能源 P_wt_e = sdpvar(1, k_win_len); % 风电出力 % 电储能 P_cha_e = sdpvar(1, k_win_len); % 充电功率 P_dis_e = sdpvar(1, k_win_len); % 放电功率 SOC_ESS = sdpvar(1, k_win_len+1); % SOC状态 % 氢系统 P_FC_e = sdpvar(1, k_win_len); % 燃料电池出力 P_e_EL = sdpvar(1, k_win_len); % 电解槽用电 SOC_HSS = sdpvar(1, k_win_len+1); % 氢SOC状态 % 网络变量 Bus_Theta = sdpvar(30, k_win_len, 'full'); % 节点相角 %% 约束初始化 Cons = []; % 机组约束 for t_wd = 1:k_win_len t_global = k_start_t + t_wd - 1; % 获取前一时刻状态 (保持列向量格式) if t_global == 1 prev_status = B_ug_z0; % N_Gens×1 列向量 prev_power = P_GT_z0; % N_Gens×1 列向量 else prev_status = B_ug_96_full(:, t_global-1); % N_Gens×1 列向量 prev_power = Results.P_GT(:, t_global-1); % N_Gens×1 列向量 end % 机组出力上下限 Cons = [Cons, B_U_GT(:, t_wd).*P_GT_Min <= P_GT_e(:, t_wd) <= B_U_GT(:, t_wd).*P_GT_Max]; % 启动状态逻辑 (统一使用列向量) if t_global == 1 Cons = [Cons, B_Start_GT(:, t_wd) == B_U_GT(:, t_wd)]; else Cons = [Cons, B_Start_GT(:, t_wd) >= B_U_GT(:, t_wd) - prev_status]; Cons = [Cons, B_Start_GT(:, t_wd) <= B_U_GT(:, t_wd)]; Cons = [Cons, B_Start_GT(:, t_wd) <= 1 - prev_status]; end % 爬坡约束 (统一使用列向量计算) if t_global >= 1 % 上爬坡约束 Cons = [Cons, P_GT_e(:, t_wd) - prev_power <= Ramp_e .* prev_status + ... (B_U_GT(:, t_wd) - prev_status) .* P_GT_Min]; % 下爬坡约束 Cons = [Cons, prev_power - P_GT_e(:, t_wd) <= Ramp_e .* B_U_GT(:, t_wd) + ... (prev_status - B_U_GT(:, t_wd)) .* P_GT_Min]; end end % 风电机组约束 Cons = [Cons, 0 <= P_wt_e <= P_wt_96(k_start_t:k_end_t)]; % 电储能约束 Cons = [Cons, SOC_ESS(1) == Results.SOC_ESS(k_start_t)]; % SOC衔接 % 添加二进制变量实现充放电互斥 bin_cha = binvar(1, k_win_len); % 充电状态指示器 bin_dis = binvar(1, k_win_len); % 放电状态指示器 for t_wd = 1:k_win_len % SOC动态更新 Cons = [Cons, SOC_ESS(t_wd+1) == SOC_ESS(t_wd) + ... K_ess_Eff * P_cha_e(t_wd) - P_dis_e(t_wd)/K_ess_Eff]; % SOC上下限 Cons = [Cons, 0 <= SOC_ESS(t_wd+1) <= C_ess_Max]; % 充放电功率限制 (使用二进制变量实现互斥) Cons = [Cons, 0 <= P_cha_e(t_wd) <= bin_cha(t_wd)*P_Cha_Max]; Cons = [Cons, 0 <= P_dis_e(t_wd) <= bin_dis(t_wd)*P_Dis_Max]; Cons = [Cons, bin_cha(t_wd) + bin_dis(t_wd) <= 1]; % 互斥约束 end % 氢系统约束 Cons = [Cons, SOC_HSS(1) == Results.SOC_HSS(k_start_t)]; % SOC衔接 for t_wd = 1:k_win_len % 功率约束 (使用调整后的上限) Cons = [Cons, 0 <= P_FC_e(t_wd) <= P_FC_Max_adj]; Cons = [Cons, 0 <= P_e_EL(t_wd) <= P_EL_Max_adj]; % 氢流量计算 (kg/15min) H2_FC = P_FC_e(t_wd) * 0.25 * K_H2P_xs; % 燃料电池耗氢量 H2_EL = P_e_EL(t_wd) * 0.25 * K_P2H_xs; % 电解槽产氢量 % 氢平衡约束 Cons = [Cons, SOC_HSS(t_wd+1) == SOC_HSS(t_wd) - H2_FC + H2_EL]; % SOC上下限 Cons = [Cons, C_HSS_Min <= SOC_HSS(t_wd+1) <= C_HSS_Max]; % 充放氢速率限制 Cons = [Cons, 0 <= H2_FC <= H2_dis_Max]; Cons = [Cons, 0 <= H2_EL <= H2_cha_Max]; end % 网络约束 for t_wd = 1:k_win_len t_global = k_start_t + t_wd - 1; % 节点功率平衡 total_injection = G_map * P_GT_e(:, t_wd) + ... W_map * P_wt_e(t_wd) + ... E_map * (P_dis_e(t_wd) - P_cha_e(t_wd)) + ... H_map * (P_FC_e(t_wd) - P_e_EL(t_wd)); % 基础负荷 load_vector = P_load_96(:, t_global); % 节点净注入功率 P_net = total_injection - load_vector; % 直流潮流约束 Cons = [Cons, Bred * Bus_Theta(non_ref, t_wd) == P_net(non_ref) / baseMVA]; end % 参考节点相角约束 Cons = [Cons, Bus_Theta(ref_node, :) == 0]; % 目标函数 Obj = 0; for t_wd = 1:k_win_len % 1. 燃料成本 for g = 1:N_Gens Obj = Obj + a2(g)*(P_GT_e(g,t_wd)^2) + a1(g)*P_GT_e(g,t_wd) + a0(g)*B_U_GT(g,t_wd); end % 2. 启动成本 for g = 1:N_Gens Obj = Obj + C_Start_qt(g) * B_Start_GT(g, t_wd); end % 3. 碳排放成本 for g = 1:N_Gens Obj = Obj + C_carbon * Carbon_Factor(g) * P_GT_e(g,t_wd) * 0.25; % 转换为15分钟成本 end end %% CPLEX求解优化问题 Opts = sdpsettings('solver', 'cplex', 'verbose', 1, 'showprogress', 1); Sol = optimize(Cons, Obj, Opts); %% 保存结果 if Sol.problem == 0 % 提取当前窗口结果 Results.P_GT(:, k_start_t:k_end_t) = value(P_GT_e); Results.P_wt(k_start_t:k_end_t) = value(P_wt_e); Results.P_cha(k_start_t:k_end_t) = value(P_cha_e); Results.P_dis(k_start_t:k_end_t) = value(P_dis_e); Results.SOC_ESS(k_start_t+1:k_end_t+1) = value(SOC_ESS(2:end)); Results.P_FC(k_start_t:k_end_t) = value(P_FC_e); Results.P_EL(k_start_t:k_end_t) = value(P_e_EL); Results.SOC_HSS(k_start_t+1:k_end_t+1) = value(SOC_HSS(2:end)); B_ug_96_full(:, k_start_t:k_end_t) = value(B_U_GT); % 保存机组状态 fprintf('>>>窗口 %d/%d 优化完成 | 成本: %.2f 元\n', win_idx, N_Windows, value(Obj)); else error('>>>窗口 %d 优化失败: %s', win_idx, Sol.info); end end %% 结果后处理与可视化 Total_Load = P_load_96; % 仅基础负荷 % 创建时间标签 (15分钟间隔) time_labels = cell(1, T); for h = 0:23 for q = 1:4 idx = h*4 + q; time_labels{idx} = sprintf(‘%02d:%02d’, h, (q-1)*15); end end %% 功率平衡图 figure(‘Name’, ‘日内调度功率平衡’, ‘Position’, [100, 100, 1200, 600]); subplot(2,1,1); % 创建堆叠区域图 components = [ sum(Results.P_GT, 1); Results.P_wt; Results.P_dis; Results.P_FC; -Results.P_cha; -Results.P_EL ]'; area(1:T, components(:,1:5), ‘FaceAlpha’, 0.7); hold on; plot(1:T, Total_Load, ‘k-’, ‘LineWidth’, 2); plot(1:T, sum(components(:,1:4), 2) - components(:,5) - components(:,6), ‘r–’, ‘LineWidth’, 1.5); % 设置图形属性 legend(‘常规机组’, ‘风电’, ‘储能放电’, ‘燃料电池’, ‘储能充电’, ‘电解槽’, ‘总负荷’, ‘净发电量’, … ‘Location’, ‘northoutside’, ‘Orientation’, ‘horizontal’); title(‘日内调度功率平衡’); xlabel(‘时间 (15分钟间隔)’); ylabel(‘功率 (MW)’); xlim([1 T]); xticks(1:4:T); xticklabels(time_labels(1:4:end)); grid on; %% 储能状态图 subplot(2,1,2); yyaxis left; plot(1:T, Results.SOC_ESS(2:end), ‘b-o’, ‘LineWidth’, 1.5); ylabel(‘电储能SOC (MWh)’); hold on; plot(1:T, Results.SOC_HSS(2:end), ‘r-s’, ‘LineWidth’, 1.5); ylabel(‘氢储能SOC (kg)’); yyaxis right; bar(1:T, Results.P_cha, ‘FaceColor’, [0 0.5 0], ‘FaceAlpha’, 0.5); hold on; bar(1:T, -Results.P_dis, ‘FaceColor’, [0.5 0 0], ‘FaceAlpha’, 0.5); ylabel(‘储能充放电功率 (MW)’); % 设置图形属性 legend(‘电储能SOC’, ‘氢储能SOC’, ‘充电功率’, ‘放电功率’, … ‘Location’, ‘northoutside’, ‘Orientation’, ‘horizontal’); title(‘储能系统状态’); xlabel(‘时间 (15分钟间隔)’); xlim([1 T]); xticks(1:4:T); xticklabels(time_labels(1:4:end)); grid on; %% 机组出力详情图 figure(‘Name’, ‘机组出力详情’, ‘Position’, [100, 100, 1200, 800]); % 机组出力堆叠图 subplot(2,1,1); plot_data = Results.P_GT’; area(1:T, plot_data, ‘FaceAlpha’, 0.7); hold on; plot(1:T, sum(Results.P_GT, 1), ‘k-’, ‘LineWidth’, 1.5); % 设置图形属性 legend(arrayfun(@(x) sprintf(‘机组%d’, x), 1:N_Gens, ‘UniformOutput’, false), … ‘总出力’, ‘Location’, ‘best’); title(‘机组出力分布’); xlabel(‘时间 (15分钟间隔)’); ylabel(‘出力 (MW)’); xlim([1 T]); xticks(1:4:T); xticklabels(time_labels(1:4:end)); grid on; % 机组启停状态热图 subplot(2,1,2); commitment = zeros(N_Gens, T); for g = 1:N_Gens commitment(g, 😃 = (Results.P_GT(g, 😃 > 0); end imagesc(commitment); colormap([0.9 0.9 0.9; 0.2 0.6 0.2]); % 灰色-停机, 绿色-运行 % 设置坐标轴 colorbar(‘Ticks’, [0.25, 0.75], ‘TickLabels’, {‘停机’, ‘运行’}); set(gca, ‘YTick’, 1:N_Gens, … ‘YTickLabel’, arrayfun(@(x) sprintf(‘机组%d’, x), 1:N_Gens, ‘UniformOutput’, false), … ‘XTick’, 1:4:T, … ‘XTickLabel’, time_labels(1:4:end)); title(‘机组启停状态’); xlabel(‘时间 (15分钟间隔)’); ylabel(‘机组编号’); %% 氢系统运行状态 figure(‘Name’, ‘氢系统运行状态’, ‘Position’, [100, 100, 1200, 600]); % 氢流量图 subplot(1,2,1); bar_data = [Results.P_EL; -Results.P_FC]'; bar(1:T, bar_data, ‘stacked’); legend(‘电解槽产氢’, ‘燃料电池耗氢’, ‘Location’, ‘best’); title(‘氢流量平衡’); xlabel(‘时间 (15分钟间隔)’); ylabel(‘氢流量 (kg/15min)’); xlim([1 T]); xticks(1:4:T); xticklabels(time_labels(1:4:end)); grid on; % 储氢状态图 subplot(1,2,2); yyaxis left; plot(1:T, Results.SOC_HSS(2:end), ‘b-o’, ‘LineWidth’, 1.5); ylabel(‘储氢量 (kg)’); yyaxis right; bar(1:T, Results.P_FC, ‘FaceColor’, [0.5 0 0.5], ‘FaceAlpha’, 0.5); hold on; bar(1:T, -Results.P_EL, ‘FaceColor’, [0 0.5 0.5], ‘FaceAlpha’, 0.5); ylabel(‘功率 (MW)’); legend(‘储氢量’, ‘燃料电池出力’, ‘电解槽用电’, ‘Location’, ‘best’); title(‘氢系统功率与储量’); xlabel(‘时间 (15分钟间隔)’); xlim([1 T]); xticks(1:4:T); xticklabels(time_labels(1:4:end)); grid on; disp(‘====== 日内滚动优化完成 ======’);
11-27
S3DIS 填写此Google 表单下载 S3DIS 数据。下载Stanford3dDataset_v1.2.zip文件并解压。 修复Area_5/office_19/Annotations/ceiling第 323474 行中的错误 (103.0�0000 => 103.000000)。 (可选)从这里下载完整的 2D-3D S3DIS 数据集(无 XYZ)以解析正常数据。 运行 S3DIS 的预处理代码如下: # S3DIS_DIR: the directory of downloaded Stanford3dDataset_v1.2 dataset. # RAW_S3DIS_DIR: the directory of Stanford2d3dDataset_noXYZ dataset. (optional, for parsing normal) # PROCESSED_S3DIS_DIR: the directory of processed S3DIS dataset (output dir). # S3DIS without aligned angle python pointcept/datasets/preprocessing/s3dis/preprocess_s3dis.py --dataset_root ${S3DIS_DIR} --output_root ${PROCESSED_S3DIS_DIR} # S3DIS with aligned angle python pointcept/datasets/preprocessing/s3dis/preprocess_s3dis.py --dataset_root ${S3DIS_DIR} --output_root ${PROCESSED_S3DIS_DIR} --align_angle # S3DIS with normal vector (recommended, normal is helpful) python pointcept/datasets/preprocessing/s3dis/preprocess_s3dis.py --dataset_root ${S3DIS_DIR} --output_root ${PROCESSED_S3DIS_DIR} --raw_root ${RAW_S3DIS_DIR} --parse_normal python pointcept/datasets/preprocessing/s3dis/preprocess_s3dis.py --dataset_root ${S3DIS_DIR} --output_root ${PROCESSED_S3DIS_DIR} --raw_root ${RAW_S3DIS_DIR} --align_angle --parse_normal (可选)我们的预处理数据也可以从 [这里]下载(带有法线矢量和对齐角度),下载前请同意官方许可。 将处理过的数据集链接到代码库。 # PROCESSED_S3DIS_DIR: the directory of processed S3DIS dataset. mkdir data ln -s ${PROCESSED_S3DIS_DIR} ${CODEBASE_DIR}/data/s3dis
09-03
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值