PSO

PSO(Particle Swarm Optimization,粒子族群最佳化)
developed in 1995 by Kennedy and Eberhart
inspired by social behavior of bird flocking
原理:一群particle分散在一個soluction space移動,每個particle會記錄自己的position和velocity

概念
1在pso中,每個particle代表一個solution
2判斷solution的好壞,是根據每個particle的fitness,該值依不同position會重新計算
3fitness的產生方式,依據解不同的問題,方法也會不同
ex:NOP(Numerical Optimization),MNO(Multiobjective Numerical Optimization),BIN(Binary Problem),clustering
4每個particle會在solution space中根據velocity不斷移動position,一直到停止條件成立
5velocity由3個值組成:1上次的velocity,2該particle本身的思考,3所有particle相互作用下的思考

pso資源
http://www.adaptivebox.net/research/bookmark/psocodes_link.html
http://www.particleswarm.info/
http://clerc.maurice.free.fr/pso/

資料結構
在D維空間下,
第i個粒子代表Xi,也等於(Xi1,Xi2,…XiD)
pbest:第i個粒子經歷過最好適應值的位置 Pi,也等於(Pi1,Pi2,…PiD)
gbest,所有第g群的粒子經歷過最好適應值的位置,Pg = (Pg1, Pg2, …, PgD),
第i個粒子的速度(移動的距離)代表Vi

在第d個維度的移動公式為
Vid_new=w*Vid+c1*rand()*(Pid-Xid)+c2*rand()*(Pgd-Xid)
ps:為了限制velocity在solution space,ex:若Vid_new > MAX_Vid,則Vid=MAX_Vid
Xid_new = Xid + Vid_new
ps:為了限制position在solution space,可檢查xid是否在該維度最大值與最小值之間

公式說明
w*Vid 粒子之前速度的慣性,值小容易local exploration有助於收斂,值大則相反
c1*rand()*(Pid-Xid)為粒子本身的思考
c2*rand()*(Pgd-Xid)為粒子相互影響下的思考
參數說明
Vid 粒子上次的速度,也就是粒子上次移動的距離
Vmax 粒子最大速度,太高會超過好解,太低則探索不足,常見的設定值是該維度最大及最小的距離
w:Inertia Weight (慣性權重)
c1,c2:acceleration constants
rand():0到1之間的隨機數

步驟大致如下
1random initial,隨機初始化Vid,Xid
2evaluate:計算每個粒子的fitneess(適應值)
3record pbest and gbest
對每個粒子,若該粒子現在位置的適應值優於pbest位置的適應值,則pbest位置=現在位置
對每個粒子,若該粒子現在位置的適應值優於gbest位置的適應值,則gbest位置=現在位置
4update Vid and Xid,根據移動公式調整粒子位置
5若達到結束條件(),否則返回步驟2
ps:結束條件設定通常為:設定固定世代數,或解一直沒有改善

演算法大致如下

//
initialize //step1
//
for t=1 to t_max do  //step5
 For each particle i do
assign zp to cluster
  calculate fitness //step2
 update pid(local best) and pgd(global best) //step3
//
 For each particle i do
  For each data vector zp
   vid = w*vid+c1*rand()*(pid-xid)+c2*Rand()*(pgd-xid) //step4
   if(vid>vmax)vid=vmax
   xid = xid+vid

 …………………………………………………….


運用在分群的做法
1
維度部份

在目標資料有D維空間下,若要分N群
粒子要維護每群的群中心,所以粒子Xi等於(Xi1,Xi2,…,XiN)
將每一群拆開後如下
(Xi11,Xi12,…,Xi1D
,Xi21,Xi22,…,Xi2D
,…
,XiN1,XiN2,…,XiND)
因此會對N*D個維度做計算
2
fitness部份

公式有以下兩種:
第一種:sum(j,sum_distance(Zp,Mj)/count(Cij) ) / N
第二種:sum(j,sum_distance(Zp,Mj))
j,表示第j群
count(Cij),第i個粒子的第j個群的資料筆數
Zp,第p個資料向量
Mj,第j個群中心
sum_distance(Zp,Mj),將組成Mj的所有Zp,拿來計算與Mj的距離並加總

…………………………………………………………………………………….
pso clustering for php example

/*********************************************
pso example

$data=psodata_fromfile(data_iris.txt);
$para['gene']=200;
$para['in']=10;
$para['jn']=3;
$para['w']=0.72;
$para['c1']=1.49;
$para['c2']=1.49;
$pso=pso($data['z'],$para);
print_r($pso['x']);
while($c=each($pso['c'])){
  while($j=each($c[1])){
    echo 'cluster'.$c[0].'=';
    while($p=each($data['full'][$j[1]])){
      echo $p[1].' : ';
    }
    echo "n";
  }
}
echo '['.$pso['fitness'].'][time:'.$pso['time_total'].']';

***********************************************/

