% fnnga7.m										% 程序名
% 采用遺傳算法訓練具有7個隸屬函數的FNN控制器權值的程序
%
clear
pausetime = 0.1;
load squarin.mat                % 系統跟蹤曲線
load plantv.mat                 % 取被控對象神經網絡的權矢量PW10, PW20, PB10, PB20

% 賦初始值
N=20;                           % 確定種群數
k1=0.000024015;
k2=1/32768;
k3=50;k4=0.5;
fn =7;  ns = fn*4+fn*fn;		% 由7個模糊標記計算個體中的所含實數編碼的總數目ns
ns
in=k2*inpi(1:100:4000);         % 取40個數進行訓練
T=in';
[mn,Q]=size(T);

A2=[zeros(N,Q)];
y=[zeros(N,Q)];

e=T-0;ec=zeros([1,Q]);P=k3*[e;ec];

SSE=[zeros(1,N)];
WB1=[zeros(N,ns)*2];
NWB=rands(N,ns)*2;

m=[zeros(1,N)];	mm=1;
min=1000000;

G =100;
for j=1:G,          % 開始直至G代的遺傳進化
  SS=0;
  for i= 1:N,

  W1=[NWB(i,1) 0;NWB(i,2) 0;NWB(i,3) 0;NWB(i,4) 0;NWB(i,5) 0;...
	0 NWB(i,6);0 NWB(i,7);0 NWB(i,8);0 NWB(i,9);0 NWB(i,10);...
	0 NWB(i,11);0 NWB(i,12);0 NWB(i,13);0 NWB(i,14)];
  B1=[NWB(i,15);NWB(i,16);NWB(i,17);NWB(i,18);NWB(i,19);NWB(i,20);...
	NWB(i,21);NWB(i,22);NWB(i,23);NWB(i,24);NWB(i,25);...
	NWB(i,26);NWB(i,27);NWB(i,28)];
  W2=[NWB(i,29) NWB(i,30) NWB(i,31) NWB(i,32) NWB(i,33) NWB(i,34)...
       NWB(i,35) NWB(i,36) NWB(i,37) NWB(i,38) NWB(i,39) NWB(i,40)...
       NWB(i,41) NWB(i,42) NWB(i,43) NWB(i,44) NWB(i,45) NWB(i,46)...
       NWB(i,47) NWB(i,48) NWB(i,49) NWB(i,50) NWB(i,51) NWB(i,52)...
       NWB(i,53) NWB(i,54) NWB(i,55) NWB(i,56) NWB(i,57) NWB(i,58)...
       NWB(i,59) NWB(i,60) NWB(i,61) NWB(i,62) NWB(i,63) NWB(i,64)...
       NWB(i,65) NWB(i,66) NWB(i,67) NWB(i,68) NWB(i,69) NWB(i,70)...
       NWB(i,71) NWB(i,72) NWB(i,73) NWB(i,74) NWB(i,75) NWB(i,76)...
       NWB(i,77)];
  B2=0;

%calculate the fnnet's output

e(1)=0.05;ec(1)=0.05;
P=k3*[e(1);ec(1)];

for k=1:Q,
 A10 = gause(W1*P,B1);          % 計算FNN神經網絡的輸出A20
 A11 = multil(A10);
 A20(k) = purelin(W2*A11);

 A21(k)=A20(k)*k4;              % 檢查所計算出的控制量是否正向飽和
if A21(k)>1. 
  cm(k)=0.7869;	              	%32768*0.000024015=0.7869
 	else if A21(k)<-1.          % 檢查所計算出的控制量是否負向飽和
      cm(k)=-0.7869;            % 飽和時輸出最大飽和值
   	else 
      cm(k)=A21(k)*0.7869;      % 未飽和則輸出計算值A21
    end
end

