Why there is an error in the line 132?

% %this program for comparison between MGA & DE

clear;

%load input-output data

tt=cputime;

load ui.dat; %load input

u=ui(:,1);

mx=2;

y(1:mx+1)=0;

for t=mx+1:500

% for i=1:500

% e(i)=(rand-0.5)*0.002;

% end

% y(t)=0.797*y(t-1)-0.237*y(t-3)+0.441*u(t-1)+0.105*y(t-4)*u(t-4)...

% +0.333*u(t-3)*u(t-5);%+e(t);

y(t)=0.3*u(t-1)*u(t-2)-0.8*y(t-2)*u(t-1)*u(t-1)-0.1*u(t-1)*u(t-1)*u(t-2);%+e(t);

end

y=y';%load output

nodata=500;

noy=2; nou=2;

nl=3; start=max(nou,noy)+1;

m=nou+noy;

mphi=[]; %form regressor for ARMAX model

for i=start:nodata

phi=[];

if nl==1

for j=1:noy

phi=[phi y(i-j)];

end

for k=0:nou

phi=[phi u(i-k)];

end

end

if nl==2

for j=1:noy

phi=[phi y(i-j)];

end

for k=0:nou

phi=[phi u(i-k)];

end

for j=1:noy

for l=j:noy

phi=[phi y(i-j)*y(i-l)];

end

for k=1:nou

phi=[phi y(i-j)*u(i-k)];

end

end

for k=1:nou

for n=k:nou

phi=[phi u(i-k)*u(i-n)];

end

end

end

if nl==3

for j=1:noy

phi=[phi y(i-j)];

end

for k=0:nou

phi=[phi u(i-k)];

end

for j=1:noy

for l=j:noy

phi=[phi y(i-j)*y(i-l)];

end

for k=1:nou

phi=[phi y(i-j)*u(i-k)];

end

end

for k=1:nou

for n=k:nou

phi=[phi u(i-k)*u(i-n)];

end

end

for j=1:noy

for l=j:noy

for p=l:noy

phi=[phi y(i-j)*y(i-l)*y(i-p)];

end

for k=1:nou

phi=[phi y(i-j)*y(i-l)*u(i-k)];

end

end

for k=1:nou

for n=k:nou

phi=[phi y(i-j)*u(i-k)*u(i-n)];

end

end

end

for k=1:nou

for n=k:nou

for q=n:nou

phi=[phi u(i-k)*u(i-n)*u(i-q)];

end

end

end

end

mphi=[mphi;phi]; %end regressor

end

[temp m]=size(mphi);

%DIFFERENTIAL EVOLUTION

%CONTROL PARAMETERS

D=m; %dimension of problem

NP=30; %Population size

F=0.9; %differentiation constant

CR=0.5; %crossover rate

GEN=100; %generations

L=0.0; %low boundary constraint 1

H=1.0; %high boundary constraint 1

% L2=4.1; %low boundary constraint 2

% H2=5.8; %high boundary constraint 2

%**********************************%

% Algorithm Variables

%**********************************%

X=zeros(D,1); %trial vctor X1

% X2=zeros(D,1); %trial vctor X2

Pop=zeros(D,NP); %population X1

% Pop2=zeros(D,NP); %population X1

Fit=zeros(1,NP); %fitness population

iBest=1; %index of best solution

r=zeros(3,1); %randomly selected indices

%**********************************%

% Create Population

%**********************************%

%Initialize random number of generator

%rand('state',sum(100*clock));

for j=1:NP %Initialize each individual

Pop(:,j)=L+(H-L)*rand(D,1); %within b.constraint 1

% Pop2(:,j)=L2+(H2-L2)*rand(D,1); %within b.constraint 2

[Fit(1,j), tetaDE]=selectmodel2(Pop(:,j),NP,D,mphi,start,nodata,y,u); %evaluate fitness 1

end

%**********************************%

% Optimization

%**********************************%

fbestDE=[];

for g=1:GEN %for each generation

for j=1:NP %for each individual

%choose 3 random individuals from population

%mutually different and different from j

r(1)=floor(rand()*NP)+1;

while r(1)==j

r(1)=floor(rand()*NP)+1;

end

r(2)=floor(rand()*NP)+1;

while (r(2)==r(1))||(r(2)==j)

r(2)=floor(rand()*NP)+1;

end

r(3)=floor(rand()*NP)+1;

while (r(3)==r(2))||(r(3)==r(1))||(r(3)==j)

r(3)=floor(rand()*NP)+1;

end

%create trial individual

%in which at least one parameter is changed

Rnd=floor(rand()*D)+1;

for i=1:D

if (rand()<CR)||(Rnd==i)

X(i)=Pop(i,r(3))+F*(Pop(i,r(1))-Pop(i,r(2)));

% X2(i)=Pop2(i,r(3))+F*(Pop2(i,r(1))-Pop2(i,r(2)));

else

