粒子群演算法c語言

『壹』 優化演算法筆記（五）粒子群演算法（3）

(已合並本篇內容至粒子群演算法（1）)

上一節中，我們看到小鳥們聚集到一個較小的范圍內後，不會再繼續集中。這是怎麼回事呢？
猜測：
1.與最大速度限制有關，權重w只是方便動態修改maxV。
2.與C1和C2有關，這兩個權重限制了鳥兒的搜索行為。
還是上一節的實驗， 。現在我們將maxV的值有5修改為50，即maxV=50，其他參數不變。參數如下

此時得到的最優位值的適應度函數值為0.25571,可以看出與maxV=5相比，結果差了很多而且小鳥們聚集的范圍更大了。
現在我們設置maxV=1，再次重復上面的實驗，實驗結果如下：

這次最終的適應度函數值為，比之前的結果都要好0.00273。從圖中我們可以看出，小鳥們在向一個點集中，但是他們飛行的速度比之前慢多了，如果問題更復雜，可能無法等到它們聚集到一個點，迭代就結束了。
為什麼maxV會影響鳥群的搜索結果呢？
我們依然以maxV=50為例，不過這次為了看的更加清晰，我們的鳥群只有2隻鳥，同時將幀數放慢5倍以便觀察。

思路一：限制鳥的最大飛行速率,由於慣性系數W的存在，使得控制最大速率和控制慣性系數的效果是等價的，取其一即可。
方案1：使慣性系數隨著迭代次數增加而降低，這里使用的是線性下降的方式，即在第1次迭代，慣性系數W=1,最後一次迭代時，慣性系數W=0，當然，也可以根據自己的意願取其他值。
實驗參數如下：

小鳥們的飛行過程如上圖，可以看到效果很好，最後甚至都聚集到了一個點。再看看最終的適應度函數值8.61666413451519E-17，這已經是一個相當精確的值了，說明這是一個可行的方案，但是由於其最後種群過於集中，有陷入局部最優的風險。
方案2：給每隻鳥一個隨機的慣性系數，那麼鳥的飛行軌跡也將不再像之前會出現周期性。每隻鳥的慣性系數W為（0,2）中的隨機數（保持W的期望為1）。
實驗參數如下：

可以看到小鳥們並沒有像上一個實驗一樣聚集於一個點，而是仍在一個較大的范圍內進行搜索。其最終的適應度函數為0.01176，比最初的0.25571稍有提升，但並不顯著。什麼原因造成了這種情況呢？我想可能是由於慣性系數成了期望為1的隨機數，那麼小鳥的飛行軌跡的期望可能仍然是繞著一個四邊形循環，只不過這個四邊形相比之前的平行四邊形更加復雜，所以其結果也稍有提升，當然對於概率演算法，得到這樣的結果可能僅僅是因為運氣不好
我們看到慣性系數W值減小，小鳥們聚攏到一處的速度明顯提升，那麼，如果我們去掉慣性系數這個參數會怎麼樣呢。
方案3：取出慣性系數，即取W=0，小鳥們只向著那兩個最優位置飛行。

可以看見鳥群們迅速聚集到了一個點，再看看得到的結果，最終的適應度函數值為2.9086886073362966E-30，明顯優於之前的所有操作。
那麼問題來了，為什麼粒子群演算法需要一個慣性速度，它的作用是什麼呢？其實很明顯，當鳥群迅速集中到了一個點之後它們就喪失了全局的搜索能力，所有的鳥會迅速向著全局最優點飛去，如果當前的全局最優解是一個局部最優點，那麼鳥群將會陷入局部最優。所以，慣性系數和慣性速度的作用是給鳥群提供跳出局部最優的可能性，獲得這個跳出局部最優能力的代價是它們的收斂速度減慢，且局部的搜索能力較弱（與當前的慣性速度有關）。
為了平衡局部搜索能力和跳出局部最優能力，我們可以人為的干預一下慣性系數W的大小，結合方案1和方案2，我們可以使每隻鳥的慣性系數以一個隨機周期，周期性下降，若小於0，則重置為初始值。

這樣結合了方案1和方案2的慣性系數，也能得到不錯的效果，大家可以自己一試。

思路二：改變小鳥們向群體最優飛行和向歷史最優飛行的權重。
方案4：讓小鳥向全局最優飛行的系數C2線性遞減。

