Question

粒子群算法（PSO）是一种进化算法。它与大名鼎鼎的遗传算法非常类似，可以用来解决一些优化问题。大家知道，一些优化问题（比如旅行商问题TSP）都属于NP问题。它们的时间复杂度可能会达到O(n!)或者O(2ⁿ)，这种在多项式时间内不可解的问题总是会让人望而生畏。而以PSO算法为代表的进化计算，往往可以将这些NP问题，转变为一个多项式问题。但这种转变是有代价的，进化算法往往都不保证你可以从结果中得到最优解。我这么说，也许就有人会问了，这个算法都不能保证得到最优解，那有什么用呢？其实，在生活中的很多场景下，并不是特别需要最优解，我们更加希望得到的是一个满意解。比如说，去水果店买西瓜，店里可能放着一大堆西瓜，每个人都想挑一个最好的。但你想拿到最好的那个西瓜必须得挨个检查过去，并且还得认真做好记录才行。我相信，没有一个人会这么买西瓜，因为成本太高了。对于大部分人来说，更倾向于在表面上挑几个顺眼的看看，如果还过得去，也就下手了。这也就是说只要这个结果不要差得太离谱就行了。

既然最优的方案很难得到，那么我们就想办法以很低的成本获得一个还算过得去的方案，也不失为一计良策。在后面给出的小案例中大家也可以看到，在很多情况下，虽然进化算法无法让你获得最优解，也无法证明它得到的解与最优解到底有多少差距，但实际中，通过进化算法搜索到的满意解很可能与最优解已经非常接近了。

7.13.1　什么是粒子群算法

粒子群优化算法（PSO）是一种进化计算技术，最早由Kenny与Eberhart于1995年提出。它源于对鸟群捕食行为的研究，与遗传算法相似，是一种基于迭代的优化算法，广泛应用于函数优化和神经网络训练等方面。与遗传算法相比，PSO算法的实现简单得多，参数配置也相对较少，对使用人员的经验要求不高，因此更加易于实际工程应用。

从日常生活的观察中可以知道，鸟类的觅食往往会表现成群体特性。如果在地上有一小撮食物，那么鸟群很可能就会聚集在这一堆食物旁边。如果其中一只小鸟发现了另外一堆更丰盛的食物，那它可能会离群飞向更丰盛的食物，而这有可能带动整个鸟群一起飞向新的地点。当然了，在整个种群中，难免会出现几只特别有“个性”的小鸟，它们不喜欢太热闹的地方，当整个种群迁移时，它们不会跟着种群走，或者自己散步，或者自行游荡。

粒子群算法正是对上述过程的模拟。在程序中，我们可以模拟大量的小鸟，小鸟的觅食点正是要求解的问题的解。解越是优秀，意味着食物越是丰盛，因此，模拟的小鸟会从自己的位置出发以一定的速度向最优点的方向移动。在移动过程中，任何一只小鸟都有可能发现更好的解，这又会进一步影响群体的行为。就这样如此反复迭代，最终，将得到一个不错的答案。

7.13.2　粒子群算法的计算过程

粒子群算法的大体步骤如下：

1. 初始化所有粒子，粒子的位置随机生成。计算每个粒子当前的适应度，并将此设为当前粒子的个体最优值（记为pBest）。

2. 所有粒子将自己的个体最优值发送给管理者Master。Master获得所有粒子的信息后，筛选出全局最优的解（记为gBest）。

3. Master将gBest通知所有粒子，所有粒子便知道全局最优点的位置。

4. 接着，所有粒子根据自己的pBest和全局gBest，更新自己的速度，在有了速度后，再更新自己的位置。
v_k+1=c₀×rand()×v_k+c₁×rand()×(pbest_k−x_k)+c₂×rand()×(gbest_k−x_k)
x_k+1=x_k+v_k+1
其中，rand()函数产生一个0，1之间的随即数。c₀=1，c₁=2，c₂=2，k表示进化的代数。v_k表示当前速度，pbest_k和gbest_k表示个体最优解和全局最优解。当然，对于每一个维度上的速度分量，我们可以为它限定一个最大值。确保“小鸟”不会飞得太快，错过了重要的信息。

5. 如果粒子产生了新的个体最优点，则发送给Master，在此，转到步骤2。

整体过程的示意图如图7.2所示。

图7-2　PSO算法示意图

从这个计算步骤中可以看到，计算过程拥有一定的随机性。但由于我们可以启用大量的例子，因此其计算效果在统计学意义上是稳定的。在这个标准的粒子群算法中，由于所有粒子都会向全局最优靠拢，因此，其跳出局部最优的能力并不算太强。因此，我们也可以想办法对标准的粒子群算法进行一些合理的改进。比如，允许各个粒子随机移动，甚至逆向移动来试图突破局部最优。在这里为简单起见，我不打算做这些复杂的实现。

