Error in user-supplied fitness function evaluation - Genetic Algorithm

Hello everyone. I am currently using the ga function in order to find the best parameters for a discrete GPC controller. However, more often than not, the execution of this algorithm ends up in error.

First, let me show the fitness function that I am currently using:

function J = simplanta(individuo)

Np = individuo(1);
Nc = individuo(2);
rw = individuo(3);

ts=1/2000;  % sampling time
t=0:ts:1;   % sim time
f0=50;      
r=1+1.5.*sin(2*pi*f0*t); % reference signal
Ldata=length(t); % # of time points

%Discrete Model
Ad =[0.4134 0.4536  0; -0.2404 0.7876 0; -0.4366 0.4224 0.7738];
Bd =[0.1331; 0.4528; 0.1273];
Cd=[0 0 1];

% Specified Gain Calc

n=length(Ad);       %order%
Omega = 2*pi*f0*ts;
gamma = ( 2*cos(Omega) ) + 1;
E=[0 1 0; 0 0 1; 1 -gamma gamma]; 
sigma = [0;0;1];

% prediction matrices
A=[Ad zeros(n,3); sigma*Cd*Ad E];
B=[Bd; sigma*Cd*Bd];
C=[zeros(1,3) sigma' ];

% MPC
nau = length(A);
Rb = rw*eye(Nc,Nc); 

% F e Phi
F=zeros(Np,nau);
for i=1:Np
   F(i,:)=C*(A^i);  
end
Phi = zeros(Np, Nc);
%1st column:
for linha=1:Np
    Phi(linha,1) = C*(A^(linha-1))*B;
end
%other column
for col =2: Nc
    Phi(col:Np, col) = Phi(1:Np-(col-1),1);
end

PTP  = Phi'*Phi;
if (Nc == 1)
   vetor10 = [1];
else
   vetor10 = [ 1 zeros(1,Nc-1)]; 
end
Hessiana=inv(PTP + Rb);
Kbs = -1*vetor10*Hessiana*(Phi')*F;

% control sim

x = zeros(n,Ldata+1); %
y = zeros(1,Ldata); %
erro =zeros(1,Ldata); %
u =zeros(1,Ldata); %
xk_1=zeros(n,1); %x(k-1)
xk_2=zeros(n,1); %x(k-2)
xk_3=zeros(n,1); %x(k-3)
uk_1=0; %u(k-1)
uk_2=0; %u(k-2)
uk_3=0; %u(k-3)
ek_1=0; %e(k-1)
ek_2=0; %e(k-2)

% closed loop sim
for k=1:Ldata
    xk = x(:,k); %
    rk = r(k);
    yk = Cd*xk; %
    ek = rk-yk;  %
    qk = xk_3+(-gamma*xk_2)+(gamma*xk_1)-xk;
    xbarra= [qk;ek_2;ek_1;ek];
    vk=Kbs*xbarra;
    uk= uk_3 + gamma*(uk_1-uk_2)-vk;
    
    % 
    x(:,k+1) = (Ad*xk)+(Bd*uk);
    % 
    erro(k) = ek;
    u(k) = uk;
    y(k) = yk;
    %
    xk_3 = xk_2;
    xk_2 = xk_1;
    xk_1 = xk;
    uk_3 = uk_2;
    uk_2 = uk_1;
    uk_1 = uk;
    ek_2 = ek_1;
    ek_1 = ek;
end

%Fit func

up = max(abs(u))/1000; %peak control law
erro_rmse = rmse(r,y); % RMSE 

% settling time 
max_e = max(erro); %
valor_acom = 0.02*max_e; % 
ind_acom = find(erro >= -valor_acom & erro <= valor_acom, 1); %
t_ae = t(ind_acom); %

J=0.2*up + 0.7*erro_rmse + 0.1*t_ae;

With this fitness function I am calling the ga function with the following parameters:

A_des = [0 0 -1; -1 1 0; 1 0 0; 0 -1 0];
b_des = [0 0 100 -2]'; 
intcon = [1 2]; os
lb = [2 2 0]; 
ub = [100 99 60]; 

%Declaração das opções de funcionamento do AG
options = optimoptions("ga","CreationFcn","gacreationuniformint",...
    "Display","iter","MutationFcn","mutationpower",...
    "SelectionFcn","selectionroulette",...
    "PlotFcn",["gaplotdistance", "gaplotselection","gaplotstopping",...
    "gaplotbestf","gaplotbestindiv"]);

[best_ind, cost] = ga(@simplanta, 3, A_des, b_des, [], [], lb, ub, [],...
    intcon, options)
param = infoPlant(best_ind)

However, sometimes (and to be quite frankly more often than not) MATLAB shows me this error during the execution of this code:

Unable to perform assignment because the size of the left side is 1-by-1 and the
size of the right side is 1-by-0.

Error in fcnvectorizer (line 19)
            y(i,:) = feval(fun,(pop(i,:)));

Error in gaminlppenaltyfcn

Error in gapenalty

Error in stepGA (line 63)
    nextScore = FitnessFcn(pop);

Error in galincon (line 86)
        [score,population,state] = stepGA(score,population,options,state,GenomeLength,FitnessFcn);

Error in gapenalty

Error in ga (line 412)
    [x,fval,exitFlag,output,population,scores] = gapenalty(FitnessFcn,nvars,...

Error in ga_implementation (line 60)
[best_ind, cost] = ga(@simplanta, 3, A_des, b_des, [], [], lb, ub, [],...

Caused by:
    Failure in user-supplied fitness function evaluation. GA cannot continue.

I am quite frankly at loss as of why that would happen, since it sometimes runs smoothly, so it seems that it isn't necessarily a code syntax error of my part. Or perhaps it IS and I am just not able to detect it.

Any help and input would be quite appreciated

[Aquatris]

For certain inputs to 'simplanta' function, the output becomes empty. That is your problem.

The reason is you cannot find the index in 'erro'' variable that is between certain value in line::

ind_acom = find(erro >= -valor_acom & erro <= valor_acom, 1); %

Since I do not know the optimization problem you are trying to solve, I cannot comment further.

One simple solution is, to check if ind_acom is empty and assign J to be really high value like 1e9 or something to indicate it is not a valid solution.

Here is one such case:

J = simplanta([37 2 4.3535])

ind_acom =

  1x0 empty double row vector


J =

  1x0 empty double row vector

function J = simplanta(individuo)

Np = individuo(1);
Nc = individuo(2);
rw = individuo(3);

ts=1/2000;  % sampling time
t=0:ts:1;   % sim time
f0=50;      
r=1+1.5.*sin(2*pi*f0*t); % reference signal
Ldata=length(t); % # of time points

%Discrete Model
Ad =[0.4134 0.4536  0; -0.2404 0.7876 0; -0.4366 0.4224 0.7738];
Bd =[0.1331; 0.4528; 0.1273];
Cd=[0 0 1];

% Specified Gain Calc

n=length(Ad);       %order%
Omega = 2*pi*f0*ts;
gamma = ( 2*cos(Omega) ) + 1;
E=[0 1 0; 0 0 1; 1 -gamma gamma]; 
sigma = [0;0;1];

% prediction matrices
A=[Ad zeros(n,3); sigma*Cd*Ad E];
B=[Bd; sigma*Cd*Bd];
C=[zeros(1,3) sigma' ];

% MPC
nau = length(A);
Rb = rw*eye(Nc,Nc); 

% F e Phi
F=zeros(Np,nau);
for i=1:Np
   F(i,:)=C*(A^i);  
end
Phi = zeros(Np, Nc);
%1st column:
for linha=1:Np
    Phi(linha,1) = C*(A^(linha-1))*B;
end
%other column
for col =2: Nc
    Phi(col:Np, col) = Phi(1:Np-(col-1),1);
end

PTP  = Phi'*Phi;
if (Nc == 1)
   vetor10 = [1];
else
   vetor10 = [ 1 zeros(1,Nc-1)]; 
end
Hessiana=inv(PTP + Rb);
Kbs = -1*vetor10*Hessiana*(Phi')*F;

% control sim

x = zeros(n,Ldata+1); %
y = zeros(1,Ldata); %
erro =zeros(1,Ldata); %
u =zeros(1,Ldata); %
xk_1=zeros(n,1); %x(k-1)
xk_2=zeros(n,1); %x(k-2)
xk_3=zeros(n,1); %x(k-3)
uk_1=0; %u(k-1)
uk_2=0; %u(k-2)
uk_3=0; %u(k-3)
ek_1=0; %e(k-1)
ek_2=0; %e(k-2)

% closed loop sim
for k=1:Ldata
    
    xk = x(:,k); %
    rk = r(k);
    yk = Cd*xk; %
    ek = rk-yk;  %
    qk = xk_3+(-gamma*xk_2)+(gamma*xk_1)-xk;
    xbarra= [qk;ek_2;ek_1;ek];
    vk=Kbs*xbarra;
    uk= uk_3 + gamma*(uk_1-uk_2)-vk;
    
    % 
    x(:,k+1) = (Ad*xk)+(Bd*uk);
    % 
    erro(k) = ek;
    u(k) = uk;
    y(k) = yk;
    %
    xk_3 = xk_2;
    xk_2 = xk_1;
    xk_1 = xk;
    uk_3 = uk_2;
    uk_2 = uk_1;
    uk_1 = uk;
    ek_2 = ek_1;
    ek_1 = ek;
    
end

%Fit func

up = max(abs(u))/1000; %peak control law
erro_rmse = sum(sqrt((r-y).^2)); % RMSE 

% settling time 
max_e = max(erro); %
valor_acom = 0.02*max_e; % 
ind_acom = find(erro >= -valor_acom & erro <= valor_acom, 1)
t_ae = t(ind_acom); %

J=0.2*up + 0.7*erro_rmse + 0.1*t_ae;
end

[Matheus Caramalac]

Thank you so much! That was exactly what I needed.

For some reason, since I selected "gaplotbestf" as one of my plots, I just assumed that every single member of my population was able to find a proper solution to the fitness function. I was so tired when I made this that I totally forgot about this possibility.

I wish I could accept your comment as the proper answer, @Aquatris!

Answer 1

[Ganesh]

Hi @Matheus Caramalac,

The issue here seems to be arising due to the following piece of code:

ind_acom = find(erro >= -valor_acom & erro <= valor_acom, 1);

In some cases, it seems that there is no index found satisfying the given conditions. And hence, it creates a 1x0 Array. Eventually during the calculations, it reflects in the end product.

I am unsure of the logic implemented in the "Genetic Algorithm", but I would suggest you to add a default value to ind_acom in case no suitable index is found.

Hope this helped you narrow down the issue!

[Matheus Caramalac]

This answer, as well as the comment made by @Aquatris really helped me out! Thank you so much, it was exactly that situation that was happening. Using a default value solved this problem.