學(xué)建模優(yōu)化難題破解：隨機(jī)搜索算法原理與應(yīng)急選址實(shí)戰(zhàn))
1. 項(xiàng)目緣起當(dāng)數(shù)學(xué)建模遇上“大海撈針”做數(shù)學(xué)建模尤其是國賽、美賽這種高強(qiáng)度的比賽最讓人頭疼的環(huán)節(jié)之一往往不是模型建立而是模型求解。你花了大量心血構(gòu)建了一個(gè)邏輯嚴(yán)密、變量眾多的非線性規(guī)劃模型感覺已經(jīng)成功了一大半。但當(dāng)你打開MATLAB或者Python準(zhǔn)備用那些教科書上的經(jīng)典算法比如梯度下降、牛頓法去求解時(shí)卻發(fā)現(xiàn)要么是初始點(diǎn)敏感動(dòng)不動(dòng)就掉進(jìn)局部最優(yōu)的“坑”里出不來要么是目標(biāo)函數(shù)或約束條件太復(fù)雜導(dǎo)數(shù)都求不出來算法直接“罷工”。這時(shí)候你需要的可能不是更精密的微調(diào)而是一種思路上的轉(zhuǎn)變——一種不那么“聰明”但足夠“魯棒”和“通用”的方法。這就是隨機(jī)搜索算法Random Search Algorithm登場(chǎng)的時(shí)刻。很多人一聽到“隨機(jī)搜索”第一反應(yīng)可能是“這不就是瞎蒙嗎”。確實(shí)它的核心思想非常樸素在變量的定義域可行域內(nèi)隨機(jī)地生成大量的候選解然后從中找出目標(biāo)函數(shù)值最好的那個(gè)。聽起來毫無技術(shù)含量甚至有些“暴力”。但正是這種“暴力”在應(yīng)對(duì)多變量、非線性、非凸、甚至“黑箱”函數(shù)的最優(yōu)化問題時(shí)展現(xiàn)出了驚人的實(shí)用性。它不依賴于函數(shù)的梯度信息對(duì)函數(shù)的形態(tài)幾乎沒有要求只要你能計(jì)算出任何一個(gè)點(diǎn)的函數(shù)值它就能工作。在數(shù)學(xué)建模競(jìng)賽中當(dāng)你面對(duì)一個(gè)結(jié)構(gòu)復(fù)雜、機(jī)理不甚清晰的實(shí)際問題需要快速得到一個(gè)“還不錯(cuò)”的、可用的解時(shí)隨機(jī)搜索往往能成為你的救命稻草。我最初深入使用隨機(jī)搜索是在一次模擬復(fù)雜供應(yīng)鏈網(wǎng)絡(luò)優(yōu)化的項(xiàng)目中。模型涉及幾十個(gè)決策變量如倉庫選址、運(yùn)輸路徑、庫存水平約束條件相互耦合目標(biāo)函數(shù)總成本是高度非線性的。嘗試了多種基于梯度的優(yōu)化器效果都不理想。最后抱著試試看的心態(tài)寫了一個(gè)簡(jiǎn)單的隨機(jī)搜索程序雖然每次運(yùn)行結(jié)果都有波動(dòng)但總能穩(wěn)定地找到一個(gè)比之前任何方法都好的解。從那以后隨機(jī)搜索就成了我解決棘手優(yōu)化問題的“標(biāo)準(zhǔn)備選方案”之一。本文將結(jié)合一個(gè)具體的建模案例拆解隨機(jī)搜索從原理到實(shí)現(xiàn)的每一個(gè)細(xì)節(jié)并分享我在實(shí)戰(zhàn)中積累的調(diào)參經(jīng)驗(yàn)和避坑指南。2. 隨機(jī)搜索的核心原理為什么“笨辦法”有時(shí)更有效要理解隨機(jī)搜索的價(jià)值我們得先看看它的對(duì)手們?yōu)槭裁磿?huì)在某些場(chǎng)景下失效。傳統(tǒng)的最優(yōu)化算法如梯度下降法其哲學(xué)是“局部尋優(yōu)”。它假設(shè)函數(shù)是光滑的、可微的并且當(dāng)前點(diǎn)的梯度方向能指引我們?nèi)ネ粋€(gè)更優(yōu)的鄰域。這就像在一個(gè)連綿起伏的山丘中你蒙著眼睛但能用腳感受地面的傾斜梯度然后朝著感覺是下坡的方向走。這個(gè)方法在“山丘”形狀良好凸函數(shù)時(shí)非常高效能快速找到谷底全局最優(yōu)。然而現(xiàn)實(shí)世界的優(yōu)化問題尤其是數(shù)學(xué)建模中抽象出來的問題其“地形”往往更像一片布滿深坑、斷崖和緩坡的復(fù)雜地貌非凸函數(shù)。梯度下降法從某個(gè)起點(diǎn)出發(fā)很容易掉進(jìn)最近的一個(gè)坑里局部最優(yōu)并且由于在坑底梯度為零或無法計(jì)算它就認(rèn)為已經(jīng)到達(dá)終點(diǎn)再也出不來了。這就是“局部最優(yōu)陷阱”。隨機(jī)搜索采取了完全不同的策略。它放棄了“局部感知”轉(zhuǎn)而進(jìn)行“全局采樣”。想象一下你不是那個(gè)蒙眼下坡的人而是指揮成千上萬個(gè)無人機(jī)在整個(gè)區(qū)域上空隨機(jī)空投傳感器。每個(gè)傳感器落地后立刻報(bào)告所在地的海拔高度函數(shù)值。最后你只需從所有報(bào)告中找出海拔最低的那個(gè)點(diǎn)。這個(gè)方法不關(guān)心地形是否連續(xù)、是否可導(dǎo)它只關(guān)心1. 我能否定義出整個(gè)搜索區(qū)域變量上下界2. 我能否對(duì)區(qū)域內(nèi)的任何一個(gè)點(diǎn)進(jìn)行評(píng)估。2.1 算法流程的標(biāo)準(zhǔn)化描述一個(gè)最基礎(chǔ)的多變量隨機(jī)搜索算法可以歸納為以下幾步定義問題明確你的目標(biāo)是最小化還是最大化一個(gè)目標(biāo)函數(shù) f(x)。確定決策變量向量 x [x1, x2, ..., xn] 的搜索空間通常由每個(gè)變量的下限 Lower Bound (LB) 和上限 Upper Bound (UB) 定義即 LB_i ≤ x_i ≤ UB_i。初始化設(shè)定算法要進(jìn)行的總迭代次數(shù)max_iter或總函數(shù)評(píng)估次數(shù)max_evaluations。初始化當(dāng)前最優(yōu)解best_x和對(duì)應(yīng)的最優(yōu)值best_f為一個(gè)很差的初始值例如對(duì)于最小化問題best_f設(shè)為無窮大。迭代搜索 a.生成候選點(diǎn)在當(dāng)前迭代k在變量的定義域內(nèi)均勻隨機(jī)地生成一個(gè)候選解candidate_x。即對(duì)于每個(gè)變量 icandidate_x[i] LB_i random() * (UB_i - LB_i)其中random()生成一個(gè)[0,1)之間的隨機(jī)數(shù)。 b.可行性檢查可選如果問題存在除變量上下界外的其他復(fù)雜約束如線性/非線性不等式約束需要檢查candidate_x是否滿足所有約束。如果不滿足則丟棄該點(diǎn)回到步驟a生成下一個(gè)隨機(jī)點(diǎn)或采取懲罰函數(shù)等策略。對(duì)于入門我們先處理無復(fù)雜約束的問題。 c.性能評(píng)估計(jì)算該候選點(diǎn)的目標(biāo)函數(shù)值current_f f(candidate_x)。 d.更新最優(yōu)解將current_f與歷史最優(yōu)值best_f比較。如果更優(yōu)對(duì)于最小化問題是更小則更新best_f current_f同時(shí)更新best_x candidate_x。終止與輸出重復(fù)步驟3直到達(dá)到預(yù)設(shè)的迭代次數(shù)或評(píng)估次數(shù)。最終輸出找到的best_x和best_f。這個(gè)流程簡(jiǎn)單到幾乎可以用任何編程語言在十幾行代碼內(nèi)實(shí)現(xiàn)。但正是這種簡(jiǎn)單性帶來了幾個(gè)關(guān)鍵優(yōu)勢(shì)全局性由于采樣是全局均勻的算法有概率采樣到整個(gè)可行域的任何角落因此理論上只要采樣點(diǎn)足夠多就有機(jī)會(huì)逼近甚至命中全局最優(yōu)解。魯棒性對(duì)目標(biāo)函數(shù) f(x) 的性質(zhì)幾乎無要求。f(x) 可以是離散的、不可微的、不連續(xù)的甚至是一個(gè)需要調(diào)用外部仿真軟件才能得到結(jié)果的“黑箱”函數(shù)。易并行性每一次隨機(jī)采樣和評(píng)估都是完全獨(dú)立的可以非常容易地分配到多個(gè)CPU核心或計(jì)算節(jié)點(diǎn)上并行執(zhí)行從而大幅縮短搜索時(shí)間。注意隨機(jī)搜索找到的“最優(yōu)解”是一個(gè)隨機(jī)變量。每次運(yùn)行的結(jié)果都可能不同。我們通常通過多次獨(dú)立運(yùn)行取最好的結(jié)果或統(tǒng)計(jì)結(jié)果的分布來評(píng)估算法的性能。3. 建模案例實(shí)戰(zhàn)應(yīng)急物資儲(chǔ)備庫的選址優(yōu)化為了讓大家更直觀地理解隨機(jī)搜索如何應(yīng)用于實(shí)際建模我們構(gòu)造一個(gè)簡(jiǎn)化但經(jīng)典的運(yùn)籌學(xué)問題多需求點(diǎn)應(yīng)急物資儲(chǔ)備庫選址優(yōu)化。3.1 問題描述與模型建立假設(shè)某地區(qū)有M個(gè)潛在的應(yīng)急物資儲(chǔ)備庫選址點(diǎn)需要服務(wù)N個(gè)已知的居民點(diǎn)需求點(diǎn)。每個(gè)居民點(diǎn)j有一個(gè)固定的物資年需求量d_j。每個(gè)候選儲(chǔ)備庫i有一個(gè)最大的建設(shè)容量C_i和一個(gè)固定的年運(yùn)營建設(shè)成本F_i只要選中建設(shè)就會(huì)產(chǎn)生此成本。從儲(chǔ)備庫i運(yùn)輸單位物資到居民點(diǎn)j的運(yùn)費(fèi)為c_{ij}。我們的決策是選址決策決定在哪些候選點(diǎn)建設(shè)儲(chǔ)備庫。用一個(gè)0-1變量y_i表示y_i 1表示在點(diǎn)i建設(shè)y_i 0表示不建設(shè)。分配決策決定每個(gè)建設(shè)的儲(chǔ)備庫i向每個(gè)居民點(diǎn)j運(yùn)輸多少物資。用一個(gè)連續(xù)變量x_{ij}表示。目標(biāo)是最小化總成本包括所有被選中的儲(chǔ)備庫的固定建設(shè)成本以及所有的運(yùn)輸成本。數(shù)學(xué)模型如下目標(biāo)函數(shù)最小化總成本Minimize Z Σ_{i1}^{M} (F_i * y_i) Σ_{i1}^{M} Σ_{j1}^{N} (c_{ij} * x_{ij})第一項(xiàng)是固定成本第二項(xiàng)是運(yùn)輸成本。約束條件需求滿足約束每個(gè)居民點(diǎn)j的需求必須被完全滿足。Σ_{i1}^{M} x_{ij} d_j, for all j 1, 2, ..., N容量約束每個(gè)儲(chǔ)備庫i發(fā)出的物資總量不能超過其建設(shè)容量。Σ_{j1}^{N} x_{ij} ≤ C_i * y_i, for all i 1, 2, ..., M 注意如果y_i 0則右側(cè)為0意味著從該點(diǎn)運(yùn)出的物資x_{ij}必須全部為0如果y_i 1則不能超過C_i。邏輯約束只有被選中的儲(chǔ)備庫才能分配物資。x_{ij} ≥ 0, for all i, jy_i ∈ {0, 1}, for all i這是一個(gè)典型的**混合整數(shù)線性規(guī)劃MILP**問題。對(duì)于小規(guī)模問題可以使用專業(yè)的優(yōu)化求解器如Gurobi, CPLEX精確求解。但在數(shù)學(xué)建模競(jìng)賽中問題規(guī)模可能較大或者環(huán)境限制無法使用商業(yè)求解器。此時(shí)隨機(jī)搜索就可以作為一個(gè)有效的近似求解方案。3.2 隨機(jī)搜索求解策略設(shè)計(jì)直接對(duì)所有的y_i和x_{ij}進(jìn)行隨機(jī)搜索效率極低因?yàn)樽兞刻嗲掖嬖趶?fù)雜的約束關(guān)系。我們需要利用問題的結(jié)構(gòu)設(shè)計(jì)更聰明的搜索策略。一個(gè)有效的策略是將問題分解外層搜索隨機(jī)搜索負(fù)責(zé)搜索y_i的0-1組合即決定“建哪些庫”。這是一個(gè)組合優(yōu)化問題。內(nèi)層求解線性規(guī)劃/運(yùn)輸問題對(duì)于外層給定的一個(gè)具體的選址方案y即確定了哪些y_i1剩下的問題是一個(gè)標(biāo)準(zhǔn)的運(yùn)輸問題在已選定的倉庫集合內(nèi)如何分配x_{ij}以滿足所有需求且不超出選定倉庫的容量并使運(yùn)輸成本最低。這個(gè)問題是線性規(guī)劃有高效算法如單純形法可以快速精確求解甚至對(duì)于簡(jiǎn)單情況可以直接用線性規(guī)劃求解器或自己編寫算法如表上作業(yè)法。算法流程調(diào)整如下編碼一個(gè)解表示為一個(gè)長(zhǎng)度為M的0-1向量代表y_1, y_2, ..., y_M。生成候選選址方案隨機(jī)生成一個(gè)0-1向量。為了增加可行性可以加入啟發(fā)式規(guī)則比如至少生成一個(gè)y_i1?？尚行赃^濾檢查隨機(jī)生成的選址方案是否容量可行。即所有被選中倉庫的總?cè)萘喀?(C_i * y_i)是否大于等于總需求Σ d_j。如果不滿足這個(gè)方案不可能滿足所有需求直接丟棄重新生成。求解子問題對(duì)于容量可行的選址方案將其y_i值固定求解內(nèi)層的運(yùn)輸問題得到最優(yōu)的x_{ij}分配和對(duì)應(yīng)的最小運(yùn)輸成本Transport_Cost(y)。計(jì)算總成本總成本Z Σ (F_i * y_i) Transport_Cost(y)。更新全局最優(yōu)比較并更新。這樣隨機(jī)搜索的核心就變成了在指數(shù)級(jí)數(shù)量的選址組合中尋找能使“固定成本對(duì)應(yīng)最優(yōu)運(yùn)輸成本”最小的那個(gè)組合。內(nèi)層運(yùn)輸問題的精確求解保證了對(duì)于任何一個(gè)選址方案我們都能得到其可能達(dá)到的最低運(yùn)輸成本從而公平地比較不同選址方案的優(yōu)劣。3.3 Python代碼實(shí)現(xiàn)與解析下面我們用Python來實(shí)現(xiàn)這個(gè)策略。我們將使用numpy生成隨機(jī)數(shù)并使用pulp一個(gè)免費(fèi)的線性規(guī)劃庫來求解內(nèi)層的運(yùn)輸問題。pulp不是標(biāo)準(zhǔn)庫需要安裝pip install pulp。import numpy as np import pulp as lp import time def solve_facility_location_with_random_search(M, N, F, C, d, c, max_iter1000, seed42): 使用隨機(jī)搜索算法求解設(shè)施選址問題。 參數(shù): M: 潛在設(shè)施數(shù)量 N: 需求點(diǎn)數(shù)量 F: list of length M, 每個(gè)設(shè)施的固定成本 C: list of length M, 每個(gè)設(shè)施的容量 d: list of length N, 每個(gè)需求點(diǎn)的需求量 c: 2D list of shape (M, N), 運(yùn)輸成本矩陣c[i][j] 從設(shè)施i到需求點(diǎn)j的成本 max_iter: 最大隨機(jī)搜索迭代次數(shù) seed: 隨機(jī)種子用于復(fù)現(xiàn)結(jié)果 返回: best_y: 最優(yōu)的選址方案 (0-1 list) best_x: 最優(yōu)的運(yùn)輸方案 (2D list) best_cost: 最優(yōu)總成本 history: 每次迭代找到的最佳成本記錄 np.random.seed(seed) total_demand sum(d) # 初始化最優(yōu)解 best_cost float(inf) best_y None best_x None history [] for iteration in range(max_iter): # 1. 隨機(jī)生成一個(gè)選址方案 y (0-1向量) y_candidate np.random.randint(0, 2, sizeM) # 2. 可行性檢查確保至少選一個(gè)設(shè)施且總?cè)萘? 總需求 if np.sum(y_candidate) 0: continue # 沒選任何設(shè)施不可行 total_capacity np.dot(C, y_candidate) if total_capacity total_demand: continue # 容量不足不可行 # 3. 構(gòu)建并求解內(nèi)層運(yùn)輸問題 # 創(chuàng)建問題實(shí)例最小化 transport_prob lp.LpProblem(Transportation_Subproblem, lp.LpMinimize) # 創(chuàng)建決策變量 x[i][j] 0 x_vars lp.LpVariable.dicts(x, ((i, j) for i in range(M) for j in range(N) if y_candidate[i] 1), lowBound0, catContinuous) # 如果設(shè)施i未被選中(y_candidate[i]0)則對(duì)應(yīng)的x變量不會(huì)被創(chuàng)建天然為0。 # 目標(biāo)函數(shù)最小化運(yùn)輸成本 transport_prob lp.lpSum(c[i][j] * x_vars[(i, j)] for (i, j) in x_vars.keys()) # 約束條件1: 每個(gè)需求點(diǎn)的需求必須滿足 for j in range(N): # 對(duì)所有選中的設(shè)施i求和其運(yùn)往j的物資量 prob lp.lpSum(x_vars.get((i, j), 0) for i in range(M) if y_candidate[i] 1) d[j] # 約束條件2: 每個(gè)選中設(shè)施的運(yùn)出量不超過其容量 for i in range(M): if y_candidate[i] 1: prob lp.lpSum(x_vars.get((i, j), 0) for j in range(N)) C[i] # 求解運(yùn)輸子問題 # pulp默認(rèn)使用CBC求解器對(duì)于線性規(guī)劃足夠 transport_prob.solve(pulp.PULP_CBC_CMD(msgFalse)) # msgFalse關(guān)閉求解器輸出 # 檢查是否求解成功 if lp.LpStatus[transport_prob.status] ! Optimal: # 如果子問題無解理論上在容量可行下應(yīng)該不會(huì)發(fā)生但保留檢查 continue # 獲取子問題最優(yōu)運(yùn)輸成本 transport_cost lp.value(transport_prob.objective) # 4. 計(jì)算總成本 fixed_cost np.dot(F, y_candidate) total_cost fixed_cost transport_cost # 5. 更新全局最優(yōu)解 if total_cost best_cost: best_cost total_cost best_y y_candidate.copy() # 提取最優(yōu)運(yùn)輸方案 best_x np.zeros((M, N)) for (i, j), var in x_vars.items(): best_x[i, j] lp.value(var) print(fIteration {iteration1}: Found new best cost {best_cost:.2f}) # 記錄歷史最佳 history.append(best_cost) return best_y, best_x, best_cost, history # 示例數(shù)據(jù)與調(diào)用 if __name__ __main__: # 設(shè)置問題規(guī)模 M 5 # 5個(gè)候選設(shè)施 N 10 # 10個(gè)需求點(diǎn) # 隨機(jī)生成數(shù)據(jù)可替換為實(shí)際數(shù)據(jù) np.random.seed(123) F np.random.randint(500, 1500, sizeM) # 固定成本 C np.random.randint(200, 500, sizeM) # 設(shè)施容量 d np.random.randint(50, 150, sizeN) # 需求量 c np.random.rand(M, N) * 10 1 # 運(yùn)輸成本1到11之間 print(問題參數(shù)摘要:) print(f固定成本 F: {F}) print(f設(shè)施容量 C: {C}) print(f總需求 sum(d): {sum(d)} 總?cè)萘?sum(C): {sum(C)}) print(f運(yùn)輸成本矩陣 c 的形狀: {c.shape}) # 運(yùn)行隨機(jī)搜索 start_time time.time() best_y, best_x, best_cost, history solve_facility_location_with_random_search( M, N, F.tolist(), C.tolist(), d.tolist(), c.tolist(), max_iter500 ) end_time time.time() print(\n 隨機(jī)搜索結(jié)果 ) print(f搜索耗時(shí): {end_time - start_time:.2f} 秒) print(f最優(yōu)總成本: {best_cost:.2f}) print(f選址方案 (y): {best_y}) print(運(yùn)輸方案 (x):) # 只打印有運(yùn)輸量的路徑 for i in range(M): for j in range(N): if best_x[i, j] 1e-6: # 忽略極小的數(shù)值浮點(diǎn)誤差 print(f 從設(shè)施 {i} 到需求點(diǎn) {j}: {best_x[i, j]:.1f} 單位)代碼關(guān)鍵點(diǎn)解析可行性檢查前置在調(diào)用線性規(guī)劃求解器之前我們先進(jìn)行快速的容量可行性檢查 (if total_capacity total_demand)。這是一個(gè)非常重要的優(yōu)化因?yàn)榍蠼庖粋€(gè)線性規(guī)劃問題比生成一個(gè)隨機(jī)向量和做點(diǎn)乘計(jì)算要昂貴得多。提前拒絕明顯不可行的方案能極大提升搜索效率。動(dòng)態(tài)創(chuàng)建變量在構(gòu)建運(yùn)輸子問題時(shí)我們只為那些被選中的設(shè)施 (y_candidate[i]1) 創(chuàng)建運(yùn)輸變量x_{ij}。這減少了問題的規(guī)模加快了求解速度。使用專業(yè)求解器內(nèi)層運(yùn)輸問題我們使用了pulp調(diào)用CBC求解器。這保證了對(duì)于任何一個(gè)可行的選址方案我們都能得到其精確的最優(yōu)運(yùn)輸成本和分配方案。這是隨機(jī)搜索能有效工作的基礎(chǔ)。結(jié)果記錄我們記錄了每次迭代后的歷史最佳成本history這可以用來繪制算法收斂曲線直觀地看到搜索進(jìn)程。運(yùn)行這段代碼你會(huì)看到算法在隨機(jī)嘗試不同的選址組合并不斷報(bào)告找到的更優(yōu)解。最終輸出的best_y告訴你應(yīng)該建設(shè)哪幾個(gè)倉庫best_x告訴你具體的物資調(diào)運(yùn)方案。4. 性能提升從“純隨機(jī)”到“智能隨機(jī)”基礎(chǔ)的均勻隨機(jī)搜索雖然有效但效率可能不高因?yàn)樗耆珱]有利用歷史搜索到的“好解”的任何信息。在實(shí)際應(yīng)用中尤其是函數(shù)評(píng)估非常耗時(shí)例如每次評(píng)估需要運(yùn)行一個(gè)復(fù)雜的仿真模型時(shí)我們需要讓隨機(jī)搜索變得更“聰明”一些。以下是幾種常用的改進(jìn)思路你可以根據(jù)具體問題選擇或組合使用。4.1 增加局部搜索兩階段混合策略思路很簡(jiǎn)單先用全局隨機(jī)搜索找到一個(gè)不錯(cuò)的“粗解”然后在這個(gè)解的附近進(jìn)行更精細(xì)的局部搜索以期找到更好的解。對(duì)于我們的選址問題局部搜索可以這樣操作擾動(dòng)Perturbation對(duì)于一個(gè)當(dāng)前最優(yōu)的選址方案best_y隨機(jī)翻轉(zhuǎn)其中少數(shù)幾個(gè)y_i的值比如1變成0或0變成1。這相當(dāng)于在當(dāng)前的解附近探索。評(píng)估對(duì)擾動(dòng)后產(chǎn)生的新方案進(jìn)行同樣的可行性檢查和運(yùn)輸問題求解。接受準(zhǔn)則如果新方案成本更低則接受它作為新的當(dāng)前最優(yōu)如果成本更高可以按一定概率接受模擬退火思想或直接拒絕最速下降思想。def local_search_around(best_y, best_cost, F, C, d, c, max_local_trials50): 在最優(yōu)解best_y附近進(jìn)行局部搜索。 M len(best_y) current_y best_y.copy() current_cost best_cost for _ in range(max_local_trials): # 隨機(jī)擾動(dòng)隨機(jī)選擇1到2個(gè)位置進(jìn)行翻轉(zhuǎn) new_y current_y.copy() num_flips np.random.randint(1, 3) flip_indices np.random.choice(M, sizenum_flips, replaceFalse) for idx in flip_indices: new_y[idx] 1 - new_y[idx] # 翻轉(zhuǎn) 0-1 # 檢查新解的可行性并計(jì)算成本 if np.sum(new_y) 0: continue total_capacity np.dot(C, new_y) if total_capacity sum(d): continue # 求解新解的運(yùn)輸子問題此處省略具體求解代碼與主函數(shù)類似 new_cost calculate_total_cost_for_y(new_y, F, C, d, c) # 假設(shè)有這個(gè)函數(shù) # 接受更優(yōu)解 if new_cost current_cost: current_y new_y current_cost new_cost print(f 局部搜索找到更優(yōu)解: cost {new_cost:.2f}) return current_y, current_cost在主隨機(jī)搜索循環(huán)結(jié)束后調(diào)用這個(gè)局部搜索函數(shù)可以對(duì)最終結(jié)果進(jìn)行一次“微調(diào)”。4.2 自適應(yīng)調(diào)整搜索區(qū)域如果發(fā)現(xiàn)隨機(jī)搜索在前期很快找到了一個(gè)較好的區(qū)域后期的隨機(jī)采樣大部分都落在較差的區(qū)域那么可以動(dòng)態(tài)調(diào)整采樣的概率分布。例如可以記錄下那些成本較低的解對(duì)應(yīng)的y_i1的概率然后在后續(xù)的隨機(jī)生成中讓每個(gè)設(shè)施被選中的概率向這個(gè)歷史經(jīng)驗(yàn)概率靠攏。這有點(diǎn)類似“交叉熵方法”或“分布估計(jì)算法”的思想。不過對(duì)于0-1組合問題實(shí)現(xiàn)起來需要更精細(xì)的設(shè)計(jì)以避免過早收斂到局部最優(yōu)。4.3 并行化計(jì)算這是提升隨機(jī)搜索效率最直接、最有效的方法尤其適合數(shù)學(xué)建模競(jìng)賽中可能使用的多核計(jì)算機(jī)。由于每次迭代生成一個(gè)候選解并評(píng)估完全獨(dú)立我們可以輕松地將max_iter次迭代分配到多個(gè)進(jìn)程上。使用Python的multiprocessing庫或concurrent.futures模塊可以方便地實(shí)現(xiàn)?；舅悸肥菍⒖偟牡螖?shù)分成若干份交給多個(gè)工作進(jìn)程同時(shí)執(zhí)行各自的隨機(jī)搜索每個(gè)進(jìn)程獨(dú)立維護(hù)自己的“當(dāng)前最優(yōu)解”。最后從所有進(jìn)程返回的結(jié)果中挑選出全局最優(yōu)的那個(gè)。from concurrent.futures import ProcessPoolExecutor, as_completed def random_search_worker(task_args): 每個(gè)工作進(jìn)程執(zhí)行的函數(shù) M, N, F, C, d, c, iterations, seed task_args # ... 執(zhí)行指定迭代次數(shù)的隨機(jī)搜索 ... return local_best_y, local_best_x, local_best_cost # 在主程序中 if __name__ __main__: num_workers 4 iterations_per_worker max_iter // num_workers with ProcessPoolExecutor(max_workersnum_workers) as executor: futures [] for i in range(num_workers): # 為每個(gè)worker分配不同的隨機(jī)種子確保獨(dú)立性 task (M, N, F, C, d, c, iterations_per_worker, 12345i) future executor.submit(random_search_worker, task) futures.append(future) # 收集結(jié)果 results [] for future in as_completed(futures): results.append(future.result()) # 從所有worker的結(jié)果中找出全局最優(yōu) global_best_cost float(inf) global_best_y None for y, x, cost in results: if cost global_best_cost: global_best_cost cost global_best_y y global_best_x x通過并行化你可以幾乎線性地減少搜索時(shí)間假設(shè)CPU核心充足。在時(shí)間緊迫的數(shù)學(xué)建模比賽中這可能是決定你能否在截止前跑出結(jié)果的關(guān)鍵。5. 實(shí)戰(zhàn)心得與避坑指南經(jīng)過多個(gè)項(xiàng)目和比賽的使用我總結(jié)了以下幾點(diǎn)關(guān)于在數(shù)學(xué)建模中應(yīng)用隨機(jī)搜索算法的經(jīng)驗(yàn)和教訓(xùn)1. 它不是萬能的要明確適用場(chǎng)景隨機(jī)搜索最適合作為“基線方法”或“最后的手段”。當(dāng)你的問題滿足以下條件時(shí)優(yōu)先考慮它目標(biāo)函數(shù)或約束條件不可微、不連續(xù)、評(píng)估代價(jià)高。問題規(guī)模中等但結(jié)構(gòu)復(fù)雜傳統(tǒng)優(yōu)化器難以建?；蚯蠼?。你對(duì)解的最優(yōu)性要求不是“絕對(duì)精確”而是“足夠好、可用”。 如果問題有明顯的數(shù)學(xué)結(jié)構(gòu)如凸性、線性或者有成熟的專用算法應(yīng)該優(yōu)先使用那些方法。2. 迭代次數(shù)與解的質(zhì)量是概率關(guān)系隨機(jī)搜索的性能嚴(yán)重依賴于采樣次數(shù)。理論上采樣點(diǎn)越多找到更好解的概率越大。你需要做一個(gè)權(quán)衡評(píng)估一次目標(biāo)函數(shù)需要多長(zhǎng)時(shí)間你總共有多少計(jì)算時(shí)間通常我會(huì)先做一個(gè)快速的“偵察跑”設(shè)置一個(gè)較小的迭代次數(shù)比如1000次看看成本下降的趨勢(shì)。如果成本在幾百次迭代后就不再顯著改善可能說明當(dāng)前的搜索空間下解的質(zhì)量已經(jīng)接近極限或者算法陷入了某個(gè)區(qū)域。如果成本還在持續(xù)緩慢下降那么增加迭代次數(shù)很可能帶來收益。3. 隨機(jī)種子的影響與統(tǒng)計(jì)評(píng)估由于算法的隨機(jī)性單次運(yùn)行的結(jié)果具有偶然性。務(wù)必多次運(yùn)行例如30次記錄每次找到的最優(yōu)解和成本。然后你可以報(bào)告最好解、最差解、平均解、解的標(biāo)準(zhǔn)差。這比只報(bào)告一次運(yùn)行的結(jié)果要嚴(yán)謹(jǐn)?shù)枚?。在論文中你可以說“我們獨(dú)立運(yùn)行算法30次最佳結(jié)果為XXX平均結(jié)果為YYY±ZZZ”這體現(xiàn)了方法的魯棒性。4. 可行性檢查是效率的關(guān)鍵如前文代碼所示在調(diào)用耗時(shí)的精確求解器或復(fù)雜仿真之前盡可能用簡(jiǎn)單、快速的條件過濾掉不可行的候選解。對(duì)于選址問題容量檢查就是這樣一個(gè)“守門員”。在其他問題中可能是變量的簡(jiǎn)單邊界檢查或者一些必須滿足的硬性邏輯約束。每過濾掉一個(gè)不可行解就節(jié)省了一次昂貴的評(píng)估。5. 與精確解或已知下界對(duì)比如果問題規(guī)模較小可以嘗試用商業(yè)求解器如Gurobi求出精確最優(yōu)解作為對(duì)比的“黃金標(biāo)準(zhǔn)”。如果求不出精確解可以嘗試計(jì)算一個(gè)問題的下界例如線性規(guī)劃松弛的解。將隨機(jī)搜索得到的最好解與精確解或下界進(jìn)行比較可以量化你的近似解的質(zhì)量。例如“我們的隨機(jī)搜索算法在5000次迭代內(nèi)找到的解與問題下界的差距在5%以內(nèi)”這是一個(gè)非常有說服力的結(jié)果。6. 可視化搜索過程將每次迭代找到的“當(dāng)前最優(yōu)成本”記錄下來并繪圖是分析算法行為的利器。你可以看到成本是如何隨著迭代下降的下降的速度如何何時(shí)趨于平穩(wěn)。這張圖放在論文的附錄或正文中能直觀地展示算法的收斂性。如果曲線下降很快然后平緩說明算法初期探索有效如果曲線一直緩慢下降說明可能需要更多迭代或改進(jìn)采樣策略。最后隨機(jī)搜索的魅力在于其簡(jiǎn)單性與強(qiáng)大通用性之間的平衡。它可能不是最優(yōu)雅、最快速的算法但在面對(duì)數(shù)學(xué)建模中那些“不講武德”的復(fù)雜現(xiàn)實(shí)問題時(shí)它往往是最忠實(shí)、最可靠的伙伴。掌握它意味著你在優(yōu)化工具箱里又多了一件應(yīng)對(duì)不確定性的利器。下次當(dāng)你面對(duì)一個(gè)看似無從下手的多變量?jī)?yōu)化模型時(shí)不妨先試試隨機(jī)搜索讓它為你照亮一片可能的解空間或許驚喜就在其中。 數(shù)學(xué)建模優(yōu)化難題破解：隨機(jī)搜索算法原理與應(yīng)急選址實(shí)戰(zhàn)
1. 項(xiàng)目緣起當(dāng)數(shù)學(xué)建模遇上“大海撈針”做數(shù)學(xué)建模尤其是國賽、美賽這種高強(qiáng)度的比賽最讓人頭疼的環(huán)節(jié)之一往往不是模型建立而是模型求解。你花了大量心血構(gòu)建了一個(gè)邏輯嚴(yán)密、變量眾多的非線性規(guī)劃模型感覺已經(jīng)成功了一大半。但當(dāng)你打開MATLAB或者Python準(zhǔn)備用那些教科書上的經(jīng)典算法比如梯度下降、牛頓法去求解時(shí)卻發(fā)現(xiàn)要么是初始點(diǎn)敏感動(dòng)不動(dòng)就掉進(jìn)局部最優(yōu)的“坑”里出不來要么是目標(biāo)函數(shù)或約束條件太復(fù)雜導(dǎo)數(shù)都求不出來算法直接“罷工”。這時(shí)候你需要的可能不是更精密的微調(diào)而是一種思路上的轉(zhuǎn)變——一種不那么“聰明”但足夠“魯棒”和“通用”的方法。這就是隨機(jī)搜索算法Random Search Algorithm登場(chǎng)的時(shí)刻。很多人一聽到“隨機(jī)搜索”第一反應(yīng)可能是“這不就是瞎蒙嗎”。確實(shí)它的核心思想非常樸素在變量的定義域可行域內(nèi)隨機(jī)地生成大量的候選解然后從中找出目標(biāo)函數(shù)值最好的那個(gè)。聽起來毫無技術(shù)含量甚至有些“暴力”。但正是這種“暴力”在應(yīng)對(duì)多變量、非線性、非凸、甚至“黑箱”函數(shù)的最優(yōu)化問題時(shí)展現(xiàn)出了驚人的實(shí)用性。它不依賴于函數(shù)的梯度信息對(duì)函數(shù)的形態(tài)幾乎沒有要求只要你能計(jì)算出任何一個(gè)點(diǎn)的函數(shù)值它就能工作。在數(shù)學(xué)建模競(jìng)賽中當(dāng)你面對(duì)一個(gè)結(jié)構(gòu)復(fù)雜、機(jī)理不甚清晰的實(shí)際問題需要快速得到一個(gè)“還不錯(cuò)”的、可用的解時(shí)隨機(jī)搜索往往能成為你的救命稻草。我最初深入使用隨機(jī)搜索是在一次模擬復(fù)雜供應(yīng)鏈網(wǎng)絡(luò)優(yōu)化的項(xiàng)目中。模型涉及幾十個(gè)決策變量如倉庫選址、運(yùn)輸路徑、庫存水平約束條件相互耦合目標(biāo)函數(shù)總成本是高度非線性的。嘗試了多種基于梯度的優(yōu)化器效果都不理想。最后抱著試試看的心態(tài)寫了一個(gè)簡(jiǎn)單的隨機(jī)搜索程序雖然每次運(yùn)行結(jié)果都有波動(dòng)但總能穩(wěn)定地找到一個(gè)比之前任何方法都好的解。從那以后隨機(jī)搜索就成了我解決棘手優(yōu)化問題的“標(biāo)準(zhǔn)備選方案”之一。本文將結(jié)合一個(gè)具體的建模案例拆解隨機(jī)搜索從原理到實(shí)現(xiàn)的每一個(gè)細(xì)節(jié)并分享我在實(shí)戰(zhàn)中積累的調(diào)參經(jīng)驗(yàn)和避坑指南。2. 隨機(jī)搜索的核心原理為什么“笨辦法”有時(shí)更有效要理解隨機(jī)搜索的價(jià)值我們得先看看它的對(duì)手們?yōu)槭裁磿?huì)在某些場(chǎng)景下失效。傳統(tǒng)的最優(yōu)化算法如梯度下降法其哲學(xué)是“局部尋優(yōu)”。它假設(shè)函數(shù)是光滑的、可微的并且當(dāng)前點(diǎn)的梯度方向能指引我們?nèi)ネ粋€(gè)更優(yōu)的鄰域。這就像在一個(gè)連綿起伏的山丘中你蒙著眼睛但能用腳感受地面的傾斜梯度然后朝著感覺是下坡的方向走。這個(gè)方法在“山丘”形狀良好凸函數(shù)時(shí)非常高效能快速找到谷底全局最優(yōu)。然而現(xiàn)實(shí)世界的優(yōu)化問題尤其是數(shù)學(xué)建模中抽象出來的問題其“地形”往往更像一片布滿深坑、斷崖和緩坡的復(fù)雜地貌非凸函數(shù)。梯度下降法從某個(gè)起點(diǎn)出發(fā)很容易掉進(jìn)最近的一個(gè)坑里局部最優(yōu)并且由于在坑底梯度為零或無法計(jì)算它就認(rèn)為已經(jīng)到達(dá)終點(diǎn)再也出不來了。這就是“局部最優(yōu)陷阱”。隨機(jī)搜索采取了完全不同的策略。它放棄了“局部感知”轉(zhuǎn)而進(jìn)行“全局采樣”。想象一下你不是那個(gè)蒙眼下坡的人而是指揮成千上萬個(gè)無人機(jī)在整個(gè)區(qū)域上空隨機(jī)空投傳感器。每個(gè)傳感器落地后立刻報(bào)告所在地的海拔高度函數(shù)值。最后你只需從所有報(bào)告中找出海拔最低的那個(gè)點(diǎn)。這個(gè)方法不關(guān)心地形是否連續(xù)、是否可導(dǎo)它只關(guān)心1. 我能否定義出整個(gè)搜索區(qū)域變量上下界2. 我能否對(duì)區(qū)域內(nèi)的任何一個(gè)點(diǎn)進(jìn)行評(píng)估。2.1 算法流程的標(biāo)準(zhǔn)化描述一個(gè)最基礎(chǔ)的多變量隨機(jī)搜索算法可以歸納為以下幾步定義問題明確你的目標(biāo)是最小化還是最大化一個(gè)目標(biāo)函數(shù) f(x)。確定決策變量向量 x [x1, x2, ..., xn] 的搜索空間通常由每個(gè)變量的下限 Lower Bound (LB) 和上限 Upper Bound (UB) 定義即 LB_i ≤ x_i ≤ UB_i。初始化設(shè)定算法要進(jìn)行的總迭代次數(shù)max_iter或總函數(shù)評(píng)估次數(shù)max_evaluations。初始化當(dāng)前最優(yōu)解best_x和對(duì)應(yīng)的最優(yōu)值best_f為一個(gè)很差的初始值例如對(duì)于最小化問題best_f設(shè)為無窮大。迭代搜索 a.生成候選點(diǎn)在當(dāng)前迭代k在變量的定義域內(nèi)均勻隨機(jī)地生成一個(gè)候選解candidate_x。即對(duì)于每個(gè)變量 icandidate_x[i] LB_i random() * (UB_i - LB_i)其中random()生成一個(gè)[0,1)之間的隨機(jī)數(shù)。 b.可行性檢查可選如果問題存在除變量上下界外的其他復(fù)雜約束如線性/非線性不等式約束需要檢查candidate_x是否滿足所有約束。如果不滿足則丟棄該點(diǎn)回到步驟a生成下一個(gè)隨機(jī)點(diǎn)或采取懲罰函數(shù)等策略。對(duì)于入門我們先處理無復(fù)雜約束的問題。 c.性能評(píng)估計(jì)算該候選點(diǎn)的目標(biāo)函數(shù)值current_f f(candidate_x)。 d.更新最優(yōu)解將current_f與歷史最優(yōu)值best_f比較。如果更優(yōu)對(duì)于最小化問題是更小則更新best_f current_f同時(shí)更新best_x candidate_x。終止與輸出重復(fù)步驟3直到達(dá)到預(yù)設(shè)的迭代次數(shù)或評(píng)估次數(shù)。最終輸出找到的best_x和best_f。這個(gè)流程簡(jiǎn)單到幾乎可以用任何編程語言在十幾行代碼內(nèi)實(shí)現(xiàn)。但正是這種簡(jiǎn)單性帶來了幾個(gè)關(guān)鍵優(yōu)勢(shì)全局性由于采樣是全局均勻的算法有概率采樣到整個(gè)可行域的任何角落因此理論上只要采樣點(diǎn)足夠多就有機(jī)會(huì)逼近甚至命中全局最優(yōu)解。魯棒性對(duì)目標(biāo)函數(shù) f(x) 的性質(zhì)幾乎無要求。f(x) 可以是離散的、不可微的、不連續(xù)的甚至是一個(gè)需要調(diào)用外部仿真軟件才能得到結(jié)果的“黑箱”函數(shù)。易并行性每一次隨機(jī)采樣和評(píng)估都是完全獨(dú)立的可以非常容易地分配到多個(gè)CPU核心或計(jì)算節(jié)點(diǎn)上并行執(zhí)行從而大幅縮短搜索時(shí)間。注意隨機(jī)搜索找到的“最優(yōu)解”是一個(gè)隨機(jī)變量。每次運(yùn)行的結(jié)果都可能不同。我們通常通過多次獨(dú)立運(yùn)行取最好的結(jié)果或統(tǒng)計(jì)結(jié)果的分布來評(píng)估算法的性能。3. 建模案例實(shí)戰(zhàn)應(yīng)急物資儲(chǔ)備庫的選址優(yōu)化為了讓大家更直觀地理解隨機(jī)搜索如何應(yīng)用于實(shí)際建模我們構(gòu)造一個(gè)簡(jiǎn)化但經(jīng)典的運(yùn)籌學(xué)問題多需求點(diǎn)應(yīng)急物資儲(chǔ)備庫選址優(yōu)化。3.1 問題描述與模型建立假設(shè)某地區(qū)有M個(gè)潛在的應(yīng)急物資儲(chǔ)備庫選址點(diǎn)需要服務(wù)N個(gè)已知的居民點(diǎn)需求點(diǎn)。每個(gè)居民點(diǎn)j有一個(gè)固定的物資年需求量d_j。每個(gè)候選儲(chǔ)備庫i有一個(gè)最大的建設(shè)容量C_i和一個(gè)固定的年運(yùn)營建設(shè)成本F_i只要選中建設(shè)就會(huì)產(chǎn)生此成本。從儲(chǔ)備庫i運(yùn)輸單位物資到居民點(diǎn)j的運(yùn)費(fèi)為c_{ij}。我們的決策是選址決策決定在哪些候選點(diǎn)建設(shè)儲(chǔ)備庫。用一個(gè)0-1變量y_i表示y_i 1表示在點(diǎn)i建設(shè)y_i 0表示不建設(shè)。分配決策決定每個(gè)建設(shè)的儲(chǔ)備庫i向每個(gè)居民點(diǎn)j運(yùn)輸多少物資。用一個(gè)連續(xù)變量x_{ij}表示。目標(biāo)是最小化總成本包括所有被選中的儲(chǔ)備庫的固定建設(shè)成本以及所有的運(yùn)輸成本。數(shù)學(xué)模型如下目標(biāo)函數(shù)最小化總成本Minimize Z Σ_{i1}^{M} (F_i * y_i) Σ_{i1}^{M} Σ_{j1}^{N} (c_{ij} * x_{ij})第一項(xiàng)是固定成本第二項(xiàng)是運(yùn)輸成本。約束條件需求滿足約束每個(gè)居民點(diǎn)j的需求必須被完全滿足。Σ_{i1}^{M} x_{ij} d_j, for all j 1, 2, ..., N容量約束每個(gè)儲(chǔ)備庫i發(fā)出的物資總量不能超過其建設(shè)容量。Σ_{j1}^{N} x_{ij} ≤ C_i * y_i, for all i 1, 2, ..., M 注意如果y_i 0則右側(cè)為0意味著從該點(diǎn)運(yùn)出的物資x_{ij}必須全部為0如果y_i 1則不能超過C_i。邏輯約束只有被選中的儲(chǔ)備庫才能分配物資。x_{ij} ≥ 0, for all i, jy_i ∈ {0, 1}, for all i這是一個(gè)典型的**混合整數(shù)線性規(guī)劃MILP**問題。對(duì)于小規(guī)模問題可以使用專業(yè)的優(yōu)化求解器如Gurobi, CPLEX精確求解。但在數(shù)學(xué)建模競(jìng)賽中問題規(guī)模可能較大或者環(huán)境限制無法使用商業(yè)求解器。此時(shí)隨機(jī)搜索就可以作為一個(gè)有效的近似求解方案。3.2 隨機(jī)搜索求解策略設(shè)計(jì)直接對(duì)所有的y_i和x_{ij}進(jìn)行隨機(jī)搜索效率極低因?yàn)樽兞刻嗲掖嬖趶?fù)雜的約束關(guān)系。我們需要利用問題的結(jié)構(gòu)設(shè)計(jì)更聰明的搜索策略。一個(gè)有效的策略是將問題分解外層搜索隨機(jī)搜索負(fù)責(zé)搜索y_i的0-1組合即決定“建哪些庫”。這是一個(gè)組合優(yōu)化問題。內(nèi)層求解線性規(guī)劃/運(yùn)輸問題對(duì)于外層給定的一個(gè)具體的選址方案y即確定了哪些y_i1剩下的問題是一個(gè)標(biāo)準(zhǔn)的運(yùn)輸問題在已選定的倉庫集合內(nèi)如何分配x_{ij}以滿足所有需求且不超出選定倉庫的容量并使運(yùn)輸成本最低。這個(gè)問題是線性規(guī)劃有高效算法如單純形法可以快速精確求解甚至對(duì)于簡(jiǎn)單情況可以直接用線性規(guī)劃求解器或自己編寫算法如表上作業(yè)法。算法流程調(diào)整如下編碼一個(gè)解表示為一個(gè)長(zhǎng)度為M的0-1向量代表y_1, y_2, ..., y_M。生成候選選址方案隨機(jī)生成一個(gè)0-1向量。為了增加可行性可以加入啟發(fā)式規(guī)則比如至少生成一個(gè)y_i1?？尚行赃^濾檢查隨機(jī)生成的選址方案是否容量可行。即所有被選中倉庫的總?cè)萘喀?(C_i * y_i)是否大于等于總需求Σ d_j。如果不滿足這個(gè)方案不可能滿足所有需求直接丟棄重新生成。求解子問題對(duì)于容量可行的選址方案將其y_i值固定求解內(nèi)層的運(yùn)輸問題得到最優(yōu)的x_{ij}分配和對(duì)應(yīng)的最小運(yùn)輸成本Transport_Cost(y)。計(jì)算總成本總成本Z Σ (F_i * y_i) Transport_Cost(y)。更新全局最優(yōu)比較并更新。這樣隨機(jī)搜索的核心就變成了在指數(shù)級(jí)數(shù)量的選址組合中尋找能使“固定成本對(duì)應(yīng)最優(yōu)運(yùn)輸成本”最小的那個(gè)組合。內(nèi)層運(yùn)輸問題的精確求解保證了對(duì)于任何一個(gè)選址方案我們都能得到其可能達(dá)到的最低運(yùn)輸成本從而公平地比較不同選址方案的優(yōu)劣。3.3 Python代碼實(shí)現(xiàn)與解析下面我們用Python來實(shí)現(xiàn)這個(gè)策略。我們將使用numpy生成隨機(jī)數(shù)并使用pulp一個(gè)免費(fèi)的線性規(guī)劃庫來求解內(nèi)層的運(yùn)輸問題。pulp不是標(biāo)準(zhǔn)庫需要安裝pip install pulp。import numpy as np import pulp as lp import time def solve_facility_location_with_random_search(M, N, F, C, d, c, max_iter1000, seed42): 使用隨機(jī)搜索算法求解設(shè)施選址問題。 參數(shù): M: 潛在設(shè)施數(shù)量 N: 需求點(diǎn)數(shù)量 F: list of length M, 每個(gè)設(shè)施的固定成本 C: list of length M, 每個(gè)設(shè)施的容量 d: list of length N, 每個(gè)需求點(diǎn)的需求量 c: 2D list of shape (M, N), 運(yùn)輸成本矩陣c[i][j] 從設(shè)施i到需求點(diǎn)j的成本 max_iter: 最大隨機(jī)搜索迭代次數(shù) seed: 隨機(jī)種子用于復(fù)現(xiàn)結(jié)果 返回: best_y: 最優(yōu)的選址方案 (0-1 list) best_x: 最優(yōu)的運(yùn)輸方案 (2D list) best_cost: 最優(yōu)總成本 history: 每次迭代找到的最佳成本記錄 np.random.seed(seed) total_demand sum(d) # 初始化最優(yōu)解 best_cost float(inf) best_y None best_x None history [] for iteration in range(max_iter): # 1. 隨機(jī)生成一個(gè)選址方案 y (0-1向量) y_candidate np.random.randint(0, 2, sizeM) # 2. 可行性檢查確保至少選一個(gè)設(shè)施且總?cè)萘? 總需求 if np.sum(y_candidate) 0: continue # 沒選任何設(shè)施不可行 total_capacity np.dot(C, y_candidate) if total_capacity total_demand: continue # 容量不足不可行 # 3. 構(gòu)建并求解內(nèi)層運(yùn)輸問題 # 創(chuàng)建問題實(shí)例最小化 transport_prob lp.LpProblem(Transportation_Subproblem, lp.LpMinimize) # 創(chuàng)建決策變量 x[i][j] 0 x_vars lp.LpVariable.dicts(x, ((i, j) for i in range(M) for j in range(N) if y_candidate[i] 1), lowBound0, catContinuous) # 如果設(shè)施i未被選中(y_candidate[i]0)則對(duì)應(yīng)的x變量不會(huì)被創(chuàng)建天然為0。 # 目標(biāo)函數(shù)最小化運(yùn)輸成本 transport_prob lp.lpSum(c[i][j] * x_vars[(i, j)] for (i, j) in x_vars.keys()) # 約束條件1: 每個(gè)需求點(diǎn)的需求必須滿足 for j in range(N): # 對(duì)所有選中的設(shè)施i求和其運(yùn)往j的物資量 prob lp.lpSum(x_vars.get((i, j), 0) for i in range(M) if y_candidate[i] 1) d[j] # 約束條件2: 每個(gè)選中設(shè)施的運(yùn)出量不超過其容量 for i in range(M): if y_candidate[i] 1: prob lp.lpSum(x_vars.get((i, j), 0) for j in range(N)) C[i] # 求解運(yùn)輸子問題 # pulp默認(rèn)使用CBC求解器對(duì)于線性規(guī)劃足夠 transport_prob.solve(pulp.PULP_CBC_CMD(msgFalse)) # msgFalse關(guān)閉求解器輸出 # 檢查是否求解成功 if lp.LpStatus[transport_prob.status] ! Optimal: # 如果子問題無解理論上在容量可行下應(yīng)該不會(huì)發(fā)生但保留檢查 continue # 獲取子問題最優(yōu)運(yùn)輸成本 transport_cost lp.value(transport_prob.objective) # 4. 計(jì)算總成本 fixed_cost np.dot(F, y_candidate) total_cost fixed_cost transport_cost # 5. 更新全局最優(yōu)解 if total_cost best_cost: best_cost total_cost best_y y_candidate.copy() # 提取最優(yōu)運(yùn)輸方案 best_x np.zeros((M, N)) for (i, j), var in x_vars.items(): best_x[i, j] lp.value(var) print(fIteration {iteration1}: Found new best cost {best_cost:.2f}) # 記錄歷史最佳 history.append(best_cost) return best_y, best_x, best_cost, history # 示例數(shù)據(jù)與調(diào)用 if __name__ __main__: # 設(shè)置問題規(guī)模 M 5 # 5個(gè)候選設(shè)施 N 10 # 10個(gè)需求點(diǎn) # 隨機(jī)生成數(shù)據(jù)可替換為實(shí)際數(shù)據(jù) np.random.seed(123) F np.random.randint(500, 1500, sizeM) # 固定成本 C np.random.randint(200, 500, sizeM) # 設(shè)施容量 d np.random.randint(50, 150, sizeN) # 需求量 c np.random.rand(M, N) * 10 1 # 運(yùn)輸成本1到11之間 print(問題參數(shù)摘要:) print(f固定成本 F: {F}) print(f設(shè)施容量 C: {C}) print(f總需求 sum(d): {sum(d)} 總?cè)萘?sum(C): {sum(C)}) print(f運(yùn)輸成本矩陣 c 的形狀: {c.shape}) # 運(yùn)行隨機(jī)搜索 start_time time.time() best_y, best_x, best_cost, history solve_facility_location_with_random_search( M, N, F.tolist(), C.tolist(), d.tolist(), c.tolist(), max_iter500 ) end_time time.time() print(\n 隨機(jī)搜索結(jié)果 ) print(f搜索耗時(shí): {end_time - start_time:.2f} 秒) print(f最優(yōu)總成本: {best_cost:.2f}) print(f選址方案 (y): {best_y}) print(運(yùn)輸方案 (x):) # 只打印有運(yùn)輸量的路徑 for i in range(M): for j in range(N): if best_x[i, j] 1e-6: # 忽略極小的數(shù)值浮點(diǎn)誤差 print(f 從設(shè)施 {i} 到需求點(diǎn) {j}: {best_x[i, j]:.1f} 單位)代碼關(guān)鍵點(diǎn)解析可行性檢查前置在調(diào)用線性規(guī)劃求解器之前我們先進(jìn)行快速的容量可行性檢查 (if total_capacity total_demand)。這是一個(gè)非常重要的優(yōu)化因?yàn)榍蠼庖粋€(gè)線性規(guī)劃問題比生成一個(gè)隨機(jī)向量和做點(diǎn)乘計(jì)算要昂貴得多。提前拒絕明顯不可行的方案能極大提升搜索效率。動(dòng)態(tài)創(chuàng)建變量在構(gòu)建運(yùn)輸子問題時(shí)我們只為那些被選中的設(shè)施 (y_candidate[i]1) 創(chuàng)建運(yùn)輸變量x_{ij}。這減少了問題的規(guī)模加快了求解速度。使用專業(yè)求解器內(nèi)層運(yùn)輸問題我們使用了pulp調(diào)用CBC求解器。這保證了對(duì)于任何一個(gè)可行的選址方案我們都能得到其精確的最優(yōu)運(yùn)輸成本和分配方案。這是隨機(jī)搜索能有效工作的基礎(chǔ)。結(jié)果記錄我們記錄了每次迭代后的歷史最佳成本history這可以用來繪制算法收斂曲線直觀地看到搜索進(jìn)程。運(yùn)行這段代碼你會(huì)看到算法在隨機(jī)嘗試不同的選址組合并不斷報(bào)告找到的更優(yōu)解。最終輸出的best_y告訴你應(yīng)該建設(shè)哪幾個(gè)倉庫best_x告訴你具體的物資調(diào)運(yùn)方案。4. 性能提升從“純隨機(jī)”到“智能隨機(jī)”基礎(chǔ)的均勻隨機(jī)搜索雖然有效但效率可能不高因?yàn)樗耆珱]有利用歷史搜索到的“好解”的任何信息。在實(shí)際應(yīng)用中尤其是函數(shù)評(píng)估非常耗時(shí)例如每次評(píng)估需要運(yùn)行一個(gè)復(fù)雜的仿真模型時(shí)我們需要讓隨機(jī)搜索變得更“聰明”一些。以下是幾種常用的改進(jìn)思路你可以根據(jù)具體問題選擇或組合使用。4.1 增加局部搜索兩階段混合策略思路很簡(jiǎn)單先用全局隨機(jī)搜索找到一個(gè)不錯(cuò)的“粗解”然后在這個(gè)解的附近進(jìn)行更精細(xì)的局部搜索以期找到更好的解。對(duì)于我們的選址問題局部搜索可以這樣操作擾動(dòng)Perturbation對(duì)于一個(gè)當(dāng)前最優(yōu)的選址方案best_y隨機(jī)翻轉(zhuǎn)其中少數(shù)幾個(gè)y_i的值比如1變成0或0變成1。這相當(dāng)于在當(dāng)前的解附近探索。評(píng)估對(duì)擾動(dòng)后產(chǎn)生的新方案進(jìn)行同樣的可行性檢查和運(yùn)輸問題求解。接受準(zhǔn)則如果新方案成本更低則接受它作為新的當(dāng)前最優(yōu)如果成本更高可以按一定概率接受模擬退火思想或直接拒絕最速下降思想。def local_search_around(best_y, best_cost, F, C, d, c, max_local_trials50): 在最優(yōu)解best_y附近進(jìn)行局部搜索。 M len(best_y) current_y best_y.copy() current_cost best_cost for _ in range(max_local_trials): # 隨機(jī)擾動(dòng)隨機(jī)選擇1到2個(gè)位置進(jìn)行翻轉(zhuǎn) new_y current_y.copy() num_flips np.random.randint(1, 3) flip_indices np.random.choice(M, sizenum_flips, replaceFalse) for idx in flip_indices: new_y[idx] 1 - new_y[idx] # 翻轉(zhuǎn) 0-1 # 檢查新解的可行性并計(jì)算成本 if np.sum(new_y) 0: continue total_capacity np.dot(C, new_y) if total_capacity sum(d): continue # 求解新解的運(yùn)輸子問題此處省略具體求解代碼與主函數(shù)類似 new_cost calculate_total_cost_for_y(new_y, F, C, d, c) # 假設(shè)有這個(gè)函數(shù) # 接受更優(yōu)解 if new_cost current_cost: current_y new_y current_cost new_cost print(f 局部搜索找到更優(yōu)解: cost {new_cost:.2f}) return current_y, current_cost在主隨機(jī)搜索循環(huán)結(jié)束后調(diào)用這個(gè)局部搜索函數(shù)可以對(duì)最終結(jié)果進(jìn)行一次“微調(diào)”。4.2 自適應(yīng)調(diào)整搜索區(qū)域如果發(fā)現(xiàn)隨機(jī)搜索在前期很快找到了一個(gè)較好的區(qū)域后期的隨機(jī)采樣大部分都落在較差的區(qū)域那么可以動(dòng)態(tài)調(diào)整采樣的概率分布。例如可以記錄下那些成本較低的解對(duì)應(yīng)的y_i1的概率然后在后續(xù)的隨機(jī)生成中讓每個(gè)設(shè)施被選中的概率向這個(gè)歷史經(jīng)驗(yàn)概率靠攏。這有點(diǎn)類似“交叉熵方法”或“分布估計(jì)算法”的思想。不過對(duì)于0-1組合問題實(shí)現(xiàn)起來需要更精細(xì)的設(shè)計(jì)以避免過早收斂到局部最優(yōu)。4.3 并行化計(jì)算這是提升隨機(jī)搜索效率最直接、最有效的方法尤其適合數(shù)學(xué)建模競(jìng)賽中可能使用的多核計(jì)算機(jī)。由于每次迭代生成一個(gè)候選解并評(píng)估完全獨(dú)立我們可以輕松地將max_iter次迭代分配到多個(gè)進(jìn)程上。使用Python的multiprocessing庫或concurrent.futures模塊可以方便地實(shí)現(xiàn)?；舅悸肥菍⒖偟牡螖?shù)分成若干份交給多個(gè)工作進(jìn)程同時(shí)執(zhí)行各自的隨機(jī)搜索每個(gè)進(jìn)程獨(dú)立維護(hù)自己的“當(dāng)前最優(yōu)解”。最后從所有進(jìn)程返回的結(jié)果中挑選出全局最優(yōu)的那個(gè)。from concurrent.futures import ProcessPoolExecutor, as_completed def random_search_worker(task_args): 每個(gè)工作進(jìn)程執(zhí)行的函數(shù) M, N, F, C, d, c, iterations, seed task_args # ... 執(zhí)行指定迭代次數(shù)的隨機(jī)搜索 ... return local_best_y, local_best_x, local_best_cost # 在主程序中 if __name__ __main__: num_workers 4 iterations_per_worker max_iter // num_workers with ProcessPoolExecutor(max_workersnum_workers) as executor: futures [] for i in range(num_workers): # 為每個(gè)worker分配不同的隨機(jī)種子確保獨(dú)立性 task (M, N, F, C, d, c, iterations_per_worker, 12345i) future executor.submit(random_search_worker, task) futures.append(future) # 收集結(jié)果 results [] for future in as_completed(futures): results.append(future.result()) # 從所有worker的結(jié)果中找出全局最優(yōu) global_best_cost float(inf) global_best_y None for y, x, cost in results: if cost global_best_cost: global_best_cost cost global_best_y y global_best_x x通過并行化你可以幾乎線性地減少搜索時(shí)間假設(shè)CPU核心充足。在時(shí)間緊迫的數(shù)學(xué)建模比賽中這可能是決定你能否在截止前跑出結(jié)果的關(guān)鍵。5. 實(shí)戰(zhàn)心得與避坑指南經(jīng)過多個(gè)項(xiàng)目和比賽的使用我總結(jié)了以下幾點(diǎn)關(guān)于在數(shù)學(xué)建模中應(yīng)用隨機(jī)搜索算法的經(jīng)驗(yàn)和教訓(xùn)1. 它不是萬能的要明確適用場(chǎng)景隨機(jī)搜索最適合作為“基線方法”或“最后的手段”。當(dāng)你的問題滿足以下條件時(shí)優(yōu)先考慮它目標(biāo)函數(shù)或約束條件不可微、不連續(xù)、評(píng)估代價(jià)高。問題規(guī)模中等但結(jié)構(gòu)復(fù)雜傳統(tǒng)優(yōu)化器難以建?；蚯蠼?。你對(duì)解的最優(yōu)性要求不是“絕對(duì)精確”而是“足夠好、可用”。 如果問題有明顯的數(shù)學(xué)結(jié)構(gòu)如凸性、線性或者有成熟的專用算法應(yīng)該優(yōu)先使用那些方法。2. 迭代次數(shù)與解的質(zhì)量是概率關(guān)系隨機(jī)搜索的性能嚴(yán)重依賴于采樣次數(shù)。理論上采樣點(diǎn)越多找到更好解的概率越大。你需要做一個(gè)權(quán)衡評(píng)估一次目標(biāo)函數(shù)需要多長(zhǎng)時(shí)間你總共有多少計(jì)算時(shí)間通常我會(huì)先做一個(gè)快速的“偵察跑”設(shè)置一個(gè)較小的迭代次數(shù)比如1000次看看成本下降的趨勢(shì)。如果成本在幾百次迭代后就不再顯著改善可能說明當(dāng)前的搜索空間下解的質(zhì)量已經(jīng)接近極限或者算法陷入了某個(gè)區(qū)域。如果成本還在持續(xù)緩慢下降那么增加迭代次數(shù)很可能帶來收益。3. 隨機(jī)種子的影響與統(tǒng)計(jì)評(píng)估由于算法的隨機(jī)性單次運(yùn)行的結(jié)果具有偶然性。務(wù)必多次運(yùn)行例如30次記錄每次找到的最優(yōu)解和成本。然后你可以報(bào)告最好解、最差解、平均解、解的標(biāo)準(zhǔn)差。這比只報(bào)告一次運(yùn)行的結(jié)果要嚴(yán)謹(jǐn)?shù)枚?。在論文中你可以說“我們獨(dú)立運(yùn)行算法30次最佳結(jié)果為XXX平均結(jié)果為YYY±ZZZ”這體現(xiàn)了方法的魯棒性。4. 可行性檢查是效率的關(guān)鍵如前文代碼所示在調(diào)用耗時(shí)的精確求解器或復(fù)雜仿真之前盡可能用簡(jiǎn)單、快速的條件過濾掉不可行的候選解。對(duì)于選址問題容量檢查就是這樣一個(gè)“守門員”。在其他問題中可能是變量的簡(jiǎn)單邊界檢查或者一些必須滿足的硬性邏輯約束。每過濾掉一個(gè)不可行解就節(jié)省了一次昂貴的評(píng)估。5. 與精確解或已知下界對(duì)比如果問題規(guī)模較小可以嘗試用商業(yè)求解器如Gurobi求出精確最優(yōu)解作為對(duì)比的“黃金標(biāo)準(zhǔn)”。如果求不出精確解可以嘗試計(jì)算一個(gè)問題的下界例如線性規(guī)劃松弛的解。將隨機(jī)搜索得到的最好解與精確解或下界進(jìn)行比較可以量化你的近似解的質(zhì)量。例如“我們的隨機(jī)搜索算法在5000次迭代內(nèi)找到的解與問題下界的差距在5%以內(nèi)”這是一個(gè)非常有說服力的結(jié)果。6. 可視化搜索過程將每次迭代找到的“當(dāng)前最優(yōu)成本”記錄下來并繪圖是分析算法行為的利器。你可以看到成本是如何隨著迭代下降的下降的速度如何何時(shí)趨于平穩(wěn)。這張圖放在論文的附錄或正文中能直觀地展示算法的收斂性。如果曲線下降很快然后平緩說明算法初期探索有效如果曲線一直緩慢下降說明可能需要更多迭代或改進(jìn)采樣策略。最后隨機(jī)搜索的魅力在于其簡(jiǎn)單性與強(qiáng)大通用性之間的平衡。它可能不是最優(yōu)雅、最快速的算法但在面對(duì)數(shù)學(xué)建模中那些“不講武德”的復(fù)雜現(xiàn)實(shí)問題時(shí)它往往是最忠實(shí)、最可靠的伙伴。掌握它意味著你在優(yōu)化工具箱里又多了一件應(yīng)對(duì)不確定性的利器。下次當(dāng)你面對(duì)一個(gè)看似無從下手的多變量?jī)?yōu)化模型時(shí)不妨先試試隨機(jī)搜索讓它為你照亮一片可能的解空間或許驚喜就在其中。