用于电力系统简单潮流计算的 MATLAB 脚本

Angus

2020 年 04 月 06 日

1009 次浏览

暂无评论

5565字数

Geek

最近在学习电力系统运行分析，需要进行潮流计算，但纯手算的话着实有点费劲，尤其是在第二次迭代的时候。于是忙里偷闲用MATLAB写了这段代码，仅供参考。

其实这段代码也没有用到什么特殊语法，只是简单的把题目计算所需要的公式转换为代码形式了，只针对本题有效（对应题目见文末）：节点1为平衡节点，节点2、3为PQ节点，雅克比矩阵为4×4形式。

进行迭代时，只需要输入初始电压U、初始相位角Theta（或输入上次迭代的结果），运行代码，即可计算出此次迭代的结果。

% Written by angustar.com -- 2020.4.6
% 定义节点导纳矩阵Y
Y=[complex(0.05,-0.0079) complex(-0.028,0.045) complex(-0.022,0.034);...,
   complex(-0.028,0.045) complex(0.051,-0.066) complex(-0.023,0.021);...,
   complex(-0.022,0.034) complex(-0.023,0.021) complex(0.045,-0.055)];
% 定义电压U
U=[115;115;115];
% 定义相位角Theta
Theta=[0;0;0];
% 定义功率P
% 由于节点1为平衡节点，因此仅考虑P2、P3、Q2、Q3，为4x1矩阵
PQ=[20;25;15;18];
% 计算雅克比矩阵
% 注意：以下计算仅适用于4x4的雅克比矩阵，其它阶数的雅克比矩阵要做相应修改
H22=U(2,1)*(U(1,1)*(real(Y(2,1))*sin(Theta(2,1)-Theta(1,1))-imag(Y(2,1))...,
    *cos(Theta(2,1)-Theta(1,1)))+U(3,1)*(real(Y(2,3))*sin(Theta(2,1)...,
    -Theta(3,1))-imag(Y(2,3))*cos(Theta(2,1)-Theta(3,1))));
H23=-U(2,1)*U(3,1)*(real(Y(2,3))*sin(Theta(2,1)-Theta(3,1))...,
    -imag(Y(2,3))*cos(Theta(2,1)-Theta(3,1)));
H32=-U(2,1)*U(3,1)*(real(Y(3,2))*sin(Theta(3,1)-Theta(2,1))...,
    -imag(Y(3,2))*cos(Theta(3,1)-Theta(2,1)));
H33=U(3,1)*(U(1,1)*(real(Y(3,1))*sin(Theta(3,1)-Theta(1,1))-imag(Y(3,1))...,
    *cos(Theta(3,1)-Theta(1,1)))+U(2,1)*(real(Y(3,2))*sin(Theta(3,1)...,
    -Theta(2,1))-imag(Y(3,2))*cos(Theta(3,1)-Theta(2,1))));
N22=-2*U(2,1)*real(Y(2,2))-U(1,1)*(real(Y(2,1))*cos(Theta(2,1)-Theta(1,1))...,
    +imag(Y(2,1))*sin(Theta(2,1)-Theta(1,1)))-U(3,1)*(real(Y(2,3))...,
    *cos(Theta(2,1)-Theta(3,1))+imag(Y(2,3))*sin(Theta(2,1)-Theta(3,1)));
N23=-U(2,1)*(real(Y(2,3))*cos(Theta(2,1)-Theta(3,1))+imag(Y(2,3))*sin(Theta(2,1)-Theta(3,1)));
N32=-U(3,1)*(real(Y(3,2))*cos(Theta(3,1)-Theta(2,1))+imag(Y(3,2))*sin(Theta(3,1)-Theta(2,1)));
N33=-2*U(3,1)*real(Y(3,3))-U(1,1)*(real(Y(3,1))*cos(Theta(3,1)-Theta(1,1))...,
    +imag(Y(3,1))*sin(Theta(3,1)-Theta(1,1)))-U(2,1)*(real(Y(3,2))...,
    *cos(Theta(3,1)-Theta(2,1))+imag(Y(3,2))*sin(Theta(3,1)-Theta(2,1)));
K22=-U(2,1)*(U(1,1)*(real(Y(2,1))*cos(Theta(2,1)-Theta(1,1))...,
    +imag(Y(2,1))*sin(Theta(2,1)-Theta(1,1)))+U(3,1)*(real(Y(2,3))...,
    *cos(Theta(2,1)-Theta(3,1))+imag(Y(2,3))*sin(Theta(2,1)-Theta(3,1))));
K23=-U(2,1)*U(3,1)*(real(Y(2,3))*cos(Theta(2,1)-Theta(3,1))+imag(Y(2,3))*sin(Theta(2,1)-Theta(3,1)));
K32=-U(3,1)*U(2,1)*(real(Y(3,2))*cos(Theta(3,1)-Theta(2,1))+imag(Y(3,2))*sin(Theta(3,1)-Theta(2,1)));
K33=-U(3,1)*(U(1,1)*(real(Y(3,1))*cos(Theta(3,1)-Theta(1,1))...,
    +imag(Y(3,1))*sin(Theta(3,1)-Theta(1,1)))+U(2,1)*(real(Y(3,2))...,
    *cos(Theta(3,1)-Theta(2,1))+imag(Y(3,2))*sin(Theta(3,1)-Theta(2,1))));
