智能算法迭代局部搜索（Iterated Local Search, ILS）详解

Posted 2020-11-20 infroad

迭代局部搜索（Iterated Local Search, ILS）

喜欢的话可以扫码关注我们的公众号哦，更多精彩尽在微信公众号【程序猿声】

00 目录

• 局部搜索算法
• 简单局部搜索
• 迭代局部搜索

01 局部搜索算法

1.1 什么是局部搜索算法？

局部搜索是解决最优化问题的一种启发式算法。因为对于很多复杂的问题，求解最优解的时间可能是极其长的。因此诞生了各种启发式算法来退而求其次寻找次优解，局部搜索就是其中一种。它是一种近似算法（Approximate algorithms）。

局部搜索算法是从爬山法改进而来的。简单来说，局部搜索算法是一种简单的贪心搜索算法，该算法每次从当前解的临近解空间中选择一个最优解作为当前解，直到达到一个局部最优解。局部搜索从一个初始解出发，然后搜索解的邻域，如有更优的解则移动至该解并继续执行搜索，否则返回当前解。

1.2 算法思想过程

局部搜索会先从一个初始解开始，通过邻域动作。产生初始解的邻居解，然后根据某种策略选择邻居解。一直重复以上过程，直到达到终止条件。

不同局部搜索算法的区别就在于：邻域动作的定义以及选择邻居解的策略。这也是决定算法好坏的关键之处。

1.3 什么又是邻域动作？

其实邻域动作就是一个函数。那么，通过这个函数，对当前的最优解s，产生s对应的邻居解的一个集合。比如：

对于一个bool型问题，其当前解为：s = 1001，当将邻域动作定义为翻转其中一个bit时，得到的邻居解的集合N(s)={0001,1101,1011,1000}，其中N(s) ∈ S。同理，当将邻域动作定义为互换相邻bit时，得到的邻居解的集合N(s)={0101,1001,1010}.

02 简单局部搜索

在开始我们的迭代局部搜索之前，还是先来给大家科普几个简单局部搜索算法。他们也是基于个体的启发式算法（Single solution）。

2.1 爬山法（HILL-CLIMBING）

干货 | 用模拟退火(SA, Simulated Annealing)算法解决旅行商问题

2.2 模拟退火（SIMULATED ANNEALING）

干货 | 用模拟退火(SA, Simulated Annealing)算法解决旅行商问题

2.3 模拟退火（SIMULATED ANNEALING）

干货|十分钟快速复习禁忌搜索(c++版)

干货 | 十分钟掌握禁忌搜索算法求解带时间窗的车辆路径问题(附C++代码和详细代码注释)

03 迭代局部搜索（Iterated Local Search, ILS）

3.1 介绍

迭代局部搜索属于探索性局部搜索方法（EXPLORATIVE LOCAL SEARCH METHODS）的一种。它在局部搜索得到的局部最优解上，加入了扰动，然后再重新进行局部搜索。

3.2 过程描述

注：下文的局部搜索(或者LocalSearch)指定都是简单局部搜索。指上文介绍的三种中的任意一种。

1. 从初始解s中进行局部搜索，找到一个局部最优解s1。
2. 扰动s1，获得新的解s2。
3. 从新解s2中进行局部搜索，再次找到一个局部最优解s3。
4. 基于判断策略，对s3好坏进行判断。选择接受s3作为新解或者回退s2。
5. 直到找到满足条件的最优解，不然跳回第二步。

其图解如下：

image

伪代码如下：

image

关于其中的接受判断准则，这里采用了模拟退火中的概率函数：

image

04 代码时间