小鳥們的飛行過程與之前好像沒什麼變化，我甚至懷疑我做了假實驗。看看最終結果，0.7267249621552874，這是到目前為止的最差結果。看來這不是一個好方案，讓全局學習因子C2遞減，勢必會降低演算法的收斂效率，而慣性系數還是那麼大，小鳥們依然會圍繞歷史最優位置打轉，畢竟這兩個最優位置是有一定關聯的。所以讓C1線性遞減的實驗也不必做了，其效果應該與方案4相差不大。
看來只要是慣性系數不變怎麼修改C1和C2都不會有太過明顯的效果。為什麼實驗都是參數遞減，卻沒有參數遞增的實驗呢？
1.慣性系數W必須遞減，因為它會影響鳥群的搜索范圍。
2.如果C1和C2遞增，那麼小鳥的慣性速度V勢必會跟著遞增，這與W遞增會產生相同的效果。

上面我們通過一些實驗及理論分析了粒子群演算法的特點及其參數的作用。粒子群作為優化演算法中模型最簡單的演算法，通過修改這幾個簡單的參數也能夠改變演算法的優化性能可以說是一個非常優秀的演算法。
上述實驗中，我們僅分析了單個參數對演算法的影響，實際使用時（創新、發明、寫論文時）也會同時動態改變多個參數，甚至是參數之間產生關聯。
實驗中，為了展現實驗效果，maxV取值較大，一般取值為搜索空間范圍的10%-20%，按上面（-100,100）的范圍maxV應該取值為20-40，在此基礎上，方案1、方案2效果應該會更好。
粒子群演算法是一種概率演算法，所以並不能使用一次實驗結果來判斷演算法的性能，我們需要進行多次實驗，然後看看這些實驗的效果最終來判斷，結果必須使用多次實驗的統計數據來說明，一般我們都會重復實驗30-50次，為了發論文去做實驗的小夥伴們不要偷懶哦。
粒子群演算法的學習目前告一段落，如果有什麼新的發現，後面繼續更新哦！
以下指標純屬個人yy,僅供參考

目錄
上一篇 優化演算法筆記（四）粒子群演算法（2）
下一篇 優化演算法筆記（六）遺傳演算法

『貳』 C語言，用粒子群演算法求函數最大值，拜託了！！

for i=1:sizepop % 隨機產生一個種群 pop(i,:)=2*rands(1,2); % 初始化粒子 V(i,:)=0.5*rands(1,2); % 初始化速度 % 計算粒子適應度值 fitness(i)=fun(pop(i,:)); end [bestfitness bestindex]=min(fitness); zbest=pop(bestindex,:); % 群體極...

『叄』 粒子群演算法簡單介紹

粒子群演算法（也稱粒子群優化演算法（particle swarm optimization, PSO）），模擬鳥群隨機搜索食物的行為。粒子群演算法中，每個優化問題的潛在解都是搜穗穗索空間中的一隻鳥，叫做「粒子」。所有的粒子都有一個由被優化的函數決定的適應值（fitness value），每個粒子還有一個速度決定它們「飛行」的方向和距離。

粒子群演算法初始化為一群隨機的粒子（隨機解），然後根據迭代找到最優解。每一次迭代中，粒子通過跟蹤兩個極敏擾值來更新自己：第1個是粒子本身所找到的最優解，這個稱為個體極值；第2個是整個種群目前找到的最優解，這個稱為全局極值。也可以不用整個種群，而是用其中的一部分作為粒子的鄰居猜拿卜，稱為局部極值。

假設在一個D維搜索空間中，有N個粒子組成一個群落，其中第i個粒子表示為一個D維的向量：

第i個粒子的速度表示為：

還要保存每個個體的已經找到的最優解 ，和一個整個群落找到的最優解 。

第i個粒子根據下面的公式更新自己的速度和位置：

其中， 是個體已知最優解， 是種群已知最優解， 為慣性權重， , 為學習因子（或加速常數 acceleration constant)， , 是[0,1]范圍內的隨機數。

式(1)由三部分組成：

『肆』 粒子群演算法

粒子群演算法（Particle Swarm Optimization），又稱鳥群覓食演算法，是由數學家J. Kennedy和R. C. Eberhart等開發出的一種新的進化演算法。它是從隨機解開始觸發，通過迭代尋找出其中的最優解。本演算法主要是通過適應度來評價解的分數，比傳統的遺傳演算法更加的簡單，它沒有傳統遺傳演算法中的「交叉」和「變異」等操作，它主要是追隨當前搜索到的最優值來尋找到全局最優值。這種演算法實現容易，精度高，收斂快等特點被廣泛運用在各個問題中。

