请检查以下修改代码：clc
clear
%% ======================= 参数定义 =======================
T = 24; % 调度时段数
Ns = 2; % 场景数量
ps = [0.3, 0.7]; % 场景概率
% 储能参数
Nerss = 1; % 储能电站总数
Nbt = 1; % 电化学储能数量
PIbt = 2.3; % 电化学储能单位时间折旧成本
Pech = -50; % 最大充电功率
Pedis = 50; % 最大放电功率
Se = 200; % 储能容量
Se0 = 100; % 初始储能容量
socmin = 0.1; % SOC下限
socmax = 0.9; % SOC上限
% 抽水蓄能参数
Pwmin = -50; % 最小功率
Pwmax = 50; % 最大功率
Vs0 = 100; % 初始蓄水量
Vpmin = 0; % 最小容量
Vpmax = 200; % 最大容量
% 常规机组
NG = 6; % 常规机组数量
NDG = 1; % DG数量
kc = 2.5; % 弃风惩罚系数
mp = [1 150 50 0.0375 20 372.5 72 2; % 机组参数矩阵
2 60 20 0.175 17.5 352.3 48 2;
3 60 15 0.625 10 316.5 30 2;
4 50 10 0.0834 32.5 329.2 30 2;
5 40 10 0.25 30 276.4 18 2;
6 45 12 0.25 30 232.2 24 2];
% 负荷与风电数据
pload = [390 382 354 340 328 409 460 460 472 511 465 458 442 456 457 471 475 503 512 469 445 438 410 362];
pw1 = [61 60.5 70 63 64 102 131 165 167 182 153 179 182 121 99 293 272 245 210 182 145 105 96 85];
pw2 = [92 150 125 162 148 132 141 96 43 19 72 96 49 132 141 180 165 192 148 196 143 81 21 45];
% 节点索引
idg = [1, 2, 5, 8, 11, 13]; % 常规机组节点
iw = 2; % 风电节点
ie = 2; % 储能节点
is = 8; % 抽水蓄能节点
% IDR成本系数
Ca = 100; Cb = 125; Cc = 150; Cd = 150;
%% ======================= 变量定义 =======================
% ------------------------ 场景1变量 ------------------------
thetaone = sdpvar(30, T, 'full'); % 节点电压相角
Pgone = sdpvar(NG, T, 'full'); % 常规机组出力
Pwone = sdpvar(1, T, 'full'); % 风电出力
upwone = binvar(1, T, 'full'); % 风电启停状态
Psone = sdpvar(1, T, 'full'); % 抽水蓄能出力
Perssone = sdpvar(1, T, 'full'); % 储能出力
uerssone = binvar(1, T, 'full'); % 储能充放电状态
Pdgone = sdpvar(1, T, 'full'); % DG出力
Pidrbone = sdpvar(30, T, 'full'); % IDR类型1
plossone = sdpvar(30, T, 'full'); % 网损
ugone = binvar(NG, T, 'full'); % 机组启停状态
Ppdrone = sdpvar(30, T, 'full'); % PDR
Pidraone = sdpvar(30, T, 'full'); % IDR类型2
% ------------------------ 场景2变量 ------------------------
thetatwo = sdpvar(30, T, 'full');
Pgtwo = sdpvar(NG, T, 'full');
Pwtwo = sdpvar(1, T, 'full');
upwtwo = binvar(1, T, 'full');
Pstwo = sdpvar(1, T, 'full');
Persstwo = sdpvar(1, T, 'full');
uersstwo = binvar(1, T, 'full');
Pdgtwo = sdpvar(1, T, 'full');
Pidrbtwo = sdpvar(30, T, 'full');
plosstwo = sdpvar(30, T, 'full');
ugtwo = binvar(NG, T, 'full');
Ppdrtwo = sdpvar(30, T, 'full');
Pidratwo = sdpvar(30, T, 'full');
%% ======================= 约束构建 =======================
con = [];
caseName = case30; % 加载case30数据
p30 = sum(caseName.bus(:,3)); % 节点总功率
pbl = pload ./ p30; % 负荷比例
pload1 = repmat(pbl, 30, 1) .* repmat(caseName.bus(:,3), 1, T); % 节点负荷
P_D = 0.85 .* pload1; % 固定负荷
% 1. 常规机组约束
con_gen = getConsGen1(Pgone, Pgtwo, mp, T, ugone, ugtwo);
con = [con, con_gen];
% 2. 风电约束
con = [con, 0 <= Pwone <= upwone .* pw1];
con = [con, 0 <= Pwtwo <= upwtwo .* pw2];
% 3. 潮流约束
baseMVA = 100;
bus = caseName.bus;
brch = caseName.branch;
f = brch(:,1);
t_brch = brch(:,2); % 为避免和变量t冲突，重命名为t_brch
x = brch(:,4);
nbus = size(bus,1);
nbrch = size(brch,1);
% 构建拓扑矩阵
Cft = sparse(nbrch, nbus);
for ii = 1:nbrch
Cft(ii, f(ii)) = 1;
Cft(ii, t_brch(ii)) = -1;
end
Bf = sparse(nbrch, nbus);
for ii = 1:nbrch
Bf(ii, f(ii)) = 1/x(ii);
Bf(ii, t_brch(ii)) = -1/x(ii);
end
Bbus = Cft' * Bf;
% 生成节点-设备关联矩阵
Mbg = getMbgMatrix(idg', bus);
Mbe = getMbdMatrix(ie, bus);
Mbw = getMbdMatrix(iw, bus);
Mbs = getMbdMatrix(is, bus);
% 场景1潮流平衡
Pbusone = [];
for ii = 1:T
Pbusone = [Pbusone, ...
Mbg * Pgone(:,ii) + Mbs * Psone(1,ii) + Mbe * Perssone(1,ii) + ...
Mbw * Pwone(1,ii) - pload1(:,ii) - Ppdrone(:,ii) - ...
Pidraone(:,ii) - Pidrbone(:,ii) - plossone(:,ii)];
end
for ii = 1:T
con = [con, Pbusone(:,ii) == Bbus * thetaone(:,ii) * baseMVA];
end
% 场景2潮流平衡
Pbustwo = [];
for ii = 1:T
Pbustwo = [Pbustwo, ...
Mbg * Pgtwo(:,ii) + Mbs * Pstwo(1,ii) + Mbe * Persstwo(1,ii) + ...
Mbw * Pwtwo(1,ii) - pload1(:,ii) - Ppdrtwo(:,ii) - ...
Pidratwo(:,ii) - Pidrbtwo(:,ii) - plosstwo(:,ii)];
end
for ii = 1:T
con = [con, Pbustwo(:,ii) == Bbus * thetatwo(:,ii) * baseMVA];
end
% 线路潮流限制
for j = 1:T
d_thetaone = thetaone(f,j) - thetaone(t_brch,j);
pfone = d_thetaone ./ x * baseMVA;
con = [con, -500 <= pfone <= 500];
d_thetatwo = thetatwo(f,j) - thetatwo(t_brch,j);
pftwo = d_thetatwo ./ x * baseMVA;
con = [con, -500 <= pftwo <= 500];
end
con = [con, -pi <= thetaone <= pi, -pi <= thetatwo <= pi];
% 4. DR约束
con = [con, -0.1.*pload1 <= Ppdrone <= 0.1.*pload1, sum(Ppdrone,2)==0];
con = [con, -0.1.*pload1 <= Ppdrtwo <= 0.1.*pload1, sum(Ppdrtwo,2)==0];
con = [con, -0.05.*pload1 <= Pidraone+Pidrbone <= 0.05.*pload1, ...
-0.05.*pload1 <= Pidrbone <= 0.05.*pload1, ...
-0.05.*pload1 <= Pidraone <= 0.05.*pload1, ...
sum(Pidraone+Pidrbone,2)==0];
con = [con, -0.05.*pload1 <= Pidratwo+Pidrbtwo <= 0.05.*pload1, ...
-0.05.*pload1 <= Pidratwo <= 0.05.*pload1, ...
-0.05.*pload1 <= Pidrbtwo <= 0.05.*pload1, ...
sum(Pidratwo+Pidrbtwo,2)==0];
% 5. 电储能约束
con = [con, Pech .* uerssone <= Perssone <= Pedis .* uerssone];
con = [con, Pech .* uersstwo <= Persstwo <= Pedis .* uersstwo];
for i = 1:T
soc1 = Se0 + sum(Perssone(:,1:i), 2);
soc2 = Se0 + sum(Persstwo(:,1:i), 2);
con = [con, Se*socmin <= soc1 <= Se*socmax];
con = [con, Se*socmin <= soc2 <= Se*socmax];
end
% 6. 抽水蓄能约束
con = [con, Pwmin <= Psone <= Pwmax];
con = [con, Pwmin <= Pstwo <= Pwmax];
for i = 1:T
vol1 = Vs0 + sum(Psone(:,1:i), 2);
vol2 = Vs0 + sum(Pstwo(:,1:i), 2);
con = [con, Vpmin <= vol1 <= Vpmax];
con = [con, Vpmin <= vol2 <= Vpmax];
end
%% ======================= 目标函数 =======================
% 常规机组成本分段线性化
tn = T; % 时段数
Pmax = mp(:,2)';
Pmin = mp(:,3)';
gn = 5; % 分段数
% 场景1线性化变量
x_pf = sdpvar(6, tn, 'full');
% 使用元胞数组存储6台机组的变量
gw = cell(1,6);
gz = cell(1,6);
for i = 1:6
gw{i} = sdpvar(gn+1, tn, 'full');
gz{i} = binvar(gn, tn, 'full');
end
% 场景2线性化变量
x_pfz = sdpvar(6, tn, 'full');
gwz = cell(1,6);
gzz = cell(1,6);
for i = 1:6
gwz{i} = sdpvar(gn+1, tn, 'full');
gzz{i} = binvar(gn, tn, 'full');
end
% 公共参数
gl1 = (Pmax - Pmin) ./ gn;
gl2 = zeros(6, gn+1);
gl2z = zeros(6, gn+1);
for i = 1:6
gl2(i,:) = linspace(Pmin(i), Pmax(i), gn+1);
gl2z(i,:) = linspace(Pmin(i), Pmax(i), gn+1);
end
% 场景1线性化约束
for i = 1:6
% 分段线性化约束
con = [con, x_pf(i,:) == (gl2(i,:).^2) * gw{i}];
con = [con, Pgone(i,:) == gl2(i,:) * gw{i}];
% 凸组合约束
con = [con, sum(gw{i}) == 1]; % 权重和为1
con = [con, gw{i} >= 0]; % 权重非负
% 特殊有序集约束（SOS2）
for j = 1:tn
con = [con, sos2(gw{i}(:,j))];
end
% 连接二进制变量
for k = 1:gn
con = [con, implies(gz{i}(k,j), [gw{i}(k,j)>=0, gw{i}(k+1,j)>=0])];
end
con = [con, sum(gz{i},1) == 1]; % 每列只有一个1，表示激活一个分段
% 确保机组出力在范围内
con = [con, Pgone(i,:) >= Pmin(i), Pgone(i,:) <= Pmax(i)];
end
% 场景2线性化约束
for i = 1:6
con = [con, x_pfz(i,:) == (gl2z(i,:).^2) * gwz{i}];
con = [con, Pgtwo(i,:) == gl2z(i,:) * gwz{i}];
con = [con, sum(gwz{i}) == 1];
con = [con, gwz{i} >= 0];
for j = 1:tn
con = [con, sos2(gwz{i}(:,j))];
end
for k = 1:gn
con = [con, implies(gzz{i}(k,j), [gwz{i}(k,j)>=0, gwz{i}(k+1,j)>=0])];
end
con = [con, sum(gzz{i},1) == 1];
con = [con, Pgtwo(i,:) >= Pmin(i), Pgtwo(i,:) <= Pmax(i)];
end
% 机组启停逻辑约束
yyone = binvar(6, 23, 'full');
yytwo = binvar(6, 23, 'full');
for i = 1:23
for j = 1:6
% 状态变化变量yyone(j,i)表示机组j在i时刻到i+1时刻状态变化
con = [con, yyone(j,i) <= ugone(j,i+1)];
con = [con, yyone(j,i) <= 1 - ugone(j,i)];
con = [con, yyone(j,i) >= ugone(j,i+1) - ugone(j,i)];
con = [con, yytwo(j,i) <= ugtwo(j,i+1)];
con = [con, yytwo(j,i) <= 1 - ugtwo(j,i)];
con = [con, yytwo(j,i) >= ugtwo(j,i+1) - ugtwo(j,i)];
end
end
% 风电启停逻辑约束
yygone = binvar(1, 23, 'full');
yygtwo = binvar(1, 23, 'full');
for i = 1:23
con = [con, yygone(1,i) <= upwone(1,i+1)];
con = [con, yygone(1,i) <= 1 - upwone(1,i)];
con = [con, yygone(1,i) >= upwone(1,i+1) - upwone(1,i)];
con = [con, yygtwo(1,i) <= upwtwo(1,i+1)];
con = [con, yygtwo(1,i) <= 1 - upwtwo(1,i)];
con = [con, yygtwo(1,i) >= upwtwo(1,i+1) - upwtwo(1,i)];
end
% 网损约束
con = [con, 0 <= plossone <= pload1, 0 <= plosstwo <= pload1];
% 目标函数组件
% 常规机组成本
fpg1 = sdpvar(6, tn); % 预分配
fpg2 = sdpvar(6, tn);
for gt = 1:tn
for gi = 1:6
fpg1(gi,gt) = mp(gi,4)*x_pf(gi,gt) + mp(gi,5)*Pgone(gi,gt) + mp(gi,6);
fpg2(gi,gt) = mp(gi,4)*x_pfz(gi,gt) + mp(gi,5)*Pgtwo(gi,gt) + mp(gi,6);
end
end
% 总常规机组成本（考虑场景概率）
fg_scene = ps(1)*sum(sum(fpg1)) + ps(2)*sum(sum(fpg2));
% 机组启停变化成本
fg_change = 0.9 * (ps(1)*sum(sum(yyone)) + ps(2)*sum(sum(yytwo)));
fg = fg_scene + fg_change;
% 储能成本
fess = ps(1)*sum(0.53*abs(Perssone) + 0.12*abs(Perssone)) + ...
ps(2)*sum(0.53*abs(Persstwo) + 0.12*abs(Persstwo)) + ...
sum(0.012*(uersstwo + uersstwo)); % 注意：这里uersstwo是二值变量，直接求和
% 风电成本（包括运行成本和弃风惩罚）
% 注意：弃风惩罚为 kc*(可用风电-实际风电)，这里kc=5.2
fdg = ps(1)*(sum(0.02*Pwone) + 5.2*sum(pw1 - Pwone)) + ...
ps(2)*(sum(0.02*Pwtwo) + 5.2*sum(pw2 - Pwtwo)) + ...
0.9*sum(ps(1)*yygone + ps(2)*yygtwo);
% IDR成本
fload = ps(1)*(sum(Ca*abs(Pidraone) + Cb*abs(Pidrbone)) + 15.9*sum(sum(plossone))) + ...
ps(2)*(sum(Ca*abs(Pidratwo) + Cb*abs(Pidrbtwo)) + 15.9*sum(sum(plosstwo)));
% 总目标函数
f1 = fg + fess + fdg + fload;
%% ======================= 模型求解 =======================
ops = sdpsettings('solver', 'cplex', 'verbose', 2, 'usex0', 0);
% 检查约束是否包含NaN
try
% 尝试求解
result = optimize(con, f1, ops);
catch ME
disp('求解过程中出错:');
disp(ME.message);
% 检查目标函数中是否有NaN
if any(isnan(value(f1)))
disp('目标函数值包含NaN');
end
% 检查约束中是否有NaN
for i = 1:length(con)
if any(isnan(double(con(i))))
disp(['约束 ', num2str(i), ' 包含NaN']);
end
end
end
%% ======================= 结果输出 =======================
% 场景1结果
thetaone = value(thetaone);
Pgone = value(Pgone);
Pwone = value(Pwone);
upwone = value(upwone);
Psone = value(Psone);
Perssone = value(Perssone);
uerssone = value(uerssone);
Pdgone = value(Pdgone);
Pidrbone = value(Pidrbone);
plossone = value(plossone);
ugone = value(ugone);
Ppdrone = value(Ppdrone);
Pidraone = value(Pidraone);
% 场景2结果
thetatwo = value(thetatwo);
Pgtwo = value(Pgtwo);
Pwtwo = value(Pwtwo);
upwtwo = value(upwtwo);
Pstwo = value(Pstwo);
Persstwo = value(Persstwo);
uersstwo = value(uersstwo);
Pdgtwo = value(Pdgtwo);
Pidrbtwo = value(Pidrbtwo);
plosstwo = value(plosstwo);
ugtwo = value(ugtwo);
Ppdrtwo = value(Ppdrtwo);
Pidratwo = value(Pidratwo);
%% ======================= 结果可视化 =======================
% 场景1出力组成
figure;
ypicb = [min(Psone,0); min(Perssone,0)]';
ypic = [max(Psone,0); max(Perssone,0); sum(Pgone,1); Pwone; sum(Pidrbone+Ppdrone+Pidraone)]';
bar(ypicb, 'stacked');
hold on;
bar(ypic, 'stacked');
legend('水电充电', '电化学储能充电', '水电放电', '电化学储能放电', '传统机组', '风电', 'DR');
xlabel('时间');
ylabel('功率');
title('场景1出力组成');
% 场景2出力组成
figure;
ypicb = [min(Pstwo,0); min(Persstwo,0)]';
ypic = [max(Pstwo,0); max(Persstwo,0); sum(Pgtwo); Pwtwo; sum(Pidrbtwo+Ppdrtwo+Pidratwo)]';
bar(ypicb, 'stacked');
hold on;
bar(ypic, 'stacked');
legend('水电充电', '电化学储能充电', '水电放电', '电化学储能放电', '传统机组', '风电', 'DR');
xlabel('时间');
ylabel('功率');
title('场景2出力组成');
%% ======================= 日内调度 =======================
% 初始化存储变量
upwonec = zeros(1,96); upwtwoc = zeros(1,96);
Perssonec = zeros(1,96); Persstwoc = zeros(1,96);
Pidrbonec = zeros(30,96); Pidrbtwoc = zeros(30,96);
% 生成日内负荷数据
pload12 = zeros(30,96); % 预分配
for i = 1:24
for j = 4*(i-1)+1:4*i
pload12(:,j) = (0.99 + 0.02*rand(30,1)) .* pload1(:,i);
end
end
% 分6个时段进行日内调度
for k = 1:6
pload2 = pload12(:, 16*(k-1)+1:16*k);
[upwone, upwtwo, Perssone, Persstwo, Pidrbone, Pidrbtwo] = ...
main2(k, pload2, ugone, ugtwo, Psone, Pstwo, Ppdrone, Ppdrtwo, Pidraone, Pidratwo);
upwonec(:,16*(k-1)+1:16*k) = upwone;
upwtwoc(:,16*(k-1)+1:16*k) = upwtwo;
Perssonec(:,16*(k-1)+1:16*k) = Perssone;
Persstwoc(:,16*(k-1)+1:16*k) = Persstwo;
Pidrbonec(:,16*(k-1)+1:16*k) = Pidrbone;
Pidrbtwoc(:,16*(k-1)+1:16*k) = Pidrbtwo;
end
%% ======================= 实时调度 =======================
% 初始化变量
upwtwo23 = zeros(1,288); Persstwo23 = zeros(1,288);
Pidrbtwo23 = zeros(30,288); pw23 = zeros(1,288);
ugtwo23 = zeros(6,288); Pstwo23 = zeros(1,288);
Pidratwo23 = zeros(30,288); Ppdrtwo23 = zeros(30,288);
pload23 = zeros(30,288); pload2s = zeros(30,288);
% 数据转换
for k = 1:96
upwtwo23(:,3*(k-1)+1:3*k) = upwtwoc(:,k);
Persstwo23(:,3*(k-1)+1:3*k) = Persstwoc(:,k);
Pidrbtwo23(:,3*(k-1)+1:3*k) = repmat(Pidrbtwoc(:,k),1,3);
end
for k = 1:24
pw23(:,12*(k-1)+1:12*k) = pw2(:,k);
ugtwo23(:,12*(k-1)+1:12*k) = repmat(ugtwo(:,k),1,12);
Pstwo23(:,12*(k-1)+1:12*k) = Pstwo(:,k);
Pidratwo23(:,12*(k-1)+1:12*k) = repmat(Pidratwo(:,k),1,12);
Ppdrtwo23(:,12*(k-1)+1:12*k) = repmat(Ppdrtwo(:,k),1,12);
end
% 生成实时负荷数据
for i = 1:96
for j = 3*(i-1)+1:3*i
pload23(:,j) = (0.995 + 0.01*rand(30,1)) .* pload12(:,i);
pload2s(:,j) = (0.995 + 0.01*rand(30,1)) .* pload12(:,i); % 实际负荷
end
end
% 执行实时调度
Pgtwo23 = []; Pidrdtwo23 = [];
for k = 1:1 % 示例仅执行第一时段
[Pgtwo23_k, Pidrdtwo23_k] = main3(upwtwo23(:,k), Persstwo23(:,k), ...
Pidrbtwo23(:,k), pw23(:,k), ugtwo23(:,k), Pstwo23(:,k), ...
Pidratwo23(:,k), Ppdrtwo23(:,k), pload23(:,k));
Pgtwo23 = [Pgtwo23, Pgtwo23_k];
Pidrdtwo23 = [Pidrdtwo23, Pidrdtwo23_k];
end
%% ======================= 负荷曲线可视化 =======================
% 数据格式转换
pload12z = zeros(30,288);
for i = 1:96
pload12z(:,3*(i-1)+1:3*i) = repmat(pload12(:,i),1,3);
end
ploadz = zeros(30,288);
for i = 1:24
ploadz(:,12*(i-1)+1:12*i) = repmat(pload1(:,i),1,12);
end
% 绘制负荷曲线
figure;
plot(sum(pload23,1), 'r:', 'LineWidth', 1); hold on;
plot(sum(pload12z,1), 'b-', 'LineWidth', 1);
plot(sum(ploadz,1), 'k--', 'LineWidth', 1);
plot(sum(pload2s,1), 'm--', 'LineWidth', 1);
legend('实时预测负荷', '日内预测负荷', '日前预测负荷', '实际负荷');
xlabel('调度时刻');
ylabel('功率/MW');
set(gca, 'XTick', 0:12:288);
set(gca, 'XTickLabel', {'0','1','2','3','4','5','6','7','8','9','10','11','12',...
'13','14','15','16','17','18','19','20','21','22','23','24'});
title('多时间尺度负荷曲线对比');
最新发布
探讨了一个复杂的数学问题，即在未知的n进制下，如何通过等式567*456=150216来确定n的具体值。通过逐步解析，将等式两边转化为十进制，并进行化简比较，最终得出n等于18的结论。
扫一扫
182
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