///////////////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////////////
function pso($z,$para){
$time_start=caclutime();

/*the structure as follow
$para['gene']=200;
$para['in']=10;
$para['jn']=3;
$para['w']=0.72;
$para['c1']=1.49;
$para['c2']=1.49;

$p[i][1-D] 粒子i記錄的pbest
$pg[1-D] 記錄gbest
d=all atribute
p=all data

$z[1-D][1-p] 資料p的各維度數據
$z[attribute_1][data_1]
$z[attribute_1][data_2]
$z[attribute_1][data_p]
$z[attribute_2][data_1]
$z[attribute_2][data_2]
$z[attribute_2][data_p]
$z[attribute_3][data_1]
$z[attribute_3][data_2]
$z[attribute_3][data_p]

$x[i][1-N][1-D] 粒子i要維護的N個群中心
$x:record center of cluster
$x[patical_1][cluster_1][attribute_1]
$x[patical_1][cluster_1][attribute_d]
$x[patical_1][cluster_j][attribute_1]
$x[patical_1][cluster_j][attribute_d]

$c[i][j][1-count(Cij)] 在粒子i下,屬於第j群的所有資料
$c:record data of cluster
$c[patical_1][cluster_1][0]=sample1
$c[patical_1][cluster_1][1]=sample3
$c[patical_1][cluster_1][n]=sample6
$c[patical_1][cluster_2][0]=sample2
$c[patical_1][cluster_2][1]=sample4
$c[patical_1][cluster_2][n]=sample5
$c[patical_1][cluster_3][0]=sample7
*/

$pn=count($z[1]);
$dn=count($z);

for($d=1;$d<=$dn;$d++){
  arsort($z[$d]);
  $zdmax=each($z[$d]);
  asort($z[$d]);
  $zdmin=each($z[$d]);
  $xida['max'][$d]=$zdmax[1];
  $xida['min'][$d]=$zdmin[1];
  $vmax[$d]=abs($zdmax[1]-$zdmin[1]);
}
for($i=1;$i<=$para['in'];$i++){
  $fitness_pbest[$i]=9999999;
}
$fitness_gbest=9999999;

//init_clustercenter
reset($z);
$x=init_clustercenter($z,$para,$dn,$pn);

/////////
while($t<$para['gene']){
  $t++;

//evaluate fitness
  reset($x);
  while($i=each($x)){
    $state=asigncluster($i[1],$z,$dn,$pn);
    $c[$i[0]]=$state['c_i'];
    $f[$i[0]]=$state['d_total'];

    $fitness=$f[$i[0]];
    $fitness_gbest_history2_tmp[$i[0]]=$fitness;
    $fitness_pbest_history2[$i[0]][$t]=$fitness;
    if($fitness_pbest[$i[0]] > $fitness){
      $pbest[$i[0]]=$x[$i[0]];
      $fitness_pbest[$i[0]]=$fitness;
      if($fitness_gbest > $fitness_pbest[$i[0]]){
        $fitness_gbest=$fitness;
        $particle_gbest=$i[0];
        $gbest=$x[$i[0]];
        $bestc=$c[$i[0]];
      }
    }
    $fitness_pbest_history[$i[0]][$t]=$fitness_pbest[$i[0]];
  }
  asort($fitness_gbest_history2_tmp);
  $tmp_localgbest=each($fitness_gbest_history2_tmp);
  $fitness_gbest_history2[$t]=$tmp_localgbest[1];
  $fitness_gbest_history[$t]=$fitness_gbest;
  for($j=1;$j<=$para['jn'];$j++){
    $bestc_count[$j][$t]=count($bestc[$j]);
  }

//update
  reset($x);
  while($i=each($x)){
    while($j=each($i[1])){
      while($d=each($j[1])){
        $pid=$pbest[$i[0]][$j[0]][$d[0]];
        $pgd=$gbest[$j[0]][$d[0]];
        $xid=$d[1];
        $c1rand=rand(1,1000)/1000;
        $c2rand=rand(1,1000)/1000;
        $vid_new[$i[0]][$d[0]]=$para['w']*$vid[$i[0]][$d[0]]+$para['c1']*$c1rand*($pid-$xid)+$para['c2']*$c2rand*($pgd-$xid);
        if($vid_new[$i[0]][$d[0]]>$vmax[$d[0]])$vid_new[$i[0]][$d[0]]=$vmax[$d[0]];
        if($vid_new[$i[0]][$d[0]]<(0-$vmax[$d[0]]))$vid_new[$i[0]][$d[0]]=(0-$vmax[$d[0]]);
        $vid_abs[$i[0]][$d[0]]=abs($vid_new[$i[0]][$d[0]]);
        $new_xid=$vid_new[$i[0]][$d[0]]+$xid;
        if($new_xid>$xida['max'][$d[0]])$new_xid=$xida['max'][$d[0]];
        if($new_xid<$xida['min'][$d[0]])$new_xid=$xida['min'][$d[0]];
        $x[$i[0]][$j[0]][$d[0]]=$new_xid;
      }
    }
    $vidsum_history[$i[0]][$t]=array_sum($vid_abs[$i[0]]);
  }
  $vid=$vid_new;

}
////////////////////end

$pso['c']=$c[$particle_gbest];
$pso['fitness']=$fitness_gbest;
$pso['fitness_history']=$fitness_gbest_history;
$pso['fitness_history2']=$fitness_gbest_history2;
$pso['fitness_pbest_history']=$fitness_pbest_history;
$pso['fitness_pbest_history2']=$fitness_pbest_history2;
$pso['x']=$gbest;
$pso['time_total']=caclutime()-$time_start;
$pso['vidsum_history']=$vidsum_history;
$pso['bestc_count']=$bestc_count;

return $pso;

}