粒子群演算法是模擬鳥群覓食的所建立起來的一種智能演算法，一開始所有的鳥都不知道食物在哪裡，它們通過找到離食物最近的鳥的周圍，再去尋找食物，這樣不斷的追蹤，大量的鳥都堆積在食物附近這樣找到食物的幾率就大大增加了。粒子群就是這樣一種模擬鳥群覓食的過程，粒子群把鳥看成一個個粒子，它們擁有兩個屬性——位置和速度，然後根據自己的這兩個屬性共享到整個集群中，其他粒子改變飛行方向去找到最近的區域，然後整個集群都聚集在最優解附近，最後最終找到最優解。

演算法中我們需要的數據結構，我們需要一個值來存儲每個粒子搜索到的最優解，用一個值來存儲整個群體在一次迭代中搜索到的最優解，這樣我們的粒子速度和位置的更新公式如下：

其中pbest是每個粒子搜索到的最優解，gbest是整個群體在一次迭代中搜索到的最優解，v[i]是代表第i個粒子的速度，w代表慣性系數是一個超參數，rang()表示的是在0到1的隨機數。Present[i]代表第i個粒子當前的位置。我們通過上面的公式不停的迭代粒子群的狀態，最終得到全局最優解

『伍』 關於粒子群演算法的問題

粒子群的版本甚多，常用的是加有慣性權重w的
v[] = w * v[] + c1 * rand() * (pbest[] - present[]) + c2 * rand() * (gbest[] - present[])
一般選擇慣性權重在迭代過程中線性下降，目的是在迭代的初期，以比較大的權重分配給粒子的原速度，而防止粒子過早的傾向於其本身的局部最優與全局最優，此時的全局搜索能力是可以的。但粒子群是基於牛頓力學的，隨著w的減小，速度v的作用會在更新中弱化，對應的是，pbest和gbest的作用得到了加強，這也就意味著，粒子會更加趨向於pbest和gbest的方向移動。這個時候粒子就特別容易陷入局部最優了。

其實陷入局部最優不只是粒子群的問題，進化類的演算法都存在這個問題，只不過有些演算法隨機性強一些，收斂速度慢一些，所以更加容易跳出局部最優（但不是絕對避免）

『陸』 粒子群演算法的演算法介紹

如前所述，PSO模擬鳥群的捕食行為。設想這樣一個場景：一群鳥在隨機搜索食物。在這個區域里只有一塊食物。所有的鳥都不知道食物在那裡。但是他們知道當前的位置離食物還有多遠。那麼找到食物的最優策略是什麼呢。最簡單有效的就是搜尋目前離食物最近的鳥的周圍區域。
PSO從這種模型中得到啟示並用於解決優化問題。PSO中，每個優化問題的解都是搜索空間中的一隻鳥。我們稱之為「粒子」。所有的粒子都有一個由被優化的函數決定的適應值(fitness value)，每個粒子還有一個速度決定他們飛翔的方向和距離。然後粒子們就追隨當前的最優粒子在解空間中搜索。
PSO 初始化為一群隨機粒子(隨機解)。然後通過迭代找到最優解。在每一次迭代中，粒子通過跟蹤兩個極值來更新自己。第一個就是粒子本身所找到的最優解，這個解叫做個體極值pBest。另一個極值是整個種群目前找到的最優解，這個極值是全局極值gBest。另外也可以不用整個種群而只是用其中一部分作為粒子的鄰居，那麼在所有鄰居中的極值就是局部極值。 在找到這兩個最優值時，粒子根據如下的公式來更新自己的速度和新的位置：
v[] = w * v[] + c1 * rand() * (pbest[] - present[]) + c2 * rand() * (gbest[] - present[]) (a)
present[] = present[] + v[] (b)
v[] 是粒子的速度, w是慣性權重,present[] 是當前粒子的位置. pbest[] and gbest[] 如前定義 rand () 是介於（0， 1）之間的隨機數. c1, c2 是學習因子. 通常 c1 = c2 = 2.
程序的偽代碼如下
For each particle
____Initialize particle
END
Do
____For each particle
________Calculate fitness value
________If the fitness value is better than the best fitness value (pBest) in history
____________set current value as the new pBest
____End
____Choose the particle with the best fitness value of all the particles as the gBest
____For each particle
________Calculate particle velocity according equation (a)
________Update particle position according equation (b)
____End
While maximum iterations or minimum error criteria is not attained
在每一維粒子的速度都會被限制在一個最大速度Vmax，如果某一維更新後的速度超過用戶設定的Vmax，那麼這一維的速度就被限定為Vmax

