开端支路分布因子原理推导及matlab程序编写电力系统老司机都懂支路开断后潮流重新分布的骚操作得靠分布因子。这玩意儿本质上是个灵敏度系数说白了就是某条线路断了之后其他线路得帮它扛多少锅。今天咱们直接上手推导顺手撸个Matlab脚本出来。先来点硬核推导热热身。假设原始网络节点导纳矩阵是B线路l-m开断后相当于在原有导纳基础上并联了个负导纳这波反向操作666。利用矩阵求逆辅助定理分布因子D_{k-l}可以写成D (Xk - Xm)(Pl^{(0)} - Pm^{(0)}) / x_lm开端支路分布因子原理推导及matlab程序编写其中X是节点阻抗矩阵对应元素P是原始潮流。不过这么搞计算量太大咱们换个聪明法子——直接修改导纳矩阵求逆。上代码先整个IEEE 9节点系统练手% 节点数据 [节点号 类型 Pd Qd V] bus [ 1 3 0 0 1.04 2 2 0 0 1.025 3 2 0 0 1.025 4 1 0 0 1 5 1 90 30 1 6 1 0 0 1 7 1 100 35 1 8 1 0 0 1 9 1 125 50 1]; % 线路数据 [发端 收端 r x b] branch [ 1 4 0 0.0576 0 4 5 0.017 0.092 0.158 5 6 0.039 0.17 0.358 3 6 0 0.0586 0 6 7 0.0119 0.1008 0.209 7 8 0.0085 0.072 0.149 8 2 0 0.0625 0 8 9 0.032 0.161 0.306 9 4 0.01 0.085 0.176];构建导纳矩阵时要注意直流潮流只考虑电抗nbus size(bus,1); B_prime zeros(nbus); for k 1:size(branch,1) i branch(k,1); j branch(k,2); x branch(k,4); B_prime(i,i) B_prime(i,i) 1/x; B_prime(j,j) B_prime(j,j) 1/x; B_prime(i,j) B_prime(i,j) - 1/x; B_prime(j,i) B_prime(j,i) - 1/x; end B_prime(1,:) []; B_prime(:,1) []; % 消去平衡节点处理线路开断的关键在这function D calc_distribution_factor(branch_idx) br branch(branch_idx,:); i br(1); j br(2); x_ij br(4); % 构造修改项 delta_B sparse([i j], [i j], [1/x_ij, 1/x_ij], nbus, nbus)... - sparse([i j], [j i], 1/x_ij, nbus, nbus); delta_B(1,:) []; delta_B(:,1) []; % 平衡节点处理 % 计算修正逆矩阵 B_new_inv inv(B_prime delta_B); % 计算分布因子 D zeros(size(branch,1),1); for k 1:size(branch,1) m branch(k,1); n branch(k,2); D(k) (B_new_inv(m-1,:) - B_new_inv(n-1,:)) * (B_prime \ P_ref); end D D / x_ij; end这里有几个坑要注意1节点编号从1开始但矩阵计算时要去掉平衡节点2sparse矩阵能大幅提升大电网计算效率3原始潮流P_ref需要提前用直流潮流算好。跑个实例验证下% 计算原始潮流 P_ref [0; 0; -90; 0; -100; 0; -125]; theta B_prime \ P_ref; % 开断线路5-6 D calc_distribution_factor(3); disp(线路5-6开断时的分布因子:); disp(D);输出结果里负数表示潮流方向反转。比如某线路分布因子0.35意味着原线路断了之后这条线路要多承担35%的潮流甩锅量。不过要注意当系统存在多个发电机时还需要考虑发电机出力调整的影响这时候得用广义分布因子——这留到下期再唠。最后给小白提个醒实际工程中别直接用inv()数值不稳定应该用LU分解或者反斜杠运算符。但咱们教学演示嘛怎么直观怎么来。完整代码已打包关注公众号回复DF领取~手动狗头