Why unifrnd(lb,ub,[nPop,D]) gives error?

Why the following code gives error?

lb=[0 0 0 0]; ub=[10 10 pi pi]; nPop=30; D=4;
unifrnd(lb,ub,[nPop,D])

Error using unifrnd

Size information is inconsistent.

[Rik]

Deleted comments can be found in the Google cache (the archived version will require modification of the attributes, as the show/hide comments toggle doesn't execute properly).

Answer 1

[Ameer Hamza]

This cannot be done in single call to unifrnd. You need a for-loop

lb = [0 0 0 0]; ub=[10 10 pi pi]; nPop=30;
D = numel(lb);

M = rand(nPop, D);
for i = 1:D
    M(:,i) = unifrnd(lb(i),ub(i),[nPop,1]);
end

[Sadiq Akbar]

Thank you very much dear Ameer Hamza. But now it gives another error. Can you help me to run the following algorithm for my fitness function given below?

function [xmin,fmin,histout] = GQPSO(fun,D,nPop,lb,ub,maxit,maxeval)

% INPUT:
%   fun     : function handle for optimization
%   D       : problem dimension (number of variables)
%   nPop    : number of particles in the swarm
%   lb      : lower bound constrain
%   ub      : upper bound constrain
%   maxit   : max number of iterations
%   maxeval : max number of function evaluations

% OUTPUT:
%   xmin    : best solution found
%   fmin    : function value at the best solution, f(xmin)
%   histout : record of function evaluations and fitness value by iteration

% EXAMPLE:

% fun = @griewankfcn;
% D = 30;
% nPop = 50;
% lb = -600;
% ub = 600;
% maxit = 1000;
% maxeval = 10000*D;

% [xmin,fmin,histout] = GQPSO(fun,D,nPop,lb,ub,maxit,maxeval);

% OR DIRECTLY:

% [xmin,fmin,histout] = GQPSO(@griewankfcn,30,50,-600,600,1000,10000*30);

% QPSO parameters:
w1 = 0.5;
w2 = 1.0;

c1 = 1.5;
c2 = 1.5;
nPop=30; % By me
% Initializing solution
% x = unifrnd(lb,ub,[nPop,D]);  %This line gave error, so asked question on
 %Mathworks site and see below
% Correction by Ameer Hamza sab from Mathworks site
x = rand(nPop, D);
for i = 1:D
    x(:,i) = unifrnd(lb(i),ub(i),[nPop,1]);
end

% Evaluate initial population
pbest = x;

histout = zeros(maxit,2);

f_x = feval(fun,x);

fval = nPop;

f_pbest = f_x;

[~,g] = min(f_pbest);
gbest = pbest(g,:);
f_gbest = f_pbest(g);

it = 1;

histout(it,1) = fval;
histout(it,2) = f_gbest;

while it < maxit && fval < maxeval

    alpha = (w2 - w1) * (maxit - it)/maxit + w1;
    mbest = sum(pbest)/nPop;

    for i = 1:nPop

        fi = abs(randn(1,D)); % Gaussian QPSO (Coelho, 2010)

        p = (c1*fi.*pbest(i,:) + c2*(1-fi).*gbest)/(c1 + c2);

        u = abs(randn(1,D));  % Gaussian QPSO (Coelho, 2010)
        
        b = alpha*abs(x(i,:) - mbest);
        v = log(1./u);

        if rand < 0.5
            x(i,:) = p + b .* v;
        else
            x(i,:) = p - b .* v;
        end
        
        % Keeping bounds
        x(i,:) = max(x(i,:),lb);
        x(i,:) = min(x(i,:),ub);

        f_x(i) = fun(x(i,:));
        
        fval = fval + 1;

        if f_x(i) < f_pbest(i)
            pbest(i,:) = x(i,:);
            f_pbest(i) = f_x(i);
        end

        if f_pbest(i) < f_gbest
            gbest = pbest(i,:);
            f_gbest = f_pbest(i);
        end

    end
    
    it = it + 1;

    histout(it,1) = fval;
    histout(it,2) = f_gbest;

end

xmin = gbest; 
fmin = f_gbest;

histout(it+1:end,:) = [];

figure
semilogy(histout(:,1),histout(:,2))
title('Gaussian Q-PSO')
xlabel('Function evaluations')
ylabel('Fitness')
grid on

end

And my fitness function is:

function out=myfitness(b)
u=[1 2 3 4];
v=[b(1) b(2) b(3) b(4)];
out1=u-v
out=sum(out1);
end

When I type the followiong in command window, it does not run

>>[xmin,fmin,histout] = GQPSO(@myfitness,4,30,[0 0 0 0],[10 10 pi pi],100,10000*4)

[Walter Roberson]

What error do you encounter?

[Sadiq Akbar]

Thank you very much for your response. When I type the above in command window, I encounter the following error:

>> [xmin,fmin,histout] = GQPSO(@myfitness,4,30,[0 0 0 0],[10 10 pi pi],100,10000*4)
Attempted to access f_pbest(2); index out of bounds because numel(f_pbest)=1.
Error in GQPSO (line 99)
        if f_x(i) < f_pbest(i) 

[Walter Roberson]

Your fitness function is expected to receive a 2D array in which the number of rows corresponds to the population size, and the columns correspond to variables, and it is expected to return a vector which is the same size as the population size. Instead your fitness function expects a vector of four elements as input and returns a scalar.

function out=myfitness(b)
  u = [1 2 3 4];
  out1 = u-b;             %needs R2016b or later
  out = sum(out1.^2,2);
end

[Walter Roberson]

No error when I tried with the GQPSO you posted in https://www.mathworks.com/matlabcentral/answers/702837-why-unifrnd-lb-ub-npop-d-gives-error#comment_1228747 and the myfitness I posted in https://www.mathworks.com/matlabcentral/answers/702837-why-unifrnd-lb-ub-npop-d-gives-error#comment_1229062

[Walter Roberson]

I see no problem in R2020b with the original code.

Could you confirm you are using R2016b and not R2016a ? The fitness function I posted needs R2016b or later.

[xmin,fmin,histout] = GQPSO(@myfitness,4,30,[0 0 0 0],[10 10 pi pi],100,10000*4)

xmin = 1×4

0.9174 1.6567 3.1416 3.1416

fmin = 0.8816

histout = 100×2

30.0000 9.6990 60.0000 5.5803 90.0000 1.2189 120.0000 1.2189 150.0000 1.2189 180.0000 1.2189 210.0000 1.2189 240.0000 1.2189 270.0000 1.2189 300.0000 1.2189

●

function out=myfitness(b)
  u = [1 2 3 4];
  out1 = u-b;             %needs R2016b or later
  out = sum(out1.^2,2);
end

function [xmin,fmin,histout] = GQPSO(fun,D,nPop,lb,ub,maxit,maxeval)

% INPUT:
%   fun     : function handle for optimization
%   D       : problem dimension (number of variables)
%   nPop    : number of particles in the swarm
%   lb      : lower bound constrain
%   ub      : upper bound constrain
%   maxit   : max number of iterations
%   maxeval : max number of function evaluations

% OUTPUT:
%   xmin    : best solution found
%   fmin    : function value at the best solution, f(xmin)
%   histout : record of function evaluations and fitness value by iteration

% EXAMPLE:

% fun = @griewankfcn;
% D = 30;
% nPop = 50;
% lb = -600;
% ub = 600;
% maxit = 1000;
% maxeval = 10000*D;

% [xmin,fmin,histout] = GQPSO(fun,D,nPop,lb,ub,maxit,maxeval);

% OR DIRECTLY:

% [xmin,fmin,histout] = GQPSO(@griewankfcn,30,50,-600,600,1000,10000*30);

% QPSO parameters:
w1 = 0.5;
w2 = 1.0;

c1 = 1.5;
c2 = 1.5;
nPop=30; % By me
% Initializing solution
% x = unifrnd(lb,ub,[nPop,D]);  %This line gave error, so asked question on
 %Mathworks site and see below
% Correction by Ameer Hamza sab from Mathworks site
x = rand(nPop, D);
for i = 1:D
    x(:,i) = unifrnd(lb(i),ub(i),[nPop,1]);
end

% Evaluate initial population
pbest = x;

histout = zeros(maxit,2);

f_x = feval(fun,x);

fval = nPop;

f_pbest = f_x;

[~,g] = min(f_pbest);
gbest = pbest(g,:);
f_gbest = f_pbest(g);

it = 1;

histout(it,1) = fval;
histout(it,2) = f_gbest;

while it < maxit && fval < maxeval

    alpha = (w2 - w1) * (maxit - it)/maxit + w1;
    mbest = sum(pbest)/nPop;

    for i = 1:nPop

        fi = abs(randn(1,D)); % Gaussian QPSO (Coelho, 2010)

        p = (c1*fi.*pbest(i,:) + c2*(1-fi).*gbest)/(c1 + c2);

        u = abs(randn(1,D));  % Gaussian QPSO (Coelho, 2010)
        
        b = alpha*abs(x(i,:) - mbest);
        v = log(1./u);

        if rand < 0.5
            x(i,:) = p + b .* v;
        else
            x(i,:) = p - b .* v;
        end
        
        % Keeping bounds
        x(i,:) = max(x(i,:),lb);
        x(i,:) = min(x(i,:),ub);

        f_x(i) = fun(x(i,:));
        
        fval = fval + 1;

        if f_x(i) < f_pbest(i)
            pbest(i,:) = x(i,:);
            f_pbest(i) = f_x(i);
        end

        if f_pbest(i) < f_gbest
            gbest = pbest(i,:);
            f_gbest = f_pbest(i);
        end

    end
    
    it = it + 1;

    histout(it,1) = fval;
    histout(it,2) = f_gbest;

end

xmin = gbest; 
fmin = f_gbest;

histout(it+1:end,:) = [];

figure
semilogy(histout(:,1),histout(:,2))
title('Gaussian Q-PSO')
xlabel('Function evaluations')
ylabel('Fitness')
grid on

end

[Walter Roberson]

With regard to your new fitness function obj : As I posted above:

Your fitness function is expected to receive a 2D array in which the number of rows corresponds to the population size, and the columns correspond to variables, and it is expected to return a vector which is the same size as the population size.

Now let us examine your code:

u=[2 7 50*th 85*th]; % angles of all four vectors

Okay, that is a vector of 4 elements.

%%%%%%%%%%%%%%%
% Swapping vector b. Taken help from Mathworks site
%%%%%%%%%%%%%%%
[~, ix] = sort(u); 

That is going to return the vector [3 4 1 2] into ix.

[~, ix1(ix)] = sort(b);  

There you are asking to sort a 2D array with 30 rows and 4 columns. The sort information for that is going to require 30 x 4 output locations, but ix1(ix) would be a vector of 4 locations. Also, sort() by default operates along the first non-singleton dimension, which would be the first dimension for a 30 x 4 array, but it seems more likely that you want to sort along the rows, which would be the second dimension, seeing as you expect an array of length 4.

For this step I would suggest instead

[~, ix1] = sort(b, 2);
ix1 = ix1(:,ix);

which would give you a 2D array, in which each row talks about the sort order of b for its corresponding row

Then

b = b(ix1);

is a problem as it does not do per-row reordering.

Now.. what really is that code trying to do??

Let us experiment:

th=pi/180
u=[2 7 50*th 85*th]
[~, ix] = sort(u)
b = [4 1 3 2]
[~, ix1] = sort(b) 
ix1 = ix1(:,ix)
new_b = b(ix1)

b_2 = [0    17   -11    -6]
[~, ix1_2] = sort(b_2) 
ix1_2 = ix1_2(:,ix)
new_b_2 = b_2(ix1)

ix comes out as [3 4 1 2]

new_b comes out as [3 4 1 2]

new_b_2 comes out as [0 17 -11 -6]

t = sort([b;b_2],2)
t(:,ix)

comes out as

     3     4     1     2
     0    17   -11    -6

so it looks to me as if it wants to sort b per-row and then re-order according to the sorting of u. In that case I would suggest

th=pi/180;
u=[2 7 50*th 85*th];
[~, ix] = sort(u);
sort_b = sort(b,2);
b = sort_b(:,ix);

[Walter Roberson]

As I posted twice before:

Your fitness function is expected to receive a 2D array in which the number of rows corresponds to the population size, and the columns correspond to variables, and it is expected to return a vector which is the same size as the population size.

Your obj is returning a scalar, not a vector of length 30.

It is not possible to change the algorithm to account for the objective function returning the wrong results.

[Walter Roberson]

            xe(1,k)=xe(1,k)+b(i)*exp(-1i*(k-1)*pi*cos(b(P+i)));

How many rows does b have? How many rows of b does that line access?

Where is your code producing one output for each row in b?

[Walter Roberson]

Read your code!!

x = rand(nPop, D);
for i = 1:D
    x(:,i) = unifrnd(lb(i),ub(i),[nPop,1]);
end

So x is the full matrix, not one row of the matrix.

f_x = feval(fun,x);

And right there you call your function on the entire x matrix.

Your objective function must return one value for each population member.

[Walter Roberson]

% Initialize the population/solutions
for i=1:n,
  Sol(i,:)=Lb+(Ub-Lb).*rand(1,d);
f_x(i)=obj(Sol(i,:));  % My fitness function
end

[Walter Roberson]

remove

f_x = feval(fun,x);

insert

for i = n
  f_x(i) = feval(fun,x(i, :)) 
end

[Walter Roberson]

for i = nPop
  f_x(i) = feval(fun,x(i, :)) 
end

[Walter Roberson]

for i = 1:nPop
  f_x(i) = feval(fun,x(i, :)) 
end

[Walter Roberson]

I am not sure on the index out of range. I will look at it after I sleep as I am not at my desk.

[Walter Roberson]

Code seems to be working for me.

Invoke obj_driver or myfit_driver as applicable.