基于MATLAB的电力系统潮流计算

为了进行潮流计算，通常可以使用像MATPOWER这样的电力系统分析工具包，在MATLAB中的MATPOWER工具包中，多个不同的系统模型被定义为测试案例（case）。这些测试案例用于模拟和分析电力系统中的潮流计算、最优潮流问题、稳态分析等。常见的MATPOWER案例，包括 case4gs.m、case5.m、case22.m、case30.m、case118.m、case2383wp.m、case9241pegase.m 和 case13659pegase.m。

这些 case*.m 文件是MATPOWER工具包中用于测试电力系统分析算法的标准测试案例。它们的主要用途包括：

潮流计算：即研究电力系统中各个节点的电压、功率流向等。
最优潮流计算：通过优化算法寻求电力系统的最佳运行方案。
稳定性分析：模拟和分析系统在不同条件下的稳定性。

Case4gs系统是一个典型的小型电力系统模型，常用于电力系统的潮流分析与测试。在该模型中，系统通常包括发电机、负荷、输电线路等元件，通过潮流计算可以确定系统中各节点的电压幅值、相角以及各支路的功率流向等参数。

为了进行潮流计算，通常可以使用像MATPOWER这样的电力系统分析工具包，它提供了高效的算法和接口来求解电力系统的潮流问题。MATPOWER中的runpf函数能够快速计算出系统的潮流分布情况。

代码中包括case4gs.m、case5.m、case22.m、case30.m、case118.m、case2383wp.m、case9241pegase.m 和 case13659pegase.m的调用函数，本文以Case4gs系统为例，下面将通过MATLAB语言，结合MATPOWER工具包，介绍如何在Case4gs系统上进行电力系统潮流计算。

潮流计算基本流程

1. 从.m文件中读取数据
2. 节点编号
3. 形成导纳矩阵
4. 给定节点电压的初值
5. 计算潮流方程不平衡量
6. 计算雅可比矩阵
7. 求解修正方程，更新状态量
8. 若不平衡量的最大值小于给定的阈值或迭代次数大于最大次数，则结束计算，否则转第5步进行迭代。

1. 准备工作

首先需要安装MATPOWER，并将其路径加入MATLAB的工作目录。安装完成后，加载MATPOWER库。

% 添加MATPOWER的路径
addpath('matpower7.1'); % 根据实际安装路径进行调整

2. 导入Case4gs系统数据

MATPOWER中已经预定义了多个典型电力系统测试案例，其中包括Case4gs系统。可以通过loadcase函数加载该系统的基础数据。

% 加载Case4gs系统数据
mpc = loadcase('case4gs'); % Case4gs通常是MATPOWER中一个默认案例

3. 执行潮流计算

接下来，使用MATPOWER提供的runpf函数执行潮流计算，得到系统中各节点的电压幅值、相角，以及功率流等结果。

% 执行潮流计算
results = runpf(mpc);% 检查潮流计算是否成功
if results.successdisp('潮流计算成功');
elsedisp('潮流计算失败');
end

4. 提取并显示计算结果

潮流计算完成后，results结构体包含了计算的各种信息。可以提取并显示关键参数，例如节点电压、发电机功率、负荷功率等。

%%  计算平衡节点功率，线路功率
Y_balance = Y(Balance,:);                                                  %获取系统的平衡节点号
S_balance = U1(Balance)*conj(Y_balance)*conj(U1)*base;                     %计算平衡节点功率
%拿出用于计算线路之间的导纳参数y_ij=y_ji（即对角线上的数值为0)
y_ij = sparse(LineData(:,1),LineData(:,2),kGij+1i*kBij,BusNumber,BusNumber);
y_ij = y_ij+sparse(LineData(:,2),LineData(:,1),kGij+1i*kBij,BusNumber,BusNumber);
%y_ij=sparse(LineData(:,1),LineData(:,2),(1./(LineData(:,3)+1i*LineData(:,4)))./LineData(:,9),BusNumber,BusNumber);        %求y_ij
%y_ij=y_ij+sparse(LineData(:,2),LineData(:,1),(1./(LineData(:,3)+1i*LineData(:,4)))./LineData(:,9),BusNumber,BusNumber);   %求y_ji
%y_i0
y_i0 = sparse(LineData(:,1),LineData(:,2),Ga1+1i.*(LineData(:,5)+Ba1),BusNumber,BusNumber);
%y_i0=sparse(LineData(:,1),LineData(:,2),((1-LineData(:,9))./LineData(:,9).^2).*(1./(LineData(:,3)+1i*LineData(:,4))),BusNumber,BusNumber);
%y_j0
y_j0 = sparse(LineData(:,2),LineData(:,1),Ga2+1i.*(LineData(:,5)+Ba2),BusNumber,BusNumber);
%y_j0=sparse(LineData(:,1),LineData(:,2),((LineData(:,9)-1)./LineData(:,9)).*(1./(LineData(:,3)+1i*LineData(:,4))),BusNumber,BusNumber);        %计算线路功率
Sij_ji=diag(U1)*(conj(y_ij).*(diag(conj(U1))*ones(BusNumber,BusNumber)-ones(BusNumber,BusNumber)*diag(conj(U1))));%先计算Zij部分的功率
Sij_ji=sparse(Sij_ji);
Sij_ji=Sij_ji+diag(U1)*(repmat(conj(U1),1,BusNumber).*conj(y_i0));%加上yi0的功率   
Sij_ji=sparse(Sij_ji);
Sij_ji=Sij_ji+diag(U1)*(repmat(conj(U1),1,BusNumber).*conj(y_j0)); %加上yj0的功率
Sij_ji=sparse(Sij_ji);%得到总的线路功率

5. 可视化潮流计算结果

为了更加直观地展示潮流计算的结果，可以使用图形化方式对各节点的电压幅值、功率流等数据进行可视化。

%%  输出结果disp('节点电压幅值');disp(full(abs(U1))); disp('节点电压角度');disp(full(angle(U1))*180/3.141592654);disp('节点电压弧度');disp(full(angle(U1))); disp('平衡节点功率');disp(S_balance); disp('线路功率');disp(Sij_ji);                                       %显示线路功率
for i = 1:BusNumberfor j = 1:iSH(j,i)=Sij_ji(i,j)+Sij_ji(j,i);end
end
disp('线路损耗');disp(SH);                                          %显示线路损耗