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])
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
