Sedikit mengingat kembali salah satu mata kuliah favorit saya yaitu Teknik Optimasi. Salah satu dari sekian teknik optimasi yang mudah untuk dipahami, dan cukup sederhana untuk dibuat programnya adalah Particle Swarm Optimization (PSO). Berdasarkan Wikipedia, PSO adalah metode komputasi untuk mencari solusi optimal dari suatu persamaan (masalah) secara iterasi dengan menggunakan sejumlah populasi (atau kandidat solusi) yang bergerak menuju titik terbaik (solusi) berdasarkan formula sederhana. Dalam tulisan kali ini, berdasarkan refensi dari WikiPedia, kita akan mencoba membuat/menuliskan algoritma PSO dalam bahasa pemrograman Octave.
 | ||
| Ilustrasi Particle Swarm |
Dalam tutorial ini, problem yang ingin diselesaikan adalah untuk mencari titik minimum $(x_\textrm{best},y_\textrm{best})$ dari fungsi di bawah ini
$f(x,y) = x \cdot \exp(-(x^2+y^2))$
dengan rentang solusi yang ingin dicari adalah
$-3.5 \leq x \leq 3.5$ and $-2.5 \leq y \leq 2.5$
Dalam tutorial ini, persamaan di atas sengaja dipilih untuk memudahkan memahami bagaimana PSO mencari solusi minimum. Langkah pertama yang perlu dilakukan adalah membuat octave function dari persamaan di atas dan membuat plot 3D-nya.
(1) Simpan kode octave di bawah ini dengan nama "f_obj.m".
function f = f_obj(x,y)
% fungsi objektif
f = x.*exp(-(x.^2+y.^2));
(2) Simpan pula kode octave berikut dengan nama "main.m" kemudian jalankan (run).
close all
clear
clc
%% PROBLEM
x1 = linspace(-3.5,3.5,53);
x2 = linspace(-2.5,2.5,53);
[X1 X2] = meshgrid(x1,x2);
F = f_obj(X1,X2);
figure(1)
pause(2)
contour3(X1,X2,F,30)
surface(X1,X2,F,'EdgeColor',[.8 .8 .8],'FaceColor','none')
colormap cool
view(348,17)
axis([-3.5 3.5 -2.5 2.5 -.5 .5])
grid off
hasilnya akan seperti pada gambar di bawah ini
Particle Swarm Optimization (PSO)
Dalam contoh ini, kita menggunakan persamaan dengan dua variabel yang ingin diselesaikan $x$ dan $y$. Dari gambar di atas, terlihat dengan jelas sekali dimana letak titik minimum. Bagaimana cara mendapatkan solusi titik tersebut dengan menggunakan PSO? berdasarkan penjelasan dari Wikipedia, ada beberapa step yang harus dilakukan yakni sebagai berikut.
(1)
Inisialisasi populasi (kandidat solusi) secara random dan berdistribusi uniform yaitu $\textbf{x}_i \sim U(\textbf{b}_{io},\textbf{b}_{up})$ dimana $\textbf{b}_{io}$ dan $\textbf{b}_{up}$ adalah batas bawah dan atas dari solusi yang ingin dicari (rentang solusi). Inisialisasi populasi dapat dilakukan dengan menambahkan kode berikut pada "main.m",
close all
clear
clc
%% PROBLEM
x1 = linspace(-3.5,3.5,53);
x2 = linspace(-2.5,2.5,53);
[X1 X2] = meshgrid(x1,x2);
F = f_obj(X1,X2);
figure(1)
pause(2)
contour3(X1,X2,F,30)
surface(X1,X2,F,'EdgeColor',[.8 .8 .8],'FaceColor','none')
colormap cool
view(348,17)
axis([-3.5 3.5 -2.5 2.5 -.5 .5])
grid off
%% ALGORITMA PSO
% PARAMETERS
% Jumlah Partikel Swarm yang digunakan untuk mencari solusi
S = 300;
% Jumlah Parameter yang akan dioptimasi
dim = 2;
% Rentang solusi yang ingin didapat
blo = [-2.5 -3.5]; % the lower boundaries
bup = [2.5 3.5]; % the upper boundaries
%(1)
% Inisialisasi lokasi partikel (solusi) secara random
xi = zeros(S,dim);
for i = 1:dim
xi(:,i) = rand(S,1)*(bup(i)-blo(i)) + blo(i);
end
hold on
plot3(xi(:,1),xi(:,2),f_obj(xi(:,1),xi(:,2)),'b.','Markersize',6)
hold off
Ketika kode di atas dijalankan, kita akan mendapatkan hasil seperti pada gambar di bawah ini,
(2)
Inisialisasi solusi-solusi yang mungkin yakni $\textbf{p}_{i} \leftarrow \textbf{x}_{i}$, kemudian dapatkan global minimum $\textbf{g} \leftarrow \textbf{p}_{i}$. Dalam octave, step ini dapat dituliskan sebagai berikut.
(1) Simpan kode octave di bawah ini dengan nama "f_obj.m".
function f = f_obj(x,y)
% fungsi objektif
f = x.*exp(-(x.^2+y.^2));
(2) Simpan pula kode octave berikut dengan nama "main.m" kemudian jalankan (run).
close all
clear
clc
%% PROBLEM
x1 = linspace(-3.5,3.5,53);
x2 = linspace(-2.5,2.5,53);
[X1 X2] = meshgrid(x1,x2);
F = f_obj(X1,X2);
figure(1)
pause(2)
contour3(X1,X2,F,30)
surface(X1,X2,F,'EdgeColor',[.8 .8 .8],'FaceColor','none')
colormap cool
view(348,17)
axis([-3.5 3.5 -2.5 2.5 -.5 .5])
grid off
hasilnya akan seperti pada gambar di bawah ini
 |
