请你检查以下代码,指出修改建议:%% 基于共享储能的综合能源系统分布式优化调度
clc
clear
%% 参数初始化
maxIter = 100;
tolerant = 0.001;
X = zeros(maxIter+1, 48);
Power_buy = zeros(maxIter, 24);
Power_sell = zeros(maxIter, 24);
% 共享储能电价参数(24小时)
u_Db = [0.4,0.4,0.4,0.4,0.4,0.4,0.79,0.79,0.79,1.2,1.2,1.2,1.2,1.2,0.79,0.79,0.79,1.2,1.2,1.2,0.79,0.79,0.4,0.4];
u_Ds = [0.35,0.35,0.35,0.35,0.35,0.35,0.68,0.68,0.68,1.12,1.12,1.12,1.12,1.12,0.68,0.68,0.68,1.12,1.12,1.12,0.79,0.79,0.35,0.35];
% 电价上下限约束
u_Pbmax = 1.5 * u_Db;
u_Pbmin = u_Db;
u_Psmax = u_Ds;
u_Psmin = u_Psmax - 0.35*ones(1,24);
% 初始价格向量(购电+售电)
ub = [u_Pbmax, u_Psmax];
lb = [u_Pbmin, u_Psmin];
X(1, :) = (ub + lb) / 2; % 初始价格取上下界中值
% 二分法历史记录
X_up = repmat(ub, maxIter+1, 1);
X_low = repmat(lb, maxIter+1, 1);
% 目标函数存储
Obj_IES = zeros(1, maxIter);
Obj_ESO = zeros(1, maxIter);
%% 主迭代循环
converged = false;
final_iter = maxIter; % 默认最大迭代
for iter = 0:maxIter
fprintf('\n=== 迭代 %d ===\n', iter);
% === 特殊初始迭代 ===
if iter == 0
% 第一次IES优化
[Power_buy(1,:), Power_sell(1,:), Obj_IES(1)] = Fun_LA(X(1,1:48));
% 第一次ESO优化 (不增加新约束)
[X(2,:), Obj_ESO(1), X_up, X_low] = Fun_ESO(...
0, Power_buy(1,:), Power_sell(1,:), ...
X(1,:), zeros(1,48), X_up, X_low, iter);
continue;
end
% === 收敛检查 ===
if iter >= 1
% 计算相对变化 (使用欧几里得范数)
diff = X(iter+1, :) - X(iter, :);
relative_change = norm(diff) / norm(X(iter, :));
fprintf('相对变化: %.4f%%, 容差: %.4f%%\n', relative_change*100, tolerant*100);
if relative_change <= tolerant
fprintf('收敛于迭代 %d\n', iter);
converged = true;
final_iter = iter;
break;
end
end
% === 二分法处理 ===
U = 1; % 默认增加新约束
if iter >= 3 && all(abs(X(iter,:) - X(iter-2,:)) < 1e-5)
fprintf('应用二分法...\n');
% 生成新的价格点(取上一次迭代和当前迭代的中间值)
X(iter+1,:) = (X(iter,:) + X(iter-1,:)) / 2;
% 记录目标函数值(沿用上一次的值,因为这一步没有实际优化)
Obj_IES(iter) = Obj_IES(iter-1);
Obj_ESO(iter) = Obj_ESO(iter-1);
% 标记为不增加新约束,并跳过后续优化步骤
U = 0;
% 注意:这里我们直接进入下一次迭代,不进行IES和ESO优化
continue;
end
% === IES优化 ===
[Power_buy(iter,:), Power_sell(iter,:), Obj_IES(iter)] = Fun_LA(X(iter,1:48));
% === ESO优化 ===
[X(iter+1,:), Obj_ESO(iter), X_up, X_low] = Fun_ESO(...
U, Power_buy(iter,:), Power_sell(iter,:), ...
X(iter,:), X(iter-1,:), X_up, X_low, iter);
end
% 截断未使用的预分配空间
if final_iter < maxIter
Obj_IES = Obj_IES(1:final_iter);
Obj_ESO = Obj_ESO(1:final_iter);
Power_buy = Power_buy(1:final_iter, :);
Power_sell = Power_sell(1:final_iter, :);
else
final_iter = maxIter;
Obj_IES = Obj_IES(1:final_iter);
Obj_ESO = Obj_ESO(1:final_iter);
end
%% 结果可视化
% 收敛过程
figure('Name', '优化收敛过程', 'Position', [100, 100, 800, 600]);
subplot(2,1,1);
plot(1:length(Obj_IES), Obj_IES, 'r-o', 'LineWidth', 1.5, 'MarkerSize', 6);
xlabel('迭代次数');
ylabel('综合能源系统收益 (元)');
title('综合能源系统收益变化');
grid on;
subplot(2,1,2);
plot(1:length(Obj_ESO), Obj_ESO, 'b-s', 'LineWidth', 1.5, 'MarkerSize', 6);
xlabel('迭代次数');
ylabel('共享储能收益 (元)');
title('共享储能收益变化');
grid on;
% 最终交易电价
final_prices = X(final_iter+1, :);
buy_prices = final_prices(1:24);
sell_prices = final_prices(25:48);
figure('Name', '交易电价策略', 'Position', [200, 200, 900, 400]);
plot(1:24, buy_prices, 'b-*', 'LineWidth', 1.5, 'MarkerSize', 8);
hold on;
plot(1:24, sell_prices, 'r-o', 'LineWidth', 1.5, 'MarkerSize', 8);
plot(u_Db, 'g--', 'LineWidth', 1.2);
plot(u_Ds, 'm--', 'LineWidth', 1.2);
xlabel('时段');
ylabel('电价 (元/kWh)');
title('共享储能定价策略');
legend('购电价', '售电价', '电网购电价', '电网售电价', 'Location', 'best');
xlim([1, 24]);
grid on;
box on;
% 功率交互结果
figure('Name', '功率交互结果', 'Position', [300, 300, 1000, 400]);
subplot(1,2,1);
bar(Power_buy(final_iter,:), 'b');
xlabel('时段');
ylabel('购电功率 (kW)');
title('综合能源系统购电功率');
subplot(1,2,2);
bar(Power_sell(final_iter,:), 'r');
xlabel('时段');
ylabel('售电功率 (kW)');
title('综合能源系统售电功率');
%% 保存结果
results = struct();
results.iterations = final_iter;
results.ies_profit = Obj_IES;
results.eso_profit = Obj_ESO;
results.final_prices = final_prices;
results.power_buy = Power_buy(final_iter,:);
results.power_sell = Power_sell(final_iter,:);
save('optimization_results.mat', 'results');
fprintf('优化完成!总迭代次数: %d, IES最终收益: %.2f元, ESO最终收益: %.2f元\n', ...
final_iter, Obj_IES(end), Obj_ESO(end));
function [P_buy, P_sell, obj_value] = Fun_LA(price_vector)
% 综合能源系统优化函数
% 输入:price_vector - 电价向量 (1×48),前24小时为购电价,后24小时为售电价
% 输出:
% P_buy - 从共享储能购电功率 (1×24)
% P_sell - 向共享储能售电功率 (1×24)
% obj_value - IES目标函数值(收益)
%% 参数初始化
T = 24; % 24小时
hourly = 1:T;
% 电价参数
u_buy = price_vector(1:T); % 购电价 (元/kWh)
u_sell = price_vector(T+1:2*T); % 售电价 (元/kWh)
% 设备参数
P_gt_max = 500; % 燃气轮机最大输出功率 (kW)
P_gt_min = 50; % 燃气轮机最小输出功率 (kW)
eta_gt = 0.35; % 燃气轮机电效率
eta_whb = 0.6; % 余热锅炉效率
eta_gb = 0.85; % 燃气锅炉效率
COP_ec = 3.5; % 电制冷机制冷系数
COP_ac = 1.2; % 吸收式制冷机制冷系数
c_gas = 3.5; % 天然气价格 (元/m³)
LHV_gas = 9.7; % 天然气低热值 (kWh/m³)
c_om_gt = 0.05; % 燃气轮机运维成本 (元/kWh)
c_om_gb = 0.02; % 燃气锅炉运维成本 (元/kWh)
% 负荷数据 (kW)
% 实际应用中应从文件读取,这里使用示例数据
electric_load = [200 180 170 160 150 160 200 300 400 450 500 550 ...
600 580 550 500 450 500 600 650 600 500 400 300];
heat_load = [150 140 130 120 110 100 120 160 200 220 250 270 ...
300 280 250 220 200 220 250 280 250 220 180 150];
cool_load = [100 90 80 70 60 50 70 100 150 180 200 220 ...
250 230 200 180 150 180 200 230 200 180 150 120];
%% 决策变量定义
yalmip('clear');
% 主要设备出力
P_gt = sdpvar(1, T); % 燃气轮机发电功率
H_whb = sdpvar(1, T); % 余热锅炉制热功率
H_gb = sdpvar(1, T); % 燃气锅炉制热功率
C_ec = sdpvar(1, T); % 电制冷机制冷功率
C_ac = sdpvar(1, T); % 吸收式制冷机制冷功率
% 交互功率
P_buy = sdpvar(1, T); % 从共享储能购电功率
P_sell = sdpvar(1, T); % 向共享储能售电功率
% 引入二元变量处理互斥约束
b_buy = binvar(1, T); % 购电标志
b_sell = binvar(1, T); % 售电标志
M = 10000; % 大M值
%% 目标函数 - 最大化IES收益
% 收益 = 售电收入 - 购电成本 + 供热收入 - 气耗成本 - 运维成本
% 售电收入
income_sell = sum(u_sell .* P_sell);
% 购电成本
cost_buy = sum(u_buy .* P_buy);
% 气耗计算 (燃气轮机 + 燃气锅炉)
gas_consumption = (P_gt ./ (eta_gt * LHV_gas)) + (H_gb ./ (eta_gb * LHV_gas));
cost_gas = sum(gas_consumption) * c_gas;
% 运维成本
cost_om = sum(c_om_gt * P_gt) + sum(c_om_gb * H_gb);
% 供热收入 (假设供热价格为0.5元/kWh)
income_heat = 0.5 * sum(heat_load);
% 总目标函数
objective = income_sell - cost_buy + income_heat - cost_gas - cost_om;
%% 约束条件
constraints = [];
% 功率非负约束
constraints = [constraints, P_gt >= 0, H_whb >= 0, H_gb >= 0, C_ec >= 0, C_ac >= 0];
constraints = [constraints, P_buy >= 0, P_sell >= 0];
% 设备出力上下限
constraints = [constraints, P_gt_min <= P_gt <= P_gt_max];
% 能量平衡约束
for t = 1:T
% 电平衡: 燃气轮机 + 购电 = 电负荷 + 电制冷机 + 售电
constraints = [constraints, ...
P_gt(t) + P_buy(t) == electric_load(t) + C_ec(t)/COP_ec + P_sell(t)];
% 热平衡: 余热锅炉 + 燃气锅炉 = 热负荷
constraints = [constraints, ...
H_whb(t) + H_gb(t) == heat_load(t)];
% 冷平衡: 电制冷机 + 吸收式制冷机 = 冷负荷
constraints = [constraints, ...
C_ec(t) + C_ac(t) == cool_load(t)];
end
% 热电耦合约束 (燃气轮机余热回收)
for t = 1:T
constraints = [constraints, ...
H_whb(t) <= eta_whb * (P_gt(t)/eta_gt) * (1 - eta_gt)];
end
% 购售电互斥约束 (同一时刻不能同时购售电) - 使用大M法
constraints = [constraints, ...
P_buy <= M * b_buy, ...
P_sell <= M * b_sell, ...
b_buy + b_sell <= 1];
%% 求解优化问题
options = sdpsettings('solver', 'cplex', 'verbose', 0, 'showprogress', 0);
sol = optimize(constraints, -objective, options); % 最大化目标函数
%% 结果提取
if sol.problem == 0
P_buy = value(P_buy);
P_sell = value(P_sell);
obj_value = value(objective);
% 调试输出
fprintf('IES优化完成! 收益: %.2f元\n', obj_value);
fprintf('总购电量: %.2fkWh, 总售电量: %.2fkWh\n', sum(P_buy), sum(P_sell));
else
warning('IES优化未收敛! 问题: %s', sol.info);
P_buy = zeros(1, T);
P_sell = zeros(1, T);
obj_value = 0;
end
end
function [X_new, obj_value, X_up_out, X_low_out] = Fun_ESO(U, P_buy, P_sell, X_cur, X_pre, X_up_in, X_low_in, iter)
% 共享储能运营商优化函数
% 输入:
% U - 约束更新标志 (0不增加新约束, 1增加新约束)
% P_buy - IES购电功率 (1×24)
% P_sell - IES售电功率 (1×24)
% X_cur - 当前电价向量 (1×48)
% X_pre - 前一次电价向量 (1×48)
% X_up_in - 电价上界历史记录
% X_low_in- 电价下界历史记录
% iter - 当前迭代次数
% 输出:
% X_new - 新电价向量 (1×48)
% obj_value- ESO目标函数值(收益)
% X_up_out - 更新后的电价上界
% X_low_out- 更新后的电价下界
%% 参数初始化
T = 24; % 24小时
hourly = 1:T;
% 共享储能系统参数
ESS_capacity = 2000; % 储能容量 (kWh)
ESS_minSOC = 0.2; % 最小SOC
ESS_maxSOC = 0.9; % 最大SOC
ESS_eff_ch = 0.95; % 充电效率
ESS_eff_dis = 0.95; % 放电效率
ESS_Pmax_ch = 500; % 最大充电功率 (kW)
ESS_Pmax_dis = 500; % 最大放电功率 (kW)
SOC_init = 0.5; % 初始SOC
% 成本参数
c_om = 0.05; % 单位运维成本 (元/kWh)
c_loss = 0.02; % 单位损耗成本 (元/kWh)
% 电价边界约束
u_Db = [0.4,0.4,0.4,0.4,0.4,0.4,0.79,0.79,0.79,1.2,1.2,1.2,1.2,1.2,0.79,0.79,0.79,1.2,1.2,1.2,0.79,0.79,0.4,0.4];
u_Ds = [0.35,0.35,0.35,0.35,0.35,0.35,0.68,0.68,0.68,1.12,1.12,1.12,1.12,1.12,0.68,0.68,0.68,1.12,1.12,1.12,0.79,0.79,0.35,0.35];
u_Pbmax = 1.5 * u_Db; % 购电价上限
u_Pbmin = u_Db; % 购电价下限
u_Psmax = u_Ds; % 售电价上限
u_Psmin = u_Psmax - 0.35*ones(1,24); % 售电价下限
% 提取当前电价
u_buy_cur = X_cur(1:T); % 当前购电价
u_sell_cur = X_cur(T+1:2*T); % 当前售电价
% 初始化边界约束
X_up = X_up_in;
X_low = X_low_in;
%% 决策变量定义
yalmip('clear');
% 电价变量 (ESO的决策变量)
u_buy = sdpvar(1, T); % 购电价 (向IES购电的价格)
u_sell = sdpvar(1, T); % 售电价 (向IES售电的价格)
% 储能运行变量
P_ch = sdpvar(1, T); % 充电功率 (从IES购电)
P_dis = sdpvar(1, T); % 放电功率 (向IES售电)
SOC = sdpvar(1, T+1); % 储能状态 (0-1)
% 引入二元变量处理充放电互斥
b_ch = binvar(1, T); % 充电标志
b_dis = binvar(1, T); % 放电标志
M = 10000; % 大M值
%% 目标函数 - 最大化ESO收益
% 收益 = 售电收入 - 购电成本 - 运维成本 - 损耗成本
% 售电收入 (向IES售电)
income_sell = sum(u_sell .* P_dis);
% 购电成本 (从IES购电)
cost_buy = sum(u_buy .* P_ch);
% 运维成本 (按充放电量计算)
cost_om = c_om * (sum(P_ch) + sum(P_dis));
% 储能损耗成本 (按循环次数折算)
cost_loss = c_loss * (sum(P_ch) + sum(P_dis)) / 2; % 近似计算
% 总目标函数 (最大化)
objective = income_sell - cost_buy - cost_om - cost_loss;
%% 约束条件
constraints = [];
% ===== 电价约束 =====
% 1. 基本边界约束
constraints = [constraints, ...
u_Pbmin <= u_buy <= u_Pbmax, ...
u_Psmin <= u_sell <= u_Psmax];
% 2. 分布式优化约束 (当U=1时增加)
if U == 1 && iter > 0
% 价格单调性约束 (新电价应使IES的购售电行为不变)
for t = 1:T
if P_buy(t) > 0 % 当前为购电状态
constraints = [constraints, u_buy(t) <= u_buy_cur(t)];
elseif P_sell(t) > 0 % 当前为售电状态
constraints = [constraints, u_sell(t) >= u_sell_cur(t)];
else % 无交易状态
% 购电价不低于当前价,售电价不高于当前价
constraints = [constraints, ...
u_buy(t) >= u_buy_cur(t), ...
u_sell(t) <= u_sell_cur(t)];
end
end
% 更新历史边界 (用于二分法)
X_up(iter+1, :) = [u_buy_cur, u_sell_cur]; % 上界为当前价格
X_low(iter+1, :) = [u_buy_cur, u_sell_cur]; % 下界为当前价格
else
% 使用历史边界
constraints = [constraints, ...
X_low_in(iter+1, 1:T) <= u_buy <= X_up_in(iter+1, 1:T), ...
X_low_in(iter+1, T+1:2*T) <= u_sell <= X_up_in(iter+1, T+1:2*T)];
end
% ===== 储能运行约束 =====
% 1. 初始SOC
constraints = [constraints, SOC(1) == SOC_init];
% 2. SOC动态更新
for t = 1:T
constraints = [constraints, ...
SOC(t+1) == SOC(t) + (ESS_eff_ch * P_ch(t) - P_dis(t)/ESS_eff_dis) / ESS_capacity];
end
% 3. SOC边界约束
constraints = [constraints, ...
ESS_minSOC <= SOC <= ESS_maxSOC, ...
SOC(end) >= SOC_init]; % 终端SOC不低于初始值
% 4. 充放电功率约束
constraints = [constraints, ...
0 <= P_ch <= ESS_Pmax_ch, ...
0 <= P_dis <= ESS_Pmax_dis];
% 5. 充放电互斥约束 - 使用大M法
constraints = [constraints, ...
P_ch <= M * b_ch, ...
P_dis <= M * b_dis, ...
b_ch + b_dis <= 1];
% ===== 功率平衡约束 =====
% ESO与IES的功率交互必须匹配
constraints = [constraints, ...
P_ch == P_sell, ... % ESO充电 <=> IES售电
P_dis == P_buy]; % ESO放电 <=> IES购电
%% 求解优化问题
options = sdpsettings('solver', 'cplex', 'verbose', 0, 'showprogress', 0);
sol = optimize(constraints, -objective, options); % 最大化目标函数
%% 结果提取
if sol.problem == 0
X_new = [value(u_buy), value(u_sell)];
obj_value = value(objective);
% 输出边界更新
X_up_out = X_up;
X_low_out = X_low;
fprintf('ESO优化完成! 收益: %.2f元\n', obj_value);
else
warning('ESO优化未收敛! 问题: %s', sol.info);
X_new = X_cur; % 保持原价格
obj_value = 0;
X_up_out = X_up_in;
X_low_out = X_low_in;
end
end
最新发布