Progressively slower performance and task constraints using matrix and loops

1 view (last 30 days)
Hello Everyone,
Hope you are all well.
I have a question regarding loops becoming progressivley slower with each iteration. Below is my script to execute Newton-Raphson method on power flow equations. I calculate the Jacobian and inverse Jacobian matrix to find the correction factors. This script runs very slow on the second iteration and will run out of memory when given long enough. The first pass runs very fast.
My question is: Why is it slowing down? The stored matrix is quite small in dimension. I even clear the numerical value after each pass. I have tested this with single variable and double variable versions as well and it is the same thing, each pass gets much slower. In fact, I sometimes need to force stop the script or crash matlab.
Does anyone have experience with this problem?
Best Wishes,
%% Newton-Raphson Method (Three Variable)
clear
clc
%sympref('FloatingPointOutput',true);
% ********************************************* %
% Enter variable names
syms x1 x2 x3 % True Variable Names: (delta1, delta2, V1)
tol = 1e-3; % required tolerance
maxIter = 2; % maximum iterations (start with 3)
% Define Functions (f(x) = x1*x2 - y = 0 [don't write = 0])
g1 = 5*sin(x2-x1) + 10*x3*sin(-x1) - 2
g2 = 5*x3*sin(x2-x1) + 2.5*sin(x2) - 1.2
g3 = 5*x3*cos(x2-x1) + 10*x3*cos(x1) - 15*x3^2 - 0.5
% Initial Parameters
x1o = -0.1; % initial x1 guess
x2o = 0.01; % initial x2 guess
x3o = 0.95; % initial x2 guess
% ********************************************* %
% Find Jacobian Matrix
J = [diff(g1,x1) diff(g1,x2) diff(g1,x3);
diff(g2,x1) diff(g2,x2) diff(g2,x3);
diff(g3,x1) diff(g3,x2) diff(g3,x3)]
% Initialize while loop
tolMet = 0; % tolerance met flag
x1 = x1o; % substitue initial guess
x2 = x2o; % substitue initial guess
x3 = x3o; % substitue initial guess
clear x10 x2o x3o
disp(' ')
disp('Initial Evaluation:')
val1 = subs(g1) % evaluate guess
val2 = subs(g2) % evaluate guess
val3 = subs(g3) % evaluate guess
count = 0;
while (~tolMet) && (count < maxIter)
% check if functions are solved to required tolerance
if (abs(val1) < tol) && (abs(val2) < tol) && (abs(val3) < tol)
tolMet = 1;
disp(' ')
disp('********************')
disp("Tolerance Met After " + num2str(count) + " Iterations")
x1
x2
x3
end
% if tolerance not met, correct values
if ~tolMet
disp(' ')
disp("***** Iteration " + num2str(count+1) + " *****")
Jval = subs(J) % evaluate Jacobian Matrix
Jinv = inv(Jval) % calculate inv Jacobian Matrix
correction = -Jinv*[val1; val2; val3] % calculate correction factor
x1 = x1 + correction(1,1) % new x1 value
x2 = x2 + correction(2,1) % new x2 value
x3 = x3 + correction(3,1) % new x3 value
count = count + 1; % update iteration counter
val1 = subs(g1) % evaluate new guess
val2 = subs(g2) % evaluate new guess
val3 = subs(g3) % evaluate new guess
clear Jval Jinv
end
disp('********************')
end
if count == maxIter
disp(' !!!! Maximum Iterations Reached !!!!')
end

Answers (1)

Alan Stevens
Alan Stevens on 13 Apr 2022
Edited: Alan Stevens on 13 Apr 2022
With these particular functions it's easy to precalculate the Jacobian function and do everything numerically (and quickly). For example:
tol = 1e-6;
maxIter = 10;
% Functions
g = @(x1,x2,x3) [5*sin(x2-x1) + 10*x3*sin(-x1) - 2;
5*x3*sin(x2-x1) + 2.5*sin(x2) - 1.2;
5*x3*cos(x2-x1) + 10*x3*cos(x1) - 15*x3^2 - 0.5];
J = @(x1,x2,x3) [-5*cos(x2-x1)-10*x3*cos(x1), 5*cos(x2-x1), -10*sin(x1);
-5*x3*cos(x2-x1), 2.5*cos(x2)+5*x3*cos(x2-x1), 5*sin(x2-x1);
5*x3*sin(x2-x1)-10*x3*sin(x1), -5*x3*sin(x2-x1), 5*cos(x2-x1)-30*x3+10*cos(x1)];
% Initial values
xo = [-0.1; 0.01; 0.95];
% Loop
tolMet = 0;
count = 0;
while (~tolMet) && (count<maxIter)
x = xo - J(xo(1),xo(2),xo(3))\g(xo(1),xo(2),xo(3)); %%%%%%%%%%%%%%%%%% Note J\g rather than Jinv*g
if abs(x(1)-xo(1))<tol && abs(x(2)-xo(2))<tol && abs(x(3)-xo(3))<tol
tolMet = 1;
end
xo = x;
count = count + 1;
end
disp(x)
-0.1047 0.0973 0.9547
disp(count)
4
Note that it is not necessary to calculate the inverse of the Jacobian explicitly - use Matlab's backslash operator instead.

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by