| Contoh Problem untuk mencari titik minimum |
Particle Swarm Optimization (PSO)
Dalam contoh ini, kita menggunakan persamaan dengan dua variabel yang ingin diselesaikan $x$ dan $y$. Dari gambar di atas, terlihat dengan jelas sekali dimana letak titik minimum. Bagaimana cara mendapatkan solusi titik tersebut dengan menggunakan PSO? berdasarkan penjelasan dari Wikipedia, ada beberapa step yang harus dilakukan yakni sebagai berikut.
(1)
Inisialisasi populasi (kandidat solusi) secara random dan berdistribusi uniform yaitu $\textbf{x}_i \sim U(\textbf{b}_{io},\textbf{b}_{up})$ dimana $\textbf{b}_{io}$ dan $\textbf{b}_{up}$ adalah batas bawah dan atas dari solusi yang ingin dicari (rentang solusi). Inisialisasi populasi dapat dilakukan dengan menambahkan kode berikut pada "main.m",
close all
clear
clc
%% PROBLEM
x1 = linspace(-3.5,3.5,53);
x2 = linspace(-2.5,2.5,53);
[X1 X2] = meshgrid(x1,x2);
F = f_obj(X1,X2);
figure(1)
pause(2)
contour3(X1,X2,F,30)
surface(X1,X2,F,'EdgeColor',[.8 .8 .8],'FaceColor','none')
colormap cool
view(348,17)
axis([-3.5 3.5 -2.5 2.5 -.5 .5])
grid off
%% ALGORITMA PSO
% PARAMETERS
% Jumlah Partikel Swarm yang digunakan untuk mencari solusi
S = 300;
% Jumlah Parameter yang akan dioptimasi
dim = 2;
% Rentang solusi yang ingin didapat
blo = [-2.5 -3.5]; % the lower boundaries
bup = [2.5 3.5]; % the upper boundaries
%(1)
% Inisialisasi lokasi partikel (solusi) secara random
xi = zeros(S,dim);
for i = 1:dim
xi(:,i) = rand(S,1)*(bup(i)-blo(i)) + blo(i);
end
hold on
plot3(xi(:,1),xi(:,2),f_obj(xi(:,1),xi(:,2)),'b.','Markersize',6)
hold off
Ketika kode di atas dijalankan, kita akan mendapatkan hasil seperti pada gambar di bawah ini,
 |
