Why I am getting "Subscript indices must either be real positive integers or logicals." error, while using PSO algorithm in MATLAB. How to fix this?

Hello Sir/Madam,

I am getting "Subscript indices must either be real positive integers or logicals." error while running pso algorithm in matlab. It is showing that the error exists in the 6th line of my objective function (i.e LoadFlow.m), but I have to make that variable negative for my objective function to work correctly. It is also showing error in line 53 of the pso.m file which I cannot understand.

This is the pso algorithm I used:

%% Problem Definition

CostFunction = @(x) LoadFlow(x); % Cost Function

nVar = 2; % Number of Decision Variables

VarSize = [1 nVar]; % Size of Decision Variables Matrix

VarMin = [2 1]; % Lower Bound of Variables

VarMax = [33 100]; % Upper Bound of Variables

%% PSO Parameters

MaxIt = 1000; % Maximum Number of Iterations

nPop = 100; % Population Size (Swarm Size)

% PSO Parameters

w = 1; % Inertia Weight

wdamp = 0.99; % Inertia Weight Damping Ratio

c1 = 1.5; % Personal Learning Coefficient

c2 = 2.0; % Global Learning Coefficient

% Velocity Limits

VelMax = 0.1*(VarMax-VarMin);

VelMin = -VelMax;

%% Initialization

empty_particle.Position = [];

empty_particle.Cost = [];

empty_particle.Velocity = [];

empty_particle.Best.Position = [];

empty_particle.Best.Cost = [];

particle = repmat(empty_particle, nPop, 1);

GlobalBest.Cost = inf;

for i = 1:nPop

% Initialize Position

particle(i).Position = unifrnd(VarMin, VarMax, VarSize);

% Initialize Velocity

particle(i).Velocity = zeros(VarSize);

% Evaluation

particle(i).Cost = CostFunction(particle(i).Position);

% Update Personal Best

particle(i).Best.Position = particle(i).Position;

particle(i).Best.Cost = particle(i).Cost;

% Update Global Best

if particle(i).Best.Cost<GlobalBest.Cost

GlobalBest = particle(i).Best;

end

end

BestCost = zeros(MaxIt, 1);

%% PSO Main Loop

for it = 1:MaxIt

for i = 1:nPop

% Update Velocity

particle(i).Velocity = w*particle(i).Velocity ...

+c1*rand(VarSize).*(particle(i).Best.Position-particle(i).Position) ...

+c2*rand(VarSize).*(GlobalBest.Position-particle(i).Position);

% Apply Velocity Limits

particle(i).Velocity = max(particle(i).Velocity, VelMin);

particle(i).Velocity = min(particle(i).Velocity, VelMax);

% Update Position

particle(i).Position = particle(i).Position + particle(i).Velocity;

% Velocity Mirror Effect

IsOutside = (particle(i).Position<VarMin | particle(i).Position>VarMax);

particle(i).Velocity(IsOutside) = -particle(i).Velocity(IsOutside);

% Apply Position Limits

particle(i).Position = max(particle(i).Position, VarMin);

particle(i).Position = min(particle(i).Position, VarMax);

% Evaluation

particle(i).Cost = CostFunction(particle(i).Position);

% Update Personal Best

if particle(i).Cost<particle(i).Best.Cost

particle(i).Best.Position = particle(i).Position;

particle(i).Best.Cost = particle(i).Cost;

% Update Global Best

if particle(i).Best.Cost<GlobalBest.Cost

GlobalBest = particle(i).Best;

end

end

end

BestCost(it) = GlobalBest.Cost;

disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);

w = w*wdamp;

end

BestSol = GlobalBest;

This is my objective function:

function y = LoadFlow(x)

x1 = x(1);

x2 = x(2);

LD = load('linedata33bus.m'); % loading the branch data

BD = load('Busdata33bus.m'); % loading the bus data

BD(x1,2) = -x2;

BD(x1,3) = 0;

Sbase = 200;

Vbase = 11;

Zbase = (Vbase^2)/Sbase;

LD(:,4:5) = LD(:,4:5)/Zbase;

BD(:,2:3) = BD(:,2:3)/(1000*Sbase);

N = max(BD(:,1));

Sload = complex(BD(:,2),BD(:,3));

V = ones(size(BD,1),1);

Z = complex(LD(:,4),LD(:,5));

Ibranch = zeros(size(LD,1),1);

Iter = 1000;

%THE ALGORITHM

for i = 1:Iter

% BACKWARD SWEEP

Iload = conj(Sload./V);

for j = size(LD,1):-1:1

[c e] = find(LD(:,2:3) == LD(j,3));

if size(c,1) == 1

Ibranch(LD(j,1)) = Iload(LD(j,3));

else

Ibranch(LD(j,1)) = Iload(LD(j,3)) + sum(Ibranch(LD(c,1))) - Ibranch(LD(j,1));

end

end

% FORWARD SWEEP

for j = 1:size(LD,1)

V(LD(j,3)) = V(LD(j,2)) - Ibranch(LD(j,1))*Z(j);

end

end

y = real_power_loss(Z,Zbase,Ibranch,Vbase,Sbase);

Please help me.

Answer 1

[Image Analyst]

x1 is probably either 0, negative, or a fractional (non-integer) number.

It's thoroughly discussed in the FAQ:

https://matlab.fandom.com/wiki/FAQ#.22Subscript_indices_must_either_be_real_positive_integers_or_logicals..22

[Baisakh Sarangi]

But actually i need x1 to be negative inside my objective function for that to work correctly.

Also i cannot understand the error in pso (line 53). Can you please explain what is it?

[Yongjian Feng]

You have in your code

BD = load('Busdata33bus.m'); % loading the bus data
BD(x1,2) = -x2;

Use load function to load a m file? Is it a mat file instead?

[Baisakh Sarangi]

@Yongjian Feng sir, Busdata33bus.m is a simple script file. Sir the file is as shown below:

[Yongjian Feng]

I see. It is more like a CSV file. If you load it as a matrix, then why do you need the index (x1 here) to be negative? Maybe you meant BD(x1, 2) is negative, instead of x1 being negative?

[Baisakh Sarangi]

@Yongjian Feng sir, actually BD(x1, 2) represents real power load in my code and sir I want to insert a distributed generation(dg) unit supplying real power i.e x2 (varying from 1KW to 1000KW) to the distribution system and as real power is being supplied to the system instead of being taken from the system, that's why x2 has to be negative i.e (-1KW to -1000KW). Let's say if x1 = 3, then BD(3,2) represents 90, i.e 90KW and by doing BD(3,2) = -10 (let's say x2 = 10KW), i am connecting a dg unit with that bus(i.e bus 3), which will be supplying 10KW of power to the system.