L22=-(U(1,1)*(real(Y(2,1))*sin(Theta(2,1)-Theta(1,1))-imag(Y(2,1))...,
    *cos(Theta(2,1)-Theta(1,1)))+U(3,1)*(real(Y(2,3))*sin(Theta(2,1)...,
    -Theta(3,1))-imag(Y(2,3))*cos(Theta(2,1)-Theta(3,1))))+2*U(2,1)*imag(Y(2,2));
L23=-U(2,1)*(real(Y(2,3))*sin(Theta(2,1)-Theta(3,1))...,
    -imag(Y(2,3))*cos(Theta(2,1)-Theta(3,1)));
L32=-U(2,1)*(real(Y(3,2))*sin(Theta(3,1)-Theta(2,1))...,
    -imag(Y(3,2))*cos(Theta(3,1)-Theta(2,1)));
L33=-(U(1,1)*(real(Y(3,1))*sin(Theta(3,1)-Theta(1,1))-imag(Y(3,1))...,
    *cos(Theta(3,1)-Theta(1,1)))+U(2,1)*(real(Y(3,2))*sin(Theta(3,1)...,
    -Theta(2,1))-imag(Y(3,2))*cos(Theta(3,1)-Theta(2,1))))+2*U(3,1)*imag(Y(3,3));
% 列出雅克比矩阵
jacobi_matrix=[H22 H23 N22 N23;H32 H33 N32 N33;K22 K23 L22 L23;K32 K33 L32 L33];
%计算功率不平衡量矩阵
Delta_P2=PQ(1,1)-U(2,1)*(U(1,1)*(real(Y(2,1))*cos(Theta(2,1)-Theta(1,1))...,
    +imag(Y(2,1))*sin(Theta(2,1)-Theta(1,1)))+U(2,1)*(real(Y(2,2))...,
    *cos(Theta(2,1)-Theta(2,1))+imag(Y(2,2))*sin(Theta(2,1)-Theta(2,1)))...,
    +U(3,1)*(real(Y(2,3))*cos(Theta(2,1)-Theta(3,1))+imag(Y(2,3))...,
    *sin(Theta(2,1)-Theta(3,1))));
Delta_P3=PQ(2,1)-U(3,1)*(U(1,1)*(real(Y(3,1))*cos(Theta(3,1)-Theta(1,1))...,
    +imag(Y(3,1))*sin(Theta(3,1)-Theta(1,1)))+U(2,1)*(real(Y(3,2))...,
    *cos(Theta(3,1)-Theta(2,1))+imag(Y(3,2))*sin(Theta(3,1)-Theta(2,1)))...,
    +U(3,1)*(real(Y(3,3))*cos(Theta(3,1)-Theta(3,1))+imag(Y(3,3))...,
    *sin(Theta(3,1)-Theta(3,1))));
Delta_Q2=PQ(3,1)-U(2,1)*(U(1,1)*(real(Y(2,1))*sin(Theta(2,1)-Theta(1,1))...,
    -imag(Y(2,1))*cos(Theta(2,1)-Theta(1,1)))+U(2,1)*(real(Y(2,2))...,
    *sin(Theta(2,1)-Theta(2,1))-imag(Y(2,2))*cos(Theta(2,1)-Theta(2,1)))...,
    +U(3,1)*(real(Y(2,3))*sin(Theta(2,1)-Theta(3,1))-imag(Y(2,3))...,
    *cos(Theta(2,1)-Theta(3,1))));
Delta_Q3=PQ(4,1)-U(3,1)*(U(1,1)*(real(Y(3,1))*sin(Theta(3,1)-Theta(1,1))...,
    -imag(Y(3,1))*cos(Theta(3,1)-Theta(1,1)))+U(2,1)*(real(Y(3,2))...,
    *sin(Theta(3,1)-Theta(2,1))-imag(Y(3,2))*cos(Theta(3,1)-Theta(2,1)))...,
    +U(3,1)*(real(Y(3,3))*sin(Theta(3,1)-Theta(3,1))-imag(Y(3,3))...,
    *cos(Theta(3,1)-Theta(3,1))));
% 列出功率不平衡量矩阵
Delta_PQ=[Delta_P2;Delta_P3;Delta_Q2;Delta_Q3];
% 计算输出结果矩阵
Delta_result=-inv(jacobi_matrix)*Delta_PQ;
% 输出计算结果
disp('雅克比矩阵为：');
disp(jacobi_matrix);
disp('功率不平衡量矩阵为：')
disp(Delta_PQ);
disp('修正方程的解为：')
disp(Delta_result);
disp('迭代后的电压为（节点1为平衡节点）：')
Adjusted_U=[U(1,1);U(2,1)+Delta_result(3,1);U(3,1)+Delta_result(4,1)];
disp(Adjusted_U);
disp('迭代后的相位角为（节点1为平衡节点）：')
Adjusted_Theta=[Theta(1,1);Theta(2,1)+Delta_result(1,1);Theta(3,1)+Delta_result(2,1)];
disp(Adjusted_Theta);
disp('此时的精度为：')
Precision=[Delta_P2;Delta_P3;Delta_result(1,1);Delta_result(2,1)];
disp(Precision);

对应题目见下图：

本文无参考文献，代码全是自己敲出来的，而且代码居然基本上一遍就能过哈哈哈，超开心，毕竟好久没写出来能一遍过的代码啦！写代码前后一共只花了不到2小时哦😀

版权属于：Angus
本文链接：https://blog.angustar.com/archives/matlab-script-for-power-flow-calculation.html
所有原创文章采用知识共享署名-非商业性使用 4.0 国际许可协议进行许可。您可以自由的转载和修改，但请务必注明文章来源并且不可用于商业目的。

用于电力系统简单潮流计算的 MATLAB 脚本