最近在学习电力系统运行分析,需要进行潮流计算,但纯手算的话着实有点费劲,尤其是在第二次迭代的时候。于是忙里偷闲用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);

对应题目见下图:

电力系统运行分析 - 作业与思考4

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

最后修改:2021 年 02 月 20 日 06 : 10 PM
如果觉得我的文章对你有用,请随意赞赏