X(i)=Pop(i,j);

% X2(i)=Pop2(i,j);

end

end

%verify b.constraint

for i=1:D

if (X(i)<L)||(X(i)>H)

X(i)=L+(H-L)*rand();

end

% if (X2(i)<L2)||(X2(i)>H2)

% X2(i)=L2+(H2-L2)*rand();

% end

end

%select the best individual

%between trial and current ones

%calculate fitness of trial vector

[f, tetaDE]=selectmodel2(X,NP,D,mphi,start,nodata,y,u);

%if trial is better or equal than current

% if f<=Fit(j)

% Pop1(:,j)=X1; %replace current by trial

% Pop2(:,j)=X2;

% Fit(j)=f;

% %if trial is better than the best

% if f<=Fit(iBest)

% iBest=j; %update the best's index

% end

% end

if f<=Fit(j)

Pop(:,j)=X; %replace current by trial

% Pop2(:,j)=X2;

Fit(j)=f;

%if trial is better than the best

if f<=Fit(iBest)

iBest=j; %update the best's index

end

end

end

fbestDE=[fbestDE Fit(iBest)];

numgenDE(g)=g;

end

% %MGA

%CONTROL PARAMETERS

evol=100;%input('# of evolution >:');

popsize=30;%input('population size >:');

pcr=0.8;%input('crossover probability=');

pmut=0.01;%input('mutation probability=');

chr=initbp(popsize,m); %initialize population

chrom=chr;

mseGA=[];

for t=1:evol %start evolution

gen1=[];

for i=1:popsize %model structure chrom's=1

mmphi=[];

krom=chrom(i,:);

for b=1:m

index=krom(:,b);

if index==1

mmphi=[mmphi mphi(:,b)];

end

end

[N, g]=size(mmphi);

yt=y(start:nodata);ut=u(start:nodata);

y2=sum(yt.*yt);

teta=pinv(mmphi'*mmphi)*mmphi'*yt; %find teta by LSE

num=0;

for c=1:g

tetahat=teta(c,:);

if tetahat<0.001

num=num+0.1; %penalty to added noise

end

end

yhat=mmphi*teta;

error=sum((yt-yhat).*(yt-yhat))+num;

gen=[krom error];

gen1=[gen1;gen];

pm=max(gen1(:,m+1));

fitness=pm-gen1(:,m+1);

% bestftn=pm-fitness-num;

gene1=[gen1 fitness];

end

generation=sortrows(gene1,[m+1]);

% mseGA=[mseGA gene1(:,m+1)];

sse=generation(:,m+1);

bestchr=0.0; secondchr=0.0;

for i=1:popsize

sume2=sse(i,:);

if sume2<1.0

bestchr=bestchr+1; %good chromosomes

else

secondchr=secondchr+1; %2nd class chromosomes

end

end

sumerror2(t)=sum(gene1(:,m+1));

fitnv=gene1(:,m+2); %assign fitness

sumfit(t)=sum(fitnv);

%newchr=selsus(fitnv,popsize);

newchr=seltour(fitnv,popsize,2);

selchrom=chrom(newchr,:);

%selchrom=generation(:,1:m);

mutchrom=mutbint(selchrom(1:bestchr,:),pmut); %mutate good chrom

xxchrom=xovsp(selchrom(bestchr+1:popsize-1,:),pcr); %xover 2nd chrom

xmchrom=mutbint(xxchrom,pmut);

newnew=initbp(1,m); %worst chrom,reproduce

chrom=[mutchrom;xmchrom;newnew];

numgen(t)=t;

end %end of evolution

% newchromosome=generation(:,1:m);

% [r f]=size(newchromosome);

% Out=[];MSE1=[];CX=[];CX2=[];MSE21=[];YHT2=[];

% for i=1:r

% selectedchrom=newchromosome(i,:);

% numbits=sum(newchromosome(i,1:m));

% cx=nou+noy+nl+numbits;

% CX=[CX;cx];

% pi=[];

% for r=1:m

% ind=selectedchrom(:,r);

% if ind==1

% pi=[pi mphi(:,r)];

% else

% pi=[pi mphi(:,r)*0];

% end

% end

% tet=pinv(pi'*pi)*pi'*yt;

% noterm=newchromosome(i,1:m);

% nterm=noterm';

% yht=pi*tet;

% msedata=mse(yt-yht);

% at=fix(log10(msedata))-1;a1=round(100*(msedata/10^at))/100;

% msedata=a1*10^at;

% MSE=msedata;

% MSE1=[MSE1;MSE];

% Out=[nterm tet];

% end

% teta=[tetaDE tetaGA]

BestMSE=[Fit(iBest) gene1(1,m+1)]

% BestTerm=Pop(:,iBest)

plot(numgenDE,fbestDE,'m-',numgen,sumfit,'g--')

xlabel('Generation');ylabel('Fitness Value');

RUN_TIME=cputime-tt;