xp(:,k)=[cm(k);T(k)];           % 構造被控對象網絡輸入矢量
AA1=tansig(PW10*xp(:,k)+PB10);
y(i,k)=purelin(PW20*AA1+PB20);  % 計算系統輸出值
E(k) = T(k)-y(i,k);                % 計算系統輸出誤差值
ec(k+1)=e(k);                   % 計算系統下一時刻的誤差值
e(k+1)=E(k)*k3;
P=[e(k+1);ec(k+1)];
end 

F(i)=sqrt(sumsqr(E));           % 求每一組適應度值

%if reach the err_goal then quit

       if F(i)<0.0001,           % 當適配值小于給定的0.00001，結束進化
           final_err = F(i); 
           save result W1 W2 B1 final_err;
           break,end
       SS=SS+F(i);               % 求N個種群適應值和
       end

% 復制操作(reproduction)

AF=SS/N;				% 求N個種群適應值的平均值
n=1;nn=N;min1=0;
  for i=1:N,				% 尋找記錄N個種群中小于平均適應值AF的個體
  dif=AF-F(i);
    if dif>=0,          
      m(n)=i;				% 記錄個數n和所處序號i
      n=n+1; 
      min1=F(i)-min;		% 尋找本代最優(最小)適應值
      if min1<=0, min = F(i); nn=i;end  % 并把其序號保存到nn, 其值保存到min中
    end
   end

plot(F(1:N)),				% 繪出本代中N個適應值各自的值
pause2(pausetime),
 minim=min				% 顯示本代中最小適應值
 nn
 AF 

  if j==mm*5,nj=j			% 每5代顯示一次網絡輸出跟蹤目標矢量的情況
  plot([T(1:k)' y(nn,:)']),pause2(pausetime),mm=mm+1;
  end

% 交換操作(crossover)

 WB1(N,:)=NWB(nn,:);				% 將本代最優個體放置N號位
   if n-1>3*N/4,s=3*N/4;				% 只取少于s = 3*N/4的個體數進行淘汰
   else s=n-1;
   end

  for i=1:n-2,		   % 將記錄中低于平均適應值的個體從N-1號位起向前順序排放
  WB1((N-i),:)=NWB(m(i),:);
  end

  for i=1:N-s,				% 用隨機值補齊被淘汰的個體
  WB1(i,:)=rands(1,ns)*2;
  end

% k=fix(rand*10);

  if j<=20, intwb=fix(WB1*10)/10;   % 20代前取權值個位及四位小數中的前一位
  else if j>20 &j<=50,intwb=fix(WB1*100)/100;  % 21~50代取權值小數中的前兩位
  else intwb=fix(WB1*1000)/1000;  %51代以后取權值小數中的前三位
   end
  end 

poiwb=WB1-intwb;				% 求權值在不同代中的小數部分
NWB(N,:)=WB1(N,:);			% 將本代最優個體直接進化到下一代
 
  for i=1:N/2,  % 將1到N/2與N-1到N/2 - 1的個體之間的權值的小數進行交換
  NWB(i,:)=intwb(N-i,:)+poiwb(i,:);
  end
  for i=1+N/2:N,
  NWB(i-1,:)=intwb(i-N/4,:)+poiwb(i,:);
  end

% 變異操作(mutation)
if j>10 & j<=30, jk=35;					% 選擇變異的頻率
 else if j>30 & j<=50, jk=30;       %從每40個個體中變異一次增長到…
      else if j>50 & j<=80, jk=25;  %每25到個個體中變異一次
          else jk=40; end
     end
end

 for i=5:jk:N-1,
 k=fix(rand*10);				% 隨機取一個9以內的整數
 mu1=num2str(NWB(i,:));   % 將所選定的個體中的權值數碼中與隨機…
 mu2=strrep(mu1,'k','k-1');   % 數k相同的數變為k - 1
 hh=sscanf(mu2,'%f')';
 [in,im]=size(hh);
 if im<ns, NWB(i,:)=[hh 0.];
 else NWB(i,:)=hh;end
 end
end

wbg7=NWB(1,:);
save fnnwb7 wbg7