Bereksperimen dengan Particle Swarm Optimization

Salah satu pepatah hidup yang pernah disampaikan ke saya adalah dalam hidup ini tidak ada yang ideal. kita hanya bisa melakukan yang terbaik hari ini dengan segala keterbatasan. Hal yang abstrak namun dapat diterjemahkan menjadi persoalan yang selalu dihadapi manusia adalah optimasi dengan harapan hari ini lebih baik dari hari kemarin dan hari esok lebih baik dari hari ini.

Suatu persoalan optimasi pada dasarnya dapat dideskripsikan sebagai pengaturan atau pengisian nilai sekumpulan variabel pada suatu fungsi sedemikian rupa sehingga keluaran fungsi tersebut mendekati nilai yang diharapkan (minimum atau maksimum). Fungsi tersebut dinamakan dengan fungsi objektif atau kriteria dan jika fungsi objektif dapat didekomposisi menjadi beberapa fungsi objektif yang lebih spesifik dan maka persoalan optimasi diberi label sebagai persoalan optimasi multiobjektif atau persoalan optimasi multikriteria.

Dari perkuliahan semester yang lalu ada salah satu mata kuliah yaitu intelejensia kolektif yang mengkaji tentang swarm intelligence untuk melakukan penyelesaian persoalan secara tersebar. Topik pembahasan kuliah ini secara umum membahas teknik optimasi yang berbasis pada populasi objek atau agen dan interaksinya dalam menyelesaikan persoalan bersama. beberapa teknik yang diceritakan mulai dari algoritma evolusioner seperti algoritma genetik dan pemrograman genetik, memetic algorithm, hingga ACO (Ant Colony Optimization) dan PSO (Particle Swarm Optimization). Kedua topik terakhir merupakan topik yang dibahas secara lebih mendalam lewat review (dan atau implementasi) makalah tentang topik tersebut yang kemudian dipresentasikan di depan kelas oleh masing-masing peserta. Topik tentang ACO sebelumnya sempat saya ceritakan secara umum sehingga sekarang saya akan bercerita tentang PSO.

Particle Swarm Optimization

PSO terinspirasi dari perilaku gerakan kawanan hewan seperti ikan (school of fish), hewan herbivor (herd), dan burung (flock) yang kemudian tiap objek hewan disederhanakan menjadi sebuah partikel. Suatu partikel dalam ruang memiliki posisi yang dikodekan sebagai vektor koordinat. Vektor posisi ini dianggap sebagai keadaan yang sedang ditempati oleh suatu partikel di ruang pencarian. Setiap posisi dalam ruang pencarian merupakan alternatif solusi yang dapat dievaluasi menggunakan fungsi objektif. Setiap partikel bergerak dengan kecepatan v.

Ciri khas dari PSO adalah pengaturan kecepatan partikel secara heuristik dan probabilistik. Jika suatu partikel memiliki kecepatan yang konstan maka jika jejak posisi suatu partikel divisualisasikan akan membentuk garis lurus. Dengan adanya faktor eksternal yang membelokkan garis tersebut yang kemudian menggerakkan partikel dalam ruang pencarian maka diharapkan partikel dapat mengarah, mendekati, dan pada akhirnya mencapai titik optimal. faktor eksternal yang dimaksud antara lain posisi terbaik yang pernah dikunjungi suatu partikel, posisi terbaik seluruh partikel (diasumsikan setiap partikel mengetahui posisi terbaik setiap partikel lainnya), serta faktor kreativitas untuk melakukan eksplorasi. Secara matematis deskripsi di atas ditampilkan sebagai berikut :

v_i(t) = u₁.k₁.v_i(t-1) + u₂.k₂.(xbest_i–x_i) + u₃.k₃.(xbest_g–x_i) + u₄.vacak

