Progressively slower performance and task constraints using matrix and loops

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

Answer 1

[Alan Stevens]

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.

[Kyle Robertson]

Hi Alan,

Thank you for taking the time for this.

Regards,

Kyle