『柒』 求粒子群演算法C#版

static void Main(string[] args)
{
PSO p = new PSO();
p.mains();
System.Console.ReadLine();
}
class PSO
{
const int S = 20; /*試驗次數 */
const int G = 2000; /*混合迭代次數*/
const int P = 40; /*個體總數*/
//#define c1 2.0 /*學習因子1*/
//#define c2 2.0 /*學習因子2*/
const int V = 30; /*個體維數*/
const double MAX = 5.12;
const double MIN = -5.12;
double D = MAX; /*蛙跳的最大值*/
int i1, i2, i3, i4;
int try_number = 0;
int try_max = 5;
double R;//0-1之間的隨機數,精度為1/10000
double Wmax = 0.9;
double Wmin = 0.4;
double PI = 3.14159265;
double W = 0.9;
double c1 = 2.0;/*學習因子1*/
double c2 = 2.0;/*學習因子2*/
double Vpso = 0.0;
double Tolerance = 0.0000001;//收斂精度
double c3 = 0.03;//擾動幅度
double e = 2.718281828459;//自然對數底數
int sm = 3;
int bz = 0;//擾動因子標志
public class Individal
{
public double[] d = new double[V];
public double fitness;
}
public class psom
{
public double[] sd = new double[V];
} public psom[] pso = new psom[P]; /*記錄每個離子各位的更新速度*/
public Individal px; /*全體中最好位置*/
public Individal[] indivial = new Individal[P]; /*全部個體*/
public Individal[] indiviala = new Individal[P]; /*全部個體——備份*/
public Individal tem; public PSO()
{
for (int i = 0; i < P; i++)
{
pso[i] = new psom();
indivial[i] = new Individal();
indiviala[i] = new Individal();
}
}
/*選擇測試函數為Sphere*/
/*選擇測試函數為Sphere*/
public double fitness(double[] a)
{
int i;
double sum = 0.0;
double sum1 = 0.0;
double s1 = 0.0, h1 = 0.0;
double[] x1 = new double[V + 1];
for (i = 0; i < V; i++) x1[i] = a[i];
for (i = 0; i < V; i++)
for (i = 0; i < V; i++)
sum = sum + (x1[i] * x1[i] - 10 * Math.Cos(2 * PI * x1[i]) + 10); return sum;
} /*對每一個個體初始化*/
public void init()
{
Random ran = new Random();
//R = ran.NextDouble(); int i, j, pmin = 0;
//srand((unsigned)time(NULL));
for (i = 0; i < P; i++)
{
for (j = 0; j < V; j++)
{
R = ran.NextDouble();
indivial[i].d[j] = R * (MAX - MIN) + MIN;
}
indivial[i].fitness = fitness(indivial[i].d);//計算初始適應值
indiviala[i] = indivial[i];//將個體復制給另一個序列
}
for (i = 0; i < P - 1; i++)
{
if (indivial[pmin].fitness > indivial[i + 1].fitness)
pmin = i + 1;//適應值小者最優
}
px = indivial[pmin];//最優者為小
}
/*按照適應度降序對全部個體進行排序和族群劃分*/ /*群組內更新*/
public void update()
{
int i, j, k, l, n;
double a;
double b;
W = Wmax - (double)(i2) * (Wmax - Wmin) / (double)(G);
for (i = 0; i < P; i++)
{ for (j = 0; j < V; j++)//更新粒子速度、位置
{//更新速度
pso[i].sd[j] = W * pso[i].sd[j] + c1 * R * (indiviala[i].d[j] - indivial[i].d[j]) + c2 * R * (px.d[j] - indivial[i].d[j]);
if (pso[i].sd[j] > D) pso[i].sd[j] = D; //D 最大速度
if (pso[i].sd[j] < -D) pso[i].sd[j] = -D;
indivial[i].d[j] = indivial[i].d[j] + pso[i].sd[j];//更新位置
}
a = fitness(indivial[i].d);//計算本次迭代的粒子適應值
// printf("old是%.16f",indivial[i].fitness);
indivial[i].fitness = a;
// printf("new是%.16f",indivial[i].fitness);
// getchar();
if (a < indiviala[i].fitness)
{
indiviala[i] = indivial[i];
if (indiviala[i].fitness < px.fitness)
px = indiviala[i];
}//比較粒子與前一次迭代的適應值 尋求最優者
}
}
public void report()
{
int i;
System.Console.WriteLine(px.fitness);

for (i = 0; i < P; i++)
{
// printf("%.16f\n",indivial[i].fitness);
// printf("%.16f\n",indiviala[i].fitness);
}
}
public void mains()
{
// int i1,i2;
//clock_t start, end;
double ave;
//FILE* f = fopen("result(SFLA).txt", "w");
//for(i4=0;i4<G+1;i4+=50)
//{
ave = 0.0;
//start = clock();
for (int i1 = 0; i1 < S; i1++)
{ init();
for (int i2 = 0; i2 < G; i2++)
{
update();
// report();
// getchar();
}
//
//
report();
ave = ave + px.fitness;
}
//end = clock();
ave = ave / S;
System.Console.WriteLine("平均極值為:");
System.Console.WriteLine(ave);
//System.Console.WriteLine("Interval=%.2fseconds\n", (double)(end - start) / ((double)CLOCKS_PER_SEC));
//printf("%d代為平均極值為%.16f\n",i4,ave);
//fprintf(f,"%d代為平均極值為%.16f\n",i4,ave);
// getchar();
//getchar();
}
}