///////////////////////////////////////////////////////////////////////////////////////////////////////
function asigncluster($x_i,$z,$dn,$pn){
$d_total=0;
for($p=1;$p<=$pn;$p++){
  reset($x_i);
  while($j=each($x_i)){
    $distance_j=0;
    $distance_pj=0;
      for($d=1;$d<=$dn;$d++){
      $data_p_dime=$z[$d][$p];
      $clustercenter_j_dime=$j[1][$d];
      $distance_pj+=pow($data_p_dime-$clustercenter_j_dime,2);
    }
  $distancepj[$j[0]]=sqrt($distance_pj);
  }
asort($distancepj);
$min=each($distancepj);
$c_i[$min[0]][]=$p;
$d_total+=$min[1];
$check[$min[0]]=1;
}
for($j=1;$j<=$para['jn'];$j++){
  if($check[$j]!=1)$c_i[$j]='';
}
$state['c_i']=$c_i;
$state['d_total']=$d_total;

return $state;
}

//////////////////////////////////////////////////////////////////////////////////////////
function init_clustercenter($z,$para,$dn,$pn){
/*
$x[1][1][1]=5.006;
$x[1][1][2]=3.418;
$x[1][1][3]=1.464;
$x[1][1][4]=0.244;
$x[1][2][1]=5.883606557377;
$x[1][2][2]=2.7409836065574;
$x[1][2][3]=4.3885245901639;
$x[1][2][4]=1.4344262295082;
$x[1][3][1]=6.8538461538462;
$x[1][3][2]=3.0769230769231;
$x[1][3][3]=5.7153846153846;
$x[1][3][4]=2.0538461538462;
*/

$len=ceil($pn/$para['jn']);
for($i=2;$i<=$para['in'];$i++){
  for($j=1;$j<=$para['jn'];$j++){
    $len_e=$len*$j;
    $len_s=$len*($j-1)+1;
    if($len_e>=$pn)$len_e=$pn;
    $center=rand($len_s,$len_e);
    for($d=1;$d<=$dn;$d++){
      $x[$i][$j][$d]=$z[$d][$center];
    }
  }
}

return $x;
}

//////////////////////////////////////////////////////////////////////////////
function caclutime(){
$time = explode( " ", microtime());
$usec = (double)$time[0];
$sec = (double)$time[1];
return $sec + $usec;
}

/////////////////////////////////////////////////////////////////////////////
function psodata_fromfile($sourcedata){
/*
data fromat description
name1,attribute1_value,attribute2_value,...attributeN_value
name2,attribute1_value,attribute2_value,...attributeN_value
...
nameN,attribute1_value,attribute2_value,...attributeN_value
ex:
Iris-setosa,5.1,3.5,1.4,0.2
Iris-setosa,4.9,3,1.4,0.2
Iris-setosa,4.7,3.2,1.3,0.2
Iris-setosa,4.6,3.1,1.5,0.2
Iris-setosa,5,3.6,1.4,0.2
ps:attribute must is continue,not discrete
*/
$data=file($sourcedata);

$dime=each($data);
$dime_a=explode(',',$dime[1]);
$dn=count($dime_a)-1;

reset($data);
while(list($key,$dot_id)=each($data)){
  $i++;
  $dot_id = trim($dot_id, " n.");
  $dot_id_a=explode(',',$dot_id);
  $fullz[$i][0]=$dot_id_a[0];
  for($d=1;$d<=$dn;$d++){
    $z[$d][$i]=$dot_id_a[$d];
    $fullz[$i][$d]=$dot_id_a[$d];
  }
}
$data['z']=$z;
$data['full']=$fullz;
return $data;
}