Continuation Power Flow does not work correctly

First of all, I am not familiar with Matlab coding, just trying learn as I go...

I want to switch the continuation parameter from lambda to voltage when the power flow fails to converge with lamba. Then once the voltage passes the point of the nose in the PV curve, I want to change the step back to lambda again. The PV curve is not turning and not in nose curve shape. Can someone help what is wrong in my code?

clc;

clear all;

close all;

% Y-bus matrix

Ybus = [-1j*20, 0, 0, 1j*20;

0, -1j*10, 1j*10, 0;

0, 1j*10, -1j*50, 1j*40;

1j*20, 0, 1j*40, -1j*60];

% Parameters

tolerance = 1e-6; % Increased precision for convergence

max_iter_corrector = 10; % Maximum iterations for corrector step

sigma_lambda = 0.1; % Initial step size for lambda continuation

sigma_voltage = 0.025; % Step size for voltage continuation

voltage_continuation = false; % Flag to indicate switch to voltage continuation

% Initial conditions

d1 = 0;

d2 = 0.0375;

d3 = -0.0125;

d4 = -0.025;

V1 = 1.05;

V2 = 1.02;

V3 = 1;

V4 = 1;

x0 = [d2; d3; d4; V3; V4];

lambda = 0; % Load parameter

lambda_values = [];

V4_values = [];

K = [0.5; 0; -1; 0; -0.8]; % Load increase direction

ek_lambda = [0, 0, 0, 0, 0, 1]; % Extra equation for lambda tracking

ek_voltage = [0, 0, 0, 0, -1, 0]; % Extra equation for voltage tracking

ek = ek_lambda; % Initially follow lambda continuation

iteration_limit = 500; % Safety limit for iterations

iteration_count = 0;

while iteration_count < iteration_limit

iteration_count = iteration_count + 1;

x = x0;

P4 = -1 * (1 + lambda); % Increase P4 by lambda

Q4 = -0.8 * (1 + lambda); % Increase Q4 by lambda

b = [0.5; 0; P4; 0; Q4];

% Define power flow functions

p2 = 10 * V2 * V3 * cos(d2 - d3 - pi/2);

p3 = 10 * V3 * V2 * cos(d3 - d2 - pi/2) + 40 * V3 * V4 * cos(d3 - d4 - pi/2);

p4 = 20 * V4 * V1 * cos(d4 - d1 - pi/2) + 40 * V4 * V3 * cos(d4 - d3 - pi/2);

q3 = 10 * V3 * V2 * sin(d3 - d2 - pi/2) + 50 * V3 * V3 * sin(pi/2) + 40 * V3 * V4 * sin(d3 - d4 - pi/2);

q4 = 20 * V4 * V1 * sin(d4 - d1 - pi/2) + 40 * V4 * V3 * sin(d4 - d3 - pi/2) + 60 * V4 * V4 * sin(pi/2);

hx = [p2; p3; p4; q3; q4];

%Jacobian matrix

% p2 = 10 * V2 * V3 * cos(d2 - d3 - pi/2);

J11 = -10 * V2 * V3 * sin(d2 - d3 - pi/2); % derivative p2/d2

J12 = 10 * V2 * V3 * sin(d2 - d3 - pi/2); % derivative p2/d3

J13 = 0; % derivative p2/d4

J14 = 10 * V2 * cos(d2 - d3 - pi/2); % derivative p2/V3

J15 = 0; % derivative p2/V4

% p3 = 10 * V3 * V2 * cos(d3 - d2 - pi/2) + 40 * V3 * V4 * cos(d3 - d4 - pi/2);

J21 = 10 * V3 * V2 * sin(d3 - d2 - pi/2); % derivative p3/d2

J22 = -10 * V3 * V2 * sin(d3 - d2 - pi/2) - 40 * V3 * V4 * sin(d3 - d4 - pi/2); % derivative p3/d3

J23 = 40 * V3 * V4 * sin(d3 - d4 - pi/2); % derivative p3/d4

J24 = 10 * V2 * cos(d3 - d2 - pi/2) + 40 * V4 * cos(d3 - d4 - pi/2); % derivative p3/V3

J25 = 40 * V3 * cos(d3 - d4 - pi/2); % derivative p3/V4

% p4 = 20 * V4 * V1 * cos(d4 - d1 - pi/2) + 40 * V4 * V3 * cos(d4 - d3 - pi/2);

J31 = 0; % derivative p4/d2

J32 = 40 * V3 * sin(d4 - d3 - pi/2); % derivative p4/d3

J33 = -20 * V4 * V1 * sin(d4 - d1 - pi/2) - 40 * V4 * V3 * sin(d4 - d3 - pi/2); % derivative p4/d4

J34 = 40 * V4 * cos(d4 - d3 - pi/2); % derivative p4/V3

