针对包含变压器分接的IEEE13节点案例研究持续潮流计算方法，并采用NR方法解决潮流问题（Matlab代码实现）

科研_G.E.M.

于 2024-09-19 14:45:56 发布

阅读量947

点赞数 29

文章标签： matlab 开发语言

本文链接：https://blog.csdn.net/qq_46009436/article/details/142360088

版权

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

⛳️赠与读者

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

⛳️赠与读者

👨‍💻做科研，涉及到一个深在的思想系统，需要科研者逻辑缜密，踏实认真，但是不能只是努力，很多时候借力比努力更重要，然后还要有仰望星空的创新点和启发点。当哲学课上老师问你什么是科学，什么是电的时候，不要觉得这些问题搞笑。哲学是科学之母，哲学就是追究终极问题，寻找那些不言自明只有小孩子会问的但是你却回答不出来的问题。建议读者按目录次序逐一浏览，免得骤然跌入幽暗的迷宫找不到来时的路，它不足为你揭示全部问题的答案，但若能让人胸中升起一朵朵疑云，也未尝不会酿成晚霞斑斓的别一番景致，万一它居然给你带来了一场精神世界的苦雨，那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。

或许，雨过云收，神驰的天地更清朗.......🔎🔎🔎

💥1 概述

一种针对包含变压器分接的IEEE13节点案例的持续潮流计算方法是使用Newton-Raphson (NR) 方法来解决潮流问题。NR方法是一种常用的电力系统潮流计算方法，通过迭代求解节点电压和支路功率的平衡方程来得到系统的潮流分布。

在NR方法中，首先需要建立系统的节点导纳矩阵和功率注入矩阵。然后根据节点电压的初始值，通过迭代计算得到各节点的电压和功率。在包含变压器分接的情况下，需要考虑变压器变比对系统的影响，可以通过将变压器的等值电路合并到系统节点导纳矩阵中来考虑变压器的影响。

具体步骤如下：
1. 建立系统的节点导纳矩阵和功率注入矩阵，包括变压器分接的等值电路。
2. 初始化节点电压的值，可以选择平衡状态下的电压作为初始值。
3. 根据节点电压的初始值，迭代计算节点电压和支路功率，直到收敛。
4. 检查收敛条件，如果满足收敛条件，则得到系统的潮流分布，否则返回第3步进行迭代计算。

通过以上步骤，可以利用NR方法解决包含变压器分接的IEEE13节点案例的潮流问题，并得到系统的潮流分布结果。这种方法可以有效地分析系统的电压稳定性和功率平衡情况，为电力系统的运行和规划提供重要参考。

📚2 运行结果

部分代码：

% bus to study
bts = 13 ;
%% setting parameters

% tolerance error for voltage and angle for consecutive values
tol_volt = 1e-3;
tol_ang = 1e-3;
tap_include = 1;
maxiters = 25;
%% getting data
[bus_data, branch_data] = data_extract();
baseMVA = 100;

%% Y bus formation
% calling function for ybus calculation
Ybus = y_bus_calculation(bus_data, branch_data, tap_include);
% getting G and B from Y bus
G = real(Ybus);
B = imag(Ybus);
% converting to polar
[Theta Y_mag]=cart2pol(G,B);

%% initializing solution

% finding types of bus
nbus = length(Ybus); % total bus number

% find indexing of PV bus
PV_bus = find(bus_data.data(:,3)==2);

% find indexing of swing bus
Swing_bus =find(bus_data.data(:,3)==3);

% find indexing of PQ bus
PQ_bus = find(bus_data.data(:,3)==0);

% flat start; initialization
% assign 1 to all voltage
Voltage= ones(nbus,1);
% fix voltage of swing bus and PV bus to given value
% column 4 of bus data contains bus voltage, you can also use column 11
Voltage(Swing_bus) = bus_data.data(Swing_bus,4);
Voltage(PV_bus) = bus_data.data(PV_bus,4);

% assign 0 to all bus angles
Delta= zeros(nbus,1);
% fix bus angle of Swing bus to given value
% column 5 of bus data contains angle
Delta(Swing_bus) = bus_data.data(Swing_bus,5) * pi/180;
Non_swing_bus = union(PQ_bus, PV_bus);

% schedule power calculation, For Pv only P, for PQ, both P and Q
% P_sch_ori is P for non-swing bus, Q_schi_ori is Q for PQ bus only
[P_sch_ori, Q_sch_ori, P_sch_all, Q_sch_all] = schedule_power_calc(bus_data, baseMVA,Swing_bus, PV_bus);

% getting K vetor
K = [P_sch_ori; Q_sch_ori];

% continuation parameters
sigma = 0.1; % for first stage
lambda = 0; % for initialization
Voltage_history1 = [];
Delta_history1 = [];
% getting voltage and angle for starting iteration
[Voltage, Delta,~] = NRPF(tol_volt,tol_ang,Voltage,Delta,Swing_bus,PQ_bus,PV_bus,nbus...
,Y_mag,Theta,bus_data,G,B,baseMVA,bts,lambda,maxiters);