7.13.3　粒子群算法能做什么

粒子群算法能为我们做些什么呢？它应用最多的场景是进行最优化计算。实际上，以粒子群算法为代表的进化计算，可以说是最优化方法中的通用方法。几乎一切最优化问题都可以通过这种随机搜索的模式解决，其成本低、难度小、效果好，因此颇受欢迎。

下面，就让我们来探讨一个典型的优化问题：

假设现在有400万资金，要求4年内使用完。若在第1年使用x万元，则可以得到效益万元（效益不能再使用），当年不用的资金可存入银行，年利率为10%。尝试制订出资金的使用规划，使4年效益之和最大。

很明显，对于这类问题，不同的方案得到的结果可能会有很大的差异。比如，若第一年把400万元全部用完，则总效益为万元；若前3年均不用而存入银行，第4年把本金和利息全部用完，则总效益为万元，显然优于第一种方案。

如果我们将此问题转为一般化的优化问题，则可以得到以下方程组，如图7.3所示。

图7.3　一般化的约束问题

其中，x1、x2、x3、x4分别表示第1、2、3、4年使用的资金。使用拉格朗日乘子法对此方程组进行求解，可以得到第一年使用86.19万元、第2年使用104.29万元，第3年使用126.19万元，第4年使用152.69万元为这个问题的最优解，此时总效益达43.09万元。

由于求解过程过于复杂，使用拉格朗日乘子法时，需要对先后12个未知数和方程进行联立求解，比较难以实现。由于求解过程与我们讨论的主题无关，所以在这里不再给出。

对于类似的优化问题，正是粒子群算法的涉猎范围。当使用粒子群算法时，我们可以先随机给出若干个满足提交的资金规划方案。接着，根据粒子群的演化公式，不断调整各个粒子的位置（粒子的每一个位置代表一套方案），逐步探索更优的方案。

7.13.4　使用Akka实现粒子群

现在，我们已经知道粒子群的原理，并且有了一个较为复杂的优化问题等待我们求解。接下来，就需要开动脑筋，使用Akka来实现一个简单的粒子群，来解决这个优化问题了。

使用Actor的模式与粒子群算法之间有着天生切合度。粒子群算法由于涉及到多个甚至是极其大量的粒子参与运算，因此它隐含着并行计算的模式。其次，从直观上我们也可以知道，粒子群算法的求解精度或者说求解的质量，与参与运算的粒子有着直接的关系。很显然，参与运算的粒子数量越多，得到的解自然也就越精确。

如果我们使用传统的多线程方式实现粒子群，一个最大的问题就是线程的数量可能是非常有限的。在当前这种应用场景中，我们希望可以拥有数万，甚至数十万的粒子，以提高计算精度，但众所周知，在一台计算机上运行数万个线程基本是不可能的，就算可以，系统的性能也会大打折扣。因此，使用多线程的模型无法很好地和粒子群的实现相融合。

但Akka的Actor的模型则不同。由于多个Actor可以复用一个线程，而Actor本身作为轻量级的并发执行单元可以有极其大量的存在。因此，我们就可以使用Actor来模拟整个粒子群计算的场景。下面就让我们仔细看一下系统的实现。

首先，我们需要两个表示pBest和gBest的消息类型，用于在多个Actor之间传递个体最优和全局最优。

01 public final class GBestMsg  {
02   final PsoValue value;
03   public GBestMsg(PsoValue v){
04     value=v;
05   }
06   public PsoValue getValue() {
07     return value;
08   }
09 }
10
11 public final class PBestMsg {
12   final PsoValue value;
13   public PBestMsg(PsoValue v){
14     value=v;
15   }
16
17   public PsoValue getValue() {
18     return value;
19   }
20
21   public String toString(){
22     return value.toString();
23   }
24 }

上述代码中，GBestMsg（代码第1行）表示携带全局最优解的消息。而PBestMsg（代码第11行）表示携带个体最优的消息。它们都使用PsoValue来表示一个可行的解。

在PsoValue中，主要包括两个信息，第一是表示投资规划的方案，即每一年分别需要投资多少钱；第二是这个投资方案的总收益：

01 public final class PsoValue {
02   final double value;
03   final List＜Double＞ x;
04   public PsoValue(double v,List＜Double＞ x){
05     value=v;
06     List＜Double＞ b=new ArrayList＜Double＞(5);
07     b.addAll(x);
08     this.x=Collections.unmodifiableList(b);
09   }
10   public double getValue(){
11     return value;
12   }
13   public List＜Double＞ getX(){
14     return x;
15   }
16
17   public String toString(){
18     StringBuffer sb=new StringBuffer();
19     sb.append("value:").append(value).append("\n")
20     .append(x.toString());
21     return sb.toString();
22   }
23 }