xbest di atas merupakan catatan khusus mengenai posisi terbaik yang pernah dikunjungi tiap partikel. Indeks g menyatakan partikel yang posisi terbaiknya merupakan posisi terbaik dibandingkan posisi partikel lainnya. vacak merupakan vektor arah acak yang diinterpretasi sebagai faktor kreativitas untuk melakukan eksplorasi. nilai skalar u1 hingga u4 merupakan variabel acak (random variate) yang terdistribusi merata (uniform distribution) sedangkan a1 hingga a4 merupakan koefisien pengaruh masing-masing suku.

Implementasi

Agar dapat bereksperimen dengan algoritma ini maka kita perlu membuat implementasinya ke dalam program. Saya telah membuat implementasi sederhana dalam delphi seperti berikut :

{bagian interface}
type
  TVector = array of real;
  TParticle = record
    Position: TVector;
    Velocity: TVector; //velocity vector
    MyBest: TVector; //personal best position
    BestCost: Real;
  end;
  TCostFunction = function( v: TVector ): real of object;

  TPSO = class
  protected
    FNumParticle: integer;
    FVectorDimension: integer;
    FBest: integer;
    FCost: TCostFunction;
  public
    FParticles: array of TParticle;
    constructor Create( NumParticle, VectorLen: integer );
    destructor Free;
    procedure Step;
    property CostFn: TCostFunction read FCost write FCost;
    property Count: integer read FNumParticle;
    property Dimension: integer read FVectorDimension;
    property BestParticle: integer read FBest;
  end;
{...}
{ bagian implementasi }
{ TPSO }

constructor TPSO.Create( NumParticle, VectorLen: integer );
var
  i                 : integer;
begin
  FNumParticle := NumParticle;
  FVectorDimension := VectorLen;
  SetLength( FParticles, NumParticle );
  for i := 0 to High( FParticles ) do begin
    SetLength( FParticles[i].Position, VectorLen );
    SetLength( FParticles[i].Velocity, VectorLen );
    SetLength( FParticles[i].MyBest, VectorLen );
  end;
  FBest := 0;
end;

destructor TPSO.Free;
var
  i                 : integer;
begin
  for i := 0 to high( FParticles ) do begin
    setlength( FParticles[i].Position, 0 );
    setlength( FParticles[i].Velocity, 0 );
    setlength( FParticles[i].MyBest, 0 );
  end;
  setlength( Fparticles, 0 );
end;

procedure TPSO.Step;
var
  i, j              : integer;
  a, b, c, d        : real;
  u1, u2, u3, u4    : real;
  cost              : real;
  rvec              : TVector;
begin
  a := 1;
  b := 1;
  c := 1;
  d := 1;
  setlength( rvec, FVectorDimension );

  for i := 0 to High( FParticles ) do begin
    u1 := random;
    u2 := random;
    u3 := random;

    u4 := random;
    for j := 0 to FVectorDimension - 1 do
      rvec[j] := random;

    //update velocity
    for j := 0 to FVectorDimension - 1 do begin
      FParticles[i].Velocity[j] :=
        u1 * a * FParticles[i].Velocity[j] //velocity inertia
        + u2 * b * ( FParticles[i].MyBest[j] - FParticles[i].Position[j] ) //personal history best
        + u3 * c * ( FParticles[FBest].MyBest[j] - FParticles[i].Position[j] ) //global best
        + u4 * d * rvec[j] //random vector for exploration
        ;
    end;
  end;

  //update position
  for i := 0 to High( FParticles ) do begin
    for j := 0 to FVectorDimension - 1 do
      FParticles[i].Position[j] := FParticles[i].Position[j] + FParticles[i].Velocity[j];

    //calculate new cost
    cost := FCost( FParticles[i].Position );

    //update mybest if necessary
    if cost < FParticles[i].BestCost then begin
      FParticles[i].BestCost := cost;
      for j := 0 to high( FParticles[i].Position ) do
        FParticles[i].MyBest[j] := FParticles[i].Position[j];
    end;
  end;

  //select global best
  for i := 0 to High( FParticles ) do
    if FParticles[i].BestCost < FParticles[FBest].BestCost then
      FBest := i;
end;