%% part I %%%% lambda as continuation parameter

Voltage_collects1 = [];
lambda_collects1 = [];
Voltage_collects1(end+1) = Voltage(bts);
lambda_collects1(end+1) = lambda;
Voltage_history1(:,end+1) = Voltage;
Delta_history1(:,end+1) = Delta;

while (1)
%%%%%%%%% Predictor %%%%%

% calculating jacobian using last voltage and angle
%listing parameters for jacobian calculation
jacobian_params.nbus = nbus;
jacobian_params.G = G;
jacobian_params.B = B;
jacobian_params.Theta = Theta;
jacobian_params.Y_mag = Y_mag;
jacobian_params.PV_bus = PV_bus;
jacobian_params.Swing_bus = Swing_bus;
jacobian_params.PQ_bus = PQ_bus;
jacobian_params.Voltage = Voltage;
jacobian_params.Delta = Delta;

% calling jacobian function
[J11,J12,J21,J22] = jacobian_calc_original(jacobian_params);
J = [J11 J12;J21 J22];
% finding ek
ek = [zeros(1,length(J)) 1];
% finding augmented jacobian
J_aug = [J K;ek];

% getting del_delta, del_voltage and del_lambda using crout based solver
del_x = crout_solver(J_aug,ek');

% getting del_delta, del_voltage and del_lambda
del_delta = del_x([1:nbus-length(Swing_bus)]);
del_voltage = del_x(nbus-length(Swing_bus)+1:end-1);
del_lambda = del_x(end);

% updating
Delta(Non_swing_bus) = Delta(Non_swing_bus) + sigma* del_delta;
Voltage(PQ_bus) = Voltage(PQ_bus) + sigma* del_voltage;
lambda = lambda + sigma*del_lambda;

%%%%%% Corrector %%%%%%%%

% updating studies load value for NR
[Voltage, Delta, iter] = NRPF(tol_volt,tol_ang,Voltage,Delta,Swing_bus,PQ_bus,PV_bus,nbus...
,Y_mag,Theta,bus_data,G,B,baseMVA,bts,lambda,maxiters);

if (iter == maxiters)
% assign previous value of V and lambda as the solution diverges
Voltage = Voltage_history1(:,end);
lambda = lambda_collects1(end);
Voltage_collects1(end+1) = Voltage(bts);
lambda_collects1(end+1) = lambda;
Delta = Delta_history1(:,end);
break;
else
Voltage_history1(:,end+1) = Voltage;
Delta_history1(:,end+1) = Delta;
Voltage_collects1(end+1) = Voltage(bts);
lambda_collects1(end+1) = lambda;
end
end

%% Part II %%%% Voltage as continuation parameter

sigma = 0.005; % for 12 and 13 use 0.0005
factor = 0.75; % factor where the continuation parameter switches to lambda for 12 and 13 use 0.75
lambda = lambda_collects1(end);
Voltage = Voltage_history1(:,end);
Delta = Delta_history1(:,end);
Voltage_history2 = [];
Delta_history2 = [];
lambda_collects2 = [lambda_collects1(end)];
Voltage_collects2 = [Voltage_collects1(end)];

%%% Predictor step %%%
lambda_2start = lambda;
iteration = 1;
while (lambda>lambda_2start*factor) && iteration < 100
iteration = iteration+1;
%finding jacobian
%listing parameters for jacobian calculation
jacobian_params.nbus = nbus;
jacobian_params.G = G;
jacobian_params.B = B;
jacobian_params.Theta = Theta;
jacobian_params.Y_mag = Y_mag;
jacobian_params.PV_bus = PV_bus;
jacobian_params.Swing_bus = Swing_bus;
jacobian_params.PQ_bus = PQ_bus;
jacobian_params.Voltage = Voltage;
jacobian_params.Delta = Delta;

% calling jacobian function
[J11,J12,J21,J22] = jacobian_calc_original(jacobian_params);
J = [J11 J12;J21 J22];

% finding ek
ek = [zeros(1,length(J)) 1];

% finding ekv i.e. last row vector of augmented jacobian for part II
ekv = [zeros(1,length(J)) 0];

% finding index of voltage continuation variable in augmented jacobian
modified_Qindex = find(PQ_bus(:)==bts);
v_index = length(Non_swing_bus) + modified_Qindex;
ekv(v_index) = -1;

% augmented jacobian
J_aug = [J K;ekv];

% getting del_delta, del_voltage and del_lambda using crout based solver
del_x = crout_solver(J_aug,ek');
% getting del_delta, del_voltage and del_lambda
del_delta = del_x([1:nbus-length(Swing_bus)]);
del_voltage = del_x(nbus-length(Swing_bus)+1:end-1);
del_lambda = del_x(end);

% updating
Delta(Non_swing_bus) = Delta(Non_swing_bus) + sigma* del_delta;
Voltage(PQ_bus) = Voltage(PQ_bus) + sigma* del_voltage;
lambda = lambda + sigma*del_lambda;

%%%% Corrector Step %%%%%%%%%%%%%%%%%%%%%%