% 复现文献：Improved Affine Arithmetic-Based Power Flow Computation for Distribution Systems Considering Uncertainties
% 使用改进仿射算法（IAA）计算IEEE33节点系统的电压不确定范围
% 主要不确定性来源：风电（风速v）和光伏（温度T、光照强度G）
% 需要安装INTLAB工具箱（区间运算工具箱）
clear all
clc
close all
% 添加INTLAB工具箱路径（请根据实际安装路径修改）
addpath(genpath('D:\Matlab2024binstall\toolbox\intlab'));
% 兼容性处理 - 在新版MATLAB中定义bitmax为flintmax
if ~exist('bitmax', 'builtin') && exist('flintmax', 'builtin')
bitmax = @() flintmax;
end
startintlab
%% 1. IEEE33节点系统参数
% 节点数据：节点编号，类型（1=PQ，2=PV，3=平衡节点），P负荷(kW)，Q负荷(kvar)，P发电(kW)，Q发电(kvar)
% 根据文献图4的改进IEEE33节点拓扑
busdata = [
1 3 0 0 0 0;
2 1 100 60 0 0;
3 1 90 40 0 0;
4 1 120 80 0 0;
5 1 60 30 0 0;
6 1 60 20 0 0;
7 1 200 100 0 0;
8 1 200 100 0 0;
9 1 60 20 0 0;
10 1 60 20 0 0;
11 1 45 30 0 0;
12 1 60 35 500 0; % 风电节点500kW
13 1 60 35 0 0;
14 1 120 80 0 0;
15 1 60 10 0 0;
16 1 60 20 500 0; % 风电节点500kW
17 1 60 20 0 0;
18 1 90 40 0 0;
19 1 90 40 0 0;
20 1 90 40 0 0;
21 1 90 40 0 0;
22 1 90 40 0 0;
23 1 90 50 0 0;
24 1 420 200 600 0; % 光伏节点600kW
25 1 420 200 0 0;
26 1 60 25 0 0;
27 1 60 25 0 0;
28 1 60 20 600 0; % 光伏节点600kW
29 1 120 70 0 0;
30 1 200 600 0 0;
31 1 150 70 500 0; % 风电节点500kW
32 1 210 100 0 0;
33 1 60 40 0 0;
];
% 线路数据：首端节点，末端节点，电阻(Ω)，电抗(Ω)
linedata = [
1 2 0.0922 0.0470;
2 3 0.4930 0.2512;
3 4 0.3660 0.1864;
4 5 0.3811 0.1941;
5 6 0.8190 0.7070;
6 7 0.1872 0.6188;
7 8 0.7114 0.2351;
8 9 1.0300 0.7400;
9 10 1.0440 0.7400;
10 11 0.1966 0.0650;
11 12 0.3744 0.1238;
12 13 1.4680 1.1550;
13 14 0.5416 0.7129;
14 15 0.5910 0.5260;
15 16 0.7463 0.5450;
16 17 1.2890 1.7210;
17 18 0.7320 0.5740;
2 19 0.1640 0.1565;
19 20 1.5042 1.3554;
20 21 0.4095 0.4784;
21 22 0.7089 0.9373;
3 23 0.4512 0.3083;
23 24 0.8980 0.7091;
24 25 0.8960 0.7011;
6 26 0.2030 0.1034;
26 27 0.2842 0.1447;
27 28 1.0590 0.9337;
28 29 0.8042 0.7006;
29 30 0.5075 0.2585;
30 31 0.9744 0.9630;
31 32 0.3105 0.3619;
32 33 0.3410 0.5302;
];
% 基准值
Vbase = 12.66e3; % 电压基准值(V)
Sbase = 10e6; % 功率基准值(VA)
Zbase = Vbase^2/Sbase; % 阻抗基准值(Ω)
% 将线路参数转换为标幺值
linedata(:,3:4) = linedata(:,3:4)/Zbase;
% 节点数
nbus = size(busdata,1);
%% 2. 构建风电和光伏的多噪声项仿射模型
% 根据文献中的公式(26)-(31)构建不确定模型
% 2.1 风电仿射模型
% 风速仿射模型：v_affine = v0 + 0.724ε1 + 0.367ε2 + 0.119ε3
v0 = 10; % 额定风速(m/s)，假设值
v_affine = affineinit(v0, [0.724, 0.367, 0.119]);
% 风电功率-风速特性曲线：Pw = 0.00422v^5 - 0.20823v^4 + 3.3675v^3 - 19.0325v^2 + 54.025v - 51.5
% 计算风电功率仿射形式
Pw_coeff = [0.00422, -0.20823, 3.3675, -19.0325, 54.025, -51.5];
Pw_affine = polyval(Pw_coeff, v_affine);
% 添加额外噪声项ε4，系数为0.37f'(v0)
% 计算导数f'(v0)
syms v_sym
Pw_sym = poly2sym(Pw_coeff, v_sym);
dPw_sym = diff(Pw_sym, v_sym);
dPw_v0 = double(subs(dPw_sym, v_sym, v0));
% 最终风电功率仿射模型
Pw_final = Pw_affine + 0.37*dPw_v0*affineinit(0, 1, 4); % ε4
% 2.2 光伏仿射模型
% 温度仿射模型：T_affine = T0 + 0.637ε1 + 0.185ε2
T0 = 25; % 额定温度(°C)
T_affine = affineinit(T0, [0.637, 0.185]);
% 光照强度仿射模型：G_affine = G0 + 0.815ε1
G0 = 800; % 额定光照强度(W/m²)
G_affine = affineinit(G0, 0.815);
% 光伏输出功率模型：Ppv = PSTC * (GT/GSTC) * [1 - K(TC - Tτ)]
PSTC = 600; % kW，最大测试功率
GSTC = 1000; % W/m²，标准测试条件光照强度
K = -0.0017; % -0.17%/°C
T = 25; % °C
% 计算光伏功率仿射形式
Ppv_affine = PSTC * (G_affine/GSTC) .* (1 - K*(T_affine - T));
%% 3. 将不确定性注入到相应节点
% 节点12,16,31：风电500kW
% 节点24,28：光伏600kW
% 创建节点功率的仿射形式
P_load = busdata(:,3)/1000; % kW转换为MW
Q_load = busdata(:,4)/1000; % kvar转换为Mvar
P_gen = zeros(nbus,1);
Q_gen = zeros(nbus,1);
% 将确定性的发电功率转换为仿射形式
for i = 1:nbus
if busdata(i,5) > 0
P_gen(i) = busdata(i,5)/1000; % kW转换为MW
end
end
% 创建仿射形式的节点注入功率
P_inj_affine = cell(nbus,1);
Q_inj_affine = cell(nbus,1);
for i = 1:nbus
% 初始化为中心值
P_inj_affine{i} = affineinit(P_gen(i) - P_load(i), zeros(1,4));
Q_inj_affine{i} = affineinit(Q_gen(i) - Q_load(i), zeros(1,4));
end
% 将不确定性注入到相应节点
% 节点12：风电
P_inj_affine{12} = P_inj_affine{12} + Pw_final/1000; % kW转换为MW
% 节点16：风电
P_inj_affine{16} = P_inj_affine{16} + Pw_final/1000; % kW转换为MW
% 节点31：风电
P_inj_affine{31} = P_inj_affine{31} + Pw_final/1000; % kW转换为MW
% 节点24：光伏
P_inj_affine{24} = P_inj_affine{24} + Ppv_affine/1000; % kW转换为MW
% 节点28：光伏
P_inj_affine{28} = P_inj_affine{28} + Ppv_affine/1000; % kW转换为MW
%% 4. 基于改进仿射算法（IAA）的前推回代潮流计算
% 初始化
V_affine = cell(nbus,1);
for i = 1:nbus
V_affine{i} = affineinit(1.0, zeros(1,4)); % 初始电压设为1.0∠0° p.u.
end
% 创建线路电流仿射形式
I_line_affine = cell(size(linedata,1),1);
for i = 1:size(linedata,1)
I_line_affine{i} = affineinit(0, zeros(1,4));
end
% 收敛精度
tolerance = 1e-10;
max_iter = 100;
% 存储每次迭代的电压
V_mag = zeros(nbus, max_iter);
V_angle = zeros(nbus, max_iter);
% 改进仿射除法函数
function result = improved_affine_division(numerator, denominator)
% 基于区间泰勒展开的改进仿射除法
% 文献中的公式(20)-(21)
% 提取分母的中心值和半径
y_c = mid(denominator);
y_r = rad(denominator);
% 避免除以零
if y_c == 0
y_c = 1e-10;
end
% 计算改进仿射除法的系数
z0 = 1/y_c + y_r^2/(4*y_c^3);
z1 = -y_r/y_c^2;
z2 = y_r^2/(4*y_c^3);
% 构建结果仿射形式
result = z0 + z1*affineinit(0,1,1) + z2*affineinit(0,1,2);
% 乘以分子
result = numerator .* result;
end
% 前推回代迭代
for iter = 1:max_iter
% 存储上一次迭代的电压
V_prev = V_affine;
% 回代过程：计算节点电流
I_inj_affine = cell(nbus,1);
for i = 1:nbus
% 使用改进仿射除法计算节点电流：I_i = conj(S_i / V_i)
S_complex = P_inj_affine{i} + 1j*Q_inj_affine{i};
I_inj_affine{i} = improved_affine_division(conj(S_complex), V_affine{i});
end
% 从末端节点开始，计算线路电流
% 首先需要确定节点的层次关系（从末端到首端）
% 这里简化处理：假设已知节点层次关系
% 初始化线路电流
for i = size(linedata,1):-1:1
from_bus = linedata(i,1);
to_bus = linedata(i,2);
% 线路电流等于下游所有节点电流之和
I_line_affine{i} = I_inj_affine{to_bus};
% 查找所有以to_bus为首端的线路
for j = i+1:size(linedata,1)
if linedata(j,1) == to_bus
I_line_affine{i} = I_line_affine{i} + I_line_affine{j};
end
end
end
% 前推过程：计算节点电压
for i = 1:size(linedata,1)
from_bus = linedata(i,1);
to_bus = linedata(i,2);
R = linedata(i,3);
X = linedata(i,4);
Z = R + 1j*X;
% 计算电压降：V_to = V_from - Z * I_line
V_affine{to_bus} = V_affine{from_bus} - Z * I_line_affine{i};
end
% 检查收敛
max_diff = 0;
for i = 1:nbus
diff_val = abs(mid(V_affine{i}) - mid(V_prev{i}));
if diff_val > max_diff
max_diff = diff_val;
end
end
% 存储当前迭代的电压幅值和相角
for i = 1:nbus
V_mag(i, iter) = abs(mid(V_affine{i}));
V_angle(i, iter) = angle(mid(V_affine{i})) * 180/pi;
end
fprintf('迭代 %d: 最大电压变化 = %.8f p.u.\n', iter, max_diff);
if max_diff < tolerance
fprintf('潮流计算在 %d 次迭代后收敛。\n', iter);
break;
end
if iter == max_iter
fprintf('达到最大迭代次数 %d，可能未收敛。\n', max_iter);
end
end
%% 5. 提取电压不确定范围
% 提取最终电压的区间范围
V_mag_interval = zeros(nbus, 2);
V_angle_interval = zeros(nbus, 2);
for i = 1:nbus
V_interval = intval(V_affine{i});
V_mag_interval(i,1) = inf(abs(V_interval));
V_mag_interval(i,2) = sup(abs(V_interval));
V_angle_interval(i,1) = inf(angle(V_interval) * 180/pi);
V_angle_interval(i,2) = sup(angle(V_interval) * 180/pi);
end
%% 6. 绘制电压不确定范围（与文献图5对比）
figure;
hold on;
grid on;
% 绘制电压幅值区间
for i = 1:nbus
plot([i, i], [V_mag_interval(i,1), V_mag_interval(i,2)], 'b-', 'LineWidth', 1.5);
end
% 绘制中心值
plot(1:nbus, mid(V_mag_interval(:,1) + V_mag_interval(:,2))/2, 'ro', 'MarkerSize', 4, 'MarkerFaceColor', 'r');
xlabel('节点编号');
ylabel('电压幅值 (p.u.)');
title('IEEE33节点系统电压幅值不确定范围（改进仿射算法）');
legend('不确定范围', '中心值', 'Location', 'best');
xlim([1, nbus]);
ylim([0.9, 1.05]);
% 绘制电压相角区间
figure;
hold on;
grid on;
for i = 1:nbus
plot([i, i], [V_angle_interval(i,1), V_angle_interval(i,2)], 'b-', 'LineWidth', 1.5);
end
plot(1:nbus, mid(V_angle_interval(:,1) + V_angle_interval(:,2))/2, 'ro', 'MarkerSize', 4, 'MarkerFaceColor', 'r');
xlabel('节点编号');
ylabel('电压相角 (度)');
title('IEEE33节点系统电压相角不确定范围（改进仿射算法）');
legend('不确定范围', '中心值', 'Location', 'best');
xlim([1, nbus]);
%% 7. 输出结果
fprintf('\n=== IEEE33节点系统电压不确定范围 ===\n');
fprintf('节点\t电压幅值下限(p.u.)\t电压幅值上限(p.u.)\t电压相角下限(度)\t电压相角上限(度)\n');
for i = 1:nbus
fprintf('%d\t%.6f\t\t%.6f\t\t%.6f\t\t%.6f\n', i, ...
V_mag_interval(i,1), V_mag_interval(i,2), ...
V_angle_interval(i,1), V_angle_interval(i,2));
end
% 保存结果
save('IEEE33_IAA_Results.mat', 'V_mag_interval', 'V_angle_interval', 'V_affine');
北京时间4月28日21:00-23:00，BitMax平台将举行STAKE代币竞拍活动，用户可使用USDT参与，竞拍区间为400-550USDT，前400名中标者将以统一价格获得1000STAKE，且无锁定期限制。
扫一扫
2500
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
