%a,b,c為KP,KI,KD的變化范圍；popsize為群體規模；lchrom為表示（KP,KI,KD）的染色體串長度；maxgen為進化的最大代數
%pcross為交叉概率；pmutation為變異概率；min為上一代最小適應度；max為最大適應度；minpp為適應度最小個體的序號；
%maxpp為適應度最大個體的序號；
function [sys,x0,str,ts]=gen1(t,x,u,flag)
global a b oldpop newpop popsize lchrom maxgen pcross pmutation fitu;
global min max minpp maxpp k lc gen s;
% %******************************參數初始化************************************
a=[0 0 0];
b=[2 2 1];
popsize=6;
lchrom=18;
lc=lchrom/3;
maxgen=100;
pcross=[0.6 0.6 0.6];
pmutation=0.02;
gen=0;
k=1;
% % %********************************群體的初始化*******************************%
for s=1:popsize
oldpop{s}.chrom=randint(1,18);
oldpop{s}.dvalue=decode(oldpop{s}.chrom);
kp=oldpop{s}.dvalue(1);
ki=oldpop{s}.dvalue(2);
kd=oldpop{s}.dvalue(3);
oldpop{s}.fitness=0.0001;                %計算每個個體的適應度%
oldpop{s}.parent1=0;
oldpop{s}.parent2=0;
oldpop{s}.xsite=0;
end
statistics(oldpop);

while(gen<maxgen)
mate1=select;
mate2=select;
%*******************判斷是否交叉***************************
for i=1:3
    pp=rand;
%************************確定交叉位置************************
   if (pp<pcross(i))
       jcross(i)=randint(1,1,[lc*(i-1)+1,lc*i]);
    else
       jcross(i)=lchrom/(4-i);
   end
%**********************進行交叉與變異操作*********************
   if(jcross(i) ~= lchrom/(4-i))
       for j=(lc*(i-1)+1):jcross(i)
     newpop{k}.chrom(j)=mutation(oldpop{mate1}.chrom(j));
     newpop{k+1}.chrom(j)=mutation(oldpop{mate2}.chrom(j));
       end
       for j=jcross(i):lc*i
     newpop{k}.chrom(j)=mutation(oldpop{mate2}.chrom(j));
     newpop{k+1}.chrom(j)=mutation(oldpop{mate1}.chrom(j));   
       end
   else
       for j=(lc*(i-1)+1):lc*i
     newpop{k}.chrom(j)=mutation(oldpop{mate1}.chrom(j));
     newpop{k+1}.chrom(j)=mutation(oldpop{mate2}.chrom(j));
       end
    end
 %*************************************************************
end
 %*************************************************************
newpop{k}.dvalue=decode(newpop{k}.chrom);
newpop{k}.fitness=objfunc(u(2));
newpop{k}.parent1=mate1;
newpop{k}.parent2=mate2;
newpop{k}.xsite=jcross;
sim('test1')
disp('kp=');
kp=newpop{k}.dvalue(1);
disp('ki=')
ki=newpop{k}.dvalue(2)
disp('kd=')
kd=newpop{k}.dvalue(3)

newpop{k+1}.dvalue=decode(newpop{k+1}.chrom);
newpop{k+1}.parent1=mate1;
newpop{k+1}.parent2=mate2;
newpop{k+1}.xsite=jcross;
newpop{k+1}.fitness=objfunc(u(2));
disp('kp=');
kp= newpop{k+1}.dvalue(1);
disp('ki=')
ki=newpop{k+1}.dvalue(2);
disp('kd=')
kd=newpop{k+1}.dvalue(3);
k=k+2;
%***************************************************************
   if mod(k,51)==0
      k=1;
     gen=gen+1
     oldmax=max;
     oldmaxpp=maxpp;
     statistics(newpop);
    if(max<oldmax)
        newpop{minpp}.chrom=oldpop{oldmaxpp}.chrom;
        newpop{minpp}.dvalue=oldpop{oldmaxpp}.dvalue;
        newpop{minpp}.fitness=oldpop{oldmaxpp}.fitness;
        statistics(newpop);
    end
    for s=1:popsize
    oldpop{s}.chrom=newpop{s}.chrom;
    end
   end
%***************************************************************   
end



%*********************譯碼*************************************
function dval=decode(x)
global a b;
m=[0 0 0];
for i=1:3
    for j=1:6
    m(i)=m(i)+x((i-1)*6+j)*(2^(6-j));
    end
end
for i=1:3
    dval(i)=a(i)+m(i)*(b(i)-a(i))/(2^6-1);
end
end
%**************************************************************


%***************************適應度函數**************************
function fit=objfunc(c)
if c~=0
    fit=1/c;
else
    fit=0;
end
end
%***************************************************************

%***************************選擇操作****************************
function rs=select
global oldpop popsize;
index(1)=randint(1,1,[1,popsize]);
index(2)=randint(1,1,[1,popsize]);
if (oldpop{index(1)}.fitness > oldpop{index(2)}.fitness)
    rs=index(1);
else
    rs=index(2);
end    
end
%***************************************************************

%*******************群體適應度計算******************************
function statistics(pop)
global newpop oldpop popsize min max maxpp minpp;
sumfitness=pop{1}.fitness;
min=pop{1}.fitness;
max=pop{1}.fitness;
maxpp=1;
minpp=1;
for j=1:popsize
sumfitness=sumfitness+pop{j}.fitness;
    if (pop{j}.fitness>max)
        max=pop{j}.fitness;
        maxpp=j;
    end
    if(pop{j}.fitness<min)
        min=pop{j}.fitness;
        minpp=j;
    end
end
avg=sumfitness/popsize;
end
%************************************************************

%*****************************變異操作***********************
function vmut=mutation(ch)
global pmutation
pp1=rand;
if (pp1<pmutation) 
   ch=~ch;
end
   vmut=ch;
end
%************************************************************