代码窗口最大化_工具推荐：Carbon 代码图片生成器

最新推荐文章于 2025-08-12 09:22:04 发布
原创 最新推荐文章于 2025-08-12 09:22:04 发布 · 238 阅读 · 0 ·
CC 4.0 BY-SA版权
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
文章标签：

#代码窗口最大化

部署运行你感兴趣的模型镜像

694aa170a5eec53998ee8d3efe9a2075.png

Carbon 是一款生成精美代码片段图片的工具。有多精美？美轮美奂。

之前在逛推特的时候，偶尔能看到一些类似下面这样的代码片段图片，每次看到都觉得特别精致。

d0e95cb459eeb1cfbab4212948bc0b30.png

有一次也特别想要这样的图片，无奈当时搜索了一下如何制作无果后，比较愚蠢的决定用我略显笨拙的 PS 技术 P 一张图，怎么 P 呢，把他们的图下载下来，然后把代码去掉，换成自己的。。。

作为一款优秀的代码图片生成器，它几乎可以操作这个图片上的所有元素！

大部分敲代码的人总归是有一款喜爱的编辑器主题的，比如我自己就特别喜欢 Material Theme 这款主题。我们可以在设置栏中的主题选择框中选择自己喜欢的一款主题。

e8a4dc4fed822b2496ff501c24bfd52b.png

接下来是不同编程语言的语法高亮，支持所有主流的编程语言，也包含一些不同格式的配置文件类型，如 YAML，TOML 等。

846ce0e3a22d9a6386254c27598e67e9.png

你可以选择纯背景色，也可以选择其他图片作为背景。不过个人觉得纯背景色应该已经能满足需求啦～

b4ff101aa66d45e811f5a5507f44032d.png

贴心的 Carbon 提供了更为细微处的控制能力，比如编辑器是圆角还是方角。详细的配置如下：

对于窗口编辑器的基本样式：

  • 编辑器的边框样式
  • 编辑器距离四周的边距
  • 是否显示编辑器的阴影
  • 编辑器阴影的范围
  • 是否显示编辑器的控制按钮（最大化/最小化/关闭）
  • 是否显示代码行号
  • 是否自动调整编辑器宽度
  • 是否添加水印

对于代码部分的基本样式：

  • 字体类型
  • 字体大小
  • 行间距

因为我最喜欢的字体（Operator Mono）不在里面，所以我倒是希望能加个自定义字体的功能～

7c201d20a714550bef589ae4ab7b6854.png

Carbon 不仅可以导出 PNG 格式的图片，也可以导出 SVG，或者在线打开。

77d26f657baf7af8440a253bdecbe962.png

作为一款小工具，其功能无疑已经足够强大了，操作上没有什么难度，都是很简洁明了的功能。如果硬要说一些自己觉得不尽如人意的地方，就是每次打开生成图片，配置项都需要重新设置一遍（解决办法也有，就是保存这一次编辑好后的网站链接，下次再打开配置就是一样的了）。

如果你也苦于自己文章中的代码排版不满意却又没办法改变，不妨试试这款小工具，对我们技术自媒体人来说挺实用的。

欢迎关注微信公众号「X1nFly」~

093f95116981317a35a27d1220e99025.png

您可能感兴趣的与本文相关的镜像

Qwen-Image

Qwen-Image

图片生成
Qwen

Qwen-Image是阿里云通义千问团队于2025年8月发布的亿参数图像生成基础模型，其最大亮点是强大的复杂文本渲染和精确图像编辑能力，能够生成包含多行、段落级中英文文本的高保真图像

