粒子群演算法求解旅行商問題TSP

寫在開頭：
最近師妹的結課作業問我，關於使用粒子群求解TSP問題的思路。我想了想，自己去年的作業用的是遺傳演算法，貌似有些關聯，索性給看了看程式碼。重新學習了一遍粒子群演算法，在這裡記錄一下，算是對知識的總結，鞏固一下。

正文部分

本文主要是使用粒子群來求解旅行商問題，即TSP問題，這裡主要講解程式碼和實現思路，原理會簡單帶過。詳細的具體原理，請讀者移步參考連結。
本文將從如下幾個方面進行描述：

• 1. 粒子群簡單介紹
• 2. TSP問題簡單介紹
• 3. JAVA程式碼實現
• 4. 執行結果展示
• 5. 原始碼下載連結

1. 粒子群簡單介紹

粒子群演算法簡稱PSO，它的基本思想是模擬鳥群的捕食行為。設想這樣一個場景：一群鳥在隨機搜尋食物。在這個區域裡只有一塊食物。所有的鳥都不知道食物在那裡。但是他們知道當前的位置離食物還有多遠。那麼找到食物的最優策略是什麼呢。最簡單有效的就是搜尋目前離食物最近的鳥的周圍區域。

PSO從這種模型中得到啟示並用於解決優化問題。PSO中，每個優化問題的解都是搜尋空間中的一隻鳥。我們稱之為“粒子”。所有的粒子都有一個由被優化的函式決定的適應值(fitness value)，每個粒子還有一個速度決定他們飛翔的方向和距離。然後粒子們就追隨當前的最優粒子在解空間中搜索。

PSO 初始化為一群隨機粒子(隨機解)。然後通過迭代找到最優解。在每一次迭代中，粒子通過跟蹤兩個”極值”來更新自己。第一個就是粒子本身所找到的最優解，這個解叫做個體極值pBest。另一個極值是整個種群目前找到的最優解，這個極值是全域性極值gBest。另外也可以不用整個種群而只是用其中一部分作為粒子的鄰居，那麼在所有鄰居中的極值就是區域性極值。

2. TSP問題簡單介紹

TSP問題（Travelling Salesman Problem）即旅行商問題，又譯為旅行推銷員問題、貨郎擔問題，是數學領域中著名問題之一。假設有一個旅行商人要拜訪n個城市，他必須選擇所要走的路徑，路徑的限制是每個城市只能拜訪一次，而且最後要回到原來出發的城市。路徑的選擇目標是要求得的路徑路程為所有路徑之中的最小值。

TSP問題是一個組合優化問題。該問題可以被證明具有NPC計算複雜性。TSP問題可以分為兩類，一類是對稱TSP問題（Symmetric TSP），另一類是非對稱問題（Asymmetric TSP）。

3. JAVA程式碼實現

首先整理一下具體實現邏輯，具體流程如下：

• 初始化種群
• 初始化慣性因子，即自身的交換序列
• 尋找每個粒子的歷代中的最優解
• 尋找全域性最優解
• 計算每個個體的自身最優解的交換序列
• 計算每個個體對全域性最優解的交換序列

• 進化——其實就是交換
  交換序列對由3部分組成：自身交換序列、自身最優解的交換序列、全域性最優解的交換序列。

• 列印結果

使用面向物件的思想，首先定義一個city.java城市實體類，有id、name、經度和緯度屬性。

public class City {
    private int mId;       //城市編號
    private String mName;  //城市名稱
    private float mLongitude; //經度
    private float mLatitude; //經度

    public City(int id, String name, float longitude, float latitude) {
        mId = id;
        mName = name;
        mLongitude = longitude;
        mLatitude = latitude;
    }

定義一個粒子實體類Unit.java，包括行走城市的路徑path和適應度fitness。並新增兩個方法，輸出城市序列和計算自己的適應度。

private int[] mPath;  //行走的城市路徑，儲存城市編號
    private int mFitness; //適應度值，為當前個體走這個路徑的總距離。越小越好。

    public Unit(int[] path) {
        mPath = path;
        mFitness = calculateFitness();
    }
        public void printPath() {
        if (mPath == null) {
            System.out.println("mPath為null,當前個體的路徑為空");
        } else {
            for (int i = 0; i < mPath.length - 1; i++) {
                System.out.print(CityLab.getInstance().getmCities().get(mPath[i]).getName() + "——》");
            }
            System.out.println(CityLab.getInstance().getmCities().get(mPath[mPath.length - 1]).getName());
        }
    }

    public void upDateFitness() {
        this.mFitness = calculateFitness();
    }