| Kandidat Solusi |
Inisialisasi solusi-solusi yang mungkin yakni $\textbf{p}_{i} \leftarrow \textbf{x}_{i}$, kemudian dapatkan global minimum $\textbf{g} \leftarrow \textbf{p}_{i}$. Dalam octave, step ini dapat dituliskan sebagai berikut.
close all
clear
clc
%% PROBLEM
x1 = linspace(-3.5,3.5,53);
x2 = linspace(-2.5,2.5,53);
[X1 X2] = meshgrid(x1,x2);
F = f_obj(X1,X2);
figure(1)
pause(2)
contour3(X1,X2,F,30)
surface(X1,X2,F,'EdgeColor',[.8 .8 .8],'FaceColor','none')
colormap cool
view(348,17)
axis([-3.5 3.5 -2.5 2.5 -.5 .5])
grid off
%% ALGORITMA PSO
% PARAMETERS
% Jumlah Partikel Swarm yang digunakan untuk mencari solusi
S = 300;
% Jumlah Parameter yang akan dioptimasi
dim = 2;
% Rentang solusi yang ingin didapat
blo = [-2.5 -3.5]; % the lower boundaries
bup = [2.5 3.5]; % the upper boundaries
%(1)
% Inisialisasi lokasi partikel (solusi) secara random
xi = zeros(S,dim);
for i = 1:dim
xi(:,i) = rand(S,1)*(bup(i)-blo(i)) + blo(i);
end
hold on
plot3(xi(:,1),xi(:,2),f_obj(xi(:,1),xi(:,2)),'b.','Markersize',6)
hold off
%(2)
% Inisialisasi solusi-solusi yang mungkin
pi = xi;
% Dapatkan global minimum
f = f_obj(xi(:,1),xi(:,2));
[min_f id_g] = min(f);
g = pi(id_g,:);
(3)
Inisialisasi perubahan kandidat solusi per satuan iterasi (kecepatan) $\textbf{v}_i \sim U(-|\textbf{b}_{up}-\textbf{b}_{lo}|,|\textbf{b}_{up}-\textbf{b}_{lo}|)$ yakni sebagai berikut.
close all
clear
clc
%% PROBLEM
x1 = linspace(-3.5,3.5,53);
x2 = linspace(-2.5,2.5,53);
[X1 X2] = meshgrid(x1,x2);
F = f_obj(X1,X2);
figure(1)
pause(2)
contour3(X1,X2,F,30)
surface(X1,X2,F,'EdgeColor',[.8 .8 .8],'FaceColor','none')
colormap cool
view(348,17)
axis([-3.5 3.5 -2.5 2.5 -.5 .5])
grid off
%% ALGORITMA PSO
% PARAMETERS
% Jumlah Partikel Swarm yang digunakan untuk mencari solusi
S = 300;
% Jumlah Parameter yang akan dioptimasi
dim = 2;
% Rentang solusi yang ingin didapat
blo = [-2.5 -3.5]; % the lower boundaries
bup = [2.5 3.5]; % the upper boundaries
%(1)
% Inisialisasi lokasi partikel (solusi) secara random
xi = zeros(S,dim);
for i = 1:dim
xi(:,i) = rand(S,1)*(bup(i)-blo(i)) + blo(i);
end
hold on
plot3(xi(:,1),xi(:,2),f_obj(xi(:,1),xi(:,2)),'b.','Markersize',6)
hold off
%(2)
% Inisialisasi solusi-solusi yang mungkin
pi = xi;
% Dapatkan global minimum
f = f_obj(xi(:,1),xi(:,2));
[min_f id_g] = min(f);
g = pi(id_g,:);
%(3)
% Inisialisasi perubahan kandidat solusi per satuan iterasi (kecepatan)
vi = zeros(S,dim);
for i = 1:dim
vi(:,i) = rand(S,1)*(2*abs(bup(i)-blo(i))) - abs(bup(i)-blo(i));
end
(4)
Tahap ini adalah tahap iterasi dimana iterasi untuk meng-update nilai $\textbf{x}_i$ akan terus berjalan hingga memenuhi kriteria yang kita tetapkan. Tahap ini terdiri dari beberapa langkah yakni sebagai berikut
- update kecepatan partikel, $\textbf{v}_{i,d} \leftarrow \omega \cdot \textbf{v}_{i,d} + \varphi_p \cdot r_p \cdot (\textbf{p}_{i,d}-\textbf{x}_{i,d}) + \varphi_g \cdot r_g \cdot (\textbf{g}_{i,d}-\textbf{x}_{i,d})$. $d$ dalam persamaan menunjukkan jumlah parameter yang digunakan (dalam contoh ini $d=2$ yakni $x$ dan $y$) Pemilihan parameter untuk proses update ini memiliki kajian tersendiri, jadi persamaan ini akan dibuat menjadi lebih sederhana seperti berikut ini: $\textbf{v}_{i,d} \leftarrow \textrm{rand} \cdot \textbf{v}_{i,d} + \textrm{rand} \cdot (\textbf{p}_{i,d}-\textbf{x}_{i,d}) + \textrm{rand} \cdot (\textbf{g}_{i,d}-\textbf{x}_{i,d})$ dimana $\textrm{rand}$ menunjukan angka random dalam rentang 0 dan 1.
- update $\textbf{x}_i \leftarrow \textbf{x}_i + \textbf{v}_i$
- update juga $\textbf{p}_{i} \leftarrow \textbf{x}_{i}$
- cek apakah global minimum yang baru lebih kecil dari global minimum yang sebelumnya, jika iya update global minimum. Jika $f(\textbf{p}_i) < f(\textbf{g})$, maka $\textbf{g} \leftarrow \textbf{p}_i$.
- langkah terakhir, jika semua solusi dalam $\textbf{p}_i$ memiliki nilai yang hampir sama, maka iterasi harus diberhentikan. Di akhir iterasi $\textbf{g}$ adalah solusi yang diperoleh oleh PSO.
close all
clear
clc
%% PROBLEM
x1 = linspace(-3.5,3.5,53);
x2 = linspace(-2.5,2.5,53);
[X1 X2] = meshgrid(x1,x2);
F = f_obj(X1,X2);
figure(1)
pause(2)
contour3(X1,X2,F,30)
surface(X1,X2,F,'EdgeColor',[.8 .8 .8],'FaceColor','none')
colormap cool
view(348,17)
axis([-3.5 3.5 -2.5 2.5 -.5 .5])
grid off
%% ALGORITMA PSO
% PARAMETERS
% Jumlah Partikel Swarm yang digunakan untuk mencari solusi
S = 300;
% Jumlah Parameter yang akan dioptimasi
dim = 2;
% Rentang solusi yang ingin didapat
blo = [-2.5 -3.5]; % the lower boundaries
bup = [2.5 3.5]; % the upper boundaries
%(1)
% Inisialisasi lokasi partikel (solusi) secara random
xi = zeros(S,dim);
for i = 1:dim
xi(:,i) = rand(S,1)*(bup(i)-blo(i)) + blo(i);
end
hold on
plot3(xi(:,1),xi(:,2),f_obj(xi(:,1),xi(:,2)),'b.','Markersize',6)
hold off
%(2)
% Inisialisasi solusi-solusi yang mungkin
pi = xi;
% Dapatkan global minimum
f = f_obj(xi(:,1),xi(:,2));
[min_f id_g] = min(f);
g = pi(id_g,:);
%(3)
% Inisialisasi perubahan kandidat solusi per satuan iterasi (kecepatan)
vi = zeros(S,dim);
for i = 1:dim
vi(:,i) = rand(S,1)*(2*abs(bup(i)-blo(i))) - abs(bup(i)-blo(i));
end
%(4)
% Lakukan iterasi hingga kriteria terpenuhi
for iterasi = 1:1111
% Update perubahan solusi per satuan iterasi (kecepatan)
vi = rand*vi + rand*(pi-xi) + rand*(g-xi);
% Update lokasi partikel (solusi)
xi = xi + vi;
% Solusi harus dalam rentang yang sudah ditentukan
for i = 1:dim
xi(xi(:,i) < blo(i),i) = blo(i);
xi(xi(:,i) > bup(i),i) = bup(i);
end
% Update solusi-solusi yang mungkin
pi = xi;
% Update solusi terbaik
fp = f_obj(pi(:,1),pi(:,2));
fg = f_obj(g(1),g(2));
[min_fp id_g] = min(fp);
if min_fp < fg
g = pi(id_g,:);
end
% Jika semua solusi dalam pi memiliki nilai yang hampir sama
% Maka, iterasi harus distop
if sum(var(pi)) < 1e-6
disp(['Total Iterasi = ' num2str(iterasi)])
break
end
contour3(X1,X2,F,30)
surface(X1,X2,F,'EdgeColor',[.8 .8 .8],'FaceColor','none')
colormap cool
view(348,17)
axis([-3.5 3.5 -2.5 2.5 -.5 .5])
grid off
hold on
plot3(pi(:,1),pi(:,2),f_obj(pi(:,1),pi(:,2)),'b.','Markersize',7)
pause(0.1)
hold off
end
Contoh perubahan tiap iterasi dapat dilihat pada gambar di bawah ini,
Akhirnya titik minimum berhasil ditemukan pada $(-0.70711, 8.2835 \times 10^{-6})$
Download
kode octave dapat didownload di link ini.
Selamat mencoba dan semoga bermanfaat.
Referensi
https://en.wikipedia.org/wiki/Particle_swarm_optimization
No comments:
Post a Comment