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

博主： Angus
发布时间：2020 年 04 月 06 日
870 次浏览
暂无评论
5565字数
分类： Geek

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

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

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

```matlab
% 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](https://blog.angustar.com/usr/uploads/2021/01/4257278044.jpg)

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

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

最后修改：2021 年 02 月 20 日

如果觉得我的文章对你有用，请随意赞赏

发表评论取消回复
本站使用 Cookie 技术以便您下次快速评论，继续评论表示您已同意该条款。请您在评论时尽量使用 QQ 邮箱，邮箱信息将严格保密。

评论 *

私密评论

名称 *

🎲

邮箱 *

地址

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

Angus • 2020 年 04 月 06 日

最近在学习电力系统运行分析，需要进行潮流计算，但纯手算的话着实有点费劲，尤其是在第二次迭代的时候。于是忙里偷闲用MATLAB写了这段代码，仅供参考。其实这段代码也没有用到什么特殊语法，只是简单的把题目计算所需要的公式转换为代码形式了，只针对本题有效（对应题目见文末）：**节点1为平衡节点，节点2、3为PQ节点，雅克比矩阵为4×4形式。**进行迭代时，只需要输入初始电压U、初始相位角Theta（或输入上次迭代的结果），运行代码，即可计算出此次迭代的结果。```matlab 
% 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](https://blog.angustar.com/usr/uploads/2021/01/4257278044.jpg)<div class="tip inlineBlock success">

*本文无参考文献，代码全是自己敲出来的，而且代码居然基本上一遍就能过哈哈哈，超开心，毕竟好久没写出来能一遍过的代码啦！写代码前后一共只花了不到2小时哦*😀
</div><hr class="content-copyright" style="margin-top:50px" /><blockquote class="content-copyright" style="font-style:normal">版权属于：Angus本文链接：<a class="content-copyright" href="https://blog.angustar.com/archives/matlab-script-for-power-flow-calculation.html">https://blog.angustar.com/archives/matlab-script-for-power-flow-calculation.html</a>所有原创文章采用<a class="no-external-link" href="http://creativecommons.org/licenses/by-nc-sa/4.0/" target="_blank">知识共享署名-非商业性使用 4.0 国际许可协议</a>进行许可。 您可以自由的转载和修改，但请务必注明文章来源并且不可用于商业目的。</blockquote>

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

发表评论取消回复
本站使用 Cookie 技术以便您下次快速评论，继续评论表示您已同意该条款。请您在评论时尽量使用 QQ 邮箱，邮箱信息将严格保密。

论文排版神器 —— LaTeX 套装推荐

在 Typecho 的评论区中插入 DPlayer 播放器

Tailscale 的 DERP 中继服务搭建与配置

搭建一个属于自己的Server酱吧！

A translation for "Flexible Camera Calibration by Viewing a Plane from Unknown Orientations"

又一个优秀的 Bing 每日壁纸展示项目

论文排版神器 —— LaTeX 套装推荐

工程领域常用希腊字母的发音

Rabbit in The Wonderland

憾别医生李文亮：他爱美食，喜欢音乐，喜欢漫威，还会做手工

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

发表评论 取消回复 本站使用 Cookie 技术以便您下次快速评论，继续评论表示您已同意该条款。请您在评论时尽量使用 QQ 邮箱，邮箱信息将严格保密。

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

发表评论取消回复
本站使用 Cookie 技术以便您下次快速评论，继续评论表示您已同意该条款。请您在评论时尽量使用 QQ 邮箱，邮箱信息将严格保密。