模擬退火(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”。
【多項式級時間復(fù)雜度: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,試圖從較大的集合中選擇一組元素時解決公平與效率的平衡,應(yīng)用廣泛然而計算難度較大。
簡單來講,就是要從一個集合中選擇一個子集合,使得子集合中某個所選元素到其他所選元素之間距離的最小和最大化。
舉個例子,假如說你一天要完成5個任務(wù),現(xiàn)在有10項任務(wù)供你選擇,煮飯30分鐘,做菜30分鐘,衣服機洗30分鐘,做作業(yè)20分鐘,燒開水10分鐘,吃飯15分鐘等,你是個聰明的無聊人,知道有些活兒可以一起干以節(jié)省時間,又想讓這5項任務(wù)占用你更多的時間以打發(fā)無聊,這便是個簡單的Max-Minsum DP。
再舉個更貼近實際的,對于網(wǎng)頁排名問題,即是標(biāo)識網(wǎng)頁等級或重要性。最早的搜索引擎采用的是分類目錄的方法,即通過人工對網(wǎng)頁進(jìn)行分類并整理出高質(zhì)量網(wǎng)站。隨著網(wǎng)頁數(shù)目的急劇增大,這種方法顯然無法實現(xiàn),Larry Page和Sergey Brin受學(xué)術(shù)界對學(xué)術(shù)論文重要性的評估方法(論文引用次數(shù))的啟發(fā),提出了PageRank算法,如果一個網(wǎng)頁被很多其它網(wǎng)頁鏈接到,說明這個網(wǎng)頁很重要,它的PageRank值也會相應(yīng)較高,如果一個PageRank值很高的網(wǎng)頁鏈接到另外某個網(wǎng)頁,那么那個網(wǎng)頁的PageRank值也會相應(yīng)地提高。所以說找重要論文、相關(guān)論文就不免涉及到Max-Minsum DP。
更多的應(yīng)用如下圖:
1.2 Max-Minsum DP的數(shù)學(xué)描述
考慮一個含有n個元素的集合N={1,2,...,n},每個元素包含著r個屬性,我們可以將一個元素用向量表示。問題在于選擇N的一個子集M,|M|為固定正整數(shù)m(m<n),即從n個元素中選出m個元素,最大化所選元素到其他元素之間距離的最小和。
這個距離有多種算法,如歐幾里得距離,曼哈頓距離等。在這里我們使用最為常用的歐幾里得距離
問題可以表達(dá)為:
Part 2 模擬退火算法(SA)再回顧
在之前的推文【算法進(jìn)階】用模擬退火(SA, Simulated Annealing)算法解決旅行商問題中,已經(jīng)對模擬退火算法有了詳細(xì)介紹并給出了偽代碼及實例。在這里,我們簡要復(fù)習(xí),詳細(xì)參見以上推文。
2.1 SA算法介紹
模擬退火算法的基礎(chǔ)是metropolis算法。metropolis算法又稱為metropolis抽樣,其核心思想是:當(dāng)能量增加的時候以一定概率接納,而非一味拒絕。
所以,當(dāng)Y(i+1)>Y(i),則無條件接受;
???當(dāng)Y(i+1)<Y(i),則以一定的概率接受,而非全然拒絕。
??以一定概率接受一個比當(dāng)前解較差的解,從而在一定程度上避免陷入局部最優(yōu)。
??然而應(yīng)當(dāng)如何計算這個概率呢?根據(jù)熱力學(xué)的原理,在溫度為T時,出現(xiàn)能量差為dE的降溫的概率為P(dE),表示為:
其中k是一個常數(shù),且dE<0(溫度總是降低的)。
1)溫度越高,出現(xiàn)一次能量差為dE的降溫的概率就越大。
2)溫度越低,則出現(xiàn)降溫的概率就越小。
3)本問題將內(nèi)能E模擬為目標(biāo)函數(shù)值 f。
具體算法介紹
通過鄰域動作產(chǎn)生新的集合M,產(chǎn)生新解及當(dāng)前解,計算即smallestDelta。若當(dāng)前解的smallestDelta小于最優(yōu)解的smallestDelta,則更新最優(yōu)解為當(dāng)前解,否則以模擬退火的那個概率接受當(dāng)前解,然后降溫。重復(fù)之前步驟,直到滿足退出條件。
現(xiàn)在拿一個小算例來操作一下:
現(xiàn)有一點集N={(0,1),(1,2),(3,4),(4,5),(6,6),(8,7)},我們要從中選出m個點構(gòu)成點集M,就取m=3吧,目標(biāo)函數(shù)是我們挑選的這3個點中的某一點到其余2個點的距離之和的最小值,而問題在于找到使目標(biāo)函數(shù)值最大的那3個點。
3.1 初始解生成
就小算例而言,我們就隨機選3個點,你不妨可以擲骰子,我擲的是5,2,1,那我們就。6,6), (1,2), (0,1)這三個點,不妨將這三個點重新標(biāo)記為,
以為中心點,則;
以為中心點,則;
以為中心點,則;
不難看出,smallestDelta為Δ2,故初始解也是當(dāng)前最優(yōu)解即為M={(6,6),(1,2),(0,1)},對應(yīng)為。
??而就本問題而言,對于初始解,我們亦是隨機產(chǎn)生,距離矩陣?yán)孟磁扑惴S機生成1-100的距離,隨機選擇m個元素構(gòu)成s1,未被選中的即為s0,為了識別M和NM,我們利用n維向量?= (?1 , ?2 , . . . , ?n ),其中,則,對應(yīng)的集合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)},對應(yīng)為。
3.3 去重優(yōu)化
對于本問題,給定鄰域解和對應(yīng)向量(?1 , ?2 , . . . , ?n ),目標(biāo)值可以在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īng)的Δ計算過程只改變了一半,這部分就是我們可以優(yōu)化的部分,因為另外一半我們就不用再重復(fù)計算了,
;
;
當(dāng)數(shù)據(jù)越來越龐大之后,這部分優(yōu)化帶來的效益就會體現(xiàn)得更加明顯,時間復(fù)雜度大幅減少。
而對于改變的點,
,基本上就是重算。
代碼分享
算例為隨機生成,具體實現(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(距離默認(rèn)為1-K)
const int max_count = 10; //當(dāng)前溫度的最大迭代次數(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,及目標(biāo)函數(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以及對應(yīng)的中心數(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) // 當(dāng)溫度低于結(jié)束溫度時,退火結(jié)束
{
for (int count = 0;count <= max_count;count++) //count記錄當(dāng)前溫度迭代次數(shù)
{
int deltav = 0; //計算deltav
//產(chǎn)生新解solution1
init_solution(&solution1);
copy_solution(&bestSolution, &solution1);
double r1 = ((double)rand()) / (RAND_M(jìn)AX + 1.0);
double r2 = ((double)rand()) / (RAND_M(jìn)AX + 1.0);
int v = (int)((N - M) * r1); //s0中交換點的位置
int u = (int)(M * r2); //s1中交換點的位置
for (int u = 0;u < M;u++) //對選中的數(shù)(s1)進(jìn)行循環(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_M(jìn)AX); //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 << "次降溫, 當(dāng)前溫度:" << T << ",當(dāng)前最優(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(); //模擬退火算法進(jìn)行搜索
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ū)
參考文獻(xiàn):
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-
文案/代碼/排版:朱正雄
指導(dǎo)學(xué)長:周航
指導(dǎo)老師:秦虎 華中科技大學(xué)管理學(xué)院
請輸入評論內(nèi)容...
請輸入評論/評論長度6~500個字
最新活動更多
-
10月31日立即下載>> 【限時免費下載】TE暖通空調(diào)系統(tǒng)高效可靠的組件解決方案
-
即日-11.13立即報名>>> 【在線會議】多物理場仿真助跑新能源汽車
-
11月28日立即報名>>> 2024工程師系列—工業(yè)電子技術(shù)在線會議
-
12月19日立即報名>> 【線下會議】OFweek 2024(第九屆)物聯(lián)網(wǎng)產(chǎn)業(yè)大會
-
即日-12.26火熱報名中>> OFweek2024中國智造CIO在線峰會
-
即日-2025.8.1立即下載>> 《2024智能制造產(chǎn)業(yè)高端化、智能化、綠色化發(fā)展藍(lán)皮書》
推薦專題
- 高級軟件工程師 廣東省/深圳市
- 自動化高級工程師 廣東省/深圳市
- 光器件研發(fā)工程師 福建省/福州市
- 銷售總監(jiān)(光器件) 北京市/海淀區(qū)
- 激光器高級銷售經(jīng)理 上海市/虹口區(qū)
- 光器件物理工程師 北京市/海淀區(qū)
- 激光研發(fā)工程師 北京市/昌平區(qū)
- 技術(shù)專家 廣東省/江門市
- 封裝工程師 北京市/海淀區(qū)
- 結(jié)構(gòu)工程師 廣東省/深圳市