How to optimise this code?

Nested loops suck...how can I optimise this code?

function [a, v, x] = model1(x_s, x_a, OB, rho_0, K_a, K_r, K_v)
  ss = 0.01; % Step Size
  tmax = 60;
  c = 0; % Row index ID (for identification only)
  k = 1:length(OB); % Vectorisation
  x = [x_s(1); x_s(2)]; % Initial position
  v = [0; 0]; % Initial velocity
  a = [0; 0]; % Initial acceleration
  dUa = [1; 1];
  while ((dUa(1)^2 + dUa(2)^2) > 1e-10) && (c <= (tmax/ss)+1)
      dUr = 0; % Reset the force
      c = c + 1; % Add one to the row count (first iteration c = 1)
      dUa = K_a*(x(:,c) - x_a); % Calculate attractive force
      % Calculate distance from each obstacle
      rho_r = sqrt( (x(1,c)-OB(1,k)).^2 + (x(2,c)-OB(2,k)).^2 );
      for numobs = 1:length(OB) % Loop for the number of obstacles created
          if rho_r(numobs) > rho_0 % If outwith the obstacle radius
              dUrobs = 0; % No force experienced
          else % Otherwise
              dUrobs = K_r * ((1/rho_r(numobs)) + (1/rho_0)) * (OB(:,numobs)-x(:,c))/(rho_r(numobs)^3);
          end
          dUr = dUr + dUrobs; % Combine force from each obstacle to total
      end
      a(:,c+1) = -dUa - dUr - K_v * v(:,c); % Calculate acceleration
      v(:,c+1) = v(:,c) + a(:,c+1) * ss; % Calculate new velocity
      x(:,c+1) = x(:,c) + v(:,c+1) * ss; % Calculate new position
  end

Answer 1

[Jan]

Pre-allocation!

n = ceil(c <= (tmax/ss) + 1);
x = zeros(2, n);
x(:, 1) = x_s(1:2); % Initial position
v = zeros(2, n); % Initial velocity
a = zeros(2, n); % Initial acceleration
...

The power operation is expensive, so do not compute it repeatedly: Omit zeros in the sum:

rho_r3 = rho_r ^ 3;
...
    for numobs = 1:length(OB) % Loop for the number of obstacles created
        if rho_r(numobs) <= rho_0 % If outwith the obstacle radius
            dUr = dUr + K_r * (1/rho_r(numobs) + (1/rho_0)) * ...
                              (OB(:,numobs)-x(:,c)) / rho_r3(numobs);
        end
    end

Crop the unused memory at the end:

x = x(:, 1:c);
v = v(:, 1:c);
a = a(:, 1:c);

[J]

Awesome! But is there a way to completely remove the nested for loop and instead replace it with a vectorised equivalent?

[Jan]

Maybe. But it is not worth to spend hours in optimizing the code, when you save some seconds of runtime only. You can try this instead of the "for numobs" loop:

K = (rho_r <= rho_0);
dUR = sum(K_r .* (1 ./ rho_r(K) +  1 ./ rho_0) .* ...
          (OB(:, K) - x(:, c)) ./ rho_r3(K));

I cannot debug this due to the missing input data. Maybe some vectors must be transposed. But there migt be further problems with the sizes. Anyway, the underlying method is cleared.