模擬退火(SA)算法求解Max-Minsum Dispersion Problem
把決策問題可以按照難易程度分為幾類:
(1)P問題:可以在多項式( polynomial )時間內(nèi)解決的問題,稱為P問題。
(2)NP問題:對于一類問題,我們可能沒有一個已知的快速的方法得到問題的答案,但給定一個解,我們可以在P時間內(nèi)檢查他正確與否的決策問題,成為NP( Non-deterministic polynomial)問題。
(3)NP-h(huán)ard問題:用一句話概括他們的特征就是“at least as hard as the hardest problems in NP Problem”。
【多項式級時間復雜度:O(1),O(log(n)),O(n^a)等。因為規(guī)模n出現(xiàn)在底數(shù)的位置!
1.1 Max-Minsum Dispersion Problem
Max-Minsum DP就是一個典型的NP-Hard問題,屬Equity-based dispersion problems,試圖從較大的集合中選擇一組元素時解決公平與效率的平衡,應用廣泛然而計算難度較大。
簡單來講,就是要從一個集合中選擇一個子集合,使得子集合中某個所選元素到其他所選元素之間距離的最小和最大化。
舉個例子,假如說你一天要完成5個任務,現(xiàn)在有10項任務供你選擇,煮飯30分鐘,做菜30分鐘,衣服機洗30分鐘,做作業(yè)20分鐘,燒開水10分鐘,吃飯15分鐘等,你是個聰明的無聊人,知道有些活兒可以一起干以節(jié)省時間,又想讓這5項任務占用你更多的時間以打發(fā)無聊,這便是個簡單的Max-Minsum DP。
再舉個更貼近實際的,對于網(wǎng)頁排名問題,即是標識網(wǎng)頁等級或重要性。最早的搜索引擎采用的是分類目錄的方法,即通過人工對網(wǎng)頁進行分類并整理出高質(zhì)量網(wǎng)站。隨著網(wǎng)頁數(shù)目的急劇增大,這種方法顯然無法實現(xiàn),Larry Page和Sergey Brin受學術(shù)界對學術(shù)論文重要性的評估方法(論文引用次數(shù))的啟發(fā),提出了PageRank算法,如果一個網(wǎng)頁被很多其它網(wǎng)頁鏈接到,說明這個網(wǎng)頁很重要,它的PageRank值也會相應較高,如果一個PageRank值很高的網(wǎng)頁鏈接到另外某個網(wǎng)頁,那么那個網(wǎng)頁的PageRank值也會相應地提高。所以說找重要論文、相關(guān)論文就不免涉及到Max-Minsum DP。
更多的應用如下圖:
1.2 Max-Minsum DP的數(shù)學描述
考慮一個含有n個元素的集合N={1,2,...,n},每個元素包含著r個屬性,我們可以將一個元素用向量表示。問題在于選擇N的一個子集M,|M|為固定正整數(shù)m(m<n),即從n個元素中選出m個元素,最大化所選元素到其他元素之間距離的最小和。
這個距離有多種算法,如歐幾里得距離,曼哈頓距離等。在這里我們使用最為常用的歐幾里得距離
問題可以表達為:
Part 2 模擬退火算法(SA)再回顧
在之前的推文【算法進階】用模擬退火(SA, Simulated Annealing)算法解決旅行商問題中,已經(jīng)對模擬退火算法有了詳細介紹并給出了偽代碼及實例。在這里,我們簡要復習,詳細參見以上推文。
2.1 SA算法介紹
模擬退火算法的基礎(chǔ)是metropolis算法。metropolis算法又稱為metropolis抽樣,其核心思想是:當能量增加的時候以一定概率接納,而非一味拒絕。
所以,當Y(i+1)>Y(i),則無條件接受;
???當Y(i+1)<Y(i),則以一定的概率接受,而非全然拒絕。
??以一定概率接受一個比當前解較差的解,從而在一定程度上避免陷入局部最優(yōu)。
??然而應當如何計算這個概率呢?根據(jù)熱力學的原理,在溫度為T時,出現(xiàn)能量差為dE的降溫的概率為P(dE),表示為:
其中k是一個常數(shù),且dE<0(溫度總是降低的)。
1)溫度越高,出現(xiàn)一次能量差為dE的降溫的概率就越大。
2)溫度越低,則出現(xiàn)降溫的概率就越小。
3)本問題將內(nèi)能E模擬為目標函數(shù)值 f。
具體算法介紹
通過鄰域動作產(chǎn)生新的集合M,產(chǎn)生新解及當前解,計算即smallestDelta。若當前解的smallestDelta小于最優(yōu)解的smallestDelta,則更新最優(yōu)解為當前解,否則以模擬退火的那個概率接受當前解,然后降溫。重復之前步驟,直到滿足退出條件。
現(xiàn)在拿一個小算例來操作一下:
現(xiàn)有一點集N={(0,1),(1,2),(3,4),(4,5),(6,6),(8,7)},我們要從中選出m個點構(gòu)成點集M,就取m=3吧,目標函數(shù)是我們挑選的這3個點中的某一點到其余2個點的距離之和的最小值,而問題在于找到使目標函數(shù)值最大的那3個點。
3.1 初始解生成
就小算例而言,我們就隨機選3個點,你不妨可以擲骰子,我擲的是5,2,1,那我們就。6,6), (1,2), (0,1)這三個點,不妨將這三個點重新標記為,
以為中心點,則;
以為中心點,則;
以為中心點,則;
不難看出,smallestDelta為Δ2,故初始解也是當前最優(yōu)解即為M={(6,6),(1,2),(0,1)},對應為。
??而就本問題而言,對于初始解,我們亦是隨機產(chǎn)生,距離矩陣利用洗牌算法隨機生成1-100的距離,隨機選擇m個元素構(gòu)成s1,未被選中的即為s0,為了識別M和NM,我們利用n維向量?= (?1 , ?2 , . . . , ?n ),其中,則,對應的集合M即為初始解,也作為最優(yōu)解。
3.2 鄰域動作
采用exchange算子:從被選擇的元素的集合中隨機選擇元素u,即u∈M,從不被選擇的元素的集合中隨機選擇元素v,即v∈NM,交換u, v。拿上文小算例N={(0,1),(1,2),(3,4),(4,5),(6,6),(8,7)}舉個例子,從、、中隨機選擇,即∈M,從、、中隨機選擇,即∈NM,交換、,此時得到新解,以三點分別為中心點,故、不變得到
以為中心點,則;
以為中心點,則;
以為中心點,則;
不難看出,smallestDelta為Δ2,故最優(yōu)解更新,變?yōu)镸={(4,5),(1,2),(0,1)},對應為。
3.3 去重優(yōu)化
對于本問題,給定鄰域解和對應向量(?1 , ?2 , . . . , ?n ),目標值可以在O(M)時間內(nèi)計算,此外,若是兩個元素u∈M,v∈NM交換,則向量?= (?1 , ?2 , . . . , ?n )可以在O(N)時間內(nèi)快速更新,具體可表示為下圖:
為了通俗易懂,接著拿上文小算例N={(0,1),(1,2),(3,4),(4,5),(6,6),(8,7)}舉例,比較3.1及3.2計算Δ過程不難看出,對于未改變的點,即以為中心點、以為中心點時,對應的Δ計算過程只改變了一半,這部分就是我們可以優(yōu)化的部分,因為另外一半我們就不用再重復計算了,
;
;
當數(shù)據(jù)越來越龐大之后,這部分優(yōu)化帶來的效益就會體現(xiàn)得更加明顯,時間復雜度大幅減少。
而對于改變的點,
,基本上就是重算。
代碼分享
算例為隨機生成,具體實現(xiàn)如下:
#include<iostream>
#include<cstdlib>
#include<cmath>
#include<string>
#include<ctime>
const int MAX = 0x7fffffff;
const int N = 1000; //最大的范圍
const int M = 500; //要選擇的集合大小
const int K = 100; //兩點間距離的最大值為K(距離默認為1-K)
const int max_count = 10; //當前溫度的最大迭代次數(shù)
const double T0 = 50000.0; //初始溫度
const double T_end = 1e-8; //退火結(jié)束溫度
const double q = 0.98; //退火系數(shù)
int* elements; //共計N個點
int** distance; //距離矩陣
clock_t start_total, end_total; //計時器,整個程序
clock_t start_delta, end_delta; //計時器,直接計算delta的步驟
struct Solution //解
{
int* s0; //未被選中的數(shù)
int* s1; //被選中的數(shù)
int* delta; //到其他s1中的數(shù)的距離和
int smallestDelta; //最大的delta,及目標函數(shù)值
int center; //核心數(shù)
}iniSolution, bestSolution,solution1;
//分配存儲空間
void init_solution(Solution* s)
{
s->s0 = new int[M];
s->s1 = new int[N - M];
s->delta = new int[M];
s->smallestDelta = 0;
}
//撤銷iniSolution, bestSolution,solution1所占存儲空間
void dispose(Solution* s)
{
delete[](s->s0);s->s0 = NULL;
delete[](s->s1);s->s1 = NULL;
delete[](s->delta);s->delta = NULL;
}
//深拷貝solution類型
void copy_solution(Solution* ini, Solution* obj)
{
for (int i = 0; i < M; i++)
obj->s1[i] = ini->s1[i];
for (int i = 0; i < M; i++)
obj->s0[i] = ini->s0[i];
for (int i = 0; i < M; i++)
obj->delta[i] = ini->delta[i];
obj->smallestDelta = ini->smallestDelta;
obj->center = ini->center;
}
//計算所有delta的值
void calculate_delta(Solution* s)
{
for (int i = 0; i < M; i++)
{
s->delta[i] = 0;
for (int j = 0; j < M; j++)
s->delta[i] += distance[s->s1[i]][s->s1[j]];
}
}
// 在所有delta中找出smallest delta以及對應的中心數(shù)
void calculate_sum(Solution* s)
{
s->smallestDelta = s->delta[0];
s->center = s->s1[0];
for (int i = 0; i < M; i++)
{
if (s->delta[i] < s->smallestDelta)
{
s->smallestDelta = s->delta[i];
s->center = s->s1[i];
}
}
}
//更新的方法算出delta的值(將s0[v]與s1[u]交換)
void update_delta(Solution* s, int u, int v, int deltav)
{
for (int i = 0; i < M; i++)
{
if (i == u) //其自身delta的改變
s->delta[u] = deltav - distance[s->s0[v]][s->s1[u]];
else //其他delta需將與s1[u]的距離轉(zhuǎn)換為與s0[v]的距離
s->delta[i] = s->delta[i] - distance[s->s1[u]][s->s1[i]] + distance[s->s0[v]][s->s1[i]];
}
}
void init()
{
//為距離矩陣隨機生成1-100的距離
for (int i = 0;i < N;i++)
for (int j = 0;j < i;j++) //因為距離矩陣是對稱的
distance[i][j] = distance[j][i] = rand() % K + 1; //距離為1-K
for (int i = 0;i < N;i++)
distance[i][i] = 0;
//隨機生成初始解
for (int i = 0; i < N; i++)
elements[i] = i;
//洗牌算法打亂
for (int i = 0; i < N; i++)
{
int index = rand() % (N - i) + i;
if (index 。 i)
{
int temp = elements[i];
elements[i] = elements[index];
elements[index] = temp;
}
}
//初始化,分配數(shù)組空間
init_solution(&iniSolution);
//前M個為s1,后面為s0
for (int i = 0;i < M;i++)
iniSolution.s1[i] = elements[i];
for (int i = M, j = 0;i < N;i++, j++)
iniSolution.s0[j] = elements[i];
//計算delta
start_delta = clock();
calculate_delta(&iniSolution);
end_delta = clock();
//計算smallest_delta
calculate_sum(&iniSolution);
//bestSolution拷貝iniSolution
init_solution(&bestSolution);
copy_solution(&iniSolution, &bestSolution);
dispose(&iniSolution);
//for (int i = 0;i < M;i++) std::cout << bestSolution.s1[i] << std::endl;
}
void SA_search() //模擬退火算法Simulated Annealing
{
srand((unsigned)time(NULL)); //初始化隨機數(shù)種子
double T = T0; //初始溫度
int count_total = 0; //記錄降溫次數(shù)
while (T > T_end) // 當溫度低于結(jié)束溫度時,退火結(jié)束
{
for (int count = 0;count <= max_count;count++) //count記錄當前溫度迭代次數(shù)
{
int deltav = 0; //計算deltav
//產(chǎn)生新解solution1
init_solution(&solution1);
copy_solution(&bestSolution, &solution1);
double r1 = ((double)rand()) / (RAND_MAX + 1.0);
double r2 = ((double)rand()) / (RAND_MAX + 1.0);
int v = (int)((N - M) * r1); //s0中交換點的位置
int u = (int)(M * r2); //s1中交換點的位置
for (int u = 0;u < M;u++) //對選中的數(shù)(s1)進行循環(huán)
deltav += distance[bestSolution.s0[v]][bestSolution.s1[u]];
update_delta(&solution1, u, v, deltav);
int temp = solution1.s0[v];
solution1.s0[v] = solution1.s1[u];
solution1.s1[u] = temp;
calculate_sum(&solution1); //計算smallest_delta
double f1, f2, df;
f1 = bestSolution.smallestDelta;
f2 = solution1.smallestDelta;
df = f2 - f1;
double r = ((double)rand()) / (RAND_MAX); //0-1之間的隨機數(shù),用來決定是否接受新解
if (df >= 0)
copy_solution(&solution1, &bestSolution);
else if (r < exp(df / T)) //若隨機數(shù)小于p,接受新解
copy_solution(&solution1, &bestSolution);
dispose(&solution1);
count++;
}
T *= q; //降溫
count_total++;
std::cout << "第" << count_total << "次降溫, 當前溫度:" << T << ",當前最優(yōu)解:" << bestSolution.smallestDelta << std::endl;
}
}
void print_info()
{
std::cout << "max-min sum answer:" << bestSolution.smallestDelta << std::endl;
std::cout << "one delta run time:" << (double)(end_delta - start_delta) / CLOCKS_PER_SEC << std::endl;
std::cout << "total run time:" << (double)(end_total - start_total) / CLOCKS_PER_SEC << std::endl;
std::cout << "模擬退火算法,初始溫度T0=" << T0 << ",降溫系數(shù)q=" << q << ",每個溫度迭代" << max_count << "次" << std::endl;
}
int main()
{
//初始化數(shù)組,分配空間
elements = new int[N];
distance = new int* [N];
for (int i = 0;i < N;i++)
distance[i] = new int[N];
init();
start_total = clock();
SA_search(); //模擬退火算法進行搜索
end_total = clock();
//打印結(jié)果
print_info();
dispose(&bestSolution);
delete[]elements;
for (int i = 0;i < N;i++)
delete[]distance[i];
delete[]distance;
system("pause");
return 0;
}
結(jié)果如圖:
欲下載本文相關(guān)代碼,請移步留言區(qū)
參考文獻:
xiangjing Lai,Dong Yue,Jin-Kao Hao,F(xiàn)red Glover "Solution-based tabu search for the maximum min-sum dispersion problem." Information Sciences 441 (2018) 79-94.
-The End-
文案/代碼/排版:朱正雄
指導學長:周航
指導老師:秦虎 華中科技大學管理學院

請輸入評論內(nèi)容...
請輸入評論/評論長度6~500個字
最新活動更多
推薦專題
- 1 UALink規(guī)范發(fā)布:挑戰(zhàn)英偉達AI統(tǒng)治的開始
- 2 北電數(shù)智主辦酒仙橋論壇,探索AI產(chǎn)業(yè)發(fā)展新路徑
- 3 降薪、加班、裁員三重暴擊,“AI四小龍”已折戟兩家
- 4 “AI寒武紀”爆發(fā)至今,五類新物種登上歷史舞臺
- 5 國產(chǎn)智駕迎戰(zhàn)特斯拉FSD,AI含量差幾何?
- 6 光計算迎來商業(yè)化突破,但落地仍需時間
- 7 東陽光:2024年扭虧、一季度凈利大增,液冷疊加具身智能打開成長空間
- 8 地平線自動駕駛方案解讀
- 9 封殺AI“照騙”,“淘寶們”終于不忍了?
- 10 優(yōu)必選:營收大增主靠小件,虧損繼續(xù)又逢關(guān)稅,能否乘機器人東風翻身?