1. 從現(xiàn)場問題說起基站覆蓋率為什么不是一個“多建站”就能解決的問題我在做LTE網(wǎng)絡(luò)規(guī)劃優(yōu)化的時候最常聽到的一句話就是“覆蓋不行就加站”。這個說法聽起來簡單但實際做起來完全不是這么回事。某片區(qū)基站數(shù)量并不算少可是因為地物遮擋、站距不合理、天線參數(shù)失配往往會出現(xiàn)一種尷尬局面——路測數(shù)據(jù)里RSRP平均值挺好看但弱覆蓋柵格依然占了一成多投訴電話還是會打到網(wǎng)優(yōu)那邊去。這里面的本質(zhì)是一個組合優(yōu)化問題在一個給定的目標區(qū)域內(nèi)從備選站址集合中挑選出若干個基站使得整體覆蓋率最大。別小看這個“挑選”它的搜索空間極其恐怖。假設(shè)備選站址有40個需要選出12個那么組合數(shù)大約是C(40,12)的量級差不多是8.5億種方案。每一種方案都要做一遍覆蓋評估。想靠窮舉遍歷或者簡單貪心策略去找全局最優(yōu)基本是不現(xiàn)實的。更麻煩的是覆蓋率的計算本身就不是一個線性問題。信號傳播要經(jīng)過路徑損耗、穿透損耗、陰影衰落再加上高大建筑的遮擋一個位置能不能被覆蓋取決于它和所有基站之間的信號強度關(guān)系。也就是說你改變?nèi)魏我粋€基站的位置或者參數(shù)整個區(qū)域的覆蓋分布都會跟著變。這種目標函數(shù)往往存在大量局部極值傳統(tǒng)的最優(yōu)化算法很難穩(wěn)定發(fā)揮。所以我在做這個模擬項目X的時候把目光放到了群體智能算法上。粒子群優(yōu)化這類算法不依賴梯度信息適合這種非常規(guī)的多峰優(yōu)化問題。但標準粒子群又有一個老毛病——收斂快也容易過早陷入局部最優(yōu)。為了在收斂速度和全局探索之間取得平衡我在量子的思路基礎(chǔ)上引入了柯西分布擾動最終實現(xiàn)了一個“柯西分布量子粒子群優(yōu)化”CQPSO的Matlab求解方案。下面把這套東西的來龍去脈、數(shù)學(xué)模型和實現(xiàn)細節(jié)完整拆一遍希望能給正在做類似“選址優(yōu)化”“資源規(guī)劃”的同學(xué)一點啟發(fā)。2. 算法選型思路從經(jīng)典PSO到量子粒子群再到柯西分布加持2.1 標準PSO的困境跑得快也容易“撞墻”說起粒子群優(yōu)化很多做優(yōu)化的朋友都不陌生。它的想法來自鳥群覓食每個粒子代表一個候選解通過不斷向自身歷史最優(yōu)位置和群體全局最優(yōu)位置學(xué)習(xí)更新自己的“飛行”方向。標準PSO的速度更新公式如下v(i,:) w * v(i,:) c1 * rand(1,d) .* (pbest(i,:) - x(i,:)) c2 * rand(1,d) .* (gbest(1,:) - x(i,:)); x(i,:) x(i,:) v(i,:);這里w是慣性權(quán)重c1、c2是學(xué)習(xí)因子rand在[0,1]之間隨機取值。這版本最大的問題在于“速度盲區(qū)”。當粒子飛向gbest的時候如果速度沒有及時衰減粒子容易沖過頭在最優(yōu)解附近震蕩而不是精細收斂反過來如果w設(shè)置偏小種群又會快速趨向同一個位置多樣性迅速流失一旦gbest本身不是全局最優(yōu)整個群體就被“鎖死”在那里。我在實驗里經(jīng)?？吹綐藴蔖SO前20輪收斂氣勢很猛但后面80輪紋絲不動基本可以斷定是陷入局部極值了。2.2 QPSO核心機制量子勢阱和平均最優(yōu)位置量子粒子群優(yōu)化QPSO和標準PSO最本質(zhì)的區(qū)別是它徹底取消了“速度”這個概念粒子的位置更新基于量子力學(xué)中粒子在勢阱中的概率分布特性。QPSO的位置更新公式核心是兩個量個體最優(yōu)的加權(quán)平均位置mbest以及每個粒子當前最優(yōu)與全局最優(yōu)之間的距離。具體公式如下mbest sum(pbest, 1) / N; % N為種群規(guī)模 fi rand(1, d) * c; % c 1.49445源自QPSO文獻常用值 p fi .* pbest(i,:) (1 - fi) .* gbest(1,:); % 個體吸引子 u rand(1, d); if rand 0.5 x(i,:) p beta * abs(mbest - x(i,:)) .* log(1 ./ u); else x(i,:) p - beta * abs(mbest - x(i,:)) .* log(1 ./ u); end這里的beta收縮擴張系數(shù)非常關(guān)鍵。它控制著粒子搜索的“烈度”實驗中通常從1.0線性遞減到0.5這樣前期探索范圍大后期逐漸收斂到最優(yōu)解附近。QPSO的優(yōu)勢在于它能以概率1覆蓋整個解空間理論上全局收斂能力優(yōu)于標準PSO而且沒有速度參數(shù)需要調(diào)的參數(shù)更少。2.3 柯西分布擾動長尾巴帶來的“跳躍力”QPSO本身已經(jīng)比標準PSO穩(wěn)健不少但我在Test函數(shù)測試中發(fā)現(xiàn)一個現(xiàn)象當問題維度升高、局部極值增多時QPSO容易出現(xiàn)階段性“卡頓”也就是種群多樣性下降得還是太快。這時候引入柯西分布正是為了提高解空間的開拓能力。柯西分布和正態(tài)分布最大的區(qū)別在于它的尾部更厚也就是說雖然大多數(shù)采樣值集中在中心附近但出現(xiàn)離群值的概率明顯更高。這種“偶爾跳一大步”的特性用在實際算法里就是一種天然的變異機制——粒子有機會脫離當前吸引域跳到另一個較遠區(qū)域重新探索。我在實現(xiàn)時沒有直接替換QPSO的位置更新分布而是把柯西擾動疊加在每次迭代產(chǎn)生的坑位位置上。具體方式是生成一個服從標準柯西分布的隨機增量以一定概率比如0.2~0.3施加到部分粒子上cauchy_rnd tan(pi * (rand(1, d) - 0.5)); % 標準柯西分布采樣 x(i,:) x(i,:) cauchy_weight * cauchy_rnd;這個擾動幅度如果設(shè)置得當可以在不破壞收斂趨勢的前提下有效避免整個群體跌入同一個陷阱。簡單打個比方正常QPSO就像一隊人在陌生山區(qū)里沿著已知的小路找營地穩(wěn)但可能會走偏柯西擾動則是每隔一陣安排一個隊員向隨機方向踢一腳路邊的石頭偶爾就能發(fā)現(xiàn)一條沒標在地圖上的岔路。2.4 CQPSO整體流程和偽代碼結(jié)合起來的CQPSO流程如下初始化種群每個粒子的位置就是一組候選基站編號的編碼向量。計算每個粒子對應(yīng)的覆蓋率適應(yīng)度。更新個體最優(yōu)pbest和全局最優(yōu)gbest。計算平均最優(yōu)位置mbest。對每個粒子執(zhí)行QPSO基礎(chǔ)位置更新。按概率施加柯西擾動并做邊界修正。重新計算適應(yīng)度更新歷史最優(yōu)。判斷是否達到最大迭代次數(shù)否則回到第3步。這套流程實現(xiàn)起來并不復(fù)雜但每一環(huán)節(jié)的細節(jié)都會影響最終效果。尤其是柯西擾動的幅度和施加概率我在實驗中反復(fù)調(diào)了幾組參數(shù)后面會單獨說。3. 覆蓋率問題的數(shù)學(xué)模型與目標函數(shù)設(shè)計3.1 信號傳播模型用Okumura-Hata估算RSRP在LTE網(wǎng)絡(luò)規(guī)劃中覆蓋率的核心判據(jù)通常看參考信號接收功率RSRP。RSRP低于一定門限比如-110dBm就視為弱覆蓋這個柵格就不算被覆蓋。RSRP的計算需要先估算路徑損耗。在城市宏站場景下Okumura-Hata模型是一個經(jīng)典選擇。我在模擬項目X中使用的參數(shù)大致如下基站發(fā)射功率46dBm約40W天線增益15dBi不含饋線損耗載波頻率2GHz接收終端天線高度1.5m基站天線高度35m路徑損耗指數(shù)和修正項直接采用Hata城市模型路徑損耗公式2GHz城市環(huán)境可以寫成L 46.3 33.9 * log10(f) - 13.82 * log10(hb) - a(hm) (44.9 - 6.55 * log10(hb)) * log10(d);其中f單位是MHzhb是基站天線高度d單位是kma(hm)是終端高度修正因子。對于中等城市a(hm) (1.1 * log10(f) - 0.7) * hm - (1.56 * log10(f) - 0.8)。然后RSRP估算為RSRP txPower antennaGain - L - penetrationLoss - shadowFading;簡單場景可以不考慮陰影衰落直接采用距離相關(guān)路徑損耗如果要更貼近實際可以加一個服從對數(shù)正態(tài)分布的陰影衰落項。我做的模擬項目X里為了讓算法效果更直觀也為了減少隨機因素干擾選擇關(guān)閉陰影衰落只考慮距離和地物穿透損耗。3.2 目標區(qū)域離散化與覆蓋率計算連續(xù)區(qū)域的覆蓋率沒法直接計算通常做法是把目標區(qū)域劃分為均勻網(wǎng)格每個網(wǎng)格點代表一個采樣位置。判斷一個網(wǎng)格點是否被覆蓋就看它的RSRP是否大于預(yù)先設(shè)置的門限我取-110dBm。覆蓋率定義為coverage_rate covered_grid_count / total_grid_count;這個值就是算法的適應(yīng)度函數(shù)。目標就是讓它最大化。這里有一個容易被忽視的點柵格劃分粒度。網(wǎng)格太細計算量指數(shù)增長網(wǎng)格太粗又可能漏掉小面積的覆蓋空洞。我的經(jīng)驗是3km乘3km的區(qū)域劃到80乘80左右就已經(jīng)能看出比較細致的覆蓋覆蓋形態(tài)了再細化到100以上覆蓋率的數(shù)值變化很小但耗時明顯增加。3.3 約束條件與編碼方式基站選址問題不僅要考慮覆蓋率還要考慮實際工程約束。比如兩個基站不能靠得太近避免干擾預(yù)算限制下最多只能建設(shè)K個基站站址只能從候選站址集合中選取。在編碼上我采用了一個簡單直接的方案每個粒子的位置是一個K維向量每一維的取值范圍是候選站址的編號索引。如果候選站址離散排布那么編號就是1到N_candidate之間的整數(shù)。為了保證粒子的位置在迭代過程中始終有物理意義我在位置更新之后做了一個取整和邊界鉗制同時檢查任意兩個基站索引是否相同如果重復(fù)則做隨機擾動。這種編碼會損失一些連續(xù)優(yōu)化算法的優(yōu)勢因為位置是離散的但勝在直觀、容易實現(xiàn)而且CQPSO本身對解編碼形式不敏感實際優(yōu)化效果依然不錯。如果想要更高精度可以采用實數(shù)編碼再映射到最近的候選站址但這個方案我沒有在模擬項目X里采用原因后面講。4. Matlab實現(xiàn)中的關(guān)鍵細節(jié)與核心代碼走讀4.1 場景建模與候選站址生成先構(gòu)建一個20乘20平方公里的區(qū)域道路和建筑物分布簡化為不同權(quán)重的柵格地物衰減。我在模擬項目X里用的是3km乘3km區(qū)域地物分為三類市區(qū)建筑密集區(qū)域、園區(qū)開闊區(qū)域、郊野低矮區(qū)域分別對應(yīng)穿透損耗25dB、15dB、5dB。柵格為30m乘30m也就是100乘100個柵格點。候選站址不是隨機撒的而是根據(jù)現(xiàn)有網(wǎng)絡(luò)布局模擬出一批“可選位置”比如現(xiàn)有桿塔站點、樓頂站點、街邊燈桿站等。我生成了30個候選點最終希望從中選出10個站址。每個站址帶上位置坐標、高度、可用天線類型這些信息。所有場景數(shù)據(jù)存放在結(jié)構(gòu)體里scenario.areaLen 3000; % 區(qū)域邊長單位米 scenario.hotspotGrid zeros(100, 100); % 地物類型權(quán)重 scenario.candidateSite [x_coords, y_coords, height, antennaGain];4.2 覆蓋率計算函數(shù)全套RSRP計算和柵格判斷這是整個項目最核心的函數(shù)。輸入是當前選中的基站坐標和參數(shù)輸出是該方案下的覆蓋率。function coverage calcCoverage(sitePos, scenario) gridX scenario.gridX; % 柵格坐標矩陣 gridY scenario.gridY; thrRSRP -110; % 覆蓋門限單位dBm covered zeros(size(gridX)); N size(sitePos, 1); rxPowerSum -Inf * ones(size(gridX)); for i 1:N d sqrt((gridX - sitePos(i,1)).^2 (gridY - sitePos(i,2)).^2); d(d 1) 1; % 防止短路 L 46.3 33.9*log10(2000) - 13.82*log10(35) - 0.98 ... (44.9 - 6.55*log10(35)) * log10(d/1000); rsrp 46 15 - L - scenario.penetrationLoss; rxPowerSum max(rxPowerSum, rsrp); % 最大RSRP原則 end covered rxPowerSum thrRSRP; coverage sum(covered(:)) / numel(covered); end這里有一個工程細節(jié)一個柵格點可能同時被多個基站覆蓋實際終端會選擇接收功率最強的那個小區(qū)駐留所以計算時用max合并多個基站的RSRP而不是做簡單疊加。這個細節(jié)曾經(jīng)讓我吃過虧一開始我用的是所有基站RSRP的線性平均值結(jié)果覆蓋率虛高后面發(fā)現(xiàn)和實際不符改成取最強信號之后才靠譜。4.3 基于柯西分布的變異擾動實現(xiàn)這一小節(jié)是整個算法的“靈魂”。在設(shè)計CQPSO時柯西擾動我打算加在QPSO位置更新之后但并不是每個粒子都加。具體實現(xiàn)里我采用的是“按粒子概率的擾動策略”% QPSO更新 x_new ... 基礎(chǔ)QPSO更新結(jié)果... % 柯西擾動 if rand cauchyProb scale cauchyWeight * (1 - iter / maxIter) 0.01; x_new x_new scale * tan(pi * (rand(1,dim) - 0.5)); end % 邊界處理 x_new max(x_new, lb); x_new min(x_new, ub);擾動權(quán)重為什么要隨時間遞減原因很簡單算法前期需要強探索能力需要長尾大跳躍來探索陌生區(qū)域后期接近收斂如果擾動依然很大粒子會在最優(yōu)解附近來回橫跳導(dǎo)致收斂精度變差。我試驗下來cauchyProb取0.25、初始cauchyWeight取0.2倍搜索范圍比較均衡。4.4 主循環(huán)與可視化輸出主循環(huán)邏輯很直觀。在每次迭代結(jié)束后我會記錄當前最優(yōu)覆蓋率并且每5代重新繪制一次覆蓋熱力圖用不同顏色標識覆蓋/未覆蓋區(qū)域。Matlab的imagesc配合自定義colormap比直接用plot散點圖直觀得多。for iter 1:maxIter beta 1.0 - 0.5 * (iter / maxIter); % 計算mbest、更新所有粒子 % ... % 柯西擾動 % ... % 適應(yīng)度評估 for i 1:popSize fitness(i) calcCoverage(decode(x(i,:)), scenario); end % 更新pbest和gbest % 記錄convergence(iter) if mod(iter, 5) 0 plotCoverage(gbest, scenario); end end這里需要額外提醒一點calcCoverage函數(shù)是整個Simulation里計算量最大的瓶頸。群體規(guī)模30、迭代200次就是6000次覆蓋率計算每次覆蓋率計算要遍歷10000個柵格點。所以函數(shù)里的矢量化操作必須做充分盡量避免逐個柵格的for循環(huán)否則Matlab跑起來會非常痛苦后面我會講具體的加速技巧。5. 實驗設(shè)計、參數(shù)調(diào)優(yōu)與結(jié)果對比分析5.1 實驗場景與參數(shù)配置我的實驗分三組算法對比標準PSO、QPSO、CQPSO共同使用相同的場景、相同的初始種群種子、相同的迭代次數(shù)以保證對比公平性。主要參數(shù)如下參數(shù)PSOQPSOCQPSO種群規(guī)模303030最大迭代次數(shù)150150150慣性權(quán)重w0.6~0.9--學(xué)習(xí)因子c1/c22.0/2.0--收縮擴張系數(shù)beta-1.0→0.51.0→0.5柯西擾動概率--0.25初始擾動尺度--0.2候選站址30個目標建站數(shù)10個。覆蓋率初始隨機方案平均在72%左右。5.2 覆蓋率與收斂行為對比三組算法各跑10次記錄最好覆蓋率、平均覆蓋率、收斂代數(shù)。從我的實驗結(jié)果看標準PSO最好覆蓋率約為86.7%平均85.9%但多次實驗的方差較大有的實驗收斂到84%就完全停滯。QPSO最好覆蓋率約為89.2%平均88.5%整體明顯優(yōu)于PSO收斂曲線平滑后期仍有緩慢爬升。CQPSO最好覆蓋率約為91.4%平均90.8%主要優(yōu)勢在于多次實驗的最差結(jié)果依然能到89%以上說明它更穩(wěn)定而且后期被“鎖死”的概率大幅降低。從收斂曲線看CQPSO初期收斂速度和QPSO接近但在迭代60代之后QPSO往往進入平臺期而CQPSO還能通過柯西擾動跳出來找到新的更優(yōu)位置。那批“跳變點”在曲線上表現(xiàn)為一個向上臺階隔幾代出現(xiàn)一次。這正是重尾分布發(fā)揮作用的表現(xiàn)。5.3 最優(yōu)站址方案的覆蓋熱力圖解讀我把最優(yōu)方案畫成熱力圖后發(fā)現(xiàn)幾個有意思的現(xiàn)象算法自動把站址往“地物衰減大”的區(qū)域適當加密而寬敞開闊的區(qū)域站間距拉得更大。區(qū)域邊緣地帶出現(xiàn)明顯的“塔下黑”緩解因為基站并非均勻分布而是和信號傳播損耗做了匹配。存在冗余覆蓋的情況比較少選出的10個站址之間距離都大于最小間隔約束沒有出現(xiàn)兩個站幾乎重疊的浪費情況。通過對比初始隨機方案和CQPSO最優(yōu)方案可以看到未覆蓋柵格從28%下降到8.6%并且這些未覆蓋柵格基本集中在場景角落和極少數(shù)被高層阻擋的區(qū)域這是純路徑損耗模型下能接受的結(jié)果。6. 實戰(zhàn)中的坑與填坑實錄6.1 柵格分辨率帶來的“假收斂”我在實驗早期把柵格設(shè)得非常粗只有20乘20算出來的覆蓋率輕松到了95%。后來改到80乘80同樣的方案掉到88%。原因很簡單粗網(wǎng)格天然忽略了小面積的弱覆蓋縫隙。所以對比算法效果時一定要固定柵格粒度不能今天用20乘20明天用80乘80。建議最終評估統(tǒng)一用100乘100的細網(wǎng)格寧可多花計算時間也要讓結(jié)果經(jīng)得起推敲。6.2 柯西擾動“尺度太大會炸”我在調(diào)參階段踩過一個很典型的坑初始cauchyWeight設(shè)成0.8倍搜索范圍結(jié)果粒子整天往邊界外面飛修邊界修到手軟覆蓋率反而比不加擾動還差。后來分析發(fā)現(xiàn)柯西分布采樣值偶爾會出現(xiàn)十幾倍于標準差的極端值如果不限制擾動幅度的上限粒子位置會嚴重脫離可行區(qū)域邊界截斷又會讓大量粒子堆積在邊界上損害種群多樣性。最終我把擾動幅度分成兩項第一是整體縮放因子隨時間遞減第二是限制單次擾動的絕對幅度不能超過解空間范圍的0.3倍。按這個邏輯重新調(diào)參后效果立刻好轉(zhuǎn)。6.3 重復(fù)站址沖突的處理離散編碼的問題是位置更新后粒子向量的各維可能收斂到同一個候選站址這樣實際建站數(shù)會少于預(yù)設(shè)值覆蓋率虛高。我一開始沒注意后來檢查輸出才發(fā)現(xiàn)的。解決辦法很簡單在評估覆蓋率之前做一個去重和補缺把重復(fù)的維度替換為當前概率最低的未使用候選站或者隨機替換一個未使用站。6.4 Matlab執(zhí)行效率大幅優(yōu)化覆蓋率函數(shù)里三重for循環(huán)是最容易拖垮Matlab的寫法。我初始版本用了嵌套循環(huán)跑一次迭代花了40秒整個人都裂開了。后來做了幾個改造所有柵格點和基站的距離計算全部矢量化利用隱式擴展處理。RSRP計算不在循環(huán)里調(diào)用log10等耗時函數(shù)而是把所有常數(shù)項預(yù)先合并成緊湊公式。利用profiler定位熱點原來是max函數(shù)在混合維度下反復(fù)創(chuàng)建臨時矩陣通過預(yù)先分配輸出矩陣單次覆蓋率計算時間壓縮到0.08秒。如果你拿到的工程規(guī)模很大還可以考慮把覆蓋率函數(shù)打包成MEX或使用并行池不過我的場景經(jīng)過矢量化之后已經(jīng)夠用就沒有再往復(fù)雜里改。6.5 隨機種子的控制如果你想嚴謹對比算法一定要在實驗?zāi)_本開頭固定隨機種子。否則一組實驗結(jié)果可能是運氣好、另一組運氣差完全失去可比性。我在測試腳本里使用rng(42, twister);每條獨立實驗前重置一次并記錄隨機種子到結(jié)果文件名里方便復(fù)查。7. 最后再分享一點個人心得做這類“智能優(yōu)化算法工程問題”的項目最容易犯的一個誤區(qū)是過度迷信算法把整篇文章的重心全放在了優(yōu)化算法本身反而把覆蓋率計算、場景建模這些基礎(chǔ)工程做得很粗糙。我實際做下來的體會是覆蓋率模型的可信度決定了整個優(yōu)化結(jié)果的上限。算法再強如果信號傳播模型脫離實際或者場景柵格劃分根本不反映真實地物那跑出來的“最優(yōu)解”只是數(shù)學(xué)上的最優(yōu)放到現(xiàn)實中一堆站址根本沒法評判。所以在動手調(diào)算法之前建議先把場景建模和覆蓋率評估函數(shù)做到自己敢拿出去復(fù)盤的程度。另一方面柯西分布引入QPSO這件事本身并不復(fù)雜但參數(shù)對結(jié)果的影響比想象中大。如果你準備在自己項目中借鑒建議按“先固定隨機種子→固定場景→跑基準QPSO→逐組調(diào)柯西概率→逐組調(diào)擾動尺度”這個順序做實驗每一步只動一個變量才不會陷入?yún)?shù)調(diào)優(yōu)的大坑。這個模擬項目X的框架后續(xù)還可以擴展不少方向比如把覆蓋率判據(jù)從單一RSRP門限改成RSRP和SINR聯(lián)合判斷或者在適應(yīng)度函數(shù)中加入建站成本約束、負載均衡因子等等。Matlab代碼整體上保持著清晰的模塊邊界替換傳播模型和評估函數(shù)都很方便拿來當算法研究和網(wǎng)絡(luò)規(guī)劃的測試骨架應(yīng)該能省下不少重復(fù)造輪子的時間。