以下代码还是用于求解TSP旅行商问题。至于什么是旅行商问题，读者可以从爬山算法那篇文章了解。

  1////////////////////////
  2//TSP问题 迭代局部搜索求解代码
  3//基于Berlin52例子求解
  4//作者：infinitor
  5//时间：2018-04-12
  6////////////////////////
  7
  8
  9#include <iostream>
 10#include <cmath>
 11#include <stdlib.h>
 12#include <time.h>
 13#include <vector>
 14#include <windows.h>
 15#include <memory.h>
 16#include <string.h>
 17#include <iomanip>
 18
 19#define DEBUG
 20
 21using namespace std;
 22
 23#define CITY_SIZE 52 //城市数量
 24
 25
 26//城市坐标
 27typedef struct candidate
 28{
 29    int x;
 30    int y;
 31}city, CITIES;
 32
 33//优化值
 34int **Delta; 
 35
 36//解决方案
 37typedef struct Solution
 38{
 39    int permutation[CITY_SIZE]; //城市排列
 40    int cost;                        //该排列对应的总路线长度
 41}SOLUTION;
 42// 计算邻域操作优化值 
 43int calc_delta(int i, int k, int *tmp, CITIES * cities);
 44
 45//计算两个城市间距离
 46int distance_2city(city c1, city c2);
 47
 48//根据产生的城市序列，计算旅游总距离
 49int cost_total(int * cities_permutation, CITIES * cities);
 50
 51//获取随机城市排列, 用于产生初始解
 52void random_permutation(int * cities_permutation);
 53
 54//颠倒数组中下标begin到end的元素位置, 用于two_opt邻域动作
 55void swap_element(int *p, int begin, int end);
 56
 57//邻域动作 反转index_i <-> index_j 间的元素
 58void two_opt_swap(int *cities_permutation, int *new_cities_permutation, int index_i, int index_j);
 59
 60//本地局部搜索，边界条件 max_no_improve
 61void local_search(SOLUTION & best, CITIES * cities, int max_no_improve);
 62
 63//判断接受准则
 64bool AcceptanceCriterion(int *cities_permutation, int *old_cities_permutation, CITIES * p_cities);
 65
 66//将城市序列分成4块，然后按块重新打乱顺序。
 67//用于扰动函数
 68void double_bridge_move(int *cities_permutation, int * new_cities_permutation);
 69
 70//扰动
 71void perturbation(CITIES * cities, SOLUTION &best_solution, SOLUTION &current_solution);
 72
 73//迭代搜索
 74void iterated_local_search(SOLUTION & best, CITIES * cities, int max_iterations, int max_no_improve);
 75
 76// 更新Delta 
 77void Update(int i, int k,  int *tmp, CITIES * cities);
 78
 79//城市排列
 80int permutation[CITY_SIZE];
 81//城市坐标数组
 82CITIES cities[CITY_SIZE];
 83
 84
 85//berlin52城市坐标，最优解7542好像
 86CITIES berlin52[CITY_SIZE] = { { 565,575 },{ 25,185 },{ 345,750 },{ 945,685 },{ 845,655 },
 87{ 880,660 },{ 25,230 },{ 525,1000 },{ 580,1175 },{ 650,1130 },{ 1605,620 },
 88{ 1220,580 },{ 1465,200 },{ 1530,5 },{ 845,680 },{ 725,370 },{ 145,665 },
 89{ 415,635 },{ 510,875 },{ 560,365 },{ 300,465 },{ 520,585 },{ 480,415 },
 90{ 835,625 },{ 975,580 },{ 1215,245 },{ 1320,315 },{ 1250,400 },{ 660,180 },
 91{ 410,250 },{ 420,555 },{ 575,665 },{ 1150,1160 },{ 700,580 },{ 685,595 },
 92{ 685,610 },{ 770,610 },{ 795,645 },{ 720,635 },{ 760,650 },{ 475,960 },
 93{ 95,260 },{ 875,920 },{ 700,500 },{ 555,815 },{ 830,485 },{ 1170,65 },
 94{ 830,610 },{ 605,625 },{ 595,360 },{ 1340,725 },{ 1740,245 } };
 95
 96int main()
 97{
 98    srand(1);
 99    int max_iterations = 600;
100    int max_no_improve = 50;
101    //初始化指针数组 
102    Delta = new int*[CITY_SIZE];
103    for (int i = 0; i < CITY_SIZE; i ++)
104        Delta[i] = new int[CITY_SIZE];
105
106    SOLUTION best_solution;
107
108    iterated_local_search(best_solution, berlin52, max_iterations, max_no_improve);
109
110    cout << endl<<endl<<"搜索完成！ 最优路线总长度 = " << best_solution.cost << endl;
111    cout << "最优访问城市序列如下：" << endl;
112    for (int i = 0; i < CITY_SIZE;i++)
113    {
114        cout << setw(4) << setiosflags(ios::left) << best_solution.permutation[i];
115    }
116
117    cout << endl << endl;
118
119    return 0;
120}
121
122
123
124//计算两个城市间距离
125int distance_2city(city c1, city c2)
126{
127    int distance = 0;
128    distance = sqrt((double)((c1.x - c2.x)*(c1.x - c2.x) + (c1.y - c2.y)*(c1.y - c2.y)));
129
130    return distance;
131}
132
133//根据产生的城市序列，计算旅游总距离
134//所谓城市序列，就是城市先后访问的顺序，比如可以先访问ABC，也可以先访问BAC等等
135//访问顺序不同，那么总路线长度也是不同的
136//p_perm 城市序列参数
137int cost_total(int * cities_permutation, CITIES * cities)
138{
139    int total_distance = 0;
140    int c1, c2;
141    //逛一圈，看看最后的总距离是多少
142    for (int i = 0; i < CITY_SIZE; i++)
143    {
144        c1 = cities_permutation[i];
145        if (i == CITY_SIZE - 1) //最后一个城市和第一个城市计算距离
146        {
147            c2 = cities_permutation[0];
148        }
149        else
150        {
151            c2 = cities_permutation[i + 1];
152        }
153        total_distance += distance_2city(cities[c1], cities[c2]);
154    }
155
156    return total_distance;
157}
158
159//获取随机城市排列
160void random_permutation(int * cities_permutation)
161{
162    int i, r, temp;
163    for (i = 0; i < CITY_SIZE; i++)
164    {
165        cities_permutation[i] = i; //初始化城市排列，初始按顺序排
166    }
167
168
169    for (i = 0; i < CITY_SIZE; i++)
170    {
171        //城市排列顺序随机打乱
172        r = rand() % (CITY_SIZE - i) + i;
173        temp = cities_permutation[i];
174        cities_permutation[i] = cities_permutation[r];
175        cities_permutation[r] = temp;
176    }
177}
178
179
180
181
182//颠倒数组中下标begin到end的元素位置
183void swap_element(int *p, int begin, int end)
184{
185    int temp;
186    while (begin < end)
187    {
188        temp = p[begin];
189        p[begin] = p[end];
190        p[end] = temp;
191        begin++;
192        end--;
193    }
194}
195
196
197//邻域动作 反转index_i <-> index_j 间的元素
198void two_opt_swap(int *cities_permutation, int *new_cities_permutation, int index_i, int index_j)
199{
200    for (int i = 0; i < CITY_SIZE; i++)
201    {
202        new_cities_permutation[i] = cities_permutation[i];
203    }
204
205    swap_element(new_cities_permutation, index_i, index_j);
206}
207
208
209
210int calc_delta(int i, int k,  int *tmp, CITIES * cities){
211    int delta = 0;
212    /*
213                以下计算说明：
214                对于每个方案，翻转以后没必要再次重新计算总距离
215                只需要在翻转的头尾做个小小处理
216
217                比如：
218                有城市序列   1-2-3-4-5 总距离 = d12 + d23 + d34 + d45 + d51 = A
219                翻转后的序列 1-4-3-2-5 总距离 = d14 + d43 + d32 + d25 + d51 = B
220                由于 dij 与 dji是一样的，所以B也可以表示成 B = A - d12 - d45 + d14 + d25
221                下面的优化就是基于这种原理
222    */
223    if (i == 0)
224    {
225        if (k == CITY_SIZE - 1)
226        {
227           delta = 0;
228        }
229        else
230        {
231            delta = 0
232                - distance_2city(cities[tmp[k]], cities[tmp[k + 1]])
233                + distance_2city(cities[tmp[i]], cities[tmp[k + 1]])
234                - distance_2city(cities[tmp[CITY_SIZE - 1]], cities[tmp[i]])
235                + distance_2city(cities[tmp[CITY_SIZE - 1]], cities[tmp[k]]);
236        }
237
238    }
239    else
240    {
241        if (k == CITY_SIZE - 1)
242        {
243            delta = 0
244                - distance_2city(cities[tmp[i - 1]], cities[tmp[i]])
245                + distance_2city(cities[tmp[i - 1]], cities[tmp[k]])
246                - distance_2city(cities[tmp[0]], cities[tmp[k]])
247                + distance_2city(cities[tmp[i]], cities[tmp[0]]);
248        }
249        else
250        {
251            delta = 0
252                - distance_2city(cities[tmp[i - 1]], cities[tmp[i]])
253                + distance_2city(cities[tmp[i - 1]], cities[tmp[k]])
254                - distance_2city(cities[tmp[k]], cities[tmp[k + 1]])
255                + distance_2city(cities[tmp[i]], cities[tmp[k + 1]]);
256        }
257    }
258
259    return delta;
260}
261
262
263/*
264    去重处理，对于Delta数组来说，对于城市序列1-2-3-4-5-6-7-8-9-10，如果对3-5应用了邻域操作2-opt ， 事实上对于
265    7-10之间的翻转是不需要重复计算的。 所以用Delta提前预处理一下。
266
267    当然由于这里的计算本身是O（1） 的，事实上并没有带来时间复杂度的减少（更新操作反而增加了复杂度） 
268    如果delta计算 是O（n）的，这种去重操作效果是明显的。 
269*/
270
271void Update(int i, int k,  int *tmp, CITIES * cities){
272    if (i && k != CITY_SIZE - 1){
273        i --; k ++;
274        for (int j = i; j <= k; j ++){
275            for (int l = j + 1; l < CITY_SIZE; l ++){
276                Delta[j][l] = calc_delta(j, l, tmp, cities);
277            }
278        }
279
280        for (int j = 0; j < k; j ++){
281            for (int l = i; l <= k; l ++){
282                if (j >= l) continue;
283                Delta[j][l] = calc_delta(j, l, tmp, cities);
284            }
285        }
286    }// 如果不是边界，更新(i-1, k + 1)之间的 
287    else{
288        for (i = 0; i < CITY_SIZE - 1; i++)
289            {
290              for (k = i + 1; k < CITY_SIZE; k++)
291            {
292                Delta[i][k] = calc_delta(i, k, tmp, cities);
293              }
294            }  
295    }// 边界要特殊更新 
296
297} 
298
299//本地局部搜索，边界条件 max_no_improve
300//best_solution最优解
301//current_solution当前解
302void local_search(SOLUTION & best_solution, CITIES * cities, int max_no_improve)
303{
304    int count = 0;
305    int i, k;
306
307    int inital_cost = best_solution.cost; //初始花费
308
309    int now_cost = 0;
310
311    SOLUTION *current_solution = new SOLUTION; //为了防止爆栈……直接new了，你懂的
312
313    for (i = 0; i < CITY_SIZE - 1; i++)
314    {
315       for (k = i + 1; k < CITY_SIZE; k++)
316        {
317            Delta[i][k] = calc_delta(i, k, best_solution.permutation, cities);
318        }
319    }
320
321    do
322    {
323        //枚举排列
324        for (i = 0; i < CITY_SIZE - 1; i++)
325        {
326            for (k = i + 1; k < CITY_SIZE; k++)
327            {
328                //邻域动作
329                two_opt_swap(best_solution.permutation, current_solution->permutation, i, k);
330                now_cost = inital_cost + Delta[i][k];
331                current_solution->cost = now_cost;
332                if (current_solution->cost < best_solution.cost)
333                {
334                    count = 0; //better cost found, so reset
335                    for (int j = 0; j < CITY_SIZE; j++)
336                    {
337                        best_solution.permutation[j] = current_solution->permutation[j];
338                    }
339                    best_solution.cost = current_solution->cost;
340                    inital_cost = best_solution.cost;
341                    Update(i, k, best_solution.permutation, cities);
342                }
343
344            }
345        }
346
347        count++;
348
349    } while (count <= max_no_improve);
350}
351//判断接受准则
352bool AcceptanceCriterion(int *cities_permutation, int *old_cities_permutation, CITIES * p_cities)
353{
354    int acceptance = 500; //接受条件,与当前最解相差不超过acceptance
355    int old_cost = cost_total(old_cities_permutation, p_cities);
356    int new_cost = cost_total(cities_permutation, p_cities);
357
358    if ((new_cost <= (old_cost + acceptance)) || (new_cost >= (old_cost - acceptance)))
359    {
360        return true;
361    }
362
363    return false;
364}
365
366//将城市序列分成4块，然后按块重新打乱顺序。
367//用于扰动函数
368void double_bridge_move(int *cities_permutation, int * new_cities_permutation)
369{
370    int temp_perm[CITY_SIZE];
371
372    int pos1 = 1 + rand() % (CITY_SIZE / 4);
373    int pos2 = pos1 + 1 + rand() % (CITY_SIZE / 4);
374    int pos3 = pos2 + 1 + rand() % (CITY_SIZE / 4);
375
376    int i;
377    vector<int> v;
378    //第一块
379    for (i = 0; i < pos1; i++)
380    {
381        v.push_back(cities_permutation[i]);
382    }
383
384    //第二块
385    for (i = pos3; i < CITY_SIZE; i++)
386    {
387        v.push_back(cities_permutation[i]);
388    }
389    //第三块
390    for (i = pos2; i < pos3; i++)
391    {
392        v.push_back(cities_permutation[i]);
393    }
394
395    //第四块
396    for (i = pos1; i < pos2; i++)
397    {
398        v.push_back(cities_permutation[i]);
399    }
400
401
402    for (i = 0; i < (int)v.size(); i++)
403    {
404        temp_perm[i] = v[i];
405    }
406    //if accept判断是否接受当前解
407    if (AcceptanceCriterion(cities_permutation, temp_perm, cities))
408    {
409        memcpy(new_cities_permutation, temp_perm, sizeof(temp_perm));//accept
410    }
411
412
413}
414
415//扰动
416void perturbation(CITIES * cities, SOLUTION &best_solution, SOLUTION &current_solution)
417{
418    double_bridge_move(best_solution.permutation, current_solution.permutation);
419    current_solution.cost = cost_total(current_solution.permutation, cities);
420}
421
422//迭代搜索
423//max_iterations用于迭代搜索次数
424//max_no_improve用于局部搜索边界条件
425void iterated_local_search(SOLUTION & best_solution, CITIES * cities, int max_iterations, int max_no_improve)
426{
427    SOLUTION *current_solution = new SOLUTION;
428
429    //获得初始随机解
430    random_permutation(best_solution.permutation);
431
432
433    best_solution.cost = cost_total(best_solution.permutation, cities);
434    local_search(best_solution, cities, max_no_improve); //初始搜索
435
436    for (int i = 0; i < max_iterations; i++)
437    {
438        perturbation(cities, best_solution, *current_solution); //扰动+判断是否接受新解
439        local_search(*current_solution, cities, max_no_improve);//继续局部搜索
440
441        //找到更优解
442        if (current_solution->cost < best_solution.cost)
443        {
444            for (int j = 0; j < CITY_SIZE; j++)
445            {
446                best_solution.permutation[j] = current_solution->permutation[j];
447            }
448            best_solution.cost = current_solution->cost;
449        }
450        cout << setw(13) << setiosflags(ios::left) <<"迭代搜索 " << i << " 次	" << "最优解 = " << best_solution.cost << " 当前解 = " << current_solution->cost << endl;
451    }
452
453}

运行结果