让看代码成为一种享受！使用Carbon生成漂亮的代码图片
05-26 1071
一、前言 当我们想展示自己写的代码给别人看的时候，希望让代码保持原有的格式，并且要美观一点，有一个强大的在线工具可以帮助我们实现。 这款工具就是：carbon，碳棒的意思。功能非常强大，操作很方便，简直好用极了！有哪些功能特点呢？如下所示： 支持市面上几乎所有的语言代码风格特点 可以自己选择生成图片的代码风格 可以选择生成什么语言的代码风格 还支持自定义选择背景色 二、实践 Carbon是开源免费的代码图片生成器，可以帮助我们很方便地生成漂亮的代码图片。操.
图片代码生成器
04-19
图片代码生成器，对刚入门HTML的初学者很有帮助
Pip:一个小工具，可以在正在运行的程序中以图片形式显示图片
05-26
点子 这只是此推文中要求的示例： : s 画中画工具，目前仅是一个原型。 它是如何工作的？ 当要镜像为小缩略图的窗口处于活动状态时，请按热键（Ctrl + Shift + P是默认设置） 出现一个小的缩略图窗口，您可以用鼠标左键移动它或通过右键单击将其关闭 您也可以在以前使用过的同一窗口上再次使用热键将其关闭 有一个配置，可以更改不透明度，也可以仅对窗口的客户区域进行镜像。 当前它同时支持5个缩略图，默认位置固定在屏幕右侧。 当未运行时，您可以在位于％APPDATA％\ Pip的pip.ini中更改热键。
在线代码图片生成工具carbon
weixin_43885532的博客
12-15 8546
在日常工作中时常需要和同事间进行代码的沟通和交流，有时只是需要讨论某一段代码的内容，因此不必将整个文件发给同事。通常可以将部分代码进行截图，或者直接将部分代码复制粘贴发送给同事。但以上方法或因为代码太长需要多次截图，破坏阅读代码时的连贯性（也可以生成长截图，但麻烦一点，也不美观），或因为复制粘贴发送过去后会显示格式混乱导致代码阅读性差。因此，为了代码分享时的格式保持不变、美观，可以使用如下代码图片生成工具carbon进行图片的生成。 carbon主要特如下： 几乎支持所有的语言代码风格 可选择主题、背
工具推荐：Carbon 代码图片生成器
weixin_34148340的博客
01-13 683
Carbon 是一款生成精美代码片段图片的工具。有多精美？美轮美奂。 之前在逛推特的时候，偶尔能看到一些类似下面这样的代码片段图片，每次看到都觉得特别精致。 有一次也特别想要这样的图片，无奈当时搜索了一下如何制作无果后，比较愚蠢的决定用我略显笨拙的 PS 技术 P 一张图，怎么 P 呢，把他们的图下载下来，然后把代码去掉，换成自己的。。。 作为一款优秀的代码图片生成器，它几乎可以操作这个图片上的...
请检查以下代码：%% 基于电-氢混合储能的电力系统优化调度模型（日内滚动调度） clc clear %% ==================== 参数设置 ==================== T = 96; % 总时段数 (15分钟间隔) W_size = 16; % 滚动窗口大小 (4小时) R_step = 4; % 滚动步长 (1小时) N_windows = floor((T - W_size)/R_step) + 1; % 滚动次数 N_Gens = 6; % 机组数量 dt = 0.25; % 时间步长 (小时) %% ==================== 数据加载 ==================== % 电负荷数据 (小时级) P_load_Hourly = [297 258 220 208 226 258 280 330 372 364 400 468 636 600 506 554 392 392 412 403 426 390 322 268]; % 风电数据 (小时级) P_wt_Hourly = [50 55 54 58 60 100 102 125 130 152 153 140 157 130 127 300 275 260 210 200 169 150 130 100]; % 转换为15分钟数据 (96时段) P_load_Min_96 = repelem(P_load_Hourly, 1, 4); % 系统总负荷 P_wt_Min_96 = repelem(P_wt_Hourly, 1, 4); % 风电出力 mpc = loadcase('case30'); 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) Carbon_Factor = Gen_Params(:,8); % 碳排放系数 (tCO2/MWh) P_GT_Min = Gen_Params(:,3); % 最小出力 P_GT_Max = Gen_Params(:,2); % 最大出力 % 启动/停机成本参数 (单位: 元) C_Start_qt = [2768; 2874; 2136; 2677; 2262; 2869]; % 启动成本 C_Stop_qt = [1500; 1600; 1200; 1400; 1300; 1500]; % 停机成本 % 电储能参数 P_Cha_Max = 45; % 最大充电功率 (MW) P_Dis_Max = 90; % 最大放电功率 (MW) K_ess_Eff = 0.92; % 充放电效率 C_ess_Max = 300; % 最大容量 (MWh) SOC_ESS_Init = 10; % 初始SOC (MWh) % 氢系统参数 P_FC_Max = 400; % 燃料电池最大功率 (MW) K_H2P_xs = 9.1; % 氢转电系数 (kg/MWh) P_EL_Max = 400; % 电解槽最大功率 (MW) K_P2H_xs = 0.22; % 电转氢系数 (kg/MWh) H_cha_Max = 100; % 最大充氢速率 (kg/h) H_dis_Max = 100; % 最大放氢速率 (kg/h) C_HSS_Max = 400; % 最大储氢容量 (kg) C_HSS_Min = 10; % 最小储氢容量 (kg) SOC_HSS_Init = 20; % 初始氢储量 (kg) % 环境成本系数 C_carbon = 150; % 碳排放成本 (元/吨CO2) % 节点位置 Bus_GEN_e = [1 2 5 8 11 13]; % 火电机组节点 Bus_WT_e = 2; % 风电节点 Bus_ESS_e = 11; % 电储能节点 Bus_HSS_h = 8; % 氢储能节点 %% ==================== 节点电负荷矩阵生成 ==================== % 电力系统基准负荷 (MW) P_Load_Base = sum(mpc.bus(:,3)); % 创建节点负荷矩阵 (30×96) P_load_96 = zeros(size(mpc.bus, 1), T); for i = 1:size(mpc.bus, 1) % 计算节点负荷比例 node_ratio = mpc.bus(i,3) / P_Load_Base; % 生成节点负荷曲线 P_load_96(i,:) = node_ratio * P_load_Min_96; end %% ==================== 电网络拓扑初始化 ==================== [Ybus, ~, ~] = makeYbus(mpc); B = imag(Ybus); % 电纳矩阵 ref_node = 1; % 参考节点 non_ref = 2:size(mpc.bus,1); % 非参考节点 Bred = B(non_ref, non_ref); % 降阶B矩阵 % 创建设备-节点映射矩阵 N = size(mpc.bus, 1); % 节点数量 G_map = sparse(Bus_GEN_e, 1:N_Gens, 1, N, N_Gens); % 常规机组 W_map = sparse(Bus_WT_e, 1, 1, N, 1); % 风电机组 E_map = sparse(Bus_ESS_e, 1, 1, N, 1); % 电储能 H_map = sparse(Bus_HSS_h, 1, 1, N, 1); % 氢系统 %% ==================== 状态初始化 ==================== % 结果存储 Results = struct(... 'P_GT', zeros(N_Gens, T), ... % 机组出力 'P_wt', zeros(1, T), ... % 风电出力 'P_cha', zeros(1, T), ... % 充电功率 'P_dis', zeros(1, T), ... % 放电功率 'SOC_ESS', SOC_ESS_Init*ones(1, T+1), ... % 电储能状态 'P_FC', zeros(1, T), ... % 燃料电池出力 'P_EL', zeros(1, T), ... % 电解槽用电 'SOC_HSS', SOC_HSS_Init*ones(1, T+1), ... % 氢储能状态 'Obj', zeros(1, N_windows) ... % 各窗口目标值 ); % 机组状态跟踪 B_ug_96 = zeros(N_Gens, T); % 机组状态历史 B_ug_96(:,1) = ones(N_Gens,1); % 初始全运行 %% ==================== 主循环滚动优化 ==================== for idx = 1:N_windows % 窗口时间范围 k_start = (idx-1)*R_step + 1; k_end = min(k_start + W_size - 1, T); win_len = k_end - k_start + 1; fprintf('\n====== 滚动窗口 %d/%d [时段 %d-%d] ======\n', idx, N_windows, k_start, k_end); % 获取初始状态 if k_start == 1 SOC_ESS_init = SOC_ESS_Init; SOC_HSS_init = SOC_HSS_Init; P_GT_prev = P_GT_Min; % 机组上一时段出力 prev_status = ones(N_Gens,1); % 机组上一时段状态 else SOC_ESS_init = Results.SOC_ESS(k_start); SOC_HSS_init = Results.SOC_HSS(k_start); P_GT_prev = Results.P_GT(:, k_start-1); prev_status = B_ug_96(:, k_start-1); end % 定义决策变量 (当前窗口) P_GT_e = sdpvar(N_Gens, win_len, 'full'); % 机组出力 B_U_G = binvar(N_Gens, win_len, 'full'); % 机组运行状态 B_U_Start = binvar(N_Gens, win_len, 'full');% 启动状态 B_U_Stop = binvar(N_Gens, win_len, 'full'); % 停机状态 P_wt_e = sdpvar(1, win_len); % 风电出力 P_cha_e = sdpvar(1, win_len); % 充电功率 P_dis_e = sdpvar(1, win_len); % 放电功率 SOC_ESS = sdpvar(1, win_len+1); % 电储能状态 P_FC_e = sdpvar(1, win_len); % 燃料电池出力 P_e_EL = sdpvar(1, win_len); % 电解槽用电 SOC_HSS = sdpvar(1, win_len+1); % 氢储能状态 Bus_Theta = sdpvar(N, win_len, 'full'); % 节点相角 % 约束条件 Cons = []; % 初始状态约束 Cons = [Cons, SOC_ESS(1) == SOC_ESS_init]; Cons = [Cons, SOC_HSS(1) == SOC_HSS_init]; % 机组约束 for t = 1:win_len % 出力上下限 Cons = [Cons, B_U_G(:,t).*P_GT_Min <= P_GT_e(:,t) <= B_U_G(:,t).*P_GT_Max]; % 爬坡约束 if t == 1 ramp_up = P_GT_e(:,1) - P_GT_prev; ramp_down = P_GT_prev - P_GT_e(:,1); else ramp_up = P_GT_e(:,t) - P_GT_e(:,t-1); ramp_down = P_GT_e(:,t-1) - P_GT_e(:,t); end Cons = [Cons, -Ramp_e <= ramp_up <= Ramp_e]; Cons = [Cons, -Ramp_e <= ramp_down <= Ramp_e]; % 启停逻辑 if t == 1 Cons = [Cons, B_U_Start(:,1) >= B_U_G(:,1) - prev_status]; Cons = [Cons, B_U_Start(:,1) <= B_U_G(:,1)]; Cons = [Cons, B_U_Start(:,1) <= 1 - prev_status]; Cons = [Cons, B_U_Stop(:,1) >= prev_status - B_U_G(:,1)]; Cons = [Cons, B_U_Stop(:,1) <= prev_status]; Cons = [Cons, B_U_Stop(:,1) <= 1 - B_U_G(:,1)]; else Cons = [Cons, B_U_Start(:,t) >= B_U_G(:,t) - B_U_G(:,t-1)]; Cons = [Cons, B_U_Start(:,t) <= B_U_G(:,t)]; Cons = [Cons, B_U_Start(:,t) <= 1 - B_U_G(:,t-1)]; Cons = [Cons, B_U_Stop(:,t) >= B_U_G(:,t-1) - B_U_G(:,t)]; Cons = [Cons, B_U_Stop(:,t) <= B_U_G(:,t-1)]; Cons = [Cons, B_U_Stop(:,t) <= 1 - B_U_G(:,t)]; end end % 风电约束 Cons = [Cons, 0 <= P_wt_e <= P_wt_Min_96(k_start:k_end)]; % 电储能约束 b_in_cha = binvar(1, win_len); % 充电状态 b_in_dis = binvar(1, win_len); % 放电状态 Cons = [Cons, b_in_cha + b_in_dis <= 1]; % 充放电互斥 for t = 1:win_len % SOC更新 (加入时间步长0.25小时) Cons = [Cons, SOC_ESS(t+1) == SOC_ESS(t) + dt*(K_ess_Eff*P_cha_e(t) - P_dis_e(t)/K_ess_Eff)]; % 容量约束 Cons = [Cons, 0 <= SOC_ESS(t+1) <= C_ess_Max]; % 功率约束 Cons = [Cons, 0 <= P_cha_e(t) <= b_in_cha(t)*P_Cha_Max]; Cons = [Cons, 0 <= P_dis_e(t) <= b_in_dis(t)*P_Dis_Max]; end % 氢系统约束 for t = 1:win_len % 产氢量和耗氢量计算 H_prod = P_e_EL(t) * dt * K_P2H_xs; % 产氢量 (kg) H_cons = P_FC_e(t) * dt * K_H2P_xs; % 耗氢量 (kg) % SOC更新 Cons = [Cons, SOC_HSS(t+1) == SOC_HSS(t) + H_prod - H_cons]; % 容量约束 Cons = [Cons, C_HSS_Min <= SOC_HSS(t+1) <= C_HSS_Max]; % 功率约束 Cons = [Cons, 0 <= P_FC_e(t) <= P_FC_Max]; Cons = [Cons, 0 <= P_e_EL(t) <= P_EL_Max]; Cons = [Cons, 0 <= H_prod <= H_cha_Max*dt]; Cons = [Cons, 0 <= H_cons <= H_dis_Max*dt]; end % 网络约束 for t = 1:win_len t_global = k_start + t - 1; % 节点注入功率 P_inj = G_map * P_GT_e(:, t) + ... % 常规机组 W_map * P_wt_e(t) + ... % 风电 E_map * (P_dis_e(t) - P_cha_e(t)) + ... % 电储能 H_map * (P_FC_e(t) - P_e_EL(t));% 氢系统 % 节点净功率 P_net = P_inj - P_load_96(:, t_global); % 直流潮流约束 Cons = [Cons, Bred * Bus_Theta(non_ref, t) == P_net(non_ref)/baseMVA]; end Cons = [Cons, Bus_Theta(ref_node, :) == 0]; % 参考节点 % 目标函数 Obj = 0; for t = 1:win_len % 燃料成本 for g = 1:N_Gens fuel_cost = a2(g)*P_GT_e(g,t)^2 + a1(g)*P_GT_e(g,t) + a0(g); Obj = Obj + dt * fuel_cost; end % 启停成本 for g = 1:N_Gens Obj = Obj + C_Start_qt(g)*B_U_Start(g,t) + C_Stop_qt(g)*B_U_Stop(g,t); end % 碳排放成本 for g = 1:N_Gens carbon_cost = C_carbon * Carbon_Factor(g) * P_GT_e(g,t); Obj = Obj + dt * carbon_cost; end end % 求解优化问题 opts = sdpsettings('solver', 'cplex', 'verbose', 1, 'showprogress', 1); sol = optimize(Cons, Obj, opts); % 结果处理 if sol.problem == 0 fprintf('优化成功 | 窗口成本: %.2f 元\n', value(Obj)); Results.Obj(idx) = value(Obj); % 计算实际可保存长度 valid_len = min(R_step, T - k_start + 1); % 保存当前窗口结果 idx_range = k_start:k_start+valid_len-1; Results.P_GT(:, idx_range) = value(P_GT_e(:,1:valid_len)); B_ug_96(:, idx_range) = value(B_U_G(:,1:valid_len)); Results.P_wt(idx_range) = value(P_wt_e(1:valid_len)); Results.P_cha(idx_range) = value(P_cha_e(1:valid_len)); Results.P_dis(idx_range) = value(P_dis_e(1:valid_len)); Results.P_FC(idx_range) = value(P_FC_e(1:valid_len)); Results.P_EL(idx_range) = value(P_e_EL(1:valid_len)); % 更新状态量 Results.SOC_ESS(k_start+1:k_end+1) = value(SOC_ESS(2:end)); Results.SOC_HSS(k_start+1:k_end+1) = value(SOC_HSS(2:end)); else fprintf('优化失败: %s\n', sol.info); % 使用前一时段值填充当前窗口 Results.P_GT(:, k_start:k_end) = repmat(Results.P_GT(:, max(1,k_start-1)), 1, win_len); Results.P_wt(k_start:k_end) = P_wt_Min_96(k_start:k_end); Results.Obj(idx) = Inf; end end
11-28
对于机组状态，代码使用B_ug_96(:,1)初始化为1，但这里注意：在第一个窗口，上一时段（t=0）的状态并没有定义，代码使用P_GT_prev = P_GT_Min（一个6x1向量）和prev_status=ones(N_Gens,1)（全为1）是合理的。...
请修改以下代码：%% 基于电-氢混合储能的电力系统优化调度模型（日内滚动优化） 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
### 修改后的完整代码 ```matlab %% 基于电-氢混合储能的电力系统优化调度模型（增强版） clc clear %% ==================== 参数设置 ==================== T = 96; % 时间周期数 (15分钟间隔) N_Gens = 6; % ...
请检查以下代码，并整理出一个修改后的完整代码：%% 基于电-氢混合储能的电力系统优化调度模型（最终修正版） 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
代码修复了数据加载、状态传递、约束条件和可视化等多个问题，并添加了详细注释： ```matlab %% 基于电-氢混合储能的电力系统优化调度模型（完整修正版） clc clear close all %% ==================== 参数设置 =...
我还希望继续优化，形成日内滚动优化调度，可以参考以下代码：%% 基于电-氢混合储能的电力系统优化调度模型（日内经济调度） clc clear %% 加载电力系统数据 mpc = loadcase('case30'); baseMVA = mpc.baseMVA; % 系统基准功率 % 原始24点数据 P_load_Hourly_24 = [297 258 220 208 226 258 280 330 372 364 400 468 636 600 506 554 392 392 412 403 426 390 322 268]; P_wt_Hourly_24 = [50 55 54 58 60 100 102 125 130 152 153 140 157 130 127 300 275 260 210 200 169 150 130 100]; %% 数据插值（24点→96点） x_24 = 1:24; x_96 = linspace(1, 24, 96); P_load_Hourly = interp1(x_24, P_load_Hourly_24, x_96, 'linear'); P_wt_Hourly = interp1(x_24, P_wt_Hourly_24, x_96, 'linear'); %% 参数设置 T = 96; % 时间周期数（15分钟间隔） dt = 0.25; % 时间步长（小时） N_Gens = 6; % 常规机组数量 % 常规机组参数 [机组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); % 爬坡率(每小时) Ramp_15min = Ramp_e * dt; % 15分钟爬坡率 Carbon_Factor = Gen_Params(:,8); % 碳排放系数 (tCO2/MWh) P_GT_Min = Gen_Params(:,3); % 最小出力 P_GT_Max = Gen_Params(:,2); % 最大出力 % 启动/停机成本参数 (单位: 元) C_Start_qt = [2768; 2874; 2136; 2677; 2262; 2869]; % 启动成本 C_Stop_qt = [1500; 1600; 1200; 1400; 1300; 1500]; % 停机成本 % 电储能参数 P_Cha_Max = 45; % 电储能最大充电功率 P_Dis_Max = 90; % 电储能最大放电功率 K_ess_Eff = 0.92; % 电储能转换效率 K_Sel_Dis = 0.98; % 自放电率(99.8%/小时) K_Sel_Dis_15min = K_Sel_Dis^dt; % 15分钟自放电率 C_ess_Max = 300; % 储能容量限制 % 氢系统参数 P_FC_Max = 400; % 燃料电池最大功率 K_H2P_xs = 9.1; % 氢转电系数 (kg/MWh) P_EL_Max = 400; % 电解槽最大功率 K_P2H_xs = 0.22; % 电转氢系数 (kg/MWh) H_cha_Max = 120; % 最大充氢速率 (kg/h) H_dis_Max = 100; % 最大放氢速率 (kg/h) C_HSS_Max = 400; % 最大储氢容量 (kg) C_HSS_Min = 10; % 最小储氢容量 (kg) % 节点位置 Bus_GEN_e = [1 2 5 8 11 13]; % 火电机组节点 Bus_WT_e = 2; % 风电节点 Bus_ESS_e = 11; % 电储能节点 Bus_HSS_h = 8; % 氢储能节点 % 环境成本系数 C_carbon = 150; % 碳排放成本 (元/吨CO2) C_FC_env = -500; % 燃料电池环境补贴 (负成本) C_EL_env = -300; % 电解槽环境补贴 (负成本) % 峰谷时段定义（15分钟点） t_peak = [12*4+1:14*4, 19*4+1:21*4]; % 高峰时段 [49:56, 77:84] t_valley = [1:20, 96]; % 低谷时段 [1:20, 96] %% 决策变量定义 Bus_Theta = sdpvar(30, T, 'full'); % 节点相角 P_GT_e = sdpvar(N_Gens, T, 'full'); % 火电机组出力 B_U_GT = binvar(N_Gens, T, 'full'); % 火电运行状态 B_Start_GT = binvar(N_Gens, T, 'full'); % 火电启动 B_Stop_GT = binvar(N_Gens, T, 'full'); % 火电停机 P_wt_e = sdpvar(1, T, 'full'); % 风电出力 % 电储能变量 P_cha_e = sdpvar(1, T, 'full'); % 充电功率 P_dis_e = sdpvar(1, T, 'full'); % 放电功率 B_U_cha = binvar(1, T, 'full'); % 充电状态 B_U_dis = binvar(1, T, 'full'); % 放电状态 SOC_ESS = sdpvar(1, T+1, 'full'); % 电储能SOC % 氢系统变量 P_FC_e = sdpvar(1, T, 'full'); % 燃料电池出力 P_e_EL = sdpvar(1, T, 'full'); % 电解槽用电 B_U_EL = binvar(1, T, 'full'); % 电解槽状态 B_U_FC = binvar(1, T, 'full'); % 燃料电池状态 P_H2_FC = sdpvar(1, T, 'full'); % 燃料电池耗氢 P_EL_H2 = sdpvar(1, T, 'full'); % 电解槽产氢 P_cha_H2 = sdpvar(1, T, 'full'); % 储氢罐充氢 P_dis_H2 = sdpvar(1, T, 'full'); % 储氢罐放氢 SOC_HSS = sdpvar(1, T+1, 'full'); % 储氢罐SOC % 节点注入功率矩阵 P_net_e = sdpvar(30, T, 'full'); %% 生成节点电负荷矩阵 P_Load_Base = sum(mpc.bus(:,3)); % 系统基准负荷 Ratio_Load = P_load_Hourly / P_Load_Base; P_load_e = mpc.bus(:,3) * Ratio_Load; % 节点负荷矩阵 (30×T, MW) %% 约束条件 Cons = []; % ========== 火电机组约束 ========== for t = 1:T % 机组出力限制 Cons = [Cons, P_GT_Min .* B_U_GT(:,t) <= P_GT_e(:,t) <= P_GT_Max .* B_U_GT(:,t)]; % 机组爬坡约束（15分钟爬坡） if t > 1 Cons = [Cons, -Ramp_15min <= P_GT_e(:,t) - P_GT_e(:,t-1) <= Ramp_15min]; end % 启停状态逻辑 if t == 1 Cons = [Cons, B_Start_GT(:,t) == B_U_GT(:,t)]; Cons = [Cons, B_Stop_GT(:,t) == 0]; else Cons = [Cons, B_Start_GT(:,t) <= 1 - B_U_GT(:,t-1)]; Cons = [Cons, B_Start_GT(:,t) <= B_U_GT(:,t)]; Cons = [Cons, B_Start_GT(:,t) >= B_U_GT(:,t) - B_U_GT(:,t-1)]; Cons = [Cons, B_Stop_GT(:,t) <= B_U_GT(:,t-1)]; Cons = [Cons, B_Stop_GT(:,t) <= 1 - B_U_GT(:,t)]; Cons = [Cons, B_Stop_GT(:,t) >= B_U_GT(:,t-1) - B_U_GT(:,t)]; end end % ========== 风电机组约束 ========== Cons = [Cons, 0 <= P_wt_e <= P_wt_Hourly]; % ========== 电储能约束 ========== Cons = [Cons, B_U_cha + B_U_dis <= 1]; % 充放电互斥 Cons = [Cons, 0 <= P_cha_e <= P_Cha_Max .* B_U_cha]; Cons = [Cons, 0 <= P_dis_e <= P_Dis_Max .* B_U_dis]; Cons = [Cons, SOC_ESS(1) == 10]; % 初始SOC Cons = [Cons, SOC_ESS(T+1) == 10]; % 终止SOC for t = 1:T Cons = [Cons, SOC_ESS(t+1) == K_Sel_Dis_15min * SOC_ESS(t) + ... K_ess_Eff * P_cha_e(t) * dt - ... % 充电能量 P_dis_e(t)/K_ess_Eff * dt]; % 放电能量 Cons = [Cons, 0 <= SOC_ESS(t+1) <= C_ess_Max]; end % ========== 氢系统约束 ========== Cons = [Cons, 0 <= P_e_EL <= P_EL_Max .* B_U_EL]; Cons = [Cons, P_EL_H2 == P_e_EL * K_P2H_xs * dt]; % 产氢量 (kg) Cons = [Cons, 0 <= P_FC_e <= P_FC_Max .* B_U_FC]; Cons = [Cons, P_FC_e == P_H2_FC / K_H2P_xs / dt]; % 耗氢功率 (kg/h) % 氢能平衡 Cons = [Cons, P_H2_FC * dt + P_cha_H2 == P_EL_H2 + P_dis_H2]; % 互斥运行 Cons = [Cons, B_U_EL + B_U_FC <= 1]; % 充放氢速率限制 (kg/h) Cons = [Cons, 0 <= P_cha_H2 <= H_cha_Max .* B_U_EL]; Cons = [Cons, 0 <= P_dis_H2 <= H_dis_Max .* B_U_FC]; % 储氢容量约束 Cons = [Cons, SOC_HSS(1) == 20]; % 初始氢储量 Cons = [Cons, SOC_HSS(T+1) == 20]; % 终止氢储量 for t = 1:T Cons = [Cons, SOC_HSS(t+1) == SOC_HSS(t) + ... (P_cha_H2(t) - P_dis_H2(t)) * dt]; Cons = [Cons, C_HSS_Min <= SOC_HSS(t+1) <= C_HSS_Max]; end % FC/EL爬坡约束 (15分钟) for t = 2:T Cons = [Cons, -7.5 <= P_FC_e(t)-P_FC_e(t-1) <= 7.5]; % 30 MW/h ÷ 4 Cons = [Cons, -5 <= P_e_EL(t)-P_e_EL(t-1) <= 5]; % 20 MW/h ÷ 4 end % ========== 节点功率平衡 ========== % 创建设备-节点映射矩阵 N = size(P_load_e, 1); G_map = sparse(Bus_GEN_e, 1:N_Gens, 1, N, N_Gens); % 常规机组 W_map = sparse(Bus_WT_e, 1, 1, N, 1); % 风电机组 E_map = sparse(Bus_ESS_e, 1, 1, N, 1); % 电储能 H_map = sparse(Bus_HSS_h, 1, 1, N, 1); % 氢系统 for t = 1:T total_injection = G_map * P_GT_e(:,t) + ... W_map * P_wt_e(t) + ... E_map * (P_dis_e(t) - P_cha_e(t)) + ... H_map * (P_FC_e(t) - P_e_EL(t)); Cons = [Cons, P_net_e(:,t) == total_injection - P_load_e(:,t)]; end % ========== 直流潮流约束 ========== [Ybus, ~, ~] = makeYbus(mpc); B = imag(Ybus); % 电纳矩阵虚部 ref_node = 1; % 参考节点 non_ref = 2:30; % 非参考节点 Bred = B(non_ref, non_ref); % 降阶B矩阵 Cons = [Cons, Bus_Theta(ref_node, :) == 0]; % 参考节点相角 for t = 1:T Cons = [Cons, Bred * Bus_Theta(non_ref, t) == P_net_e(non_ref, t) / baseMVA]; end Cons = [Cons, -pi/2 <= Bus_Theta <= pi/2]; % 相角范围 % ========== 全局功率平衡 ========== Total_Gen_e = sum(P_GT_e, 1) + P_wt_e + P_dis_e - P_cha_e + P_FC_e - P_e_EL; Total_Load_e = sum(P_load_e, 1); Cons = [Cons, Total_Gen_e == Total_Load_e]; % 严格功率平衡 %% 目标函数 (使用循环计算) % 1. 燃料成本 (乘以时间步长) F1_cost = 0; for i = 1:N_Gens for t = 1:T F1_cost = F1_cost + (a2(i) * P_GT_e(i,t)^2 + a1(i) * P_GT_e(i,t) + a0(i) * B_U_GT(i,t)) * dt; end end % 2. 启停成本 (不变) F2_cost = 0; for i = 1:N_Gens for t = 1:T F2_cost = F2_cost + C_Start_qt(i) * B_Start_GT(i,t) + C_Stop_qt(i) * B_Stop_GT(i,t); end end % 3. 环境成本 Carbon_Emission = 0; for i = 1:N_Gens for t = 1:T Carbon_Emission = Carbon_Emission + Carbon_Factor(i) * P_GT_e(i,t) * dt; end end F_carbon = C_carbon * Carbon_Emission; % 氢储能环境补贴 (乘以时间步长) F_env = 0; for t = 1:T % 燃料电池补贴 if ismember(t, t_peak) F_env = F_env + 1.3 * C_FC_env * P_FC_e(t) * dt; else F_env = F_env + C_FC_env * P_FC_e(t) * dt; end % 电解槽补贴 if ismember(t, t_valley) F_env = F_env + 1.5 * C_EL_env * P_e_EL(t) * dt; else F_env = F_env + C_EL_env * P_e_EL(t) * dt; end end F3_cost = F_carbon + F_env; F_Obj_Total = F1_cost + F2_cost + F3_cost; %% 模型优化求解 Opt = sdpsettings('solver', 'cplex', 'verbose', 1, 'showprogress', 1); Sol = optimize(Cons, F_Obj_Total, Opt); % 结果处理与分析 if Sol.problem == 0 disp('======= 优化成功 ======='); fprintf('>>>总成本: %.2f 元\n', value(F_Obj_Total)); fprintf('-燃料成本: %.2f 元\n', value(F1_cost)); fprintf('-启停成本: %.2f 元\n', value(F2_cost)); fprintf('-环境成本: %.2f 元\n', value(F3_cost)); % 提取结果 P_GT_e = value(P_GT_e); P_wt_e = value(P_wt_e); P_cha_e = value(P_cha_e); P_dis_e = value(P_dis_e); P_FC_e = value(P_FC_e); P_e_EL = value(P_e_EL); SOC_ESS = value(SOC_ESS); SOC_HSS = value(SOC_HSS); % 输出氢储能使用情况 disp('==== 氢储能使用情况 ===='); fprintf('电解槽总用电量: %.2f MWh\n', sum(P_e_EL)*dt); fprintf('燃料电池总发电量: %.2f MWh\n', sum(P_FC_e)*dt); else error('>>>求解失败: %s', Sol.info); end %% 结果可视化 (96点) figure(1) components = [sum(P_GT_e, 1)', ... % 常规机组 P_wt_e', ... % 风电机组 P_dis_e', ... % 储能放电 P_FC_e', ... % 燃料电池 -P_cha_e', ... % 储能充电（负值） - P_e_EL']; % 电解槽机（负值） he1 = bar(components, 'stacked'); hold on; total_load = sum(P_load_e, 1); plot(1:T, total_load, 'g-*', 'LineWidth', 1.5); % 设置颜色和透明度 colors = lines(6); for i = 1:6 he1(i).FaceColor = colors(i,:); he1(i).FaceAlpha = 0.9; end legend('常规机组', '风电机组', '储能放电', '燃料电池', '储能充电', '电解槽机', '总电负荷'); grid on; xlabel('时间 (15分钟间隔)'); ylabel('功率 (MW)'); title('电力系统功率平衡 (日内调度)'); xlim([0 T+1]); figure(2) Plot_Ele_Net = zeros(T, N_Gens); for t = 1:T for i = 1:N_Gens Plot_Ele_Net(t, i) = P_GT_e(i, t); end end % 绘制堆叠条形图 bar(Plot_Ele_Net, 'stacked'); hold on % 添加负荷曲线 plot(1:T, sum(P_GT_e,1), 'g-*', 'LineWidth', 1.5); grid on legend('机组1','机组2','机组3','机组4','机组5','机组6','常规机组'); xlabel('时间 (15分钟间隔)'); ylabel('功率/MW'); title('常规机组出力分布 (日内调度)'); xlim([0 T+1]); figure(3) commitment = double(B_U_GT); % 逻辑矩阵转数值 [num_units, T] = size(commitment); % 创建增强型颜色映射 cmap = [0.95 0.95 0.95; % 停机 - 浅灰色 0.25 0.65 0.35]; % 运行 - 绿色 % 创建热图 imagesc(commitment); colormap(cmap); % 设置坐标轴 set(gca, 'YTick', 1:num_units, ... 'YTickLabel', arrayfun(@(x) sprintf('机组 %d', x), 1:num_units, 'UniformOutput', false), ... 'TickDir', 'out', ... 'FontSize', 10); % 仅在可读性好的情况下添加1/0 for i = 1:num_units for t = 1:T if num_units <= 30 && T <= 96 % 仅在矩阵尺寸可读时添加文本 if commitment(i,t) > 0.5 text(t, i, '1', 'HorizontalAlignment', 'center', ... 'Color', [0 0.3 0], 'FontWeight', 'bold', 'FontSize', 9); else text(t, i, '0', 'HorizontalAlignment', 'center', ... 'Color', [0.4 0.4 0.4], 'FontSize', 9); end end end end grid on; % 添加颜色条和图例 c = colorbar('Ticks', [0.25, 0.75], 'TickLabels', {'停机', '运行'}); c.Label.String = '状态'; xlabel('时间 (15分钟间隔)'); ylabel('机组编号'); title('机组启停状态分析 (日内调度)'); xlim([0 T+1]); figure(4); hourly_carbon = value(Carbon_Factor' * P_GT_e * dt); bar(hourly_carbon', 'FaceColor', [0.25 0.65 0.35], 'EdgeColor', 'none'); hold on; plot(1:T, value(P_FC_e), 'b-*', 'LineWidth', 1.5); plot(1:T, value(P_e_EL), 'r-*', 'LineWidth', 1.5); xlabel('时间 (15分钟间隔)'); ylabel('碳排放 (吨CO₂) / 功率 (MW)'); title('碳排放量与氢能系统运行 (日内调度)'); legend('机组碳排放量', '燃料电池出力', '电解槽器用电', 'Location', 'best'); grid on; xlim([0 T+1]); ylim([-50 500]); figure(5) yyaxis left; bar(P_cha_e', 'FaceColor', [1 0 0]); hold on bar(-P_dis_e', 'FaceColor', [0 0 1]); ylabel('充放电功率 (MW)'); ylim([-100 100]); yyaxis right; plot(0:T, SOC_ESS, 'g-*', 'LineWidth', 1.5); ylabel('电储能容量 (MWh)'); title('电储能运行状态 (日内调度)'); legend( '充电功率', '放电功率','电储能SOC', 'Location', 'best'); grid on; xlim([0 T+1]); figure(6) yyaxis left; bar(P_cha_H2', 'FaceColor', [1 0 0]); hold on bar(-P_dis_H2', 'FaceColor', [0 0 1]); ylabel('氢流量 (kg/h)'); ylim([-150 150]); yyaxis right; plot(0:T, SOC_HSS, 'g-*', 'LineWidth', 1.5); ylabel('储氢量 (kg)'); title('氢储能状态 (日内调度)'); legend( '储能产氢', '储能耗氢','储氢量SOC', 'Location', 'best'); grid on; xlim([0 T+1]); figure(7) Plot_H2_Net = zeros(T, 4); for t = 1:T Plot_H2_Net(t, 1) = value(P_EL_H2(1, t)); Plot_H2_Net(t, 2) = -value(P_H2_FC(1, t)); Plot_H2_Net(t, 3) = -value(P_cha_H2(1, t)); Plot_H2_Net(t, 4) = value(P_dis_H2(1, t)); end bar(Plot_H2_Net, 'stacked'); hold on grid on legend('电解槽器产氢', '燃料电池耗氢', '储氢罐充氢', '储氢罐放氢'); xlabel('时间 (15分钟间隔)'); ylabel('储氢量 (kg)'); title('供需氢功率平衡状态 (日内调度)'); xlim([0 T+1]); ylim([-150 200]); %% 基于电-氢混合储能的电力系统优化调度模型（日内滚动调度） clc clear close all %% ==================== 参数设置 ==================== T = 96; % 总时段数 (15分钟间隔) W_size = 16; % 滚动窗口大小 (4小时) R_step = 4; % 滚动步长 (1小时) N_windows = floor((T - W_size)/R_step) + 1; % 滚动次数 N_Gens = 6; % 机组数量 dt = 0.25; % 时间步长 (小时) %% ==================== 数据加载与预处理 ==================== load pload_24.mat p_load_1 load pwt_24.mat p_wt load ug_24.mat ugone % 扩展24小时数据到96小时 P_wt_96 = repelem(p_wt, 1, 4); % 风电预测 P_load_96 = zeros(30, T); % 节点负荷 B_ug_96 = zeros(6, T); % 机组状态 for i = 1:30 P_load_96(i, :) = repelem(p_load_1(i, :), 1, 4); end for i = 1:6 B_ug_96(i, :) = repelem(ugone(i, :), 1, 4); end %% ==================== 系统参数设置 ==================== mpc = loadcase('case30'); baseMVA = mpc.baseMVA; % 系统基准功率 % 常规机组参数 Gen_Params = [1 250 50 0.0375 20 372.5 72; 2 80 20 0.175 17.5 352.3 48; 3 50 15 0.625 10 316.5 30; 4 35 10 0.0834 32.5 329.2 21; 5 30 10 0.25 30 276.4 18; 6 40 12 0.25 30 232.2 24]; % 提取发电机参数 a2 = Gen_Params(:,4); % 二次项系数 a1 = Gen_Params(:,5); % 一次项系数 a0 = Gen_Params(:,6); % 常数项 Ramp = Gen_Params(:,7); % 爬坡率 (MW/15min) Pg_min = Gen_Params(:,3); % 最小出力 Pg_max = Gen_Params(:,2); % 最大出力 % 启动成本参数 StartupCost = [2768; 2874; 2136; 2677; 2262; 2869]; % 环境成本系数 (新增) C_carbon = 150; % 碳排放成本 (元/吨CO2) C_FC_env = -500; % 燃料电池环境补贴 (负成本) 元/MW C_EL_env = -300; % 电解槽环境补贴 (负成本) 元/MW Carbon_Factor = [0.85; 0.92; 1.05; 0.78; 0.95; 0.88]; % 碳排放系数 (吨CO2/MWh) % 峰谷时段定义 (15分钟分辨率) (新增) t_peak_hour = [12:14, 19:21]; % 高峰时段（小时） t_valley_hour = [1:5, 24]; % 低谷时段（小时） t_peak = []; for h = t_peak_hour t_peak = [t_peak, (h-1)*4+1 : h*4]; end t_valley = []; for h = t_valley_hour t_valley = [t_valley, (h-1)*4+1 : h*4]; end % 电储能参数 P_Cha_Max = 45; % 最大充电功率 (MW) P_Dis_Max = 90; % 最大放电功率 (MW) K_ess_Eff = 0.92; % 充放电效率 C_ess_Max = 300; % 最大容量 (MWh) SOC_ESS_Init = 10; % 初始SOC (MWh) % 氢系统参数 (更新单位) P_FC_Max = 400; % 燃料电池最大功率 (MW) K_H2P_xs = 9.1; % 氢转电系数 (kg/MWh) (原名K_H2P) P_EL_Max = 400; % 电解槽最大功率 (MW) K_P2H_xs = 0.22; % 电转氢系数 (kg/MWh) (原名K_P2H) H_cha_Max = 100; % 最大充氢速率 (kg/h) H_dis_Max = 100; % 最大放氢速率 (kg/h) C_HSS_Max = 400; % 最大储氢容量 (kg) C_HSS_Min = 0; % 最小储氢容量 (kg) SOC_HSS_Init = 20; % 初始氢储量 (kg) % 成本系数 kc = 0.5; % 弃风惩罚因子 O_M_ESS = [0.56, 0.76]; % 电储能运维成本 [充电, 放电] O_M_H2 = [0.65, 0.83]; % 氢系统运维成本 [电解槽, 燃料电池] % 节点位置 Bus_GEN = [1 2 5 8 11 13]; % 火电机组节点 Bus_WT = 2; % 风电节点 Bus_ESS = 11; % 电储能节点 Bus_HSS = 8; % 氢储能节点 %% ==================== 网络拓扑初始化 ==================== [Ybus, ~, ~] = makeYbus(mpc); Gbus = real(Ybus); % 电导矩阵 Bbus = imag(Ybus); % 电纳矩阵 % 创建节点-设备映射 N = size(mpc.bus, 1); % 节点数量 Mbg = sparse(Bus_GEN, 1:N_Gens, 1, N, N_Gens); % 常规机组 Mbw = sparse(Bus_WT, 1, 1, N, 1); % 风电 Mbe = sparse(Bus_ESS, 1, 1, N, 1); % 电储能 Mbs = sparse(Bus_HSS, 1, 1, N, 1); % 氢系统 %% ==================== 状态初始化 ==================== % 结果存储 Results = struct(... 'P_GT', zeros(N_Gens, T), ... % 机组出力 'P_wt', zeros(1, T), ... % 风电出力 'P_cha', zeros(1, T), ... % 充电功率 'P_dis', zeros(1, T), ... % 放电功率 'SOC_ESS', zeros(1, T+1), ... % 电储能状态 'P_FC', zeros(1, T), ... % 燃料电池出力 'P_EL', zeros(1, T), ... % 电解槽用电 'SOC_HSS', zeros(1, T+1), ... % 氢储能状态 'Cost', zeros(1, N_windows) ... % 各窗口成本 ); % 初始状态 Results.SOC_ESS(1) = SOC_ESS_Init; Results.SOC_HSS(1) = SOC_HSS_Init; %% ==================== 主循环滚动优化 ==================== for idx = 1:N_windows % 窗口时间范围 k_start = (idx-1)*R_step + 1; k_end = min(k_start + W_size - 1, T); win_len = k_end - k_start + 1; fprintf('\n====== 滚动窗口 %d/%d [时段 %d-%d] ======\n', idx, N_windows, k_start, k_end); % 状态初始化 SOC_ESS_init = Results.SOC_ESS(k_start); SOC_HSS_init = Results.SOC_HSS(k_start); % ========== 定义决策变量 ========== P_GT = sdpvar(N_Gens, win_len, 'full'); % 机组出力 P_wt = sdpvar(1, win_len); % 风电出力 P_cha = sdpvar(1, win_len); % 充电功率 P_dis = sdpvar(1, win_len); % 放电功率 P_FC = sdpvar(1, win_len); % 燃料电池出力 P_EL = sdpvar(1, win_len); % 电解槽用电 SOC_ESS = sdpvar(1, win_len+1); % 电储能状态 SOC_HSS = sdpvar(1, win_len+1); % 氢储能状态 theta = sdpvar(N, win_len, 'full'); % 节点相角 V = sdpvar(N, win_len, 'full'); % 节点电压幅值 % 二进制变量 b_cha = binvar(1, win_len); % 充电状态 b_dis = binvar(1, win_len); % 放电状态 z_FC = binvar(1, win_len); % 燃料电池状态 z_EL = binvar(1, win_len); % 电解槽状态 % ========== 约束条件 ========== constraints = []; % 初始状态约束 constraints = [constraints, ... SOC_ESS(1) == SOC_ESS_init, ... SOC_HSS(1) == SOC_HSS_init]; % 机组约束 for t = 1:win_len for i = 1:N_Gens % 出力上下限 constraints = [constraints, ... B_ug_96(i, k_start+t-1)*Pg_min(i) <= P_GT(i,t) <= B_ug_96(i, k_start+t-1)*Pg_max(i)]; % 爬坡约束 if t == 1 && k_start > 1 constraints = [constraints, ... Results.P_GT(i, k_start-1) - P_GT(i,t) <= Ramp(i), ... P_GT(i,t) - Results.P_GT(i, k_start-1) <= Ramp(i)]; elseif t > 1 constraints = [constraints, ... P_GT(i,t-1) - P_GT(i,t) <= Ramp(i), ... P_GT(i,t) - P_GT(i,t-1) <= Ramp(i)]; end end % 风电约束 constraints = [constraints, 0 <= P_wt(t) <= P_wt_96(k_start+t-1)]; % 储能约束 constraints = [constraints, ... b_cha(t) + b_dis(t) <= 1, ... % 充放电互斥 0 <= P_cha(t) <= b_cha(t)*P_Cha_Max, ... 0 <= P_dis(t) <= b_dis(t)*P_Dis_Max]; % 氢系统约束 (修改为代码一的结构) % 互斥运行约束 constraints = [constraints, z_FC(t) + z_EL(t) <= 1]; % 功率约束 constraints = [constraints, ... 0 <= P_FC(t) <= P_FC_Max * z_FC(t), ... 0 <= P_EL(t) <= P_EL_Max * z_EL(t)]; % 充放氢速率计算 (kg/h) P_cha_H2 = P_EL(t) * K_P2H_xs; % 电解槽产氢率 P_dis_H2 = P_FC(t) / K_H2P_xs; % 燃料电池耗氢率 % 充放氢速率约束 constraints = [constraints, ... 0 <= P_cha_H2 <= H_cha_Max * z_EL(t), ... 0 <= P_dis_H2 <= H_dis_Max * z_FC(t)]; % 氢储能状态更新 (kg) constraints = [constraints, ... SOC_HSS(t+1) == SOC_HSS(t) + (P_cha_H2 - P_dis_H2) * dt]; end % 容量约束 constraints = [constraints, ... 0 <= SOC_ESS(2:end) <= C_ess_Max, ... C_HSS_Min <= SOC_HSS(2:end) <= C_HSS_Max]; % 电储能状态更新 for t = 1:win_len constraints = [constraints, ... SOC_ESS(t+1) == SOC_ESS(t) + K_ess_Eff*P_cha(t)*dt - P_dis(t)/K_ess_Eff*dt]; end % 潮流约束 for t = 1:win_len % 节点功率平衡 P_inj = Mbg * P_GT(:,t) + ... % 常规机组 Mbw * P_wt(t) + ... % 风电 Mbe * (P_dis(t) - P_cha(t)) + ... % 电储能 Mbs * (P_FC(t) - P_EL(t)); % 氢系统 P_net = P_inj - P_load_96(:, k_start+t-1); % 交流潮流约束 P_calc = Gbus * V(:,t) * baseMVA - (Bbus + diag(sum(Bbus,2))) * theta(:,t) * baseMVA; constraints = [constraints, P_calc == P_net]; % 线路功率约束 for k = 1:size(mpc.branch, 1) i = mpc.branch(k, 1); j = mpc.branch(k, 2); P_flow = -Gbus(i,j)*(V(i,t) - V(j,t)) + Bbus(i,j)*(theta(i,t) - theta(j,t)); constraints = [constraints, -500 <= P_flow*baseMVA <= 500]; end % 电压和相角约束 constraints = [constraints, ... 0.95 <= V(:,t) <= 1.05, ... -pi <= theta(:,t) <= pi]; end % ========== 目标函数 (加入环境成本) ========== objective = 0; for t = 1:win_len % 全局时段索引 t_global = k_start + t - 1; % 机组成本 for i = 1:N_Gens fuel_cost = a2(i)*(baseMVA*P_GT(i,t))^2 + ... a1(i)*(baseMVA*P_GT(i,t)) + a0(i); objective = objective + fuel_cost; % 启动成本 (固定成本) objective = objective + B_ug_96(i, k_start+t-1)*StartupCost(i); % 碳排放成本 (新增) carbon_emission = Carbon_Factor(i) * baseMVA * P_GT(i,t) * dt; % 吨CO2 objective = objective + C_carbon * carbon_emission; end % 弃风惩罚 objective = objective + kc*(P_wt_96(t_global) - P_wt(t)); % 储能运维成本 objective = objective + O_M_ESS(1)*P_cha(t) + O_M_ESS(2)*P_dis(t); % 氢系统运维成本 objective = objective + O_M_H2(1)*P_EL(t) + O_M_H2(2)*P_FC(t); % 氢系统环境补贴 (新增) % 燃料电池补贴 (高峰时段1.3倍) if ismember(t_global, t_peak) objective = objective + 1.3 * C_FC_env * P_FC(t) * dt; else objective = objective + C_FC_env * P_FC(t) * dt; end % 电解槽补贴 (低谷时段1.5倍) if ismember(t_global, t_valley) objective = objective + 1.5 * C_EL_env * P_EL(t) * dt; else objective = objective + C_EL_env * P_EL(t) * dt; end end % ========== 求解优化问题 ========== options = sdpsettings('solver', 'cplex', 'verbose', 1); sol = optimize(constraints, objective, options); % ========== 结果处理 ========== if sol.problem == 0 fprintf('优化成功! 窗口成本: %.2f\n', value(objective)); Results.Cost(idx) = value(objective); % 保存当前窗口结果 valid_len = min(R_step, T - k_start + 1); idx_save = k_start:(k_start + valid_len - 1); Results.P_GT(:, idx_save) = value(P_GT(:, 1:valid_len)); Results.P_wt(idx_save) = value(P_wt(1:valid_len)); Results.P_cha(idx_save) = value(P_cha(1:valid_len)); Results.P_dis(idx_save) = value(P_dis(1:valid_len)); Results.P_FC(idx_save) = value(P_FC(1:valid_len)); Results.P_EL(idx_save) = value(P_EL(1:valid_len)); % 更新状态量 Results.SOC_ESS(k_start+1:k_start+win_len) = value(SOC_ESS(2:end)); Results.SOC_HSS(k_start+1:k_start+win_len) = value(SOC_HSS(2:end)); else fprintf('优化失败: %s\n', sol.info); % 使用前一时段值填充当前窗口 valid_len = min(R_step, T - k_start + 1); idx_save = k_start:(k_start + valid_len - 1); Results.P_GT(:, idx_save) = repmat(Results.P_GT(:, max(1,k_start-1)), 1, valid_len); Results.P_wt(idx_save) = P_wt_96(idx_save); Results.Cost(idx) = Inf; end end %% ==================== 结果可视化 ==================== time_vec = (1:T)*0.25; % 时间向量 (小时) % 1. 功率平衡图 figure('Position', [100, 100, 1200, 600]); total_gen = sum(Results.P_GT, 1) + Results.P_wt + Results.P_dis - Results.P_cha + ... Results.P_FC - Results.P_EL; total_load = sum(P_load_96, 1); plot(time_vec, total_gen, 'b-', 'LineWidth', 1.5); hold on; plot(time_vec, total_load, 'r--', 'LineWidth', 1.5); title('系统功率平衡'); xlabel('时间 (小时)'); ylabel('功率 (MW)'); legend('总发电量', '总负荷', 'Location', 'best'); grid on; xlim([0 24]); xticks(0:2:24); % 2. 机组出力曲线 figure('Position', [100, 100, 1200, 800]); subplot(3,1,1); plot(time_vec, Results.P_GT', 'LineWidth', 1.5); title('常规机组出力'); xlabel('时间 (小时)'); ylabel('出力 (MW)'); legend('G1','G2','G3','G4','G5','G6'); grid on; xlim([0 24]); % 3. 新能源与储能 subplot(3,1,2); plot(time_vec, Results.P_wt, 'g-', 'LineWidth', 1.5); hold on; plot(time_vec, P_wt_96, 'g:'); title('风电出力'); xlabel('时间 (小时)'); ylabel('功率 (MW)'); legend('实际出力', '预测出力'); grid on; xlim([0 24]); subplot(3,1,3); plot(time_vec, Results.P_cha, 'b', time_vec, Results.P_dis, 'r', ... time_vec, Results.P_FC, 'm', time_vec, Results.P_EL, 'c', 'LineWidth', 1.5); title('储能系统功率'); xlabel('时间 (小时)'); ylabel('功率 (MW)'); legend('充电', '放电', '燃料电池', '电解槽'); grid on; xlim([0 24]); % 4. 储能状态 figure('Position', [100, 100, 1200, 400]); plot(time_vec, Results.SOC_ESS(1:T), 'b-', 'LineWidth', 1.5); hold on; plot(time_vec, Results.SOC_HSS(1:T), 'r-', 'LineWidth', 1.5); title('储能系统状态'); xlabel('时间 (小时)'); ylabel('状态'); legend('电储能 (MWh)', '氢储能 (kg)'); grid on; xlim([0 24]); %% ==================== 结果分析 ==================== total_cost = sum(Results.Cost(Results.Cost < Inf)); fprintf('总优化成本: %.2f\n', total_cost); % 保存结果 save('optimization_results.mat', 'Results'); disp('优化结果已保存'); %% 基于电-氢混合储能的电力系统优化调度模型（日内滚动调度） clc clear close all %% ==================== 参数设置 ==================== T = 96; % 总时段数 (15分钟间隔) W_size = 16; % 滚动窗口大小 (4小时) R_step = 4; % 滚动步长 (1小时) N_windows = floor((T - W_size)/R_step) + 1; % 滚动次数 N_Gens = 6; % 机组数量 dt = 0.25; % 时间步长 (小时) %% ==================== 数据加载与预处理 ==================== load pload_24.mat p_load_1 load pwt_24.mat p_wt load ug_24.mat ugone % 扩展24小时数据到96小时 P_wt_96 = repelem(p_wt, 1, 4); % 风电预测 P_load_96 = zeros(30, T); % 节点负荷 B_ug_96 = zeros(6, T); % 机组状态 for i = 1:30 P_load_96(i, :) = repelem(p_load_1(i, :), 1, 4); end for i = 1:6 B_ug_96(i, :) = repelem(ugone(i, :), 1, 4); end %% ==================== 系统参数设置 ==================== mpc = loadcase('case30'); baseMVA = mpc.baseMVA; % 系统基准功率 % 常规机组参数 Gen_Params = [1 250 50 0.0375 20 372.5 72; 2 80 20 0.175 17.5 352.3 48; 3 50 15 0.625 10 316.5 30; 4 35 10 0.0834 32.5 329.2 21; 5 30 10 0.25 30 276.4 18; 6 40 12 0.25 30 232.2 24]; % 提取发电机参数 a2 = Gen_Params(:,4); % 二次项系数 a1 = Gen_Params(:,5); % 一次项系数 a0 = Gen_Params(:,6); % 常数项 Ramp = Gen_Params(:,7); % 爬坡率 (MW/15min) Pg_min = Gen_Params(:,3); % 最小出力 Pg_max = Gen_Params(:,2); % 最大出力 % 启动成本参数 StartupCost = [2768; 2874; 2136; 2677; 2262; 2869]; % 环境成本系数 (新增) C_carbon = 150; % 碳排放成本 (元/吨CO2) C_FC_env = -500; % 燃料电池环境补贴 (负成本) 元/MW C_EL_env = -300; % 电解槽环境补贴 (负成本) 元/MW Carbon_Factor = [0.85; 0.92; 1.05; 0.78; 0.95; 0.88]; % 碳排放系数 (吨CO2/MWh) % 峰谷时段定义 (15分钟分辨率) (新增) t_peak_hour = [12:14, 19:21]; % 高峰时段（小时） t_valley_hour = [1:5, 24]; % 低谷时段（小时） t_peak = []; for h = t_peak_hour t_peak = [t_peak, (h-1)*4+1 : h*4]; end t_valley = []; for h = t_valley_hour t_valley = [t_valley, (h-1)*4+1 : h*4]; end % 电储能参数 P_Cha_Max = 45; % 最大充电功率 (MW) P_Dis_Max = 90; % 最大放电功率 (MW) K_ess_Eff = 0.92; % 充放电效率 C_ess_Max = 300; % 最大容量 (MWh) SOC_ESS_Init = 10; % 初始SOC (MWh) % 氢系统参数 (更新单位) P_FC_Max = 400; % 燃料电池最大功率 (MW) K_H2P_xs = 9.1; % 氢转电系数 (kg/MWh) (原名K_H2P) P_EL_Max = 400; % 电解槽最大功率 (MW) K_P2H_xs = 0.22; % 电转氢系数 (kg/MWh) (原名K_P2H) H_cha_Max = 100; % 最大充氢速率 (kg/h) H_dis_Max = 100; % 最大放氢速率 (kg/h) C_HSS_Max = 400; % 最大储氢容量 (kg) C_HSS_Min = 0; % 最小储氢容量 (kg) SOC_HSS_Init = 20; % 初始氢储量 (kg) % 成本系数 kc = 0.5; % 弃风惩罚因子 O_M_ESS = [0.56, 0.76]; % 电储能运维成本 [充电, 放电] O_M_H2 = [0.65, 0.83]; % 氢系统运维成本 [电解槽, 燃料电池] % 节点位置 Bus_GEN = [1 2 5 8 11 13]; % 火电机组节点 Bus_WT = 2; % 风电节点 Bus_ESS = 11; % 电储能节点 Bus_HSS = 8; % 氢储能节点 %% ==================== 网络拓扑初始化 ==================== [Ybus, ~, ~] = makeYbus(mpc); Gbus = real(Ybus); % 电导矩阵 Bbus = imag(Ybus); % 电纳矩阵 % 创建节点-设备映射 N = size(mpc.bus, 1); % 节点数量 Mbg = sparse(Bus_GEN, 1:N_Gens, 1, N, N_Gens); % 常规机组 Mbw = sparse(Bus_WT, 1, 1, N, 1); % 风电 Mbe = sparse(Bus_ESS, 1, 1, N, 1); % 电储能 Mbs = sparse(Bus_HSS, 1, 1, N, 1); % 氢系统 %% ==================== 状态初始化 ==================== % 结果存储 Results = struct(... 'P_GT', zeros(N_Gens, T), ... % 机组出力 'P_wt', zeros(1, T), ... % 风电出力 'P_cha', zeros(1, T), ... % 充电功率 'P_dis', zeros(1, T), ... % 放电功率 'SOC_ESS', zeros(1, T+1), ... % 电储能状态 'P_FC', zeros(1, T), ... % 燃料电池出力 'P_EL', zeros(1, T), ... % 电解槽用电 'SOC_HSS', zeros(1, T+1), ... % 氢储能状态 'Cost', zeros(1, N_windows) ... % 各窗口成本 ); % 初始状态 Results.SOC_ESS(1) = SOC_ESS_Init; Results.SOC_HSS(1) = SOC_HSS_Init; %% ==================== 主循环滚动优化 ==================== for idx = 1:N_windows % 窗口时间范围 k_start = (idx-1)*R_step + 1; k_end = min(k_start + W_size - 1, T); win_len = k_end - k_start + 1; fprintf('\n====== 滚动窗口 %d/%d [时段 %d-%d] ======\n', idx, N_windows, k_start, k_end); % 状态初始化 SOC_ESS_init = Results.SOC_ESS(k_start); SOC_HSS_init = Results.SOC_HSS(k_start); % ========== 定义决策变量 ========== P_GT = sdpvar(N_Gens, win_len, 'full'); % 机组出力 P_wt = sdpvar(1, win_len); % 风电出力 P_cha = sdpvar(1, win_len); % 充电功率 P_dis = sdpvar(1, win_len); % 放电功率 P_FC = sdpvar(1, win_len); % 燃料电池出力 P_EL = sdpvar(1, win_len); % 电解槽用电 SOC_ESS = sdpvar(1, win_len+1); % 电储能状态 SOC_HSS = sdpvar(1, win_len+1); % 氢储能状态 theta = sdpvar(N, win_len, 'full'); % 节点相角 V = sdpvar(N, win_len, 'full'); % 节点电压幅值 % 二进制变量 b_cha = binvar(1, win_len); % 充电状态 b_dis = binvar(1, win_len); % 放电状态 z_FC = binvar(1, win_len); % 燃料电池状态 z_EL = binvar(1, win_len); % 电解槽状态 % ========== 约束条件 ========== constraints = []; % 初始状态约束 constraints = [constraints, ... SOC_ESS(1) == SOC_ESS_init, ... SOC_HSS(1) == SOC_HSS_init]; % 机组约束 for t = 1:win_len for i = 1:N_Gens % 出力上下限 constraints = [constraints, ... B_ug_96(i, k_start+t-1)*Pg_min(i) <= P_GT(i,t) <= B_ug_96(i, k_start+t-1)*Pg_max(i)]; % 爬坡约束 if t == 1 && k_start > 1 constraints = [constraints, ... Results.P_GT(i, k_start-1) - P_GT(i,t) <= Ramp(i), ... P_GT(i,t) - Results.P_GT(i, k_start-1) <= Ramp(i)]; elseif t > 1 constraints = [constraints, ... P_GT(i,t-1) - P_GT(i,t) <= Ramp(i), ... P_GT(i,t) - P_GT(i,t-1) <= Ramp(i)]; end end % 风电约束 constraints = [constraints, 0 <= P_wt(t) <= P_wt_96(k_start+t-1)]; % 储能约束 constraints = [constraints, ... b_cha(t) + b_dis(t) <= 1, ... % 充放电互斥 0 <= P_cha(t) <= b_cha(t)*P_Cha_Max, ... 0 <= P_dis(t) <= b_dis(t)*P_Dis_Max]; % 氢系统约束 (修改为代码一的结构) % 互斥运行约束 constraints = [constraints, z_FC(t) + z_EL(t) <= 1]; % 功率约束 constraints = [constraints, ... 0 <= P_FC(t) <= P_FC_Max * z_FC(t), ... 0 <= P_EL(t) <= P_EL_Max * z_EL(t)]; % 充放氢速率计算 (kg/h) P_cha_H2 = P_EL(t) * K_P2H_xs; % 电解槽产氢率 P_dis_H2 = P_FC(t) / K_H2P_xs; % 燃料电池耗氢率 % 充放氢速率约束 constraints = [constraints, ... 0 <= P_cha_H2 <= H_cha_Max * z_EL(t), ... 0 <= P_dis_H2 <= H_dis_Max * z_FC(t)]; % 氢储能状态更新 (kg) constraints = [constraints, ... SOC_HSS(t+1) == SOC_HSS(t) + (P_cha_H2 - P_dis_H2) * dt]; end % 容量约束 constraints = [constraints, ... 0 <= SOC_ESS(2:end) <= C_ess_Max, ... C_HSS_Min <= SOC_HSS(2:end) <= C_HSS_Max]; % 电储能状态更新 for t = 1:win_len constraints = [constraints, ... SOC_ESS(t+1) == SOC_ESS(t) + K_ess_Eff*P_cha(t)*dt - P_dis(t)/K_ess_Eff*dt]; end % 潮流约束 for t = 1:win_len % 节点功率平衡 P_inj = Mbg * P_GT(:,t) + ... % 常规机组 Mbw * P_wt(t) + ... % 风电 Mbe * (P_dis(t) - P_cha(t)) + ... % 电储能 Mbs * (P_FC(t) - P_EL(t)); % 氢系统 P_net = P_inj - P_load_96(:, k_start+t-1); % 交流潮流约束 P_calc = Gbus * V(:,t) * baseMVA - (Bbus + diag(sum(Bbus,2))) * theta(:,t) * baseMVA; constraints = [constraints, P_calc == P_net]; % 线路功率约束 for k = 1:size(mpc.branch, 1) i = mpc.branch(k, 1); j = mpc.branch(k, 2); P_flow = -Gbus(i,j)*(V(i,t) - V(j,t)) + Bbus(i,j)*(theta(i,t) - theta(j,t)); constraints = [constraints, -500 <= P_flow*baseMVA <= 500]; end % 电压和相角约束 constraints = [constraints, ... 0.95 <= V(:,t) <= 1.05, ... -pi <= theta(:,t) <= pi]; end % ========== 目标函数 (加入环境成本) ========== objective = 0; for t = 1:win_len % 全局时段索引 t_global = k_start + t - 1; % 机组成本 for i = 1:N_Gens fuel_cost = a2(i)*(baseMVA*P_GT(i,t))^2 + ... a1(i)*(baseMVA*P_GT(i,t)) + a0(i); objective = objective + fuel_cost; % 启动成本 (固定成本) objective = objective + B_ug_96(i, k_start+t-1)*StartupCost(i); % 碳排放成本 (新增) carbon_emission = Carbon_Factor(i) * baseMVA * P_GT(i,t) * dt; % 吨CO2 objective = objective + C_carbon * carbon_emission; end % 弃风惩罚 objective = objective + kc*(P_wt_96(t_global) - P_wt(t)); % 储能运维成本 objective = objective + O_M_ESS(1)*P_cha(t) + O_M_ESS(2)*P_dis(t); % 氢系统运维成本 objective = objective + O_M_H2(1)*P_EL(t) + O_M_H2(2)*P_FC(t); % 氢系统环境补贴 (新增) % 燃料电池补贴 (高峰时段1.3倍) if ismember(t_global, t_peak) objective = objective + 1.3 * C_FC_env * P_FC(t) * dt; else objective = objective + C_FC_env * P_FC(t) * dt; end % 电解槽补贴 (低谷时段1.5倍) if ismember(t_global, t_valley) objective = objective + 1.5 * C_EL_env * P_EL(t) * dt; else objective = objective + C_EL_env * P_EL(t) * dt; end end % ========== 求解优化问题 ========== options = sdpsettings('solver', 'cplex', 'verbose', 1); sol = optimize(constraints, objective, options); % ========== 结果处理 ========== if sol.problem == 0 fprintf('优化成功! 窗口成本: %.2f\n', value(objective)); Results.Cost(idx) = value(objective); % 保存当前窗口结果 valid_len = min(R_step, T - k_start + 1); idx_save = k_start:(k_start + valid_len - 1); Results.P_GT(:, idx_save) = value(P_GT(:, 1:valid_len)); Results.P_wt(idx_save) = value(P_wt(1:valid_len)); Results.P_cha(idx_save) = value(P_cha(1:valid_len)); Results.P_dis(idx_save) = value(P_dis(1:valid_len)); Results.P_FC(idx_save) = value(P_FC(1:valid_len)); Results.P_EL(idx_save) = value(P_EL(1:valid_len)); % 更新状态量 Results.SOC_ESS(k_start+1:k_start+win_len) = value(SOC_ESS(2:end)); Results.SOC_HSS(k_start+1:k_start+win_len) = value(SOC_HSS(2:end)); else fprintf('优化失败: %s\n', sol.info); % 使用前一时段值填充当前窗口 valid_len = min(R_step, T - k_start + 1); idx_save = k_start:(k_start + valid_len - 1); Results.P_GT(:, idx_save) = repmat(Results.P_GT(:, max(1,k_start-1)), 1, valid_len); Results.P_wt(idx_save) = P_wt_96(idx_save); Results.Cost(idx) = Inf; end end %% ==================== 结果可视化 ==================== time_vec = (1:T)*0.25; % 时间向量 (小时) % 1. 功率平衡图 figure('Position', [100, 100, 1200, 600]); total_gen = sum(Results.P_GT, 1) + Results.P_wt + Results.P_dis - Results.P_cha + ... Results.P_FC - Results.P_EL; total_load = sum(P_load_96, 1); plot(time_vec, total_gen, 'b-', 'LineWidth', 1.5); hold on; plot(time_vec, total_load, 'r--', 'LineWidth', 1.5); title('系统功率平衡'); xlabel('时间 (小时)'); ylabel('功率 (MW)'); legend('总发电量', '总负荷', 'Location', 'best'); grid on; xlim([0 24]); xticks(0:2:24); % 2. 机组出力曲线 figure('Position', [100, 100, 1200, 800]); subplot(3,1,1); plot(time_vec, Results.P_GT', 'LineWidth', 1.5); title('常规机组出力'); xlabel('时间 (小时)'); ylabel('出力 (MW)'); legend('G1','G2','G3','G4','G5','G6'); grid on; xlim([0 24]); % 3. 新能源与储能 subplot(3,1,2); plot(time_vec, Results.P_wt, 'g-', 'LineWidth', 1.5); hold on; plot(time_vec, P_wt_96, 'g:'); title('风电出力'); xlabel('时间 (小时)'); ylabel('功率 (MW)'); legend('实际出力', '预测出力'); grid on; xlim([0 24]); subplot(3,1,3); plot(time_vec, Results.P_cha, 'b', time_vec, Results.P_dis, 'r', ... time_vec, Results.P_FC, 'm', time_vec, Results.P_EL, 'c', 'LineWidth', 1.5); title('储能系统功率'); xlabel('时间 (小时)'); ylabel('功率 (MW)'); legend('充电', '放电', '燃料电池', '电解槽'); grid on; xlim([0 24]); % 4. 储能状态 figure('Position', [100, 100, 1200, 400]); plot(time_vec, Results.SOC_ESS(1:T), 'b-', 'LineWidth', 1.5); hold on; plot(time_vec, Results.SOC_HSS(1:T), 'r-', 'LineWidth', 1.5); title('储能系统状态'); xlabel('时间 (小时)'); ylabel('状态'); legend('电储能 (MWh)', '氢储能 (kg)'); grid on; xlim([0 24]); %% ==================== 结果分析 ==================== total_cost = sum(Results.Cost(Results.Cost < Inf)); fprintf('总优化成本: %.2f\n', total_cost); % 保存结果 save('optimization_results.mat', 'Results'); disp('优化结果已保存'); %% 基于电-氢混合储能的电力系统优化调度模型（日内滚动调度） clc clear close all %% ==================== 参数设置 ==================== T = 96; % 总时段数 (15分钟间隔) W_size = 16; % 滚动窗口大小 (4小时) R_step = 4; % 滚动步长 (1小时) N_windows = floor((T - W_size)/R_step) + 1; % 滚动次数 N_Gens = 6; % 机组数量 dt = 0.25; % 时间步长 (小时) %% ==================== 数据加载 ==================== % 电负荷数据 (小时级) P_load_Hourly = 0.25*[297 258 220 208 226 258 280 330 372 364 400 468 636 600 506 554 392 392 412 403 426 390 322 268]; % 风电数据 (小时级) P_wt_Hourly = [50 55 54 58 60 100 102 125 130 152 153 140 157 130 127 300 275 260 210 200 169 150 130 100]; % 转换为15分钟数据 (96时段) P_load_Min_96 = repelem(P_load_Hourly, 1, 4); % 系统总负荷 P_wt_Min_96 = repelem(P_wt_Hourly, 1, 4); % 风电出力 % 加载电力系统案例 mpc = loadcase('case30'); 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/小时) Ramp_15min = Ramp_e * dt; % 15分钟爬坡率 (MW/时段) Carbon_Factor = Gen_Params(:,8); % 碳排放系数 (tCO2/MWh) P_GT_Min = Gen_Params(:,3); % 最小出力 P_GT_Max = Gen_Params(:,2); % 最大出力 % 启动/停机成本参数 (单位: 元) C_Start_qt = [2768; 2874; 2136; 2677; 2262; 2869]; % 启动成本 C_Stop_qt = [1500; 1600; 1200; 1400; 1300; 1500]; % 停机成本 % 电储能参数 P_Cha_Max = 45; % 最大充电功率 (MW) P_Dis_Max = 90; % 最大放电功率 (MW) K_ess_Eff = 0.92; % 充放电效率 C_ess_Max = 300; % 最大容量 (MWh) SOC_ESS_Init = 10; % 初始SOC (MWh) % 氢系统参数 P_FC_Max = 400; % 燃料电池最大功率 (MW) K_H2P_xs = 9.1; % 氢转电系数 (kg/MWh) P_EL_Max = 400; % 电解槽最大功率 (MW) K_P2H_xs = 0.22; % 电转氢系数 (kg/MWh) H_cha_Max = 100; % 最大充氢速率 (kg/h) H_dis_Max = 100; % 最大放氢速率 (kg/h) C_HSS_Max = 400; % 最大储氢容量 (kg) C_HSS_Min = 10; % 最小储氢容量 (kg) SOC_HSS_Init = 20; % 初始氢储量 (kg) % 环境成本系数 C_carbon = 150; % 碳排放成本 (元/吨CO2) % 节点位置 Bus_GEN_e = [1 2 5 8 11 13]; % 火电机组节点 Bus_WT_e = 2; % 风电节点 Bus_ESS_e = 11; % 电储能节点 Bus_HSS_h = 8; % 氢储能节点 %% ==================== 节点电负荷矩阵生成 ==================== % 电力系统基准负荷 (MW) P_Load_Base = sum(mpc.bus(:,3)); % 创建节点负荷矩阵 (30×96) P_load_96 = zeros(size(mpc.bus, 1), T); for i = 1:size(mpc.bus, 1) % 计算节点负荷比例 node_ratio = mpc.bus(i,3) / P_Load_Base; % 生成节点负荷曲线 P_load_96(i,:) = node_ratio * P_load_Min_96; end %% ==================== 电网络拓扑初始化 ==================== [Ybus, ~, ~] = makeYbus(mpc); B = imag(Ybus); % 电纳矩阵 ref_node = 1; % 参考节点 non_ref = 2:size(mpc.bus,1); % 非参考节点 % 创建设备-节点映射矩阵 N = size(mpc.bus, 1); % 节点数量 G_map = sparse(Bus_GEN_e, 1:N_Gens, 1, N, N_Gens); % 常规机组 W_map = sparse(Bus_WT_e, 1, 1, N, 1); % 风电机组 E_map = sparse(Bus_ESS_e, 1, 1, N, 1); % 电储能 H_map = sparse(Bus_HSS_h, 1, 1, N, 1); % 氢系统 %% ==================== 状态初始化 ==================== % 结果存储 Results = struct(... 'P_GT', zeros(N_Gens, T), ... % 机组出力 'P_wt', zeros(1, T), ... % 风电出力 'P_cha', zeros(1, T), ... % 充电功率 'P_dis', zeros(1, T), ... % 放电功率 'SOC_ESS', SOC_ESS_Init*ones(1, T+1), ... % 电储能状态 'P_FC', zeros(1, T), ... % 燃料电池出力 'P_EL', zeros(1, T), ... % 电解槽用电 'SOC_HSS', SOC_HSS_Init*ones(1, T+1), ... % 氢储能状态 'Obj', zeros(1, N_windows), ... % 各窗口目标值 'Ramp_Violation', zeros(1, T)... % 爬坡约束违反检测 ); % 机组状态跟踪 B_ug_96 = zeros(N_Gens, T); % 机组状态历史 B_ug_96(:,1) = ones(N_Gens,1); % 初始全运行 %% ==================== 主循环滚动优化 ==================== % 生成随机扰动种子 (确保可重复性) rng(2023, 'twister'); for idx = 1:N_windows % 窗口时间范围 k_start = (idx-1)*R_step + 1; k_end = min(k_start + W_size - 1, T); win_len = k_end - k_start + 1; fprintf('\n====== 滚动窗口 %d/%d [时段 %d-%d] ======\n', idx, N_windows, k_start, k_end); % 状态初始化 if k_start == 1 SOC_ESS_init = SOC_ESS_Init; SOC_HSS_init = SOC_HSS_Init; P_GT_prev = P_GT_Min(:); % 初始状态为最小出力 prev_status = ones(N_Gens,1); % 初始状态全运行 else % 使用上一个窗口结束时的状态 SOC_ESS_init = Results.SOC_ESS(k_start); SOC_HSS_init = Results.SOC_HSS(k_start); P_GT_prev = Results.P_GT(:, k_start); % 取k_start时刻状态 prev_status = B_ug_96(:, k_start); % 取k_start时刻状态 end % 定义决策变量 (当前窗口) P_GT_e = sdpvar(N_Gens, win_len, 'full'); % 机组出力 B_U_G = binvar(N_Gens, win_len, 'full'); % 机组运行状态 B_U_Start = binvar(N_Gens, win_len, 'full');% 启动状态 B_U_Stop = binvar(N_Gens, win_len, 'full'); % 停机状态 P_wt_e = sdpvar(1, win_len); % 风电出力 P_cha_e = sdpvar(1, win_len); % 充电功率 P_dis_e = sdpvar(1, win_len); % 放电功率 SOC_ESS = sdpvar(1, win_len+1); % 电储能状态 P_FC_e = sdpvar(1, win_len); % 燃料电池出力 P_e_EL = sdpvar(1, win_len); % 电解槽用电 SOC_HSS = sdpvar(1, win_len+1); % 氢储能状态 Bus_Theta = sdpvar(N, win_len, 'full'); % 节点相角 % 二进制变量 (充放电状态) b_in_cha = binvar(1, win_len); b_in_dis = binvar(1, win_len); % 约束条件 Cons = []; %% ======== 初始状态约束 ======== Cons = [Cons, SOC_ESS(1) == SOC_ESS_init]; Cons = [Cons, SOC_HSS(1) == SOC_HSS_init]; %% ======== 机组约束 ======== for t = 1:win_len % 出力上下限 Cons = [Cons, B_U_G(:,t).*P_GT_Min <= P_GT_e(:,t) <= B_U_G(:,t).*P_GT_Max]; % 爬坡约束 (动态参考状态机制) if t == 1 ref_state = P_GT_prev; % 窗口起始状态 else ref_state = P_GT_e(:,t-1); % 前一时段状态 end % 爬坡约束计算 delta_P = P_GT_e(:,t) - ref_state; Cons = [Cons, -Ramp_15min <= delta_P <= Ramp_15min]; % 启停逻辑 if t == 1 Cons = [Cons, B_U_Start(:,1) >= B_U_G(:,1) - prev_status]; Cons = [Cons, B_U_Start(:,1) <= B_U_G(:,1)]; Cons = [Cons, B_U_Start(:,1) <= 1 - prev_status]; Cons = [Cons, B_U_Stop(:,1) >= prev_status - B_U_G(:,1)]; Cons = [Cons, B_U_Stop(:,1) <= prev_status]; Cons = [Cons, B_U_Stop(:,1) <= 1 - B_U_G(:,1)]; else Cons = [Cons, B_U_Start(:,t) >= B_U_G(:,t) - B_U_G(:,t-1)]; Cons = [Cons, B_U_Start(:,t) <= B_U_G(:,t)]; Cons = [Cons, B_U_Start(:,t) <= 1 - B_U_G(:,t-1)]; Cons = [Cons, B_U_Stop(:,t) >= B_U_G(:,t-1) - B_U_G(:,t)]; Cons = [Cons, B_U_Stop(:,t) <= B_U_G(:,t-1)]; Cons = [Cons, B_U_Stop(:,t) <= 1 - B_U_G(:,t)]; end end %% ======== 风电约束 ======== Cons = [Cons, 0 <= P_wt_e <= P_wt_Min_96(k_start:k_end)]; %% ======== 电储能约束 ======== Cons = [Cons, b_in_cha + b_in_dis <= 1]; % 充放电互斥 for t = 1:win_len % SOC更新 (加入时间步长0.25小时) Cons = [Cons, SOC_ESS(t+1) == SOC_ESS(t) + dt*(K_ess_Eff*P_cha_e(t) - P_dis_e(t)/K_ess_Eff)]; % 容量约束 Cons = [Cons, 0 <= SOC_ESS(t+1) <= C_ess_Max]; % 功率约束 Cons = [Cons, 0 <= P_cha_e(t) <= b_in_cha(t)*P_Cha_Max]; Cons = [Cons, 0 <= P_dis_e(t) <= b_in_dis(t)*P_Dis_Max]; end %% ======== 氢系统约束 ======== for t = 1:win_len % 产氢量和耗氢量计算 H_prod = P_e_EL(t) * dt * K_P2H_xs; % 产氢量 (kg) H_cons = P_FC_e(t) * dt * K_H2P_xs; % 耗氢量 (kg) % SOC更新 Cons = [Cons, SOC_HSS(t+1) == SOC_HSS(t) + H_prod - H_cons]; % 容量约束 Cons = [Cons, C_HSS_Min <= SOC_HSS(t+1) <= C_HSS_Max]; % 功率约束 Cons = [Cons, 0 <= P_FC_e(t) <= P_FC_Max]; Cons = [Cons, 0 <= P_e_EL(t) <= P_EL_Max]; Cons = [Cons, 0 <= H_prod <= H_cha_Max*dt]; Cons = [Cons, 0 <= H_cons <= H_dis_Max*dt]; end %% ======== 网络约束 ======== % 线路约束参数 branch = mpc.branch; rateA = branch(:, 6); % 线路容量限制 (MW) rateA(rateA == 0) = max(rateA(rateA>0))*1.5; % 处理零值容量限制 for t = 1:win_len t_global = k_start + t - 1; % 节点注入功率 P_inj = G_map * P_GT_e(:, t) + ... % 常规机组 W_map * P_wt_e(t) + ... % 风电 E_map * (P_dis_e(t) - P_cha_e(t)) + ... % 电储能 H_map * (P_FC_e(t) - P_e_EL(t));% 氢系统 % 节点净功率 P_net = P_inj - P_load_96(:, t_global); % 直流潮流约束 Cons = [Cons, B(non_ref, :) * Bus_Theta(:, t) == P_net(non_ref)/baseMVA]; % 线路潮流约束 for br = 1:size(branch, 1) from = branch(br, 1); to = branch(br, 2); x = branch(br, 4); % 线路电抗 P_flow = (Bus_Theta(from,t) - Bus_Theta(to,t)) / x * baseMVA; Cons = [Cons, -rateA(br) <= P_flow <= rateA(br)]; end end Cons = [Cons, Bus_Theta(ref_node, :) == 0]; % 参考节点 %% ======== 目标函数 ======== Obj = 0; for t = 1:win_len % 燃料成本 for g = 1:N_Gens fuel_cost = a2(g)*P_GT_e(g,t)^2 + a1(g)*P_GT_e(g,t) + a0(g); Obj = Obj + dt * fuel_cost; end % 启停成本 for g = 1:N_Gens Obj = Obj + C_Start_qt(g)*B_U_Start(g,t) + C_Stop_qt(g)*B_U_Stop(g,t); end % 碳排放成本 for g = 1:N_Gens carbon_cost = C_carbon * Carbon_Factor(g) * P_GT_e(g,t)/1000; % 元/吨 Obj = Obj + dt * carbon_cost; end end %% ======== 求解优化问题 ======== opts = sdpsettings('solver', 'cplex', ... 'verbose', 1, ... 'showprogress', 1, ... 'cplex.timelimit', 120); % 增加求解时限 sol = optimize(Cons, Obj, opts); %% ======== 结果处理 ======== if sol.problem == 0 fprintf('优化成功 | 窗口成本: %.2f 元\n', value(Obj)); Results.Obj(idx) = value(Obj); % 计算实际可保存长度 valid_len = min(R_step, T - k_start + 1); idx_range = k_start:k_start+valid_len-1; % 保存当前窗口结果 Results.P_GT(:, idx_range) = value(P_GT_e(:,1:valid_len)); B_ug_96(:, idx_range) = value(B_U_G(:,1:valid_len)); Results.P_wt(idx_range) = value(P_wt_e(1:valid_len)); Results.P_cha(idx_range) = value(P_cha_e(1:valid_len)); Results.P_dis(idx_range) = value(P_dis_e(1:valid_len)); Results.P_FC(idx_range) = value(P_FC_e(1:valid_len)); Results.P_EL(idx_range) = value(P_e_EL(1:valid_len)); % 更新状态量 Results.SOC_ESS(k_start+1:k_end+1) = value(SOC_ESS(2:end)); Results.SOC_HSS(k_start+1:k_end+1) = value(SOC_HSS(2:end)); % 爬坡约束验证 for t_idx = 1:valid_len t_global = k_start + t_idx - 1; if t_global == 1 delta_P = Results.P_GT(:,1) - P_GT_Min; else delta_P = Results.P_GT(:,t_global) - Results.P_GT(:,t_global-1); end % 检测爬坡约束违反 ramp_violation = max(abs(delta_P) - Ramp_15min); Results.Ramp_Violation(t_global) = any(ramp_violation > 1e-4); end else fprintf('优化失败: %s\n', sol.info); valid_len = min(R_step, T - k_start + 1); idx_range = k_start:k_start+valid_len-1; % 失败处理（使用上一个窗口结束状态） if k_start == 1 fill_value = repmat(P_GT_Min(:), 1, valid_len); wt_fill = P_wt_Min_96(idx_range); soc_ess_fill = SOC_ESS_Init * ones(1, valid_len); soc_hss_fill = SOC_HSS_Init * ones(1, valid_len); else fill_value = repmat(Results.P_GT(:, k_start), 1, valid_len); wt_fill = Results.P_wt(k_start) * ones(1, valid_len); soc_ess_fill = Results.SOC_ESS(k_start) * ones(1, valid_len); soc_hss_fill = Results.SOC_HSS(k_start) * ones(1, valid_len); end Results.P_GT(:, idx_range) = fill_value; Results.P_wt(idx_range) = wt_fill; Results.P_cha(idx_range) = 0; Results.P_dis(idx_range) = 0; Results.P_FC(idx_range) = 0; Results.P_EL(idx_range) = 0; Results.SOC_ESS(idx_range+1) = soc_ess_fill; Results.SOC_HSS(idx_range+1) = soc_hss_fill; Results.Obj(idx) = Inf; end end %% ==================== 状态转移验证 ==================== % 精确验证爬坡约束 for t = 1:T if t == 1 delta_P = Results.P_GT(:,1) - P_GT_Min; else delta_P = Results.P_GT(:,t) - Results.P_GT(:,t-1); end % 检查爬坡约束 violation = any(abs(delta_P) > Ramp_15min + 1e-4); if violation max_viol = max(abs(delta_P) - Ramp_15min); fprintf('时段 %d: 最大爬坡违反 %.2f MW (限制: %.2f MW)\n',... t, max_viol, max(Ramp_15min)); end end %% ==================== 结果可视化 ==================== time_vec = (1:T)*dt; % 时间向量 (小时) % 1. 机组出力曲线 figure('Name', '机组出力调度', 'Position', [100, 100, 1200, 800]); subplot(4,1,1); plot(time_vec, Results.P_GT', 'LineWidth', 1.5); title('常规机组出力'); xlabel('时间 (小时)'); ylabel('出力 (MW)'); legend('机组1', '机组2', '机组3', '机组4', '机组5', '机组6'); grid on; xlim([0 24]); % 2. 风电出力 subplot(4,1,2); plot(time_vec, Results.P_wt, 'b-', time_vec, P_wt_Min_96, 'r--', 'LineWidth', 1.5); title('风电出力'); xlabel('时间 (小时)'); ylabel('出力 (MW)'); legend('实际出力', '预测出力', 'Location', 'best'); grid on; xlim([0 24]); % 3. 储能系统状态 subplot(4,1,3); plot(time_vec, Results.SOC_ESS(1:T), 'b-', time_vec, Results.SOC_HSS(1:T), 'r-', 'LineWidth', 1.5); title('储能系统状态'); xlabel('时间 (小时)'); ylabel('状态'); legend('电储能 (MWh)', '氢储能 (kg)', 'Location', 'best'); grid on; xlim([0 24]); % 4. 爬坡约束违反检测 subplot(4,1,4); bar(time_vec, Results.Ramp_Violation, 'r'); title('爬坡约束违反检测'); xlabel('时间 (小时)'); ylabel('违反 (0/1)'); ylim([0 1.5]); grid on; xlim([0 24]); % 5. 储能系统功率 figure('Name', '储能系统功率', 'Position', [100, 100, 1200, 600]); plot(time_vec, Results.P_cha, 'g', time_vec, Results.P_dis, 'm', ... time_vec, Results.P_FC, 'b', time_vec, Results.P_EL, 'r', 'LineWidth', 1.5); title('储能系统功率'); xlabel('时间 (小时)'); ylabel('功率 (MW)'); legend('电储能充电', '电储能放电', '燃料电池', '电解槽', 'Location', 'best'); grid on; xlim([0 24]); % 6. 优化成本分析 figure('Name', '优化成本分析', 'Position', [100, 100, 800, 400]); plot(1:N_windows, Results.Obj, 'bo-', 'LineWidth', 1.5); xlabel('滚动窗口编号'); ylabel('优化成本 (元)'); title('各窗口优化成本变化'); grid on; % 7. 机组爬坡率使用热力图 figure('Name', '机组爬坡率使用分析', 'Position', [100, 100, 1200, 600]); ramp_ratio = zeros(N_Gens, T-1); for t = 1:T-1 delta_P = abs(Results.P_GT(:,t+1) - Results.P_GT(:,t)); ramp_ratio(:,t) = delta_P ./ Ramp_15min; end imagesc(1:T-1, 1:N_Gens, ramp_ratio); colorbar; caxis([0 1.2]); % 设置颜色范围 title('机组爬坡率使用率 (|ΔP|/Ramp_{max})'); xlabel('时段过渡'); ylabel('机组编号'); % 标记超限区域 hold on; [gen_idx, time_idx] = find(ramp_ratio > 1); plot(time_idx, gen_idx, 'rx', 'MarkerSize', 10, 'LineWidth', 1.5); %% ==================== 结果分析 ==================== % 计算总成本 total_cost = sum(Results.Obj(isfinite(Results.Obj))); fprintf('\n====== 优化结果汇总 ======\n'); fprintf('总优化成本: %.2f 元\n', total_cost); % 计算弃风率 potential_wind = sum(P_wt_Min_96); actual_wind = sum(Results.P_wt); curtailment_rate = (potential_wind - actual_wind)/potential_wind * 100; fprintf('弃风率: %.2f%%\n', curtailment_rate); % 计算碳排放量 carbon_emission = sum(sum(dt * Carbon_Factor' .* Results.P_GT)); fprintf('总碳排放量: %.2f 吨CO2\n', carbon_emission); % 保存结果 save('optimization_results.mat', 'Results', 'P_load_96', 'P_wt_Min_96'); disp('优化结果已保存');
最新发布
11-28
每个窗口的优化模型与第一个代码类似，但时间长度变为窗口大小（16个时段），并且初始状态由上一个窗口的结果提供。 5. 保存当前窗口第一时段（或前R_step个时段）的结果，并更新系统状态。 由于滚动优化需要多...
请根据代码一的日前优化模型修改代码二的日内调度模型，代码一【%% 基于电-氢混合储能的电力系统优化调度模型（日前经济调度） clc clear %% 加载电力系统数据 mpc = loadcase(‘case30’); baseMVA = mpc.baseMVA; % 系统基准功率 % 电负荷和风电数据 P_load_Hourly = [297 258 220 208 226 258 280 330 372 364 400 468 636 600 506 554 392 392 412 403 426 390 322 268]; P_wt_Hourly = [50 55 54 58 60 100 102 125 130 152 153 140 157 130 127 300 275 260 210 200 169 150 130 100]; %% 参数设置 T = 24; % 时间周期数 N_Gens = 6; % 常规机组数量 % 常规机组参数 [机组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); % 爬坡率 Carbon_Factor = Gen_Params(:,8); % 碳排放系数 (tCO2/MWh) P_GT_Min = Gen_Params(:,3); % 最小出力 P_GT_Max = Gen_Params(:,2); % 最大出力 % 启动/停机成本参数 (单位: 元) C_Start_qt = [2768; 2874; 2136; 2677; 2262; 2869]; % 启动成本 C_Stop_qt = [1500; 1600; 1200; 1400; 1300; 1500]; % 停机成本 % 电储能参数 P_Cha_Max = 45; % 电储能最大充电功率 P_Dis_Max = 90; % 电储能最大放电功率 K_ess_Eff = 0.92; % 电储能转换效率 K_Sel_Dis = 0.98; % 自放电率(99.8%/小时) C_ess_Max = 300; % 储能容量限制 % 氢系统参数 P_FC_Max = 400; % 燃料电池最大功率 K_H2P_xs = 9.1; % 氢转电系数 (kg/MWh) P_EL_Max = 400; % 电解槽最大功率 K_P2H_xs = 0.22; % 电转氢系数 (kg/MWh) H_cha_Max = 120; % 最大充氢速率 (kg/h) H_dis_Max = 100; % 最大放氢速率 (kg/h) C_HSS_Max = 400; % 最大储氢容量 (kg) C_HSS_Min = 10; % 最小储氢容量 (kg) % 节点位置 Bus_GEN_e = [1 2 5 8 11 13]; % 火电机组节点 Bus_WT_e = 2; % 风电节点 Bus_ESS_e = 11; % 电储能节点 Bus_HSS_h = 8; % 氢储能节点 % 环境成本系数 C_carbon = 150; % 碳排放成本 (元/吨CO2) C_FC_env = -500; % 燃料电池环境补贴 (负成本) C_EL_env = -300; % 电解槽环境补贴 (负成本) % 峰谷时段定义 t_peak = [12:14, 19:21]; % 高峰时段 (中午+晚高峰) t_valley = [1:5, 24]; % 低谷时段 (凌晨+午夜) %% 决策变量定义 Bus_Theta = sdpvar(30, T, ‘full’); % 节点相角 P_GT_e = sdpvar(N_Gens, T, ‘full’); % 火电机组出力 B_U_GT = binvar(N_Gens, T, ‘full’); % 火电运行状态 B_Start_GT = binvar(N_Gens, T, ‘full’); % 火电启动 B_Stop_GT = binvar(N_Gens, T, ‘full’); % 火电停机 P_wt_e = sdpvar(1, T, ‘full’); % 风电出力 % 电储能变量 P_cha_e = sdpvar(1, T, ‘full’); % 充电功率 P_dis_e = sdpvar(1, T, ‘full’); % 放电功率 B_U_cha = binvar(1, T, ‘full’); % 充电状态 B_U_dis = binvar(1, T, ‘full’); % 放电状态 SOC_ESS = sdpvar(1, T+1, ‘full’); % 电储能SOC % 氢系统变量 P_FC_e = sdpvar(1, T, ‘full’); % 燃料电池出力 P_e_EL = sdpvar(1, T, ‘full’); % 电解槽用电 B_U_EL = binvar(1, T, ‘full’); % 电解槽状态 B_U_FC = binvar(1, T, ‘full’); % 燃料电池状态 P_H2_FC = sdpvar(1, T, ‘full’); % 燃料电池耗氢 P_EL_H2 = sdpvar(1, T, ‘full’); % 电解槽产氢 P_cha_H2 = sdpvar(1, T, ‘full’); % 储氢罐充氢 P_dis_H2 = sdpvar(1, T, ‘full’); % 储氢罐放氢 SOC_HSS = sdpvar(1, T+1, ‘full’); % 储氢罐SOC % 节点注入功率矩阵 P_net_e = sdpvar(30, T, ‘full’); %% 生成节点电负荷矩阵 P_Load_Base = sum(mpc.bus(:,3)); % 系统基准负荷 Ratio_Load = P_load_Hourly / P_Load_Base; P_load_e = mpc.bus(:,3) * Ratio_Load; % 节点负荷矩阵 (30×T, MW) %% 约束条件 Cons = []; % ========== 火电机组约束 ========== for t = 1:T % 机组出力限制 Cons = [Cons, P_GT_Min .* B_U_GT(:,t) <= P_GT_e(:,t) <= P_GT_Max .* B_U_GT(:,t)]; % 机组爬坡约束 if t > 1 Cons = [Cons, -Ramp_e <= P_GT_e(:,t) - P_GT_e(:,t-1) <= Ramp_e]; end % 启停状态逻辑 if t == 1 Cons = [Cons, B_Start_GT(:,t) == B_U_GT(:,t)]; Cons = [Cons, B_Stop_GT(:,t) == 0]; else Cons = [Cons, B_Start_GT(:,t) <= 1 - B_U_GT(:,t-1)]; Cons = [Cons, B_Start_GT(:,t) <= B_U_GT(:,t)]; Cons = [Cons, B_Start_GT(:,t) >= B_U_GT(:,t) - B_U_GT(:,t-1)]; Cons = [Cons, B_Stop_GT(:,t) <= B_U_GT(:,t-1)]; Cons = [Cons, B_Stop_GT(:,t) <= 1 - B_U_GT(:,t)]; Cons = [Cons, B_Stop_GT(:,t) >= B_U_GT(:,t-1) - B_U_GT(:,t)]; end end % ========== 风电机组约束 ========== Cons = [Cons, 0 <= P_wt_e <= P_wt_Hourly]; % ========== 电储能约束 ========== Cons = [Cons, B_U_cha + B_U_dis <= 1]; % 充放电互斥 Cons = [Cons, 0 <= P_cha_e <= P_Cha_Max .* B_U_cha]; Cons = [Cons, 0 <= P_dis_e <= P_Dis_Max .* B_U_dis]; Cons = [Cons, SOC_ESS(1) == 10]; % 初始SOC Cons = [Cons, SOC_ESS(T+1) == 10]; % 终止SOC for t = 1:T Cons = [Cons, SOC_ESS(t+1) == K_Sel_Dis * SOC_ESS(t) + K_ess_Eff * P_cha_e(t) - P_dis_e(t) / K_ess_Eff]; Cons = [Cons, 0 <= SOC_ESS(t+1) <= C_ess_Max]; end % ========== 氢系统约束 ========== Cons = [Cons, 0 <= P_e_EL <= P_EL_Max .* B_U_EL]; Cons = [Cons, P_EL_H2 == P_e_EL * K_P2H_xs]; Cons = [Cons, 0 <= P_FC_e <= P_FC_Max .* B_U_FC]; Cons = [Cons, P_FC_e == P_H2_FC / K_H2P_xs]; % 氢能平衡 Cons = [Cons, P_H2_FC + P_cha_H2 == P_EL_H2 + P_dis_H2]; % 互斥运行 Cons = [Cons, B_U_EL + B_U_FC <= 1]; % 充放氢速率限制 Cons = [Cons, 0 <= P_cha_H2 <= H_cha_Max .* B_U_EL]; Cons = [Cons, 0 <= P_dis_H2 <= H_dis_Max .* B_U_FC]; % 储氢容量约束 Cons = [Cons, SOC_HSS(1) == 20]; % 初始氢储量 Cons = [Cons, SOC_HSS(T+1) == 20]; % 终止氢储量 for t = 1:T Cons = [Cons, SOC_HSS(t+1) == SOC_HSS(t) + P_cha_H2(t) - P_dis_H2(t)]; Cons = [Cons, C_HSS_Min <= SOC_HSS(t+1) <= C_HSS_Max]; end % FC/EL爬坡约束 for t = 2:T Cons = [Cons, -30 <= P_FC_e(t)-P_FC_e(t-1) <= 30]; Cons = [Cons, -20 <= P_e_EL(t)-P_e_EL(t-1) <= 20]; end % ========== 节点功率平衡 ========== % 创建设备-节点映射矩阵 N = size(P_load_e, 1); G_map = sparse(Bus_GEN_e, 1:N_Gens, 1, N, N_Gens); % 常规机组 W_map = sparse(Bus_WT_e, 1, 1, N, 1); % 风电机组 E_map = sparse(Bus_ESS_e, 1, 1, N, 1); % 电储能 H_map = sparse(Bus_HSS_h, 1, 1, N, 1); % 氢系统 for t = 1:T total_injection = G_map * P_GT_e(:,t) + … W_map * P_wt_e(t) + … E_map * (P_dis_e(t) - P_cha_e(t)) + … H_map * (P_FC_e(t) - P_e_EL(t)); Cons = [Cons, P_net_e(:,t) == total_injection - P_load_e(:,t)]; end % ========== 直流潮流约束 ========== [Ybus, ~, ~] = makeYbus(mpc); B = imag(Ybus); % 电纳矩阵虚部 ref_node = 1; % 参考节点 non_ref = 2:30; % 非参考节点 Bred = B(non_ref, non_ref); % 降阶B矩阵 Cons = [Cons, Bus_Theta(ref_node, 😃 == 0]; % 参考节点相角 for t = 1:T Cons = [Cons, Bred * Bus_Theta(non_ref, t) == P_net_e(non_ref, t) / baseMVA]; end Cons = [Cons, -pi/2 <= Bus_Theta <= pi/2]; % 相角范围 % ========== 全局功率平衡 ========== Total_Gen_e = sum(P_GT_e, 1) + P_wt_e + P_dis_e - P_cha_e + P_FC_e - P_e_EL; Total_Load_e = sum(P_load_e, 1); Cons = [Cons, Total_Gen_e == Total_Load_e]; % 严格功率平衡 %% 目标函数 (使用循环计算) % 1. 燃料成本 F1_cost = 0; for i = 1:N_Gens for t = 1:T F1_cost = F1_cost + a2(i) * P_GT_e(i,t)^2 + a1(i) * P_GT_e(i,t) + a0(i) * B_U_GT(i,t); end end % 2. 启停成本 F2_cost = 0; for i = 1:N_Gens for t = 1:T F2_cost = F2_cost + C_Start_qt(i) * B_Start_GT(i,t) + C_Stop_qt(i) * B_Stop_GT(i,t); end end % 3. 环境成本 Carbon_Emission = 0; for i = 1:N_Gens for t = 1:T Carbon_Emission = Carbon_Emission + Carbon_Factor(i) * P_GT_e(i,t); end end F_carbon = C_carbon * Carbon_Emission; % 氢储能环境补贴 F_env = 0; for t = 1:T % 燃料电池补贴 if ismember(t, t_peak) F_env = F_env + 1.3 * C_FC_env * P_FC_e(t); else F_env = F_env + C_FC_env * P_FC_e(t); end % 电解槽补贴 if ismember(t, t_valley) F_env = F_env + 1.5 * C_EL_env * P_e_EL(t); else F_env = F_env + C_EL_env * P_e_EL(t); end end F3_cost = F_carbon + F_env; F_Obj_Total = F1_cost + F2_cost + F3_cost; %% 模型优化求解 Opt = sdpsettings(‘solver’, ‘cplex’, ‘verbose’, 1, ‘showprogress’, 1); Sol = optimize(Cons, F_Obj_Total, Opt); % 结果处理与分析 if Sol.problem == 0 disp(‘======= 优化成功 =======’); fprintf(‘>>>总成本: %.2f 元\n’, value(F_Obj_Total)); fprintf(‘-燃料成本: %.2f 元\n’, value(F1_cost)); fprintf(‘-启停成本: %.2f 元\n’, value(F2_cost)); fprintf(‘-环境成本: %.2f 元\n’, value(F_carbon)); % 提取结果 P_GT_e = value(P_GT_e); P_wt_e = value(P_wt_e); P_cha_e = value(P_cha_e); P_dis_e = value(P_dis_e); P_FC_e = value(P_FC_e); P_e_EL = value(P_e_EL); SOC_ESS = value(SOC_ESS); SOC_HSS = value(SOC_HSS); % 输出氢储能使用情况 disp('==== 氢储能使用情况 ===='); fprintf('电解槽总用电量: %.2f MWh\n', sum(P_e_EL)); fprintf('燃料电池总发电量: %.2f MWh\n', sum(P_FC_e)); else error(‘>>>求解失败: %s’, Sol.info); end %% 结果可视化 figure(1) components = [sum(P_GT_e, 1)‘, … % 常规机组 P_wt_e’, … % 风电机组 P_dis_e’, … % 储能放电 P_FC_e’, … % 燃料电池 -P_cha_e’, … % 储能充电（负值） - P_e_EL’]; % 电解槽机（负值） he1 = bar(components, ‘stacked’); hold on; total_load = sum(P_load_e, 1); plot(1:T, total_load, ‘g-*’); % 设置颜色和透明度 colors = lines(6); for i = 1:6 he1(i).FaceColor = colors(i,:); he1(i).FaceAlpha = 0.9; end legend(‘常规机组’, ‘风电机组’, ‘储能放电’, ‘燃料电池’, ‘储能充电’, ‘电解槽机’, ‘总电负荷’); grid on; xlabel(‘时间 (小时)’); ylabel(‘功率 (MW)’); title(‘电力系统功率平衡’); figure(2) Plot_Ele_Net = zeros(24, 6); for t = 1:24 Plot_Ele_Net(t, 1) = value(P_GT_e(1, t)); Plot_Ele_Net(t, 2) = value(P_GT_e(2, t)); Plot_Ele_Net(t, 3) = value(P_GT_e(3, t)); Plot_Ele_Net(t, 4) = value(P_GT_e(4, t)); Plot_Ele_Net(t, 5) = value(P_GT_e(5, t)); Plot_Ele_Net(t, 6) = value(P_GT_e(6, t)); end % 绘制堆叠条形图 bar(Plot_Ele_Net, ‘stacked’); hold on % 添加负荷曲线 plot(1:24, sum(P_GT_e,1), ‘g-*’); grid on legend(‘机组1’,‘机组2’,‘机组3’,‘机组4’,‘机组5’,‘机组6’,‘常规机组’); xlabel(‘时间/h’); ylabel(‘功率/MW’); title(‘常规机组出力分布’); figure(3) commitment = double(B_U_GT); % 逻辑矩阵转数值 [num_units, T] = size(commitment); % 创建增强型颜色映射 cmap = [0.95 0.95 0.95; % 停机 - 浅灰色 0.25 0.65 0.35]; % 运行 - 绿色 % 创建热图 imagesc(commitment); colormap(cmap); % 设置坐标轴 set(gca, ‘YTick’, 1:num_units, … ‘YTickLabel’, arrayfun(@(x) sprintf(‘机组 %d’, x), 1:num_units, ‘UniformOutput’, false), … ‘TickDir’, ‘out’, … ‘FontSize’, 10); % 仅在可读性好的情况下添加1/0 for i = 1:num_units for t = 1:T if num_units <= 30 && T <= 48 % 仅在矩阵尺寸可读时添加文本 if commitment(i,t) > 0.5 text(t, i, ‘1’, ‘HorizontalAlignment’, ‘center’, … ‘Color’, [0 0.3 0], ‘FontWeight’, ‘bold’, ‘FontSize’, 9); else text(t, i, ‘0’, ‘HorizontalAlignment’, ‘center’, … ‘Color’, [0.4 0.4 0.4], ‘FontSize’, 9); end end end end grid on; % 添加颜色条和图例 c = colorbar(‘Ticks’, [0.25, 0.75], ‘TickLabels’, {‘停机’, ‘运行’}); c.Label.String = ‘状态’; xlabel(‘时间/h’); ylabel(‘机组编号’); title(‘机组启停状态分析’); xlim([0 25]); figure(4); hourly_carbon = value(Carbon_Factor’ * P_GT_e); bar(hourly_carbon’, ‘FaceColor’, [0.25 0.65 0.35], ‘EdgeColor’, ‘none’); hold on; plot(1:24, value(P_FC_e), ‘b-'); plot(1:24, value(P_e_EL), 'r-’); xlabel(‘时间/h’); ylabel(‘碳排放 (吨CO₂) / 功率 (MW)’); title(‘碳排放量与氢能系统运行’); legend(‘机组碳排放量’, ‘燃料电池出力’, ‘电解槽器用电’, ‘Location’, ‘best’); grid on; xlim([0 25]); ylim([-50 500]); figure(5) yyaxis left; bar(P_cha_e’, ‘FaceColor’, [1 0 0]); hold on bar(-P_dis_e’, ‘FaceColor’, [0 0 1]); ylabel(‘充放电功率 (MW)’); ylim([-100 100]); yyaxis right; plot(0:T, SOC_ESS, ‘g-*’); ylabel(‘电储能容量 (MWh)’); title(‘电储能运行状态’); legend( ‘充电功率’, ‘放电功率’,‘电储能SOC’, ‘Location’, ‘best’); grid on; xlim([0 25]); figure(6) yyaxis left; bar(P_cha_H2’, ‘FaceColor’, [1 0 0]); hold on bar(-P_dis_H2’, ‘FaceColor’, [0 0 1]); ylabel(‘氢流量 (kg/h)’); ylim([-150 150]); yyaxis right; plot(0:T, SOC_HSS, ‘g-*’); ylabel(‘储氢量 (kg)’); title(‘氢储能状态’); legend( ‘储能产氢’, ‘储能耗氢’,‘储氢量SOC’, ‘Location’, ‘best’); grid on; xlim([0 25]); figure(7) Plot_H2_Net = zeros(24, 4); for t = 1:24 Plot_H2_Net(t, 1) = value(P_EL_H2(1, t)); Plot_H2_Net(t, 2) = -value(P_H2_FC(1, t)); Plot_H2_Net(t, 3) = -value(P_cha_H2(1, t)); Plot_H2_Net(t, 4) = value(P_dis_H2(1, t)); end bar(Plot_H2_Net, ‘stacked’); hold on grid on legend(‘电解槽器产氢’, ‘燃料电池耗氢’, ‘储氢罐充氢’, ‘储氢罐放氢’); xlabel(‘时间/h’); ylabel(‘储氢量 (kg)’); title(‘供需氢功率平衡状态’); xlim([0 25]); ylim([-150 200]); 】 ；代码二【%% 基于电-氢混合储能的电力系统优化调度模型（日内滚动调度） clc clear close all %% ==================== 参数设置 ==================== T = 96; % 总时段数 (15分钟间隔) W_size = 16; % 滚动窗口大小 (4小时) R_step = 4; % 滚动步长 (1小时) N_windows = floor((T - W_size)/R_step) + 1; % 滚动次数 N_Gens = 6; % 机组数量 dt = 0.25; % 时间步长 (小时) %% ==================== 数据加载 ==================== % 电负荷数据 (小时级) P_load_Hourly = 0.25*[297 258 220 208 226 258 280 330 372 364 400 468 636 600 506 554 392 392 412 403 426 390 322 268]; % 风电数据 (小时级) P_wt_Hourly = [50 55 54 58 60 100 102 125 130 152 153 140 157 130 127 300 275 260 210 200 169 150 130 100]; % 转换为15分钟数据 (96时段) P_load_Min_96 = repelem(P_load_Hourly, 1, 4); % 系统总负荷 P_wt_Min_96 = repelem(P_wt_Hourly, 1, 4); % 风电出力 % 加载电力系统案例 mpc = loadcase(‘case30’); 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/小时) Ramp_15min = Ramp_e * dt; % 15分钟爬坡率 (MW/时段) Carbon_Factor = Gen_Params(:,8); % 碳排放系数 (tCO2/MWh) P_GT_Min = Gen_Params(:,3); % 最小出力 P_GT_Max = Gen_Params(:,2); % 最大出力 % 启动/停机成本参数 (单位: 元) C_Start_qt = [2768; 2874; 2136; 2677; 2262; 2869]; % 启动成本 C_Stop_qt = [1500; 1600; 1200; 1400; 1300; 1500]; % 停机成本 % 电储能参数 P_Cha_Max = 45; % 最大充电功率 (MW) P_Dis_Max = 90; % 最大放电功率 (MW) K_ess_Eff = 0.92; % 充放电效率 C_ess_Max = 300; % 最大容量 (MWh) SOC_ESS_Init = 10; % 初始SOC (MWh) % 氢系统参数 P_FC_Max = 400; % 燃料电池最大功率 (MW) K_H2P_xs = 9.1; % 氢转电系数 (kg/MWh) P_EL_Max = 400; % 电解槽最大功率 (MW) K_P2H_xs = 0.22; % 电转氢系数 (kg/MWh) H_cha_Max = 100; % 最大充氢速率 (kg/h) H_dis_Max = 100; % 最大放氢速率 (kg/h) C_HSS_Max = 400; % 最大储氢容量 (kg) C_HSS_Min = 10; % 最小储氢容量 (kg) SOC_HSS_Init = 20; % 初始氢储量 (kg) % 环境成本系数 C_carbon = 150; % 碳排放成本 (元/吨CO2) % 节点位置 Bus_GEN_e = [1 2 5 8 11 13]; % 火电机组节点 Bus_WT_e = 2; % 风电节点 Bus_ESS_e = 11; % 电储能节点 Bus_HSS_h = 8; % 氢储能节点 %% ==================== 节点电负荷矩阵生成 ==================== % 电力系统基准负荷 (MW) P_Load_Base = sum(mpc.bus(:,3)); % 创建节点负荷矩阵 (30×96) P_load_96 = zeros(size(mpc.bus, 1), T); for i = 1:size(mpc.bus, 1) % 计算节点负荷比例 node_ratio = mpc.bus(i,3) / P_Load_Base; % 生成节点负荷曲线 P_load_96(i,:) = node_ratio * P_load_Min_96; end %% ==================== 电网络拓扑初始化 ==================== [Ybus, ~, ~] = makeYbus(mpc); B = imag(Ybus); % 电纳矩阵 ref_node = 1; % 参考节点 non_ref = 2:size(mpc.bus,1); % 非参考节点 % 创建设备-节点映射矩阵 N = size(mpc.bus, 1); % 节点数量 G_map = sparse(Bus_GEN_e, 1:N_Gens, 1, N, N_Gens); % 常规机组 W_map = sparse(Bus_WT_e, 1, 1, N, 1); % 风电机组 E_map = sparse(Bus_ESS_e, 1, 1, N, 1); % 电储能 H_map = sparse(Bus_HSS_h, 1, 1, N, 1); % 氢系统 %% ==================== 状态初始化 ==================== % 结果存储 Results = struct(… ‘P_GT’, zeros(N_Gens, T), … % 机组出力 ‘P_wt’, zeros(1, T), … % 风电出力 ‘P_cha’, zeros(1, T), … % 充电功率 ‘P_dis’, zeros(1, T), … % 放电功率 ‘SOC_ESS’, SOC_ESS_Initones(1, T+1), … % 电储能状态 ‘P_FC’, zeros(1, T), … % 燃料电池出力 ‘P_EL’, zeros(1, T), … % 电解槽用电 ‘SOC_HSS’, SOC_HSS_Initones(1, T+1), … % 氢储能状态 ‘Obj’, zeros(1, N_windows), … % 各窗口目标值 ‘Ramp_Violation’, zeros(1, T)… % 爬坡约束违反检测 ); % 机组状态跟踪 B_ug_96 = zeros(N_Gens, T); % 机组状态历史 B_ug_96(:,1) = ones(N_Gens,1); % 初始全运行 %% ==================== 主循环滚动优化 ==================== % 生成随机扰动种子 (确保可重复性) rng(2023, ‘twister’); for idx = 1:N_windows % 窗口时间范围 k_start = (idx-1)*R_step + 1; k_end = min(k_start + W_size - 1, T); win_len = k_end - k_start + 1; fprintf('\n====== 滚动窗口 %d/%d [时段 %d-%d] ======\n', idx, N_windows, k_start, k_end); % 状态初始化 if k_start == 1 SOC_ESS_init = SOC_ESS_Init; SOC_HSS_init = SOC_HSS_Init; P_GT_prev = P_GT_Min(:); % 初始状态为最小出力 prev_status = ones(N_Gens,1); % 初始状态全运行 else % 使用上一个窗口结束时的状态 SOC_ESS_init = Results.SOC_ESS(k_start); SOC_HSS_init = Results.SOC_HSS(k_start); P_GT_prev = Results.P_GT(:, k_start); % 取k_start时刻状态 prev_status = B_ug_96(:, k_start); % 取k_start时刻状态 end % 定义决策变量 (当前窗口) P_GT_e = sdpvar(N_Gens, win_len, 'full'); % 机组出力 B_U_G = binvar(N_Gens, win_len, 'full'); % 机组运行状态 B_U_Start = binvar(N_Gens, win_len, 'full');% 启动状态 B_U_Stop = binvar(N_Gens, win_len, 'full'); % 停机状态 P_wt_e = sdpvar(1, win_len); % 风电出力 P_cha_e = sdpvar(1, win_len); % 充电功率 P_dis_e = sdpvar(1, win_len); % 放电功率 SOC_ESS = sdpvar(1, win_len+1); % 电储能状态 P_FC_e = sdpvar(1, win_len); % 燃料电池出力 P_e_EL = sdpvar(1, win_len); % 电解槽用电 SOC_HSS = sdpvar(1, win_len+1); % 氢储能状态 Bus_Theta = sdpvar(N, win_len, 'full'); % 节点相角 % 二进制变量 (充放电状态) b_in_cha = binvar(1, win_len); b_in_dis = binvar(1, win_len); % 约束条件 Cons = []; %% ======== 初始状态约束 ======== Cons = [Cons, SOC_ESS(1) == SOC_ESS_init]; Cons = [Cons, SOC_HSS(1) == SOC_HSS_init]; %% ======== 机组约束 ======== for t = 1:win_len % 出力上下限 Cons = [Cons, B_U_G(:,t).*P_GT_Min <= P_GT_e(:,t) <= B_U_G(:,t).*P_GT_Max]; % 爬坡约束 (动态参考状态机制) if t == 1 ref_state = P_GT_prev; % 窗口起始状态 else ref_state = P_GT_e(:,t-1); % 前一时段状态 end % 爬坡约束计算 delta_P = P_GT_e(:,t) - ref_state; Cons = [Cons, -Ramp_15min <= delta_P <= Ramp_15min]; % 启停逻辑 if t == 1 Cons = [Cons, B_U_Start(:,1) >= B_U_G(:,1) - prev_status]; Cons = [Cons, B_U_Start(:,1) <= B_U_G(:,1)]; Cons = [Cons, B_U_Start(:,1) <= 1 - prev_status]; Cons = [Cons, B_U_Stop(:,1) >= prev_status - B_U_G(:,1)]; Cons = [Cons, B_U_Stop(:,1) <= prev_status]; Cons = [Cons, B_U_Stop(:,1) <= 1 - B_U_G(:,1)]; else Cons = [Cons, B_U_Start(:,t) >= B_U_G(:,t) - B_U_G(:,t-1)]; Cons = [Cons, B_U_Start(:,t) <= B_U_G(:,t)]; Cons = [Cons, B_U_Start(:,t) <= 1 - B_U_G(:,t-1)]; Cons = [Cons, B_U_Stop(:,t) >= B_U_G(:,t-1) - B_U_G(:,t)]; Cons = [Cons, B_U_Stop(:,t) <= B_U_G(:,t-1)]; Cons = [Cons, B_U_Stop(:,t) <= 1 - B_U_G(:,t)]; end end %% ======== 风电约束 ======== Cons = [Cons, 0 <= P_wt_e <= P_wt_Min_96(k_start:k_end)]; %% ======== 电储能约束 ======== Cons = [Cons, b_in_cha + b_in_dis <= 1]; % 充放电互斥 for t = 1:win_len % SOC更新 (加入时间步长0.25小时) Cons = [Cons, SOC_ESS(t+1) == SOC_ESS(t) + dt*(K_ess_Eff*P_cha_e(t) - P_dis_e(t)/K_ess_Eff)]; % 容量约束 Cons = [Cons, 0 <= SOC_ESS(t+1) <= C_ess_Max]; % 功率约束 Cons = [Cons, 0 <= P_cha_e(t) <= b_in_cha(t)*P_Cha_Max]; Cons = [Cons, 0 <= P_dis_e(t) <= b_in_dis(t)*P_Dis_Max]; end %% ======== 氢系统约束 ======== for t = 1:win_len % 产氢量和耗氢量计算 H_prod = P_e_EL(t) * dt * K_P2H_xs; % 产氢量 (kg) H_cons = P_FC_e(t) * dt * K_H2P_xs; % 耗氢量 (kg) % SOC更新 Cons = [Cons, SOC_HSS(t+1) == SOC_HSS(t) + H_prod - H_cons]; % 容量约束 Cons = [Cons, C_HSS_Min <= SOC_HSS(t+1) <= C_HSS_Max]; % 功率约束 Cons = [Cons, 0 <= P_FC_e(t) <= P_FC_Max]; Cons = [Cons, 0 <= P_e_EL(t) <= P_EL_Max]; Cons = [Cons, 0 <= H_prod <= H_cha_Max*dt]; Cons = [Cons, 0 <= H_cons <= H_dis_Max*dt]; end %% ======== 网络约束 ======== % 线路约束参数 branch = mpc.branch; rateA = branch(:, 6); % 线路容量限制 (MW) rateA(rateA == 0) = max(rateA(rateA>0))*1.5; % 处理零值容量限制 for t = 1:win_len t_global = k_start + t - 1; % 节点注入功率 P_inj = G_map * P_GT_e(:, t) + ... % 常规机组 W_map * P_wt_e(t) + ... % 风电 E_map * (P_dis_e(t) - P_cha_e(t)) + ... % 电储能 H_map * (P_FC_e(t) - P_e_EL(t));% 氢系统 % 节点净功率 P_net = P_inj - P_load_96(:, t_global); % 直流潮流约束 Cons = [Cons, B(non_ref, :) * Bus_Theta(:, t) == P_net(non_ref)/baseMVA]; % 线路潮流约束 for br = 1:size(branch, 1) from = branch(br, 1); to = branch(br, 2); x = branch(br, 4); % 线路电抗 P_flow = (Bus_Theta(from,t) - Bus_Theta(to,t)) / x * baseMVA; Cons = [Cons, -rateA(br) <= P_flow <= rateA(br)]; end end Cons = [Cons, Bus_Theta(ref_node, :) == 0]; % 参考节点 %% ======== 目标函数 ======== Obj = 0; for t = 1:win_len % 燃料成本 for g = 1:N_Gens fuel_cost = a2(g)*P_GT_e(g,t)^2 + a1(g)*P_GT_e(g,t) + a0(g); Obj = Obj + dt * fuel_cost; end % 启停成本 for g = 1:N_Gens Obj = Obj + C_Start_qt(g)*B_U_Start(g,t) + C_Stop_qt(g)*B_U_Stop(g,t); end % 碳排放成本 for g = 1:N_Gens carbon_cost = C_carbon * Carbon_Factor(g) * P_GT_e(g,t)/1000; % 元/吨 Obj = Obj + dt * carbon_cost; end end %% ======== 求解优化问题 ======== opts = sdpsettings('solver', 'cplex', ... 'verbose', 1, ... 'showprogress', 1, ... 'cplex.timelimit', 120); % 增加求解时限 sol = optimize(Cons, Obj, opts); %% ======== 结果处理 ======== if sol.problem == 0 fprintf('优化成功 | 窗口成本: %.2f 元\n', value(Obj)); Results.Obj(idx) = value(Obj); % 计算实际可保存长度 valid_len = min(R_step, T - k_start + 1); idx_range = k_start:k_start+valid_len-1; % 保存当前窗口结果 Results.P_GT(:, idx_range) = value(P_GT_e(:,1:valid_len)); B_ug_96(:, idx_range) = value(B_U_G(:,1:valid_len)); Results.P_wt(idx_range) = value(P_wt_e(1:valid_len)); Results.P_cha(idx_range) = value(P_cha_e(1:valid_len)); Results.P_dis(idx_range) = value(P_dis_e(1:valid_len)); Results.P_FC(idx_range) = value(P_FC_e(1:valid_len)); Results.P_EL(idx_range) = value(P_e_EL(1:valid_len)); % 更新状态量 Results.SOC_ESS(k_start+1:k_end+1) = value(SOC_ESS(2:end)); Results.SOC_HSS(k_start+1:k_end+1) = value(SOC_HSS(2:end)); % 爬坡约束验证 for t_idx = 1:valid_len t_global = k_start + t_idx - 1; if t_global == 1 delta_P = Results.P_GT(:,1) - P_GT_Min; else delta_P = Results.P_GT(:,t_global) - Results.P_GT(:,t_global-1); end % 检测爬坡约束违反 ramp_violation = max(abs(delta_P) - Ramp_15min); Results.Ramp_Violation(t_global) = any(ramp_violation > 1e-4); end else fprintf('优化失败: %s\n', sol.info); valid_len = min(R_step, T - k_start + 1); idx_range = k_start:k_start+valid_len-1; % 失败处理（使用上一个窗口结束状态） if k_start == 1 fill_value = repmat(P_GT_Min(:), 1, valid_len); wt_fill = P_wt_Min_96(idx_range); soc_ess_fill = SOC_ESS_Init * ones(1, valid_len); soc_hss_fill = SOC_HSS_Init * ones(1, valid_len); else fill_value = repmat(Results.P_GT(:, k_start), 1, valid_len); wt_fill = Results.P_wt(k_start) * ones(1, valid_len); soc_ess_fill = Results.SOC_ESS(k_start) * ones(1, valid_len); soc_hss_fill = Results.SOC_HSS(k_start) * ones(1, valid_len); end Results.P_GT(:, idx_range) = fill_value; Results.P_wt(idx_range) = wt_fill; Results.P_cha(idx_range) = 0; Results.P_dis(idx_range) = 0; Results.P_FC(idx_range) = 0; Results.P_EL(idx_range) = 0; Results.SOC_ESS(idx_range+1) = soc_ess_fill; Results.SOC_HSS(idx_range+1) = soc_hss_fill; Results.Obj(idx) = Inf; end end %% ==================== 状态转移验证 ==================== % 精确验证爬坡约束 for t = 1:T if t == 1 delta_P = Results.P_GT(:,1) - P_GT_Min; else delta_P = Results.P_GT(:,t) - Results.P_GT(:,t-1); end % 检查爬坡约束 violation = any(abs(delta_P) > Ramp_15min + 1e-4); if violation max_viol = max(abs(delta_P) - Ramp_15min); fprintf('时段 %d: 最大爬坡违反 %.2f MW (限制: %.2f MW)\n',... t, max_viol, max(Ramp_15min)); end end %% ==================== 结果可视化 ==================== time_vec = (1:T)*dt; % 时间向量 (小时) % 1. 机组出力曲线 figure(‘Name’, ‘机组出力调度’, ‘Position’, [100, 100, 1200, 800]); subplot(4,1,1); plot(time_vec, Results.P_GT’, ‘LineWidth’, 1.5); title(‘常规机组出力’); xlabel(‘时间 (小时)’); ylabel(‘出力 (MW)’); legend(‘机组1’, ‘机组2’, ‘机组3’, ‘机组4’, ‘机组5’, ‘机组6’); grid on; xlim([0 24]); % 2. 风电出力 subplot(4,1,2); plot(time_vec, Results.P_wt, ‘b-’, time_vec, P_wt_Min_96, ‘r–’, ‘LineWidth’, 1.5); title(‘风电出力’); xlabel(‘时间 (小时)’); ylabel(‘出力 (MW)’); legend(‘实际出力’, ‘预测出力’, ‘Location’, ‘best’); grid on; xlim([0 24]); % 3. 储能系统状态 subplot(4,1,3); plot(time_vec, Results.SOC_ESS(1:T), ‘b-’, time_vec, Results.SOC_HSS(1:T), ‘r-’, ‘LineWidth’, 1.5); title(‘储能系统状态’); xlabel(‘时间 (小时)’); ylabel(‘状态’); legend(‘电储能 (MWh)’, ‘氢储能 (kg)’, ‘Location’, ‘best’); grid on; xlim([0 24]); % 4. 爬坡约束违反检测 subplot(4,1,4); bar(time_vec, Results.Ramp_Violation, ‘r’); title(‘爬坡约束违反检测’); xlabel(‘时间 (小时)’); ylabel(‘违反 (0/1)’); ylim([0 1.5]); grid on; xlim([0 24]); % 5. 储能系统功率 figure(‘Name’, ‘储能系统功率’, ‘Position’, [100, 100, 1200, 600]); plot(time_vec, Results.P_cha, ‘g’, time_vec, Results.P_dis, ‘m’, … time_vec, Results.P_FC, ‘b’, time_vec, Results.P_EL, ‘r’, ‘LineWidth’, 1.5); title(‘储能系统功率’); xlabel(‘时间 (小时)’); ylabel(‘功率 (MW)’); legend(‘电储能充电’, ‘电储能放电’, ‘燃料电池’, ‘电解槽’, ‘Location’, ‘best’); grid on; xlim([0 24]); % 6. 优化成本分析 figure(‘Name’, ‘优化成本分析’, ‘Position’, [100, 100, 800, 400]); plot(1:N_windows, Results.Obj, ‘bo-’, ‘LineWidth’, 1.5); xlabel(‘滚动窗口编号’); ylabel(‘优化成本 (元)’); title(‘各窗口优化成本变化’); grid on; % 7. 机组爬坡率使用热力图 figure(‘Name’, ‘机组爬坡率使用分析’, ‘Position’, [100, 100, 1200, 600]); ramp_ratio = zeros(N_Gens, T-1); for t = 1:T-1 delta_P = abs(Results.P_GT(:,t+1) - Results.P_GT(:,t)); ramp_ratio(:,t) = delta_P ./ Ramp_15min; end imagesc(1:T-1, 1:N_Gens, ramp_ratio); colorbar; caxis([0 1.2]); % 设置颜色范围 title(‘机组爬坡率使用率 (|ΔP|/Ramp_{max})’); xlabel(‘时段过渡’); ylabel(‘机组编号’); % 标记超限区域 hold on; [gen_idx, time_idx] = find(ramp_ratio > 1); plot(time_idx, gen_idx, ‘rx’, ‘MarkerSize’, 10, ‘LineWidth’, 1.5); %% ==================== 结果分析 ==================== % 计算总成本 total_cost = sum(Results.Obj(isfinite(Results.Obj))); fprintf(‘\n====== 优化结果汇总 ======\n’); fprintf(‘总优化成本: %.2f 元\n’, total_cost); % 计算弃风率 potential_wind = sum(P_wt_Min_96); actual_wind = sum(Results.P_wt); curtailment_rate = (potential_wind - actual_wind)/potential_wind * 100; fprintf(‘弃风率: %.2f%%\n’, curtailment_rate); % 计算碳排放量 carbon_emission = sum(sum(dt * Carbon_Factor’ .* Results.P_GT)); fprintf(‘总碳排放量: %.2f 吨CO2\n’, carbon_emission); % 保存结果 save(‘optimization_results.mat’, ‘Results’, ‘P_load_96’, ‘P_wt_Min_96’); disp(‘优化结果已保存’); clc clear %% ==================== 参数设置 ==================== T = 96; % 总时段数 (15分钟间隔) W_size = 16; % 滚动窗口大小 (4小时) R_step = 4; % 滚动步长 (1小时) N_windows = floor((T - W_size)/R_step) + 1; % 滚动次数 N_Gens = 6; % 机组数量 dt = 0.25; % 时间步长 (小时) %% ==================== 数据加载与预处理 ==================== load pload_24.mat p_load_1 load pwt_24.mat p_wt load ug_24.mat ugone % 扩展24小时数据到96小时 P_wt_96 = repelem(p_wt, 1, 4); % 风电预测 P_load_96 = zeros(30, T); % 节点负荷 B_ug_96 = zeros(6, T); % 机组状态 for i = 1:30 P_load_96(i, 😃 = repelem(p_load_1(i, 😃, 1, 4); end for i = 1:6 B_ug_96(i, 😃 = repelem(ugone(i, 😃, 1, 4); end %% ==================== 系统参数设置 ==================== mpc = loadcase(‘case30’); baseMVA = mpc.baseMVA; % 系统基准功率 % 常规机组参数 Gen_Params = [1 250 50 0.0375 20 372.5 72; 2 80 20 0.175 17.5 352.3 48; 3 50 15 0.625 10 316.5 30; 4 35 10 0.0834 32.5 329.2 21; 5 30 10 0.25 30 276.4 18; 6 40 12 0.25 30 232.2 24]; % 提取发电机参数 a2 = Gen_Params(:,4); % 二次项系数 a1 = Gen_Params(:,5); % 一次项系数 a0 = Gen_Params(:,6); % 常数项 Ramp = Gen_Params(:,7); % 爬坡率 (MW/15min) Pg_min = Gen_Params(:,3); % 最小出力 Pg_max = Gen_Params(:,2); % 最大出力 % 启动成本参数 StartupCost = [2768; 2874; 2136; 2677; 2262; 2869]; % 电储能参数 P_Cha_Max = 45; % 最大充电功率 (MW) P_Dis_Max = 90; % 最大放电功率 (MW) K_ess_Eff = 0.92; % 充放电效率 C_ess_Max = 300; % 最大容量 (MWh) SOC_ESS_Init = 10; % 初始SOC (MWh) % 氢系统参数 P_FC_Max = 400; % 燃料电池最大功率 (MW) K_H2P = 9.09; % 氢转电系数 (kg/MWh) n_FC = 0.4; % 燃料电池效率 P_EL_Max = 400; % 电解槽最大功率 (MW) K_P2H = 0.22; % 电转氢系数 (kg/MWh) n_EL = 0.72; % 电解槽效率 H_in_Max = 100; % 最大充氢速率 (kg/h) H_out_Max = 100; % 最大放氢速率 (kg/h) C_HSS_Max = 400; % 最大储氢容量 (kg) C_HSS_Min = 0; % 最小储氢容量 (kg) SOC_HSS_Init = 20; % 初始氢储量 (kg) % 成本系数 kc = 0.5; % 弃风惩罚因子 O_M_ESS = [0.56, 0.76]; % 电储能运维成本 [充电, 放电] O_M_H2 = [0.65, 0.83]; % 氢系统运维成本 [电解槽, 燃料电池] % 节点位置 Bus_GEN = [1 2 5 8 11 13]; % 火电机组节点 Bus_WT = 2; % 风电节点 Bus_ESS = 11; % 电储能节点 Bus_HSS = 8; % 氢储能节点 %% ==================== 网络拓扑初始化 ==================== [Ybus, ~, ~] = makeYbus(mpc); Gbus = real(Ybus); % 电导矩阵 Bbus = imag(Ybus); % 电纳矩阵 % 创建节点-设备映射 N = size(mpc.bus, 1); % 节点数量 Mbg = sparse(Bus_GEN, 1:N_Gens, 1, N, N_Gens); % 常规机组 Mbw = sparse(Bus_WT, 1, 1, N, 1); % 风电 Mbe = sparse(Bus_ESS, 1, 1, N, 1); % 电储能 Mbs = sparse(Bus_HSS, 1, 1, N, 1); % 氢系统 %% ==================== 状态初始化 ==================== % 结果存储 Results = struct(… ‘P_GT’, zeros(N_Gens, T), … % 机组出力 ‘P_wt’, zeros(1, T), … % 风电出力 ‘P_cha’, zeros(1, T), … % 充电功率 ‘P_dis’, zeros(1, T), … % 放电功率 ‘SOC_ESS’, zeros(1, T+1), … % 电储能状态 ‘P_FC’, zeros(1, T), … % 燃料电池出力 ‘P_EL’, zeros(1, T), … % 电解槽用电 ‘SOC_HSS’, zeros(1, T+1), … % 氢储能状态 ‘Cost’, zeros(1, N_windows) … % 各窗口成本 ); % 初始状态 Results.SOC_ESS(1) = SOC_ESS_Init; Results.SOC_HSS(1) = SOC_HSS_Init; %% ==================== 主循环滚动优化 ==================== for idx = 1:N_windows % 窗口时间范围 k_start = (idx-1)*R_step + 1; k_end = min(k_start + W_size - 1, T); win_len = k_end - k_start + 1; fprintf('\n====== 滚动窗口 %d/%d [时段 %d-%d] ======\n', idx, N_windows, k_start, k_end); % 状态初始化 SOC_ESS_init = Results.SOC_ESS(k_start); SOC_HSS_init = Results.SOC_HSS(k_start); % ========== 定义决策变量 ========== P_GT = sdpvar(N_Gens, win_len, 'full'); % 机组出力 P_wt = sdpvar(1, win_len); % 风电出力 P_cha = sdpvar(1, win_len); % 充电功率 P_dis = sdpvar(1, win_len); % 放电功率 P_FC = sdpvar(1, win_len); % 燃料电池出力 P_EL = sdpvar(1, win_len); % 电解槽用电 SOC_ESS = sdpvar(1, win_len+1); % 电储能状态 SOC_HSS = sdpvar(1, win_len+1); % 氢储能状态 theta = sdpvar(N, win_len, 'full'); % 节点相角 V = sdpvar(N, win_len, 'full'); % 节点电压幅值 % 二进制变量 b_cha = binvar(1, win_len); % 充电状态 b_dis = binvar(1, win_len); % 放电状态 z_FC = binvar(1, win_len); % 燃料电池状态 z_EL = binvar(1, win_len); % 电解槽状态 % ========== 约束条件 ========== constraints = []; % 初始状态约束 constraints = [constraints, ... SOC_ESS(1) == SOC_ESS_init, ... SOC_HSS(1) == SOC_HSS_init]; % 机组约束 for t = 1:win_len for i = 1:N_Gens % 出力上下限 constraints = [constraints, ... B_ug_96(i, k_start+t-1)*Pg_min(i) <= P_GT(i,t) <= B_ug_96(i, k_start+t-1)*Pg_max(i)]; % 爬坡约束 if t == 1 && k_start > 1 % 窗口间爬坡约束 constraints = [constraints, ... Results.P_GT(i, k_start-1) - P_GT(i,t) <= Ramp(i), ... P_GT(i,t) - Results.P_GT(i, k_start-1) <= Ramp(i)]; elseif t > 1 % 窗口内爬坡约束 constraints = [constraints, ... P_GT(i,t-1) - P_GT(i,t) <= Ramp(i), ... P_GT(i,t) - P_GT(i,t-1) <= Ramp(i)]; end end % 风电约束 constraints = [constraints, 0 <= P_wt(t) <= P_wt_96(k_start+t-1)]; % 储能约束 constraints = [constraints, ... b_cha(t) + b_dis(t) <= 1, ... % 充放电互斥 0 <= P_cha(t) <= b_cha(t)*P_Cha_Max, ... 0 <= P_dis(t) <= b_dis(t)*P_Dis_Max]; % 氢系统约束 constraints = [constraints, ... 0 <= P_FC(t) <= P_FC_Max * z_FC(t), ... 0 <= P_EL(t) <= P_EL_Max * z_EL(t), ... z_FC(t) + z_EL(t) <= 1]; % 运行状态互斥 end % 状态更新约束 for t = 1:win_len % 电储能状态更新 constraints = [constraints, ... SOC_ESS(t+1) == SOC_ESS(t) + K_ess_Eff*P_cha(t)*dt - P_dis(t)/K_ess_Eff*dt]; % 氢系统状态更新 H_prod = P_EL(t) * n_EL * K_P2H * dt; % 产氢量 H_cons = P_FC(t) / (n_FC * K_H2P) * dt; % 耗氢量 constraints = [constraints, ... SOC_HSS(t+1) == SOC_HSS(t) + H_prod - H_cons]; end % 容量约束 constraints = [constraints, ... 0 <= SOC_ESS(2:end) <= C_ess_Max, ... C_HSS_Min <= SOC_HSS(2:end) <= C_HSS_Max]; % 潮流约束 for t = 1:win_len % 节点功率平衡 P_inj = Mbg * P_GT(:,t) + ... % 常规机组 Mbw * P_wt(t) + ... % 风电 Mbe * (P_dis(t) - P_cha(t)) + ... % 电储能 Mbs * (P_FC(t) - P_EL(t)); % 氢系统 P_net = P_inj - P_load_96(:, k_start+t-1); % 交流潮流约束 P_calc = Gbus * V(:,t) * baseMVA - (Bbus + diag(sum(Bbus,2))) * theta(:,t) * baseMVA; constraints = [constraints, P_calc == P_net]; % 线路功率约束 for k = 1:size(mpc.branch, 1) i = mpc.branch(k, 1); j = mpc.branch(k, 2); P_flow = -Gbus(i,j)*(V(i,t) - V(j,t)) + Bbus(i,j)*(theta(i,t) - theta(j,t)); constraints = [constraints, -500 <= P_flow*baseMVA <= 500]; end % 电压和相角约束 constraints = [constraints, ... 0.95 <= V(:,t) <= 1.05, ... -pi <= theta(:,t) <= pi]; end % ========== 目标函数 ========== objective = 0; for t = 1:win_len % 机组成本 for i = 1:N_Gens objective = objective + a2(i)*(baseMVA*P_GT(i,t))^2 + ... a1(i)*(baseMVA*P_GT(i,t)) + a0(i); objective = objective + B_ug_96(i, k_start+t-1)*StartupCost(i); end % 储能运维成本 objective = objective + O_M_ESS(1)*P_cha(t) + O_M_ESS(2)*P_dis(t); % 氢系统运维成本 objective = objective + O_M_H2(1)*P_EL(t) + O_M_H2(2)*P_FC(t); % 弃风惩罚 objective = objective + kc*(P_wt_96(k_start+t-1) - P_wt(t)); end % ========== 求解优化问题 ========== options = sdpsettings('solver', 'cplex', 'verbose', 1); sol = optimize(constraints, objective, options); % ========== 结果处理 ========== if sol.problem == 0 fprintf('优化成功! 窗口成本: %.2f\n', value(objective)); Results.Cost(idx) = value(objective); % 保存当前窗口结果 valid_len = min(R_step, T - k_start + 1); idx_save = k_start:(k_start + valid_len - 1); Results.P_GT(:, idx_save) = value(P_GT(:, 1:valid_len)); Results.P_wt(idx_save) = value(P_wt(1:valid_len)); Results.P_cha(idx_save) = value(P_cha(1:valid_len)); Results.P_dis(idx_save) = value(P_dis(1:valid_len)); Results.P_FC(idx_save) = value(P_FC(1:valid_len)); Results.P_EL(idx_save) = value(P_EL(1:valid_len)); % 更新状态量 Results.SOC_ESS(k_start+1:k_start+win_len) = value(SOC_ESS(2:end)); Results.SOC_HSS(k_start+1:k_start+win_len) = value(SOC_HSS(2:end)); else fprintf('优化失败: %s\n', sol.info); % 使用前一时段值填充当前窗口 valid_len = min(R_step, T - k_start + 1); idx_save = k_start:(k_start + valid_len - 1); Results.P_GT(:, idx_save) = repmat(Results.P_GT(:, max(1,k_start-1)), 1, valid_len); Results.P_wt(idx_save) = P_wt_96(idx_save); Results.Cost(idx) = Inf; end end %% ==================== 结果可视化 ==================== time_vec = (1:T)*0.25; % 时间向量 (小时) % 1. 功率平衡图 figure(‘Position’, [100, 100, 1200, 600]); total_gen = sum(Results.P_GT, 1) + Results.P_wt + Results.P_dis - Results.P_cha + … Results.P_FC - Results.P_EL; total_load = sum(P_load_96, 1); plot(time_vec, total_gen, ‘b-’, ‘LineWidth’, 1.5); hold on; plot(time_vec, total_load, ‘r–’, ‘LineWidth’, 1.5); title(‘系统功率平衡’); xlabel(‘时间 (小时)’); ylabel(‘功率 (MW)’); legend(‘总发电量’, ‘总负荷’, ‘Location’, ‘best’); grid on; xlim([0 24]); xticks(0:2:24); % 2. 机组出力曲线 figure(‘Position’, [100, 100, 1200, 800]); subplot(3,1,1); plot(time_vec, Results.P_GT’, ‘LineWidth’, 1.5); title(‘常规机组出力’); xlabel(‘时间 (小时)’); ylabel(‘出力 (MW)’); legend(‘G1’,‘G2’,‘G3’,‘G4’,‘G5’,‘G6’); grid on; xlim([0 24]); % 3. 新能源与储能 subplot(3,1,2); plot(time_vec, Results.P_wt, ‘g-’, ‘LineWidth’, 1.5); hold on; plot(time_vec, P_wt_96, ‘g:’); title(‘风电出力’); xlabel(‘时间 (小时)’); ylabel(‘功率 (MW)’); legend(‘实际出力’, ‘预测出力’); grid on; xlim([0 24]); subplot(3,1,3); plot(time_vec, Results.P_cha, ‘b’, time_vec, Results.P_dis, ‘r’, … time_vec, Results.P_FC, ‘m’, time_vec, Results.P_EL, ‘c’, ‘LineWidth’, 1.5); title(‘储能系统功率’); xlabel(‘时间 (小时)’); ylabel(‘功率 (MW)’); legend(‘充电’, ‘放电’, ‘燃料电池’, ‘电解槽’); grid on; xlim([0 24]); % 4. 储能状态 figure(‘Position’, [100, 100, 1200, 400]); plot(time_vec, Results.SOC_ESS(1:T), ‘b-’, ‘LineWidth’, 1.5); hold on; plot(time_vec, Results.SOC_HSS(1:T), ‘r-’, ‘LineWidth’, 1.5); title(‘储能系统状态’); xlabel(‘时间 (小时)’); ylabel(‘状态’); legend(‘电储能 (MWh)’, ‘氢储能 (kg)’); grid on; xlim([0 24]); %% ==================== 结果分析 ==================== total_cost = sum(Results.Cost(Results.Cost < Inf)); fprintf(‘总优化成本: %.2f\n’, total_cost); % 保存结果 save(‘optimization_results.mat’, ‘Results’); disp(‘优化结果已保存’); 】
11-28
由于代码二是滚动调度，每个窗口内的时间是连续的，但整个调度周期是96个时段（15分钟间隔）。因此，我们需要将峰谷时段定义转换为15分钟分辨率。 步骤： 1. 在参数设置部分，定义峰谷时段（按15分钟间隔）。...
vscode-carbon_now_sh:用于在carbon.now.sh中打开当前编辑器内容的代码包
05-19
现在的碳 VS Code扩展，用于打开的当前编辑器内容。 只需突出显示要捕获的代码，然后运行Carbon扩展即可。 您将使用选择的代码填充textarea重定向到 。 运行扩展 有三种方法可以运行carbon-now-sh扩展。 Alt+Cmd+A或（ Alt+Windows+A上为Alt+Windows+A ）-快捷键，可立即启动扩展程序 打开命令面板（在Windows上为Cmd+Shift+P或Ctrl+Shift+P ），然后键入Carbon。 从编辑器的上下文菜单中。 替代碳选择 在VSCode设置中，您可以使用以下设置覆盖碳默认值： { " carbon.backgroundColor " : " rgba(0,0,0,0) " , " carbon.theme " : " seti " , " carbon.dropShadow " : true
截图工具（源代码）非常实用
04-19
截图工具（源代码） 截图工具（源代码） 截图工具（源代码）
carbon
03-12
录用测试 概括 创建一个应用程序，该应用程序从National Grid ESO的碳强度API加载并显示碳强度数据。 在开始之前，请先阅读完整的说明。 详细要求 细分为必填项-常规，必填项-功能，奖金和津贴 必填-一般 使用React。 我们已经为您提供了一个基本的应用样板来上手。 使用CSS（如果想花哨的话，可以使用预处理器）。 添加一些设计天赋。 从两个奖励类别中至少尝试一项（并解释您的选择）。 使用最佳实践进行源代码管理（提交粒度，提交消息策略等）。 必需-功能 主页（概述）页面 显示今天的总体碳强度。 按地区显示今天的碳强度（例如，表格）。 用户应能够按碳强度对区域进行排序。 单击一个区域应在侧面板中打开该区域的亮点。 侧面板应显示过去一周和过去一个月的区域碳强度。 侧面板应具有“转到区域”链接/按钮，可将用户带到区域详细信息页面。 区域（详细信息）页面 用户应该
Carbon
10-25
Carbon
carbon- 一个能生成漂亮的代码分享图的在线工具
little_spice的博客
05-15 4660
本篇转载自https://juejin.im/post/5a82c27f5188257a76632cc0 carbon是一个在线工具，它负责生成更加漂亮不同风格的代码图片，绝对比截图工具更加搞逼格。 先来一张效果图 是不是很帅气？ 接下来介绍下这个在线工具 地址为：https://carbon.now.sh 1 carbon支持不同风格不同语言，你可以选择你需要的，比如当我需要P...
Carbon 工具包使用指南
gitblog_00390的博客
09-11 705
Carbon 工具包使用指南 项目介绍 Carbon 是一个强大的命令行工具集合，专为开发者设计，旨在简化日常开发任务。它源自对效率的追求，提供了多种功能，包括但不限于代码片段分享、文件转换、以及多种编程语言的便捷操作。通过集成一系列实用工具，Carbon 让开发者能够在终端环境中更加高效地工作，提升编码体验。 项目快速启动 首先，确保你的系统中安装了 Node.js。然后，可以通过以下步骤来安装...
一键生成代码截图的神器，让你的代码美如画！
前端开发博客
05-18 1648
大家好，我是宝哥。今天要给大家安利几款超实用的代码截图工具，它们不仅功能强大，而且使用起来也非常便捷。Chalk.isthttps://chalk.ist/由@Idered打造，这是一款易于使用的工具，拥有漂亮的背景，能让你的代码片段脱颖而出。它包含了一系列现成的主题，并且允许你创建自己的主题。它是开源的，并且是用Vue 3构建的，所以速度飞快。想要真正体验它的功能，我建议你亲自试用一下。截图中使...
想要把代码片段转成图片？试试这几款在线代码截图神器
程序员成长指北
03-19 4434
大厂技术高级前端Node进阶点击上方程序员成长指北，关注公众号 回复1，加入高级Node交流群如何将代码片段转化为好看的图片，本文带你认识一些不错的在线代码截图工具。写在前面：文章中推荐的所有网站地址，帮大家全部整理好了，公众号后台回复【 代码片段】，即可领取在数字时代，代码不仅仅是程序员的语言，它也可以成为视觉艺术的一部分。如果你想要将你的代码片段转化为引人注目的图片，这里有几款在线工...
探索 ScreenshotBOF: 一款强大的代码截图工具
gitblog_00003的博客
03-05 430
探索 ScreenshotBOF: 一款强大的代码截图工具 在日常工作中，我们需要与他人分享代码片段或屏幕截图，而传统的截图工具无法满足我们的需求。因此，我们很荣幸向您推荐 ScreenshotBOF——一个强大、易用的代码截图工具。 什么是 ScreenshotBOF？ ScreenshotBOF 是一个开源的代码截图工具，可以帮助开发者快速生成包含代码片段的高质量屏幕截图。它支持多种编程语言...
【开发工具】Carbon使用全攻略，轻松打造吸睛代码展示
大雨的博客
08-12 852
摘要：Carbon是一款将源代码转化为精美图片的开源工具，帮助开发者在社交媒体、技术博客和演示中更优雅地展示代码。文章详细介绍了Carbon的安装步骤（需Node.js环境）、基础使用方法（包括界面介绍和代码图片生成流程），以及高级自定义技巧（主题搭配、字体调整和背景设置）。同时展示了Carbon在社交媒体分享、技术写作和演讲演示中的实际应用场景，并讲解了与VSCode和Markdown工具的集成方法。该工具通过提升代码的可视化效果，显著增强了技术内容的吸引力和专业度，成为开发者展示代码的利器。
评论
成就一亿技术人!
拼手气红包6.0元
还能输入1000个字符  | 博主筛选后可见
 
红包 添加红包
表情包 插入表情
表情包 代码片
 条评论被折叠 查看
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值