Seiring dengan kode tersebut saya juga membuat program sederhana untuk bereksperimen dengan parameter PSO. Program tersebut memiliki fungsi objektif sederhana yaitu mencari satu titik di ruang berdimensi N sehingga fungsi evaluasinya pun cukup berupa fungsi jarak euclidean antara posisi tiap partikel dengan posisi titik tujuan.

program eksperimen PSO

potongan kode untuk aplikasi di atas adalah sbb:

{ di bagian definisi form }
    target: TVector;
    function dist( v: TVector ): real;
{ bagian implementasi }
procedure TForm1.FormCreate( Sender: TObject );
begin
  randomize;
end;

function TForm1.dist( v: TVector ): real;
var
  i                 : integer;
  tmp               : real;
begin
  result := 0;
  for i := 0 to high( v ) do begin
    tmp := ( target[i] - v[i] );
    tmp := tmp * tmp;
    result := result + tmp;
  end;
  result := sqrt( result );
end;

procedure TForm1.Button1Click( Sender: TObject );
var
  pso               : TPSO;
  i, j              : integer;
  Dim               : integer;
  thr               : real;

  procedure display;
  const
    w               = 4;
  var
    ptgt            : TPOint;
    i               : integer;
  begin
    with Image1.Canvas do begin
      FillRect( Image1.ClientRect );
      ptgt.X := round( target[0] );
      ptgt.Y := round( target[1] );
      Pen.Color := clRed;
      Rectangle( Rect( ptgt.X - w, ptgt.Y - w, ptgt.X + w, ptgt.Y + w ) );
      Pen.Color := clBlue;
      for i := 0 to pso.Count - 1 do begin
        ptgt.X := round( pso.FParticles[i].Position[0] );
        ptgt.Y := round( pso.FParticles[i].Position[1] );
        Rectangle( Rect( ptgt.X - w, ptgt.Y - w, ptgt.X + w, ptgt.Y + w ) );
      end;
    end;
  end;

begin
  thr := strtofloat( maskedit1.Text ); //MaskEdit1 : Error Threshold
  dim := spinedit2.Value; //SpinEdit2 : Vector Dimension

  setlength( target, dim );
  target[0] := random( image1.Width );
  target[1] := random( image1.Height );
  for j := 2 to high( target ) do
    target[j] := random( image1.Width );

  pso := TPSO.Create( spinedit1.Value, dim ); //SpinEdit1 : Number of Particle
  pso.CostFn := dist;
  for i := 0 to pso.Count - 1 do begin
    for j := 0 to pso.Dimension - 1 do begin
      pso.FParticles[i].Position[j] := random( image1.Width );
      pso.FParticles[i].Velocity[j] := 0;
      pso.FParticles[i].MyBest[j] := pso.FParticles[i].Position[j];
      pso.FParticles[i].BestCost := MAXINT;
    end;
  end;
  display;

  for i := 0 to SpinEdit3.Value do begin //SpinEdit3 : Max Iteration
    pso.Step;
    if pso.FParticles[pso.BestParticle].BestCost < thr then break;
    display;
    Image1.Canvas.TextOut( 0, 0, format( 'iteration : %d', [i] ) );
    Application.ProcessMessages;
  end;

  memo1.Lines.Add( format( 'best cost : %.8f', [pso.FParticles[pso.BestParticle].BestCost] ) );
  pso.Free;
end;

Program di atas akan membangkitkan posisi acak pada titik tujuan dan posisi partikel sebagai vektor berdimensi N dan menampilkan proyeksi dua sumbu pertama dari partikel (untuk alasan ini, jangan membuat vektor berdimensi satu ya!).

Pembahasan

Pada dasarnya rumus perubahan kecepatan di atas bisa bervariasi, namun dari eksperimen saya menggunakan program di atas saya menyimpulkan bahwa kombinasi keempat elemen di atas menghasilkan waktu konvergensi yang paling cepat dibandingkan jika salah satu suku dihilangkan untuk persoalan berdimensi > 3 (saya mencoba di 8 dimensi). Bagaimana menurut anda?