『捌』 粒子群演算法（一）：粒子群演算法概述

  本系列文章主要針對粒子群演算法進行介紹和運用，並給出粒子群演算法的經典案例，從而進一步加深對粒子群演算法的了解與運用（預計在一周內完成本系列文章）。主要包括四個部分：

  粒子群演算法也稱粒子群優化演算法（Particle Swarm Optimization, PSO），屬於群體智能優化演算法，是近年來發展起來的一種新的進化演算法（Evolutionary Algorithm, EA）。 群體智能優化演算法主要模擬了昆蟲、獸群、鳥群和魚群的群集行為，這些群體按照一種合作的方式尋找食物，群體中的每個成員通過學習它自身的經驗和其他成員的經驗來不斷地改變搜索的方向。 群體智能優化演算法的突出特點就是利用了種群的群體智慧進行協同搜索，從而在解空間內找到最優解。
  PSO 演算法和模擬退火演算法相比，也是 從隨機解出發，通過迭代尋找最優解 。它是通過適應度來評價解的品質，但比遺傳演算法規則更為簡單，沒有遺傳演算法的「交叉」和「變異」，它通過追隨當前搜索到的最大適應度來尋找全局最優。這種演算法以其 容易實現、精度高、收斂快 等優點引起了學術界的重視，並在解決實際問題中展示了其優越性。

  在粒子群演算法中，每個優化問題的解被看作搜索空間的一隻鳥，即「粒子」。演算法開始時首先生成初始解，即在可行解空間中隨機初始化 粒子組成的種群 ，其中每個粒子所處的位置 ，都表示問題的一個解，並依據目標函數計算搜索新解。在每次迭代時，粒子將跟蹤兩個「極值」來更新自己， 一個是粒子本身搜索到的最好解 ，另一個是整個種群目前搜索到的最優解 。 此外每個粒子都有一個速度 ，當兩個最優解都找到後，每個粒子根據如下迭代式更新：

  其中參數 稱為是 PSO 的 慣性權重(inertia weight) ，它的取值介於[0,1]區間；參數 和 稱為是 學習因子(learn factor) ；而 和 為介於[0,1]之間的隨機概率值。
  實踐證明沒有絕對最優的參數，針對不同的問題選取合適的參數才能獲得更好的收斂速度和魯棒性，一般情況下 , 取 1.4961 ，而 採用 自適應的取值方法 ，即一開始令 ， 使得 PSO 全局優化能力較強 ；隨著迭代的深入，遞減至 ， 從而使得PSO具有較強的局部優化能力 。

  參數 之所以被稱之為慣性權重，是因為 實際 反映了粒子過去的運動狀態對當前行為的影響，就像是我們物理中提到的慣性。 如果 ，從前的運動狀態很少能影響當前的行為，粒子的速度會很快的改變；相反， 較大，雖然會有很大的搜索空間，但是粒子很難改變其運動方向，很難向較優位置收斂，由於演算法速度的因素，在實際運用中很少這樣設置。也就是說， 較高的 設置促進全局搜索，較低的 設置促進快速的局部搜索。