    /**
     * 計算當前路徑的適應值，即為路徑長度
     */
    public int calculateFitness() {
        //根據經緯度計算距離
        //近似計算：0.00001度，距離相差約1米；0.01，距離相差1000米.1度，距離相差100km
        int distance = 0;  //單位千米(km)
        int n = mPath.length;
        for (int i = 1; i < n; i++) {
            City c1 = CityLab.getInstance().getmCities().get(mPath[i - 1]);
            City c2 = CityLab.getInstance().getmCities().get(mPath[i]);
            distance += Math.sqrt(Math.pow(100 * (c1.getLatitude() - c2.getLatitude()), 2) + Math.pow(100 * (c1.getLongitude() - c2.getLongitude()), 2));
        }
        distance += Math.sqrt(Math.pow(100 * (CityLab.getInstance().getmCities().get(mPath[0]).getLatitude() - CityLab.getInstance().getmCities().get(mPath[n - 1]).getLatitude()), 2) + Math.pow(100 * (CityLab.getInstance().getmCities().get(mPath[0]).getLongitude() - CityLab.getInstance().getmCities().get(mPath[n - 1]).getLongitude()), 2));
        return distance;
    }

接下來，定義交換對實體類，SO.java，包含兩個屬性x和y，用來做交換的。

public class SO {
    private int x;
    private int y;

    public SO(int x, int y) {
        this.x = x;
        this.y = y;
    }

最後使用單例模式來管理所有的城市列表CityLab.java

/**
 * 城市單例
 */
public class CityLab {
    private static CityLab mCityLab = null;
    private ArrayList<City> mCities = new ArrayList<>();

    private CityLab() {
        //讀取檔案，新增城市資訊到mcities中
        String filename = "d://city.csv";
        String strbuff;
        BufferedReader data = null;
        try {
            data = new BufferedReader(new InputStreamReader(
                    new FileInputStream(filename)));
            int id = 1;
            while (true) {
                //讀取一行資料:北京,116.41667,39.91667
                strbuff = data.readLine();
                if (strbuff == null) {
                    break;
                }
                String[] arr = strbuff.split(",");
                City c = new City(id, arr[1], valueOf(arr[2]), valueOf(arr[3]));
                id++;
                mCities.add(c);
            }
        } catch (FileNotFoundException e) {
            e.printStackTrace();
        } catch (IOException e) {
            e.printStackTrace();
        }

    }

    public static CityLab getInstance() {
        if (mCityLab == null) {
            mCityLab = new CityLab();
        }
        return mCityLab;
    }

    public ArrayList<City> getmCities() {
        return mCities;
    }
}

實體類我們定義好了，接下來看具體實現粒子群求解TSP問題的過程。
具體PSO.java類如下：

/**
 * 粒子群求解TSP旅行商問題
 * 更新公式：Vii=wVi+ra(Pid-Xid)+rb(Pgd-Xid)
 */
public class PSO {
    private int scale; //種群規模
    private ArrayList<Unit> mUnits = new ArrayList<>();   //粒子群

    private int MAX_GEN;// 迭代次數
    private float w;    //慣性權重
    private int cityNum = 25;    //城市數量

    private HashMap<Integer, Unit> Pd = new HashMap<>();    //一顆粒子歷代中出現最好的解
    private Unit Pgd;   // 整個粒子群經歷過的的最好的解，每個粒子都能記住自己搜尋到的最好解
    private int bestT;// 最佳出現代數

    private ArrayList<ArrayList<SO>> listV = new ArrayList<>();    //自身交換序列，即所謂的慣性因子

    Random random = new Random();

    /**
     * 構造方法
     *
     * @param scale
     * @param MAX_GEN
     * @param w
     */
    public PSO(int scale, int MAX_GEN, float w) {
        this.scale = scale;
        this.MAX_GEN = MAX_GEN;
        this.w = w;
    }

    /**
     * 初始化引數配置
     */
    private void init() {
        cityNum = CityLab.getInstance().getmCities().size();
        random = new Random(System.currentTimeMillis());

    }

    /**
     * 初始化種群
     */
    private void initGroup() {
        for (int k = 0; k < scale; k++) {
            int[] path = new int[cityNum];
            for (int i = 0; i < cityNum; ) {
                //隨機生成一個城市路徑
                //int s = random.nextInt(max)%(max-min+1) + min;
                int s = random.nextInt(65535) % cityNum;
                int j;
                for (j = 0; j < i; j++) {
                    if (s == path[j]) {
                        break;
                    }
                }
                if (i == j) {
                    path[i] = s;
                    i++;
                }
            }
            Unit unit = new Unit(path);
            mUnits.add(unit);
        }
    }

    /**
     * 初始化自身的交換序列即慣性因子
     */
    private void initListV() {
        for (int i = 0; i < scale; i++) {
            ArrayList<SO> list = new ArrayList<>();
            int n = random.nextInt(cityNum - 1) % (cityNum);    //隨機生成一個數，表示當前粒子需要交換的對數
            for (int j = 0; j < n; j++) {
                //生成兩個不相等的城市編號x,y
                int x = random.nextInt(cityNum - 1) % (cityNum);
                int y = random.nextInt(cityNum - 1) % (cityNum);
                while (x == y) {
                    y = random.nextInt(cityNum - 1) % (cityNum);
                }

                //x不等於y
                SO so = new SO(x, y);
                list.add(so);
            }
            listV.add(list);
        }
    }

