Script Algoritma Genetika Optimasi Pencarian Nilai Konstanta Sliding Mode Controller/Pengendali Modus Luncur Oleh: Ahmad Riyad Firdaus Teknik Elektro Politeknik Batam I. Program Utama clear; clc; format long; rand('state',0); %% %Simulasi Perancangan Model Pengendalian Motor DC %========================================================================== % Parameter Algoritma Genetika %========================================================================== pcross=0.8; %Probabilitas crossover pmut=0.05; %Probabilitas Mutasi L=12; %panjang bit kromosom (12 bit x2) bits=[L L L]; maxgen=100; %maksimum iterasi (menghasilkan generasi baru) siz=30; %Jumlah populasi tiap satu generasi sea=3; %Ada 2 parameter yang akan dicari bounds=[0 15;1 10;1 10]; %Batas range variabel yang akan dicari nilainya (nilai k) %========================================================================== %Parameter motor %========================================================================== Rf=120; %Rf=field resistance (ohm) Lf=20; J=1; f=0.5; K2=68.5; %========================================================================== %Ruang Status sistem %========================================================================== A=[0 1;-(Rf*f/Lf) -(Rf*J/Lf+f)]; B=[0;K2/Lf]; C=[1 0]; D=0; %========================================================================== %Inisial status dan status referensi %========================================================================== x0=[0 0]'; xr=[1 0]'; %% %========================================================================== %Mengenerate nilai k dan Sf secara random dengan populasi size %========================================================================== k_Sf=init_ga(bounds,siz,sea); k=k_Sf(:,1); c1=k_Sf(:,2); c2=k_Sf(:,3); Sf=[c1 c2]; %========================================================================== %Gain Ekivalen %========================================================================== for r=1:siz K(r,:)=-inv(Sf(r,:)*B)*(Sf(r,:)*A); end %========================================================================== %Mengubah variable phenotif dari float menjadi biner (dinamakan gen)
%========================================================================== for m=1:siz gen(m,:) = f2b(k_Sf(m,:),bounds,bits); end; %========================================================================== %Mengevaluasi semua kemungkinan nilai k yang digenerate %========================================================================== [fitness,object]=evaluasi_smc_dcmotor(K,k,A,B,Sf,x0,xr,siz); [best_fit(1),index]=max(fitness); best_gen=gen(index,:); best_k_Sf=k_Sf(index,:); [worst_fit(1),index1]=min(fitness); worst_cur_gen=gen(index1,:); worst_cur_k_Sf=k_Sf(index1,:); avg_fit(1)=0; for m=1:siz avg_fit(1)=avg_fit(1)+fitness(m); end; avg_fit(1)=avg_fit(1)/siz; fprintf(1,'%f ',best_k_Sf); fprintf('\n'); fprintf(1,'Terbaik: %f |Terburuk: %f Rata-rata:%f \n fprintf('\n');
',best_fit(1), worst_fit(1), avg_fit(1));
%% %========================================================================= %Melakukan seleksi, reproduski, crossover, mutasi dan evaluasi dari %generasi baru. Object yang memiliki nilai maksimum pada setiap generasi %akan disimpan pada best_cur_obj dan setiap k yang dipilih pada setiap %generasi akan disimpan pada best_cur_gain. Best_cur_obj pada setiap %generasi akan disimpan pada object_array %========================================================================= for i=1:maxgen gen=reproduksi(gen,fitness); gen=mate(gen); gen=xover(gen,pcross); gen=mutate(gen,pmut); for n=1:siz k_Sf(n,:)= b2f(gen(n,:),bounds,bits); end; k=k_Sf(:,1); c1=k_Sf(:,2); c2=k_Sf(:,3); Sf=[c1 c2]; for r=1:siz K(r,:)=-inv(Sf(r,:)*B)*(Sf(r,:)*A); end [fitness,object]=evaluasi_smc_dcmotor(K,k,A,B,Sf,x0,xr,siz); [best_cur_fit,index]=max(fitness); best_cur_gen=gen(index,:); best_cur_k_Sf=k_Sf(index,:); [worst_fit(i+1),index1]=min(fitness); worst_cur_gen=gen(index1,:); worst_cur_k_Sf=k_Sf(index1,:); avg_fit(i+1)=0; for m=1:siz avg_fit(i+1)=avg_fit(i+1)+fitness(m); end; avg_fit(i+1)=avg_fit(i+1)/siz; if(best_cur_fit < best_fit(i)) k_Sf(index1,:)=best_k_Sf; gen(index1,:)=best_gen; object(index1)=best_fit(i);
best_fit(i+1)=best_fit(i); elseif(best_cur_fit>=best_fit(i)) best_k_Sf=best_cur_k_Sf; best_gen=best_cur_gen; best_fit(i+1)=best_cur_fit; end; k(i+1)=best_k_Sf(1); Sf((i+1),:)=[best_k_Sf(2) best_k_Sf(3)]; fprintf(1,'Iterasi ke -%3.0f \n',i); fprintf(1,'k: %f Sf: [%f %f] ',k(i+1),Sf((i+1),:)); fprintf(1,'Best Fitness ------>%f\n ',best_fit(i+1)); fprintf('\n'); fprintf(1,'Terbaik: %f |Terburuk: %f |Rata-rata: %f\n ',best_fit(i+1),worst_fit(i+1),avg_fit(i+1)); fprintf('\n'); fprintf('\n'); end %% xx=1:maxgen+1; figure(1); plot(xx,best_fit,'k');%,xx,worst_fit,'r');%,xx,avg_fit,'g'); title('Histori Nilai Fitness Terbaik'); grid; k=best_k_Sf(1); c1=best_k_Sf(2); c2=best_k_Sf(3); Sf=[c1 c2]; K=-inv(Sf*B)*(Sf*A); %% %========================================================================== %Pengujian nilai k dan Sf apakah sesuai dengan kriteria yang diinginkan %dengan hitting time tertentu %========================================================================== u=0; x=x0; dx=0; dt=0.01; for i=1:1000 x=x+dx; xs(:,i)=x; xe=x-xr; us(i)=u; dx=(A*x+B*u)*dt; s=Sf*x-Sf*xr; ss(i)=s; %Control Law Ueq=K*x; Un=-k*inv(Sf*B)*sign(s); u=Ueq+Un; end; %% %========================================================================== %Menghitung kembali nilai hitting time %========================================================================== %a=xr(1); %Nilai referensi (merupakan batas hitting time) a=xr(1)-0.01*(xr(1)-x0(1)); %toleransi 1% dari nilai referensi t=1; stop=0; while stop==0 if xs(1,t)
end; if t==length(xs) stop=1; end; t=t+1; end; ts=t*dt; %ts pada kondisi mulai steady state %% if ts<4 %menguji apakah hitting time sesuai dengan yang diinginkan t=t+200; for m=1:100 x_steady_array(m)=xs(1,t); %nilai state dari mulai ts sampai ts10 t=t+1; ts10=t*dt; end; x_array=x_steady_array; %array steady state samapi ts10 x_steady_min=min(x_steady_array); %nilai minimum state pada kondisi steady state x_steady_max=max(x_steady_array); %nilai maximum state pada kondisi steady state CH=x_steady_max-x_steady_min; end; %========================================================================== %% disp('********************************************************************'); disp('********* Data-data Hasil Optimasi Algoritma Genetika **************') disp('********************************************************************'); HT=ts; %Hitting Time fprintf(1,'Nilai Settling Time : %f secon\n',HT); fprintf(1,'Nilai k terbaik : %f\n',best_k_Sf(1)); fprintf(1,'Nilai konstanta Sliding Surface terbaik : [c1=%f c2=%f]\n',Sf); u_rata=mean(us); ue=us-u_rata; fprintf(1,'MSE Input Kendali : %f\n',mse(ue)); xe1=xs(1,:)-xr(1); fprintf(1,'MSE Keluaran State (x1) : %f\n',mse(xe1)); xe2=xs(2,:)-xr(2); fprintf(1,'MSE state x2 : %f\n',mse(xe2)); fprintf(1,'Error Steady State (ess) : %f\n',CH); disp('********************************************************************'); %% %========================================================================== %Menggambarkan respon state terhadap waktu %-------------------------------------------------------------------------figure(2) for j=1:length(xs) t(j)=j*dt; end; plot(t,xs(1,:),'k'); grid on;hold on; o=1; plot(t,o,'k--');xlabel('time (secon)');legend('Kecepatan aktual','Kecepatan referensi') %-------------------------------------------------------------------------figure(3) plot(t,xs(2,:),'k');xlabel('time (secon)');legend('Percepatan'); grid on; %-------------------------------------------------------------------------figure(4); plot(t,us,'k'); xlabel('time (secon)'); ylabel('u'); legend('Masukan Kendali'); grid on; %-------------------------------------------------------------------------figure(5) m=-0.5:0.02:(xr(1)+0.5); for n=1:length(m) y(n)=-Sf(1)*(m(n)-xr(1)); end;
plot(m,y,'r'); grid on; hold on; plot(xs(1,:),xs(2,:),'k'); xlabel('x1 (speed)'); ylabel('x2 (acceleration'); hold off; %%
II. Program Inisialisasi Algoritma Genetika function phen=init_ga(bounds,siz,sea) for i=1:siz phen(i,:)=(bounds(:,2)-bounds(:,1))'.*rand(1,sea)+bounds(:,1)'; end
III. Program Evaluasi function [fitness,object]=evaluasi_smc_dcmotor(K,k,A,B,Sf,x0,xr,siz) for p=1:siz u=0; x0; xr; x=x0; dx=0; dt=0.01; %% for i=1:1000 x=x+dx; xs(:,i)=x; us(i)=u; dx=(A*x+B*u)*dt; s(p,:)=Sf(p,:)*x-Sf(p,:)*xr; ss(i,:)=s(p,:); %Control Law Ueq=K(p,:)*x; Un=-k(p,:)*inv(Sf(p,:)*B)*sign(s(p,:)); u=Ueq+Un; end; %% %menentukan nilai hitting time a=xr(1); % Jika state tidak sampai pada batas a, maka x_steady_array akan indices exceeds matrices dimension karena t nya % sudah sampai maksimum 1000, kemudian pada for m=1:10 % x_seady_state nya sudah tidak ada lagi. t=1; stop=0; while stop==0 if xs(1,t)
HT(p,:)=ts(p); %% u_rata=mean(us); ue=us-u_rata; u_err(p,:)=mse(ue);
%Hitting Time
xe1=xs(1,:)-xr(1); x_err(p,:)=mse(xe1);
%bobot untuk optimasi u_const=10; xe_const=1e7; c_ht=1; %Fungsi Object dan Fungsi Fitness object(p)=c_ht*(HT(p,:)^2)+u_const*(u_err(p,:)^2)+xe_const*(x_err(p,:)^2); %object(p)=c_ht*(HT(p,:)^2)+xe_const*(x_err(p,:)^2); fitness(p)=1/(object(p)+1); end
IV. Program Konversi Biner ke Float function [fval] = b2f(bval,bounds,bits) % function [fval] = b2f(bval,bounds,bits) % % Return the float number corresponing to the binary representation of bval. % % fval - the float representation of the number % bval - the binary representation of the number % bounds - the bounds on the variables % bits - the number of bits to represent each variable % % % % % % % % % % % % % % % % % %
Binary and Real-Valued Simulation Evolution for Matlab Copyright (C) 1996 C.R. Houck, J.A. Joines, M.G. Kay C.R. Houck, J.Joines, and M.Kay. A genetic algorithm for function optimization: A Matlab implementation. ACM Transactions on Mathmatical Software, Submitted 1996. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 1, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. A copy of the GNU General Public License can be obtained from the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
scale=(bounds(:,2)-bounds(:,1))'./(2.^bits-1); %The range of the variables numV=size(bounds,1); cs=[0 cumsum(bits)]; for i=1:numV a=bval((cs(i)+1):cs(i+1)); fval(i)=sum(2.^(size(a,2)-1:-1:0).*a)*scale(i)+bounds(i,1); end
V. Program Konversi Float ke Biner function [bval] = f2b(fval,bounds,bits) % function [bval] = f2b(fval,bounds,bits) % % Return the binary representation of the float number fval.
% % % % % % % % % % % % % % % % % % % % % % %
fval bval bounds bits
-
the the the the
float representation of the number binary representation of the number bounds on the variables number of bits to represent each variable
Binary and Real-Valued Simulation Evolution for Matlab Copyright (C) 1996 C.R. Houck, J.A. Joines, M.G. Kay C.R. Houck, J.Joines, and M.Kay. A genetic algorithm for function optimization: A Matlab implementation. ACM Transactions on Mathmatical Software, Submitted 1996. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 1, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. A copy of the GNU General Public License can be obtained from the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
scale=(2.^bits-1)./ (bounds(:,2)-bounds(:,1))'; %The range of the variables numV=size(bounds,1); cs=[0 cumsum(bits)]; bval=[]; for i=1:numV fval(i)=(fval(i)-bounds(i,1)) * scale(i); bval=[bval rem(floor(fval(i)*pow2(1-bits(i):0)),2)]; end
VI. Program Mate function [new_gen,mating]=mate(old_gen) [junk,mating]=sort(rand(size(old_gen,1),1)); new_gen=old_gen(mating,:);
VII. Program Mutate function [new_gen,mutated]=mutate(old_gen,pmut) mutated=find(rand(size(old_gen))
VIII. Program Reproduksi function [new_gen,selected]=reproduksi(old_gen,fitness) norm_fit=fitness/sum(fitness); selected=rand(size(fitness)); sum_fit=0;
%probabilitas fitness %membangkitkan secara random nilai fitness
%melakukan seleksi menggunakan rolute wheel for i=1:length(fitness) sum_fit=sum_fit+norm_fit(i); %nilai komulatif dari probailitas fitness index=find(selected<sum_fit); %jika nilai yang dibangkitkan lebih kecil dari sum_fit(i) maka akan diambil sum_fit(i) sebagai induk selected(index)=i*ones(size(index));%induk yang dipilih adalah index ke-i dari sum_fit(i) yang diambil end new_gen=old_gen(selected,:); %phenotif/chromosome baru hasil proses seleksi
IX. Program XOVER function [new_gen,sites]=xover(old_gen,pcross) lchrom=size(old_gen,2);
sites=ceil(rand(size(old_gen,1)/2,1)*(lchrom-1)); sites=sites.*(rand(size(sites))