上述代码中，数组x中，x[1]、x[2]、x[3]、x[4]分别表示第1年、第2年、第3年和第4年的投资额。这里为了方便起见，我忽略了x[0]（它在我们的程序中是没有作用的）。成员变量value表示这组投资方案的收益值。

因此，根据需求x与value之间的关系如下代码所示：

1 public class Fitness {
2   public static double fitness(List＜Double＞ x){
3     double sum=0;
4     for(int i=1;i＜x.size();i++){
5       sum+=Math.sqrt(x.get(i));
6     }
7     return sum;
8   }
9 }

上述代码定义的fitness()函数返回了给定投资方案的适应度。适应度也就是投资的收益，我们自然应该更倾向于选择适应度更高的投资方案。在这里适应度=。

有了这些基础工具，我们就可以来实现简单的粒子（这里我把它叫作Bird）了。

对于基本粒子，我们需要定义以下成员变量：

1 public class Bird extends UntypedActor {
2   private final LoggingAdapter log = Logging.getLogger(getContext().system(), this);
3   private PsoValue pBest=null;
4   private PsoValue gBest=null;
5   private List＜Double＞ velocity =new ArrayList＜Double＞(5);
6   private List＜Double＞ x =new ArrayList＜Double＞(5);
7   private Random r = new Random();

上述代码中，pBest和gBest分别表示个体最优和全局最优，velocity表示粒子在各个维度上的速度（在当前案例中，每一年的投资额就可以认为是一个维度，因此系统有4个维度）。x表示投资方案，即每一年的投资额。由于在粒子群算法中，需要使用随机数，因此，这里定义了r。

当一个粒子被创建时，我们需要初始化粒子的当前位置。粒子的每一个位置都代表一个投资方案，下面的代码展示了粒子的初始化逻辑：

01 @Override
02 public void preStart(){
03   for(int i=0;i＜5;i++){
04     velocity.add(Double.NEGATIVE_INFINITY);
05     x.add(Double.NEGATIVE_INFINITY);
06   }
07   //x1＜=400
08   x.set(1, (double)r.nextInt(401));
09
10   //x2＜=440-1.1*x1
11   double max=400-1.1*x.get(1);
12   if(max＜0)max=0;
13   x.set(2, r.nextDouble()*max);
14
15   //x3＜=484-1.21*x1-1.1*x2
16   max=484-1.21*x.get(1)-1.1*x.get(2);
17   if(max＜=0)max=0;
18   x.set(3, r.nextDouble()*max);
19
20   //x4＜= 532.4-1.331*x1-1.21*x2-1.1*x3
21   max=532.4-1.331*x.get(1)-1.21*x.get(2)-1.1*x.get(3);
22   if(max＜=0)max=0;
23   x.set(4, r.nextDouble()*max);
24
25   double newFit=Fitness.fitness(x);
26   pBest=new PsoValue(newFit,x);
27   PBestMsg pBestMsg=new PBestMsg(pBest);
28   ActorSelection selection = getContext().actorSelection("/user/masterbird");
29   selection.tell(pBestMsg, getSelf());
30 }

由于在当前案例中，每一年的投资额度是有条件约束的，比如第一年的投资额不能超过400万（第7～8行），而第2年的投资上限是440万（假设第一年全部存银行，代码第10～13行），依此类推。粒子初始化时，随机生成一组满足基本约束条件的投资组合，并计算它的适应度（第25行）。初始的投资方案自然也就作为当前的个体最优，并发送给Master（第29行）。

当Master计算出当前全局最优后，会将全局最优发送给每一个粒子，粒子根据全局最优更新自己的运行速度，并更新自己的速度以及当前位置。

01 @Override
02 public void onReceive(Object msg) {
03   if (msg instanceof GBestMsg) {
04     gBest=((GBestMsg) msg).getValue();
05     //更新速度
06     for(int i=1;i＜velocity.size();i++){
07       updateVelocity(i);
08     }
09     //更新位置
10     for(int i=1;i＜x.size();i++){
11       updateX(i);
12     }
13     validateX();
14     double newFit=Fitness.fitness(x);
15     if(newFit＞pBest.value){
16       pBest=new PsoValue(newFit,x);
17       PBestMsg pBestMsg=new PBestMsg(pBest);
18       getSender().tell(pBestMsg, getSelf());
19     }
20   }
21   else{
22     unhandled(msg);
23   }
24 }

上述代码中，粒子接收到了全局最优（代码第4行），接着根据粒子群的标准公式更新自己的速度（第6～8行）。接着，根据速度，更新自己的位置（第10～12行）。由于当前问题是有约束的，也就是说解空间并不是随意的。粒子很可能在更新位置后，跑出了合理的范围之外，因此，还有必要进行有效性检查（第13行）。

在更新完成后，就可以计算新位置的适应度，如果产生了新的个体最优，就将其发送给Master（第15～19行）。

在当前案例中，速度和位置的更新是依据标准的粒子群实现，如下：

01 public double updateVelocity(int i){
02   double v= Math.random()*velocity.get(i)
03       +2*Math.random()*(pBest.getX().get(i)-x.get(i))
04       +2*Math.random()*(gBest.getX().get(i)-x.get(i));
05   v=v＞0? Math.min(v, 5): Math.max(v, -5);
06   velocity.set(i, v);
07   return v;
08 }
09
10 public double updateX(int i){
11   double newX=x.get(i)+velocity.get(i);
12   x.set(i, newX);
13   return newX;
14 }

上述代码中updateVelocity()和updateX()分别更新了粒子的速度和位置。位置的更新依赖于当前的速度（第11行）。

由于每一年的投资都是有限额的，因此，要避免粒子跑到合理空间之外，下面的代码强制将粒子约束中合理的区间中。

01 public void validateX(){
02   if(x.get(1)＞400){
03     x.set(1, (double)r.nextInt(401));
04   }
05
06   //x2
07   double max=400-1.1*x.get(1);
08   if(x.get(2)＞max || x.get(2)＜0){
09     x.set(2, r.nextDouble()*max);
10   }
11   //x3
12   max=484-1.21*x.get(1)-1.1*x.get(2);
13   if(x.get(3)＞max || x.get(3)＜0){
14     x.set(3, r.nextDouble()*max);
15   }
16   //x4
17   max=532.4-1.331*x.get(1)-1.21*x.get(2)-1.1*x.get(3);
18   if(x.get(4)＞max || x.get(4)＜0){
19     x.set(4, r.nextDouble()*max);
20   }
21 }

上述代码分别对x1、x2、x3、x4进行约束，一旦发现粒子跑出了定义范围就将它进行随机化。

此外，我们还需要一只MasterBird，用于管理和通知全局最优。

01 public class MasterBird extends UntypedActor {
02   private final LoggingAdapter log = Logging.getLogger(getContext().system(), this);
03   private PsoValue gBest=null;
04
05   @Override
06   public void onReceive(Object msg) {
07     if (msg instanceof PBestMsg) {
08       PsoValue pBest = ((PBestMsg) msg).getValue();
09       if(gBest==null || gBest.value ＜ pBest.value){
10         //更新全局最优，通知所有粒子
11         System.out.println(msg+"\n");
12         gBest=pBest;
13         ActorSelection selection = getContext().actorSelection("/user/bird_*");
14         selection.tell(new GBestMsg(gBest), getSelf());
15       }
16     }
17     else{
18       unhandled(msg);
19     }
20   }
21 }

上述代码定义了MasterBird。当它收到一个个体最优的解时，会将其与全局最优进行比较，如果产生了新的全局最优，就更新这个全局最优并通知所有的粒子（第12～14行）。

好了，现在万事俱备只欠东风。下面就是主函数：

01 public class PSOMain {
02   public static final int BIRD_COUNT = 100000;
03   public static void main(String[] args) {
04     ActorSystem system = ActorSystem
05         .create("psoSystem", ConfigFactory.load("samplehello.conf"));
06     system.actorOf(Props.create(MasterBird.class), "masterbird");
07     for (int i = 0; i ＜ BIRD_COUNT; i++) {
08       system.actorOf(Props.create(Bird.class), "bird_" + i);
09     }
10   }
11 }

上述代码定义了粒子总数，这里是10万个粒子。接着创建一个MasterBird Actor（第6行），和10万个bird（第7～9行）。

执行上述代码，运行一小段时间，你就可以得到如下输出（截取部分）：

value:36.15412875487459
[-Infinity, 168.0, 18.786423873345715, 102.1742923174793, 76.5657638235272]
value:37.88452477135976
[-Infinity, 64.0, 87.66774733441137, 37.976681047619195, 206.17791816445362]
....
value:42.240797528048176
[-Infinity, 113.0, 42.37168995110633, 141.70570102409184, 174.16812834843475]
....
value:43.01934824083668
[-Infinity, 76.0, 112.89557345993592, 133.29270155682005, 147.16289926594942]

上述输出表示，当粒子群随机初始化时，最优解为36.15万元，但随着粒子的搜索，这个投资方案被逐步优化，由37.88万一直上升到43.02万元。根据我们前面的求解，我们知道这个投资方案的最优结果是43.09万元，可以看到，粒子群的搜索结果和全局最优已经非常接近了。

当然了，由于粒子群算法的随机性，每次执行结果可能并不一样，这意味着有时候，你可能会求得更好的解，或者得到一个稍差一些的解，但其偏差不会相差太远。