    public void solve() {
        initGroup();

        initListV();

        //挑選最好的個體
        for (int i = 0; i < scale; i++) {
            Pd.put(i, mUnits.get(i));
        }
        Pgd = Pd.get(0);
        for (int i = 0; i < scale; i++) {
            if (Pgd.getFitness() > Pd.get(i).getFitness()) {
                Pgd = Pd.get(i);
            }
        }
        System.out.println("初始化最好結果為：" + Pgd.getFitness());
        Pgd.printPath();

        // 進化
        evolution();

        // 列印
        System.out.println("==================最後粒子群=====================");

        System.out.println("最佳長度出現代數：");
        System.out.println(bestT);
        System.out.println("最佳長度");
        System.out.println(Pgd.getFitness());
        System.out.println("最佳路徑：");
        Pgd.printPath();

    }

    /**
     * 進化
     */
    private void evolution() {
        for (int t = 0; t < MAX_GEN; t++) {
            for (int k = 0; k < scale; k++) {
                ArrayList<SO> vii = new ArrayList<>();
                //更新公式：Vii=wVi+ra(Pid-Xid)+rb(Pgd-Xid)
                //第一部分，自身交換對
                int len = (int) (w * listV.get(k).size());
                for (int i = 0; i < len; i++) {
                    vii.add(listV.get(k).get(i));
                }

                //第二部分，和當前粒子中出現最好的結果比較，得出交換序列
                //ra(Pid-Xid)
                ArrayList<SO> a = minus(mUnits.get(k).getPath(), Pd.get(k).getPath());
                float ra = random.nextFloat();
                len = (int) (ra * a.size());
                for (int i = 0; i < len; i++) {
                    vii.add(a.get(i));
                }

                //第三部分，和全域性最優的結果比較，得出交換序列
                //rb(Pgd-Xid)
                ArrayList<SO> b = minus(mUnits.get(k).getPath(), Pgd.getPath());
                float rb = random.nextFloat();
                len = (int) (rb * b.size());
                for (int i = 0; i < len; i++) {
                    vii.add(b.get(i));
                }

                listV.remove(0);    //移除當前，加入新的
                listV.add(vii);

                //執行交換，生成下一個粒子
                exchange(mUnits.get(k).getPath(), vii);
            }

            //更新適應度的值，並挑選最好的個體
            for (int i = 0; i < scale; i++) {
                mUnits.get(i).upDateFitness();
                if (Pd.get(i).getFitness() > mUnits.get(i).getFitness()) {
                    Pd.put(i, mUnits.get(i));
                }
                if (Pgd.getFitness() > Pd.get(i).getFitness()) {
                    Pgd = Pd.get(i);
                    bestT = t;
                }
            }
            //列印當前代的結果
            if (t % 100 == 0) {
                // 列印
                System.out.println("--------第"+t+"代的最佳結果為-----------");
                System.out.println(Pgd.getFitness());
                System.out.println("最佳路徑：");
                Pgd.printPath();
            }
        }
    }

    /**
     * 執行交換，更新粒子
     *
     * @param path
     * @param vii  儲存的是需要交換的下標對
     */
    private void exchange(int[] path, ArrayList<SO> vii) {
        int tmp;
        for (SO so : vii) {
            tmp = path[so.getX()];
            path[so.getX()] = path[so.getY()];
            path[so.getY()] = tmp;
        }
    }


    /**
     * 生成交換對，把a變成和b一樣，返回需要交換的下標對列表
     *
     * @param a
     * @param b
     * @return
     */
    private ArrayList<SO> minus(int[] a, int[] b) {
        int[] tmp = a.clone();
        ArrayList<SO> list = new ArrayList<>();
        int index = 0;
        for (int i = 0; i < b.length; i++) {
            if (tmp[i] != b[i]) {
                //在tmp中找到和b[i]相等的值，將下標儲存起來
                for (int j = i + 1; j < tmp.length; j++) {
                    if (tmp[j] == b[i]) {
                        index = j;
                        break;
                    }
                }
                SO so = new SO(i, index);
                list.add(so);
            }
        }
        return list;
    }

    public static void main(String[] args) {
        PSO pso = new PSO(30, 5000, 0.5f);
        pso.init();
        pso.solve();

    }
}

註釋寫得很詳細，這裡不再贅述，各位可以直接執行程式碼，通過debug的方式，這樣我覺得很容易理解的。看原始碼是學習程式設計最快的方式

最後感謝大家的觀看，如果有問題請評論區留言，本人將盡力解答。

4. 執行結果展示

5. 原始碼下載連結

宣告
在參考下面連結的基礎上，對程式碼進行重構，使用java面向物件的思想，是程式碼可讀性大大增強