J35 = 20 * V1 * cos(d4 - d1 - pi/2) + 40 * V3 * cos(d4 - d3 - pi/2); % derivative p4/V4

% q3 = 10 * V3 * V2 * sin(d3 - d2 - pi/2) + 50 * V3 * V3 * sin(pi/2) + 40 * V3 * V4 * sin(d3 - d4 - pi/2);

J41 = -10 * V3 * V2 * cos(d3 - d2 - pi/2); % derivative q3/d2

J42 = 10 * V3 * V2 * cos(d3 - d2 - pi/2) + 40 * V3 * V4 * cos(d3 - d4 - pi/2); % derivative q3/d3

J43 = -40 * V3 * V4 * cos(d3 - d4 - pi/2); % derivative q3/d4

J44 = 10 * V2 * sin(d3 - d2 - pi/2) + 100 * V3 * sin(pi/2) + 40 * V4 * sin(d3 - d4 - pi/2); % derivative q3/V3

J45 = 40 * V3 * sin(d3 - d4 - pi/2); % derivative q3/V4

% q4 = 20 * V4 * V1 * sin(d4 - d1 - pi/2) + 40 * V4 * V3 * sin(d4 - d3 - pi/2) + 60 * V4 * V4 * sin(pi/2);

J51 = 0; % derivative q4/d2

J52 = -40 * V4 * V3 * cos(d4 - d3 - pi/2); % derivative q4/d3

J53 = 20 * V4 * V1 * cos(d4 - d1 - pi/2) + 40 * V4 * V3 * cos(d4 - d3 - pi/2); % derivative q4/d4

J54 = 40 * V4 * sin(d4 - d3 - pi/2); % derivative q4/V3

J55 = 20 * V1 * sin(d4 - d1 - pi/2) + 40 * V3 * sin(d4 - d3 - pi/2) + 120 * V4 * sin(pi/2); % derivative q4/V4

J = [J11, J12, J13, J14, J15;

J21, J22, J23, J24, J25;

J31, J32, J33, J34, J35;

J41, J42, J43, J44, J45;

J51, J52, J53, J54, J55];

% Predictor Step

try

x_predicted = [x; lambda] + sigma_lambda * inv([J, K; ek]) * [0; 0; 0; 0; 0; 1];

catch

if ~voltage_continuation

disp('Power flow failed to converge. Switching to voltage continuation.');

voltage_continuation = true;

ek = ek_voltage; % Switch ek for voltage continuation

sigma_lambda = sigma_voltage; % Reduce step size

continue;

else

break;

end

end

% Corrector Step using Newton-Raphson

x_corrected = x_predicted(1:end-1);

lambda_corrected = x_predicted(end);

V_predicted = x_corrected(end);

for iter_corrector = 1:max_iter_corrector

P4_corrected = -1 * (1 + lambda_corrected);

Q4_corrected = -0.8 * (1 + lambda_corrected);

b_corrected = [0.5; 0; P4_corrected; 0; Q4_corrected];

mismatch_corrected = b_corrected - hx;

if voltage_continuation

mismatch_corrected = [mismatch_corrected; V_predicted - x_corrected(end)];

end

if max(abs(mismatch_corrected)) < tolerance

break;

end

correction = [J, K; ek] \ [mismatch_corrected; 0];

x_corrected = x_corrected + correction(1:end-1);

lambda_corrected = lambda_corrected + correction(end);

end

x0 = x_corrected;

lambda = lambda_corrected;

lambda_values = [lambda_values; lambda];

V4_values = [V4_values; x_corrected(end)];

if voltage_continuation && V4_values(end) < 0.55

disp('Switching back to lambda continuation.');

voltage_continuation = false;

ek = ek_lambda; % Switch ek back for lambda continuation

end

end

figure;

plot(lambda_values, V4_values, '-o');

xlabel('Power (pu)');

ylabel('Voltage V4 (pu)');

title('PV Curve for Bus 4');

grid on;

Answer 1

[MULI]

Hi @Hyunmin,

I understand that you are trying to implement a continuation power flow method that switches between lambda and voltage as continuation parameters.

The following issues might be affecting your results:

• The power flow equations (hx) and Jacobian matrix (J) are not updated inside the correction loop.They should be recalculated in each Newton-Raphson iteration using the current values.
• The condition to switch back to lambda continuation is based on a fixed voltage threshold (V4 < 0.55). Instead, you can consider switching back when the voltage starts increasing, which typically indicates the PV curve’s turning point.
• The PV curve appears flat or incomplete because the system equations are not dynamically updated during corrections.

The following corrections may help you achieve the desired results:

• Move the hx and J calculations inside the corrector loop.
• Use a condition that checks for increasing voltage instead of a fixed threshold.
• You can also add error handling for matrix inversion to avoid crashes.