童問(wèn)題：庫(kù)存決策優(yōu)化與不確定性建模實(shí)踐)
1. 項(xiàng)目概述從報(bào)童到庫(kù)存決策的經(jīng)典模型報(bào)童問(wèn)題這個(gè)名字聽(tīng)起來(lái)有點(diǎn)懷舊但它絕不是只存在于歷史課本里的故事。我第一次接觸這個(gè)模型是在研究生階段的一門(mén)運(yùn)籌學(xué)課上當(dāng)時(shí)覺(jué)得它不過(guò)是個(gè)簡(jiǎn)單的概率計(jì)算練習(xí)。直到后來(lái)在電商公司的供應(yīng)鏈部門(mén)實(shí)習(xí)親眼看到每天凌晨算法是如何決定向各個(gè)倉(cāng)庫(kù)補(bǔ)多少貨而第二天又有多少商品因?yàn)槿必浕驕N被標(biāo)記處理時(shí)我才恍然大悟——那個(gè)“賣(mài)報(bào)紙的小孩”面對(duì)的困境正是現(xiàn)代商業(yè)庫(kù)存管理的核心縮影。簡(jiǎn)單來(lái)說(shuō)報(bào)童問(wèn)題描述的是這樣一個(gè)場(chǎng)景一個(gè)報(bào)童每天早晨需要決定從報(bào)社批發(fā)多少份報(bào)紙來(lái)賣(mài)。報(bào)紙的需求量是隨機(jī)的他只知道一個(gè)大概的概率分布。如果批發(fā)多了賣(mài)不完的報(bào)紙到了晚上就一文不值會(huì)造成損失如果批發(fā)少了沒(méi)買(mǎi)到的顧客就走了會(huì)損失潛在的利潤(rùn)。他的目標(biāo)就是找到一個(gè)最優(yōu)的訂購(gòu)量讓他的長(zhǎng)期平均利潤(rùn)最大化或者說(shuō)期望損失最小化。這個(gè)模型的核心就是在不確定性的環(huán)境下做單周期的庫(kù)存決策。今天我們不用真的去賣(mài)報(bào)紙而是用 MATLAB 這個(gè)強(qiáng)大的工具來(lái)親手搭建一個(gè)報(bào)童問(wèn)題的仿真環(huán)境。仿真的意義在于它允許我們?cè)谟?jì)算機(jī)里創(chuàng)造一個(gè)“虛擬世界”在這個(gè)世界里我們可以設(shè)定不同的需求分布、成本參數(shù)然后讓“報(bào)童”按照我們?cè)O(shè)定的策略去運(yùn)營(yíng)成千上萬(wàn)天快速、低成本地觀察不同決策帶來(lái)的長(zhǎng)期結(jié)果。這對(duì)于驗(yàn)證理論公式、比較不同補(bǔ)貨策略、或者處理那些理論模型難以解決的復(fù)雜情況比如需求分布未知、存在缺貨懲罰等來(lái)說(shuō)是極其有效的方法。無(wú)論你是學(xué)習(xí)運(yùn)籌學(xué)、供應(yīng)鏈管理的學(xué)生還是對(duì)數(shù)據(jù)分析和決策優(yōu)化感興趣的從業(yè)者這個(gè)仿真項(xiàng)目都能幫你直觀地理解不確定性決策的精髓。2. 問(wèn)題拆解與數(shù)學(xué)模型建立在動(dòng)手寫(xiě)代碼之前我們必須先把問(wèn)題用數(shù)學(xué)語(yǔ)言清晰地定義出來(lái)。這是所有仿真和分析的基石含糊不得。2.1 核心參數(shù)與變量定義首先我們需要明確幾個(gè)關(guān)鍵的經(jīng)濟(jì)參數(shù)這些是驅(qū)動(dòng)整個(gè)模型的“輸入”單位成本 (c)報(bào)童從報(bào)社批發(fā)一份報(bào)紙需要支付的價(jià)格。這是他的成本。單位售價(jià) (p)報(bào)童將一份報(bào)紙賣(mài)給顧客的價(jià)格。這是他的收入來(lái)源。單位殘值 (s)當(dāng)天結(jié)束時(shí)一份沒(méi)有賣(mài)出去的報(bào)紙的剩余價(jià)值。通常s c可能為零完全報(bào)廢也可能是個(gè)很小的正數(shù)回收價(jià)。缺貨懲罰 (g)這是一個(gè)可選但很實(shí)際的參數(shù)。它表示當(dāng)顧客需要報(bào)紙而報(bào)童缺貨時(shí)所造成的額外損失。這不僅包括失去本次銷售的利潤(rùn) (p - c)還可能包括商譽(yù)損失、顧客流失等隱性成本。在基礎(chǔ)模型中常設(shè)為0但加上它會(huì)讓模型更貼近現(xiàn)實(shí)。接下來(lái)是決策變量和隨機(jī)變量訂購(gòu)量 (Q)這是報(bào)童需要做出的決策也就是我們通過(guò)仿真要尋找的最優(yōu)解。它是一個(gè)非負(fù)整數(shù)。需求量 (D)這是一個(gè)隨機(jī)變量。我們假設(shè)它服從某種已知的概率分布比如正態(tài)分布、泊松分布或均勻分布。仿真的核心之一就是生成符合這個(gè)分布的隨機(jī)需求序列。最后是基于以上變量計(jì)算出的結(jié)果實(shí)際銷量 (Sales)這取決于訂購(gòu)量和需求量中較小的那個(gè)即Sales min(Q, D)。你只能賣(mài)掉你有的和顧客需要的兩者中較少的那部分。剩余庫(kù)存 (Leftover)當(dāng)天結(jié)束時(shí)沒(méi)賣(mài)出去的報(bào)紙即Leftover max(0, Q - D)。缺貨量 (Shortage)當(dāng)天未能滿足的顧客需求即Shortage max(0, D - Q)。2.2 利潤(rùn)函數(shù)與期望利潤(rùn)最大化有了這些定義一天的利潤(rùn)Π(Q, D)就可以寫(xiě)出來(lái)了Π(Q, D) p * min(Q, D) s * max(0, Q - D) - c * Q - g * max(0, D - Q)這個(gè)公式拆開(kāi)看很直觀p * min(Q, D)銷售收入。s * max(0, Q - D)剩余庫(kù)存的殘值回收收入。c * Q批發(fā)報(bào)紙的總成本。g * max(0, D - Q)缺貨造成的懲罰成本。由于需求量D是隨機(jī)的單日的利潤(rùn)也是隨機(jī)的。因此報(bào)童關(guān)心的是長(zhǎng)期平均利潤(rùn)也就是利潤(rùn)的期望值E[Π(Q)]。我們的優(yōu)化目標(biāo)是找到一個(gè)最優(yōu)訂購(gòu)量Q*使得期望利潤(rùn)最大化Q* argmax_{Q≥0} E[Π(Q)]在理論上對(duì)于某些特定的分布如正態(tài)分布存在一個(gè)著名的臨界分位數(shù) (Critical Fractile) 公式來(lái)求解Q*F(Q*) (p - c g) / (p - s g)其中F(·)是需求量D的累積分布函數(shù) (CDF)。這個(gè)公式的意義在于最優(yōu)庫(kù)存水平應(yīng)該設(shè)置在這樣一個(gè)位置需求不超過(guò)該水平的概率恰好等于“單位欠儲(chǔ)成本”與“單位欠儲(chǔ)成本加單位超儲(chǔ)成本”之比。這里(p - c g)可以理解為少進(jìn)一份報(bào)紙?jiān)斐傻倪呺H損失即欠儲(chǔ)成本(p - s g)可以理解為決策的總體邊際影響。注意這個(gè)理論解非常優(yōu)美但它依賴于我們知道準(zhǔn)確的需求分布F(·)。在現(xiàn)實(shí)中分布可能未知、可能隨時(shí)間變化、或者問(wèn)題本身更復(fù)雜如多產(chǎn)品、多周期。這時(shí)仿真 Monte Carlo Simulation 的價(jià)值就凸顯出來(lái)了——我們不需要知道F(·)的解析形式只需要能根據(jù)歷史數(shù)據(jù)或假設(shè)生成隨機(jī)需求樣本就能通過(guò)模擬來(lái)評(píng)估任何給定Q的性能甚至用搜索算法來(lái)尋找近似的Q*。3. MATLAB仿真環(huán)境搭建與核心代碼解析理論鋪墊完畢現(xiàn)在進(jìn)入實(shí)戰(zhàn)環(huán)節(jié)。我們將用 MATLAB 一步步構(gòu)建這個(gè)仿真系統(tǒng)。我個(gè)人的習(xí)慣是先搭建一個(gè)清晰、模塊化的框架這樣調(diào)試和擴(kuò)展都會(huì)很方便。3.1 參數(shù)初始化與需求數(shù)據(jù)生成首先我們創(chuàng)建一個(gè)腳本文件比如叫newsvendor_simulation.m。開(kāi)頭先定義所有基礎(chǔ)參數(shù)。%% 1. 參數(shù)設(shè)置 clear; clc; close all; % 清空環(huán)境好習(xí)慣 % 經(jīng)濟(jì)參數(shù) unit_cost 2; % c: 每份報(bào)紙批發(fā)成本元 unit_price 5; % p: 每份報(bào)紙零售價(jià)格元 unit_salvage 0.5; % s: 每份未售出報(bào)紙的殘值元 penalty_cost 1; % g: 每份缺貨的懲罰成本元可選設(shè)為0則為經(jīng)典模型 % 需求分布參數(shù) - 這里假設(shè)需求服從正態(tài)分布 demand_mean 100; % 平均日需求 demand_std 20; % 日需求標(biāo)準(zhǔn)差 % 仿真參數(shù) num_days 10000; % 模擬的天數(shù)天數(shù)越多結(jié)果越穩(wěn)定 order_quantity 90; % Q: 我們要測(cè)試的訂購(gòu)量可以先設(shè)一個(gè)值跑跑看接下來(lái)是生成隨機(jī)需求。MATLAB 的統(tǒng)計(jì)工具箱提供了豐富的隨機(jī)數(shù)生成器。%% 2. 生成隨機(jī)需求序列 % 使用正態(tài)分布生成需求。注意需求應(yīng)為非負(fù)整數(shù)所以需要取整和取最大值。 daily_demand max(round(normrnd(demand_mean, demand_std, num_days, 1)), 0); % normrnd生成正態(tài)分布隨機(jī)數(shù)round四舍五入取整max(...,0)確保非負(fù)。 % 可視化一下需求分布可選但強(qiáng)烈推薦 figure; subplot(2,1,1); histogram(daily_demand, Normalization, probability); xlabel(日需求量); ylabel(頻率); title(模擬日需求分布直方圖); grid on; subplot(2,1,2); cdfplot(daily_demand); % 繪制經(jīng)驗(yàn)累積分布函數(shù) xlabel(日需求量); ylabel(F(x)); title(需求的經(jīng)驗(yàn)CDF); grid on;實(shí)操心得生成需求時(shí)round和max(...,0)這兩個(gè)處理很重要。現(xiàn)實(shí)中需求是整數(shù)四舍五入更合理。雖然正態(tài)分布理論上可能產(chǎn)生負(fù)數(shù)但我們的demand_mean100,demand_std20產(chǎn)生負(fù)數(shù)的概率極低max(...,0)是一個(gè)安全的保護(hù)措施。如果你模擬的需求均值很小比如接近0則需要考慮使用嚴(yán)格非負(fù)的分布如泊松分布poissrnd(lambda, num_days, 1)。3.2 單周期利潤(rùn)計(jì)算與仿真循環(huán)核心的計(jì)算邏輯封裝成一個(gè)函數(shù)會(huì)非常清晰。我們先寫(xiě)一個(gè)計(jì)算單日利潤(rùn)的函數(shù)。function profit calculate_daily_profit(Q, D, p, c, s, g) % 計(jì)算報(bào)童模型單日利潤(rùn) % 輸入: Q - 訂購(gòu)量, D - 當(dāng)日實(shí)際需求, p,c,s,g - 經(jīng)濟(jì)參數(shù) % 輸出: profit - 當(dāng)日利潤(rùn) sales min(Q, D); % 實(shí)際銷量 leftover max(0, Q - D); % 剩余庫(kù)存 shortage max(0, D - Q); % 缺貨量 revenue p * sales; % 銷售收入 salvage_income s * leftover; % 殘值收入 procurement_cost c * Q; % 采購(gòu)成本 shortage_penalty g * shortage; % 缺貨懲罰 profit revenue salvage_income - procurement_cost - shortage_penalty; end然后在主腳本中我們進(jìn)行仿真循環(huán)計(jì)算長(zhǎng)期平均利潤(rùn)。%% 3. 仿真計(jì)算 daily_profits zeros(num_days, 1); % 預(yù)分配數(shù)組提升效率 for day 1:num_days current_demand daily_demand(day); daily_profits(day) calculate_daily_profit(order_quantity, ... current_demand, ... unit_price, ... unit_cost, ... unit_salvage, ... penalty_cost); end % 計(jì)算關(guān)鍵績(jī)效指標(biāo) (KPIs) average_daily_profit mean(daily_profits); profit_std std(daily_profits); service_level sum(daily_demand order_quantity) / num_days; % 需求滿足率庫(kù)存覆蓋概率 fprintf(仿真結(jié)果訂購(gòu)量 Q%d\n, order_quantity); fprintf( 平均日利潤(rùn): %.2f 元\n, average_daily_profit); fprintf( 利潤(rùn)標(biāo)準(zhǔn)差: %.2f 元\n, profit_std); % 衡量風(fēng)險(xiǎn) fprintf( 服務(wù)水平需求滿足率: %.2f%%\n, service_level * 100);3.3 結(jié)果可視化與分析數(shù)字有了但圖表更能說(shuō)明問(wèn)題。我們來(lái)繪制利潤(rùn)的分布和收斂情況。%% 4. 結(jié)果可視化 figure; % 子圖1日利潤(rùn)分布 subplot(2,2,1); histogram(daily_profits, 50, FaceColor, [0.2 0.6 0.8]); xlabel(日利潤(rùn)元); ylabel(頻數(shù)); title(sprintf(日利潤(rùn)分布 (Q%d), order_quantity)); grid on; hold on; % 標(biāo)記平均利潤(rùn)線 yl ylim; plot([average_daily_profit, average_daily_profit], [yl(1), yl(2)], r--, LineWidth, 2); legend(利潤(rùn)分布, 平均利潤(rùn), Location, best); hold off; % 子圖2累積平均利潤(rùn)看仿真收斂性 subplot(2,2,2); cumulative_avg_profit cumsum(daily_profits) ./ (1:num_days); plot(1:num_days, cumulative_avg_profit, b-, LineWidth, 1.5); xlabel(模擬天數(shù)); ylabel(累積平均利潤(rùn)元); title(平均利潤(rùn)隨仿真天數(shù)的收斂過(guò)程); grid on; hold on; plot([1, num_days], [average_daily_profit, average_daily_profit], r--); legend(累積平均, 最終平均, Location, southeast); hold off; % 子圖3利潤(rùn)與需求的關(guān)系散點(diǎn)圖 subplot(2,2,3); scatter(daily_demand, daily_profits, 10, filled, MarkerFaceAlpha, 0.6); xlabel(日需求量); ylabel(日利潤(rùn)); title(需求與利潤(rùn)關(guān)系散點(diǎn)圖); grid on; % 可以添加趨勢(shì)線或分界線 hold on; plot([order_quantity, order_quantity], ylim, k--, LineWidth, 1.5); % 標(biāo)記訂購(gòu)量 hold off; % 子圖4不同需求下的利潤(rùn)構(gòu)成示例取一天 subplot(2,2,4); sample_day find(daily_demand round(demand_mean), 1); % 找一個(gè)需求接近均值的天 if isempty(sample_day) sample_day 1; end sample_demand daily_demand(sample_day); [sales, leftover, shortage] deal(min(order_quantity, sample_demand), ... max(0, order_quantity - sample_demand), ... max(0, sample_demand - order_quantity)); profit_breakdown [unit_price*sales, unit_salvage*leftover, -unit_cost*order_quantity, -penalty_cost*shortage]; labels {銷售收入, 殘值收入, 采購(gòu)成本, 缺貨懲罰}; bar(profit_breakdown); set(gca, XTickLabel, labels); ylabel(金額元); title(sprintf(第%d天利潤(rùn)構(gòu)成 (需求%d), sample_day, sample_demand)); grid on;運(yùn)行這段代碼你就能得到一個(gè)完整的單點(diǎn)仿真結(jié)果。但我們的目標(biāo)是找到最優(yōu)的Q*所以下一步是進(jìn)行敏感性分析。4. 尋找最優(yōu)訂購(gòu)量仿真與理論對(duì)比現(xiàn)在我們讓Q動(dòng)起來(lái)觀察平均利潤(rùn)如何隨Q變化并嘗試找到那個(gè)最高點(diǎn)。4.1 遍歷搜索與利潤(rùn)曲線繪制我們?cè)O(shè)定一個(gè)Q的搜索范圍比如從demand_mean - 3*demand_std到demand_mean 3*demand_std覆蓋需求的絕大部分可能區(qū)間。%% 5. 尋找最優(yōu)訂購(gòu)量 Q* % 定義搜索范圍 Q_range floor(demand_mean - 3*demand_std) : ceil(demand_mean 3*demand_std); Q_range Q_range(Q_range 0); % 確保非負(fù) num_Q length(Q_range); avg_profit_list zeros(num_Q, 1); service_level_list zeros(num_Q, 1); fprintf(開(kāi)始掃描 %d 個(gè)不同的Q值...\n, num_Q); % 對(duì)每個(gè)Q進(jìn)行仿真。注意這里為了速度復(fù)用之前生成的需求序列。 % 如果追求絕對(duì)準(zhǔn)確應(yīng)對(duì)每個(gè)Q重新生成獨(dú)立的需求序列但計(jì)算量會(huì)大很多。 % 在Q值掃描中使用同一組需求序列是標(biāo)準(zhǔn)做法保證了比較的公平性。 for i 1:num_Q current_Q Q_range(i); temp_profits zeros(num_days, 1); for day 1:num_days temp_profits(day) calculate_daily_profit(current_Q, daily_demand(day), ... unit_price, unit_cost, ... unit_salvage, penalty_cost); end avg_profit_list(i) mean(temp_profits); service_level_list(i) sum(daily_demand current_Q) / num_days; end % 找到仿真下的最優(yōu)Q [sim_max_profit, sim_opt_idx] max(avg_profit_list); sim_opt_Q Q_range(sim_opt_idx); fprintf(【仿真結(jié)果】最優(yōu)訂購(gòu)量 Q* %d對(duì)應(yīng)平均日利潤(rùn) %.2f 元服務(wù)水平 %.2f%%\n, ... sim_opt_Q, sim_max_profit, service_level_list(sim_opt_idx)*100);繪制利潤(rùn)-訂購(gòu)量曲線。% 可視化利潤(rùn)曲線 figure; yyaxis left; plot(Q_range, avg_profit_list, b-o, LineWidth, 1.5, MarkerSize, 4); hold on; plot(sim_opt_Q, sim_max_profit, r*, MarkerSize, 15, LineWidth, 2); xlabel(訂購(gòu)量 Q); ylabel(平均日利潤(rùn)元); yyaxis right; plot(Q_range, service_level_list*100, g--s, LineWidth, 1.5, MarkerSize, 4); ylabel(服務(wù)水平 (%)); title(平均利潤(rùn)與服務(wù)水平隨訂購(gòu)量變化曲線); grid on; legend(平均利潤(rùn), sprintf(最優(yōu)點(diǎn) (Q%d), sim_opt_Q), 服務(wù)水平, ... Location, best);你會(huì)看到一條經(jīng)典的凹曲線利潤(rùn)先隨Q增加而上升因?yàn)槟茏プ「噤N售機(jī)會(huì)達(dá)到一個(gè)頂峰后開(kāi)始下降因?yàn)闇N損失開(kāi)始超過(guò)新增銷售的收益。那個(gè)頂峰對(duì)應(yīng)的Q就是我們的仿真最優(yōu)解。4.2 理論解計(jì)算與對(duì)比現(xiàn)在我們用前面提到的臨界分位數(shù)公式來(lái)計(jì)算理論最優(yōu)解并與仿真結(jié)果對(duì)比。%% 6. 理論解計(jì)算與對(duì)比 % 計(jì)算臨界分位數(shù) critical_ratio (unit_price - unit_cost penalty_cost) / ... (unit_price - unit_salvage penalty_cost); fprintf(臨界分位數(shù) (p - c g) / (p - s g) %.4f\n, critical_ratio); % 由于我們假設(shè)需求服從正態(tài)分布 N(mu, sigma^2) % 理論最優(yōu)Q*是滿足 F(Q*) critical_ratio 的值即逆CDF % 使用 norminv 函數(shù) theory_opt_Q norminv(critical_ratio, demand_mean, demand_std); theory_opt_Q round(theory_opt_Q); % 取整因?yàn)镼是整數(shù) fprintf(【理論解】最優(yōu)訂購(gòu)量 Q*_theory %.2f (取整后為 %d)\n, ... norminv(critical_ratio, demand_mean, demand_std), theory_opt_Q); % 計(jì)算理論解對(duì)應(yīng)的仿真利潤(rùn)用同一組需求數(shù)據(jù)評(píng)估 theory_profits zeros(num_days, 1); for day 1:num_days theory_profits(day) calculate_daily_profit(theory_opt_Q, daily_demand(day), ... unit_price, unit_cost, ... unit_salvage, penalty_cost); end theory_avg_profit mean(theory_profits); theory_service_level sum(daily_demand theory_opt_Q) / num_days; fprintf(理論解Q%d對(duì)應(yīng)的仿真評(píng)估平均利潤(rùn)%.2f服務(wù)水平%.2f%%\n, ... theory_opt_Q, theory_avg_profit, theory_service_level*100); % 對(duì)比分析 comparison_table table([sim_opt_Q; theory_opt_Q], ... [sim_max_profit; theory_avg_profit], ... [service_level_list(sim_opt_idx); theory_service_level]*100, ... VariableNames, {最優(yōu)訂購(gòu)量Q, 平均日利潤(rùn), 服務(wù)水平_百分比}, ... RowNames, {仿真搜索, 理論公式}); disp(comparison_table);正常情況下仿真搜索得到的Q*和理論公式計(jì)算的Q*應(yīng)該非常接近。如果差異較大可能的原因有1) 仿真天數(shù)num_days不夠多結(jié)果有波動(dòng)2) 需求分布不是完美的正態(tài)分布因?yàn)槲覀冏隽巳≌腿》秦?fù)處理3) 搜索的步長(zhǎng)不夠精細(xì)。增加num_days和縮小Q_range的步長(zhǎng)例如以1為步進(jìn)可以改善。注意事項(xiàng)norminv函數(shù)要求critical_ratio在 (0,1) 開(kāi)區(qū)間內(nèi)。如果您的成本參數(shù)設(shè)置導(dǎo)致critical_ratio非常接近0或1例如售價(jià)遠(yuǎn)低于成本norminv可能會(huì)返回-Inf或Inf。在實(shí)際業(yè)務(wù)中這通常意味著最優(yōu)策略是“不訂購(gòu)”或“訂購(gòu)極大數(shù)量”需要在實(shí)際代碼中加入邊界判斷。5. 深入分析與擴(kuò)展應(yīng)用場(chǎng)景基礎(chǔ)仿真跑通后我們可以玩點(diǎn)更花的讓模型更貼近復(fù)雜的現(xiàn)實(shí)情況。5.1 敏感性分析參數(shù)如何影響決策最優(yōu)訂購(gòu)量Q*對(duì)成本參數(shù)非常敏感。我們可以系統(tǒng)地改變一個(gè)參數(shù)比如單位成本c觀察Q*和最大利潤(rùn)的變化。%% 7. 敏感性分析示例單位成本c的影響 cost_range 1.5:0.1:2.5; % 單位成本從1.5元到2.5元變化 num_costs length(cost_range); opt_Q_vs_cost zeros(num_costs, 1); max_profit_vs_cost zeros(num_costs, 1); % 固定其他參數(shù)和需求序列 for i 1:num_costs current_cost cost_range(i); % 計(jì)算當(dāng)前成本下的臨界分位數(shù)和理論Q* current_cr (unit_price - current_cost penalty_cost) / ... (unit_price - unit_salvage penalty_cost); % 防止cr超出(0,1)范圍 current_cr max(min(current_cr, 0.999), 0.001); current_opt_Q round(norminv(current_cr, demand_mean, demand_std)); opt_Q_vs_cost(i) current_opt_Q; % 評(píng)估該Q下的仿真利潤(rùn) temp_profits zeros(num_days, 1); for day 1:num_days temp_profits(day) calculate_daily_profit(current_opt_Q, daily_demand(day), ... unit_price, current_cost, ... unit_salvage, penalty_cost); end max_profit_vs_cost(i) mean(temp_profits); end figure; subplot(2,1,1); plot(cost_range, opt_Q_vs_cost, b-o, LineWidth, 1.5); xlabel(單位成本 c (元)); ylabel(最優(yōu)訂購(gòu)量 Q*); title(最優(yōu)訂購(gòu)量隨單位成本變化); grid on; subplot(2,1,2); plot(cost_range, max_profit_vs_cost, r-s, LineWidth, 1.5); xlabel(單位成本 c (元)); ylabel(最大期望利潤(rùn) (元)); title(最大期望利潤(rùn)隨單位成本變化); grid on;你可以清晰地看到隨著批發(fā)成本c上升最優(yōu)訂購(gòu)量Q*會(huì)下降因?yàn)槊糠莘e壓的損失風(fēng)險(xiǎn)變大同時(shí)最大期望利潤(rùn)也會(huì)下降。類似的你可以分析售價(jià)p、殘值s或需求波動(dòng)demand_std的影響。5.2 需求分布誤判的風(fēng)險(xiǎn)現(xiàn)實(shí)中我們可能錯(cuò)誤地估計(jì)了需求分布。假設(shè)真實(shí)需求是泊松分布但我們誤以為是正態(tài)分布并據(jù)此制定了訂購(gòu)策略結(jié)果會(huì)怎樣%% 8. 需求分布誤判的風(fēng)險(xiǎn)分析 % 假設(shè)真實(shí)需求服從泊松分布均值 lambda 100 lambda_true 100; true_demand poissrnd(lambda_true, num_days, 1); % 決策者誤以為需求是正態(tài)分布并用歷史數(shù)據(jù)擬合了參數(shù)這里假設(shè)擬合出的均值和標(biāo)準(zhǔn)差恰好也是100和sqrt(100)10 demand_mean_wrong 100; demand_std_wrong sqrt(100); % 泊松分布方差等于均值 % 基于錯(cuò)誤的正態(tài)分布假設(shè)計(jì)算“理論最優(yōu)Q” critical_ratio (unit_price - unit_cost penalty_cost) / ... (unit_price - unit_salvage penalty_cost); Q_decision_wrong round(norminv(critical_ratio, demand_mean_wrong, demand_std_wrong)); % 基于真實(shí)的泊松分布計(jì)算真正的最優(yōu)Q通過(guò)仿真搜索 Q_range_poisson floor(lambda_true - 3*sqrt(lambda_true)) : ceil(lambda_true 3*sqrt(lambda_true)); Q_range_poisson Q_range_poisson(Q_range_poisson 0); profit_poisson zeros(length(Q_range_poisson), 1); for i 1:length(Q_range_poisson) temp_profits zeros(num_days, 1); for day 1:num_days temp_profits(day) calculate_daily_profit(Q_range_poisson(i), true_demand(day), ... unit_price, unit_cost, ... unit_salvage, penalty_cost); end profit_poisson(i) mean(temp_profits); end [true_max_profit, true_opt_idx] max(profit_poisson); Q_decision_true Q_range_poisson(true_opt_idx); % 評(píng)估錯(cuò)誤決策在真實(shí)世界中的表現(xiàn) profits_wrong zeros(num_days, 1); for day 1:num_days profits_wrong(day) calculate_daily_profit(Q_decision_wrong, true_demand(day), ... unit_price, unit_cost, ... unit_salvage, penalty_cost); end avg_profit_wrong mean(profits_wrong); fprintf(\n 需求分布誤判分析 \n); fprintf(真實(shí)需求分布泊松(λ%d)\n, lambda_true); fprintf(決策者誤認(rèn)為正態(tài)(μ%.1f, σ%.1f)\n, demand_mean_wrong, demand_std_wrong); fprintf(基于錯(cuò)誤模型決策的訂購(gòu)量 Q_wrong %d\n, Q_decision_wrong); fprintf(基于真實(shí)模型的最優(yōu)訂購(gòu)量 Q_true %d\n, Q_decision_true); fprintf(錯(cuò)誤決策在真實(shí)環(huán)境下的平均利潤(rùn)%.2f 元\n, avg_profit_wrong); fprintf(正確決策可達(dá)到的最大平均利潤(rùn)%.2f 元\n, true_max_profit); fprintf(因模型誤判導(dǎo)致的利潤(rùn)損失%.2f 元/天 (損失率 %.2f%%)\n, ... true_max_profit - avg_profit_wrong, ... (true_max_profit - avg_profit_wrong)/true_max_profit*100);這個(gè)分析能讓你直觀地感受到錯(cuò)誤的需求模型會(huì)帶來(lái)真金白銀的損失。這也說(shuō)明了在現(xiàn)實(shí)中使用更魯棒的預(yù)測(cè)方法或采用數(shù)據(jù)驅(qū)動(dòng)的仿真優(yōu)化而非依賴強(qiáng)分布假設(shè)的重要性。5.3 擴(kuò)展到多周期與動(dòng)態(tài)規(guī)劃思想經(jīng)典的報(bào)童問(wèn)題是單周期的。但現(xiàn)實(shí)中庫(kù)存可以跨期持有。我們可以做一個(gè)簡(jiǎn)單的兩周期擴(kuò)展思考今天沒(méi)賣(mài)完的報(bào)紙可以留到明天賣(mài)但可能貶值或完全報(bào)廢而明天的需求又是隨機(jī)的。這就變成了一個(gè)動(dòng)態(tài)規(guī)劃問(wèn)題。雖然用MATLAB實(shí)現(xiàn)完整的動(dòng)態(tài)規(guī)劃求解稍復(fù)雜但我們可以用仿真來(lái)近似評(píng)估一個(gè)簡(jiǎn)單的(s, S)策略當(dāng)庫(kù)存低于s時(shí)補(bǔ)貨到S。%% 9. 簡(jiǎn)單多周期仿真思路兩周期帶庫(kù)存結(jié)轉(zhuǎn) % 假設(shè)當(dāng)天未售出報(bào)紙可以以更低的殘值 s2 s 留到第二天銷售。 % 第二天報(bào)紙的批發(fā)價(jià)和售價(jià)不變。 num_periods 2; initial_inventory 0; % 期初庫(kù)存 holding_cost 0.1; % 每份報(bào)紙每周期持有成本如倉(cāng)儲(chǔ)費(fèi) salvage_period2 0.2; % 第二周期末的殘值比第一周期末s更低 % 策略每周期初如果庫(kù)存低于 reorder_point則訂購(gòu)到 order_up_to_level reorder_point 20; order_up_to_level 100; total_profit_multi 0; current_inv initial_inventory; for period 1:num_periods % 本期決策是否補(bǔ)貨補(bǔ)多少 if current_inv reorder_point order_qty order_up_to_level - current_inv; current_inv current_inv order_qty; procurement_cost_this_period unit_cost * order_qty; else order_qty 0; procurement_cost_this_period 0; end % 生成本期需求 period_demand max(round(normrnd(demand_mean, demand_std)), 0); % 計(jì)算本期銷售、剩余等 sales min(current_inv, period_demand); leftover max(0, current_inv - period_demand); shortage max(0, period_demand - current_inv); revenue unit_price * sales; shortage_penalty penalty_cost * shortage; % 本期利潤(rùn)不考慮期末庫(kù)存價(jià)值 period_profit revenue - procurement_cost_this_period - shortage_penalty; total_profit_multi total_profit_multi period_profit; % 庫(kù)存結(jié)轉(zhuǎn)剩余庫(kù)存進(jìn)入下一期但產(chǎn)生持有成本并可能貶值 if period num_periods holding_cost_this holding_cost * leftover; total_profit_multi total_profit_multi - holding_cost_this; current_inv leftover; % 庫(kù)存結(jié)轉(zhuǎn)到下期 else % 最后一期計(jì)算期末殘值 salvage_income salvage_period2 * leftover; total_profit_multi total_profit_multi salvage_income; end end fprintf(\n 簡(jiǎn)單兩周期(s,S)策略仿真 \n); fprintf(策略(s%d, S%d)\n, reorder_point, order_up_to_level); fprintf(兩周期總利潤(rùn)%.2f 元\n, total_profit_multi);這個(gè)簡(jiǎn)單的多周期仿真框架可以很容易地?cái)U(kuò)展到更多周期并用于評(píng)估不同的庫(kù)存策略參數(shù)(s, S)通過(guò)網(wǎng)格搜索或優(yōu)化算法來(lái)尋找長(zhǎng)期最優(yōu)策略。6. 常見(jiàn)問(wèn)題、調(diào)試技巧與性能優(yōu)化在仿真過(guò)程中你可能會(huì)遇到各種問(wèn)題。這里分享一些我踩過(guò)的坑和總結(jié)的技巧。6.1 仿真結(jié)果不穩(wěn)定或與理論值偏差大問(wèn)題每次運(yùn)行程序找到的仿真最優(yōu)Q*都不一樣或者與理論解差距較大。排查與解決增加仿真天數(shù) (num_days)這是最直接有效的方法。大數(shù)定律要求樣本足夠多才能收斂到期望值。對(duì)于報(bào)童問(wèn)題我建議至少num_days10000對(duì)于更精細(xì)的分析可以增加到100000甚至更多。檢查隨機(jī)數(shù)種子在調(diào)試階段為了結(jié)果可復(fù)現(xiàn)可以在腳本開(kāi)頭固定隨機(jī)數(shù)種子rng(12345); % 設(shè)置隨機(jī)種子。這樣每次運(yùn)行都會(huì)生成相同的隨機(jī)需求序列。驗(yàn)證需求分布繪制生成的需求數(shù)據(jù)的直方圖并與你假設(shè)的理論分布概率密度函數(shù)PDF進(jìn)行對(duì)比。使用histfit函數(shù)或ksdensity函數(shù)。figure; histfit(daily_demand, 50, normal); % 擬合正態(tài)分布 title(生成的需求數(shù)據(jù)與正態(tài)分布擬合對(duì)比);細(xì)化搜索步長(zhǎng)在尋找最優(yōu)Q*時(shí)確保Q_range的步長(zhǎng)是1整數(shù)。如果步長(zhǎng)太大可能會(huì)錯(cuò)過(guò)真正的峰值。6.2 代碼運(yùn)行速度慢當(dāng)num_days很大或者需要掃描很多Q值時(shí)循環(huán)嵌套會(huì)導(dǎo)致運(yùn)行變慢。優(yōu)化技巧向量化操作這是 MATLAB 性能提升的關(guān)鍵。避免在循環(huán)內(nèi)進(jìn)行逐元素計(jì)算。例如計(jì)算所有天數(shù)利潤(rùn)的循環(huán)可以改寫(xiě)為% 向量化計(jì)算針對(duì)固定的Q sales_vec min(order_quantity, daily_demand); % 向量與標(biāo)量的min生成向量 leftover_vec max(0, order_quantity - daily_demand); shortage_vec max(0, daily_demand - order_quantity); profit_vec unit_price * sales_vec unit_salvage * leftover_vec ... - unit_cost * order_quantity - penalty_cost * shortage_vec; average_daily_profit mean(profit_vec);這種方法比f(wàn)or循環(huán)快一個(gè)數(shù)量級(jí)。預(yù)分配數(shù)組在循環(huán)前使用zeros()預(yù)分配存儲(chǔ)結(jié)果的大數(shù)組避免數(shù)組在循環(huán)中動(dòng)態(tài)增長(zhǎng)這能顯著提升速度。我們的代碼中已經(jīng)這樣做了。并行計(jì)算如果掃描多個(gè)Q值可以使用parfor循環(huán)需要 Parallel Computing Toolbox。注意并行循環(huán)內(nèi)部的操作需要是獨(dú)立的。avg_profit_list zeros(num_Q, 1); parfor i 1:num_Q % 將 for 改為 parfor current_Q Q_range(i); % ... 計(jì)算 temp_profits ... avg_profit_list(i) mean(temp_profits); end6.3 理論公式計(jì)算報(bào)錯(cuò)NaN或Inf問(wèn)題使用norminv(critical_ratio, mu, sigma)時(shí)返回NaN或Inf。原因與解決norminv函數(shù)的第一個(gè)參數(shù)必須在 (0,1) 開(kāi)區(qū)間內(nèi)。檢查critical_ratio的計(jì)算公式是否正確。確保(p - s g)不為零分母為零意味著模型無(wú)意義。在計(jì)算前對(duì)critical_ratio進(jìn)行鉗制critical_ratio max(min(critical_ratio, 0.9999), 0.0001);這能保證數(shù)值穩(wěn)定性。如果critical_ratio被鉗制到極端值說(shuō)明你的成本參數(shù)設(shè)置導(dǎo)致最優(yōu)策略是“永不訂購(gòu)”或“無(wú)限訂購(gòu)”需要重新審視業(yè)務(wù)參數(shù)。6.4 如何將模型應(yīng)用于實(shí)際數(shù)據(jù)仿真模型的強(qiáng)大之處在于能處理實(shí)際數(shù)據(jù)。假設(shè)你有一份歷史日銷量數(shù)據(jù)historical_sales.csv。數(shù)據(jù)導(dǎo)入與處理data readtable(historical_sales.csv); demand_data data.SalesQuantity; % 假設(shè)列名為SalesQuantity % 注意歷史銷量可能受庫(kù)存限制存在缺貨并非真實(shí)需求。 % 更嚴(yán)謹(jǐn)?shù)淖龇ㄐ枰褂眯枨蠊烙?jì)技術(shù)來(lái)還原未觀測(cè)到的需求。經(jīng)驗(yàn)分布替代理論分布不再假設(shè)正態(tài)分布直接用歷史數(shù)據(jù)的經(jīng)驗(yàn)分布來(lái)生成隨機(jī)需求。% 方法1自助法 (Bootstrap) - 有放回地隨機(jī)抽取歷史數(shù)據(jù) num_days_sim 10000; bootstrap_demand datasample(demand_data, num_days_sim); % 方法2使用經(jīng)驗(yàn)累積分布函數(shù) (ecdf) 和逆變換采樣 [f, x] ecdf(demand_data); % f是累積概率x是對(duì)應(yīng)的需求值 % 生成均勻分布隨機(jī)數(shù)然后插值得到需求 u rand(num_days_sim, 1); ecdf_demand interp1(f, x, u, linear, extrap); ecdf_demand max(round(ecdf_demand), 0); % 取整并確保非負(fù)然后用bootstrap_demand或ecdf_demand替代之前代碼中normrnd生成的需求序列進(jìn)行仿真。這種方法完全由數(shù)據(jù)驅(qū)動(dòng)避免了錯(cuò)誤指定理論分布的風(fēng)險(xiǎn)。通過(guò)這個(gè)從理論到實(shí)踐、從基礎(chǔ)到擴(kuò)展的完整仿真流程你不僅掌握了用 MATLAB 解決報(bào)童問(wèn)題的方法更獲得了一套處理不確定性庫(kù)存決策的建模與分析框架。這個(gè)框架的核心——定義參數(shù)、建立利潤(rùn)模型、生成隨機(jī)場(chǎng)景、評(píng)估策略、優(yōu)化搜索——可以遷移到無(wú)數(shù)類似的運(yùn)營(yíng)決策問(wèn)題中去比如航空公司的超售決策、零售商的季節(jié)性商品采購(gòu)、甚至金融領(lǐng)域的風(fēng)險(xiǎn)管理。真正理解了這個(gè)簡(jiǎn)單的“報(bào)童”你就拿到了打開(kāi)運(yùn)籌優(yōu)化世界大門(mén)的一把鑰匙。 MATLAB仿真報(bào)童問(wèn)題：庫(kù)存決策優(yōu)化與不確定性建模實(shí)踐
1. 項(xiàng)目概述從報(bào)童到庫(kù)存決策的經(jīng)典模型報(bào)童問(wèn)題這個(gè)名字聽(tīng)起來(lái)有點(diǎn)懷舊但它絕不是只存在于歷史課本里的故事。我第一次接觸這個(gè)模型是在研究生階段的一門(mén)運(yùn)籌學(xué)課上當(dāng)時(shí)覺(jué)得它不過(guò)是個(gè)簡(jiǎn)單的概率計(jì)算練習(xí)。直到后來(lái)在電商公司的供應(yīng)鏈部門(mén)實(shí)習(xí)親眼看到每天凌晨算法是如何決定向各個(gè)倉(cāng)庫(kù)補(bǔ)多少貨而第二天又有多少商品因?yàn)槿必浕驕N被標(biāo)記處理時(shí)我才恍然大悟——那個(gè)“賣(mài)報(bào)紙的小孩”面對(duì)的困境正是現(xiàn)代商業(yè)庫(kù)存管理的核心縮影。簡(jiǎn)單來(lái)說(shuō)報(bào)童問(wèn)題描述的是這樣一個(gè)場(chǎng)景一個(gè)報(bào)童每天早晨需要決定從報(bào)社批發(fā)多少份報(bào)紙來(lái)賣(mài)。報(bào)紙的需求量是隨機(jī)的他只知道一個(gè)大概的概率分布。如果批發(fā)多了賣(mài)不完的報(bào)紙到了晚上就一文不值會(huì)造成損失如果批發(fā)少了沒(méi)買(mǎi)到的顧客就走了會(huì)損失潛在的利潤(rùn)。他的目標(biāo)就是找到一個(gè)最優(yōu)的訂購(gòu)量讓他的長(zhǎng)期平均利潤(rùn)最大化或者說(shuō)期望損失最小化。這個(gè)模型的核心就是在不確定性的環(huán)境下做單周期的庫(kù)存決策。今天我們不用真的去賣(mài)報(bào)紙而是用 MATLAB 這個(gè)強(qiáng)大的工具來(lái)親手搭建一個(gè)報(bào)童問(wèn)題的仿真環(huán)境。仿真的意義在于它允許我們?cè)谟?jì)算機(jī)里創(chuàng)造一個(gè)“虛擬世界”在這個(gè)世界里我們可以設(shè)定不同的需求分布、成本參數(shù)然后讓“報(bào)童”按照我們?cè)O(shè)定的策略去運(yùn)營(yíng)成千上萬(wàn)天快速、低成本地觀察不同決策帶來(lái)的長(zhǎng)期結(jié)果。這對(duì)于驗(yàn)證理論公式、比較不同補(bǔ)貨策略、或者處理那些理論模型難以解決的復(fù)雜情況比如需求分布未知、存在缺貨懲罰等來(lái)說(shuō)是極其有效的方法。無(wú)論你是學(xué)習(xí)運(yùn)籌學(xué)、供應(yīng)鏈管理的學(xué)生還是對(duì)數(shù)據(jù)分析和決策優(yōu)化感興趣的從業(yè)者這個(gè)仿真項(xiàng)目都能幫你直觀地理解不確定性決策的精髓。2. 問(wèn)題拆解與數(shù)學(xué)模型建立在動(dòng)手寫(xiě)代碼之前我們必須先把問(wèn)題用數(shù)學(xué)語(yǔ)言清晰地定義出來(lái)。這是所有仿真和分析的基石含糊不得。2.1 核心參數(shù)與變量定義首先我們需要明確幾個(gè)關(guān)鍵的經(jīng)濟(jì)參數(shù)這些是驅(qū)動(dòng)整個(gè)模型的“輸入”單位成本 (c)報(bào)童從報(bào)社批發(fā)一份報(bào)紙需要支付的價(jià)格。這是他的成本。單位售價(jià) (p)報(bào)童將一份報(bào)紙賣(mài)給顧客的價(jià)格。這是他的收入來(lái)源。單位殘值 (s)當(dāng)天結(jié)束時(shí)一份沒(méi)有賣(mài)出去的報(bào)紙的剩余價(jià)值。通常s c可能為零完全報(bào)廢也可能是個(gè)很小的正數(shù)回收價(jià)。缺貨懲罰 (g)這是一個(gè)可選但很實(shí)際的參數(shù)。它表示當(dāng)顧客需要報(bào)紙而報(bào)童缺貨時(shí)所造成的額外損失。這不僅包括失去本次銷售的利潤(rùn) (p - c)還可能包括商譽(yù)損失、顧客流失等隱性成本。在基礎(chǔ)模型中常設(shè)為0但加上它會(huì)讓模型更貼近現(xiàn)實(shí)。接下來(lái)是決策變量和隨機(jī)變量訂購(gòu)量 (Q)這是報(bào)童需要做出的決策也就是我們通過(guò)仿真要尋找的最優(yōu)解。它是一個(gè)非負(fù)整數(shù)。需求量 (D)這是一個(gè)隨機(jī)變量。我們假設(shè)它服從某種已知的概率分布比如正態(tài)分布、泊松分布或均勻分布。仿真的核心之一就是生成符合這個(gè)分布的隨機(jī)需求序列。最后是基于以上變量計(jì)算出的結(jié)果實(shí)際銷量 (Sales)這取決于訂購(gòu)量和需求量中較小的那個(gè)即Sales min(Q, D)。你只能賣(mài)掉你有的和顧客需要的兩者中較少的那部分。剩余庫(kù)存 (Leftover)當(dāng)天結(jié)束時(shí)沒(méi)賣(mài)出去的報(bào)紙即Leftover max(0, Q - D)。缺貨量 (Shortage)當(dāng)天未能滿足的顧客需求即Shortage max(0, D - Q)。2.2 利潤(rùn)函數(shù)與期望利潤(rùn)最大化有了這些定義一天的利潤(rùn)Π(Q, D)就可以寫(xiě)出來(lái)了Π(Q, D) p * min(Q, D) s * max(0, Q - D) - c * Q - g * max(0, D - Q)這個(gè)公式拆開(kāi)看很直觀p * min(Q, D)銷售收入。s * max(0, Q - D)剩余庫(kù)存的殘值回收收入。c * Q批發(fā)報(bào)紙的總成本。g * max(0, D - Q)缺貨造成的懲罰成本。由于需求量D是隨機(jī)的單日的利潤(rùn)也是隨機(jī)的。因此報(bào)童關(guān)心的是長(zhǎng)期平均利潤(rùn)也就是利潤(rùn)的期望值E[Π(Q)]。我們的優(yōu)化目標(biāo)是找到一個(gè)最優(yōu)訂購(gòu)量Q*使得期望利潤(rùn)最大化Q* argmax_{Q≥0} E[Π(Q)]在理論上對(duì)于某些特定的分布如正態(tài)分布存在一個(gè)著名的臨界分位數(shù) (Critical Fractile) 公式來(lái)求解Q*F(Q*) (p - c g) / (p - s g)其中F(·)是需求量D的累積分布函數(shù) (CDF)。這個(gè)公式的意義在于最優(yōu)庫(kù)存水平應(yīng)該設(shè)置在這樣一個(gè)位置需求不超過(guò)該水平的概率恰好等于“單位欠儲(chǔ)成本”與“單位欠儲(chǔ)成本加單位超儲(chǔ)成本”之比。這里(p - c g)可以理解為少進(jìn)一份報(bào)紙?jiān)斐傻倪呺H損失即欠儲(chǔ)成本(p - s g)可以理解為決策的總體邊際影響。注意這個(gè)理論解非常優(yōu)美但它依賴于我們知道準(zhǔn)確的需求分布F(·)。在現(xiàn)實(shí)中分布可能未知、可能隨時(shí)間變化、或者問(wèn)題本身更復(fù)雜如多產(chǎn)品、多周期。這時(shí)仿真 Monte Carlo Simulation 的價(jià)值就凸顯出來(lái)了——我們不需要知道F(·)的解析形式只需要能根據(jù)歷史數(shù)據(jù)或假設(shè)生成隨機(jī)需求樣本就能通過(guò)模擬來(lái)評(píng)估任何給定Q的性能甚至用搜索算法來(lái)尋找近似的Q*。3. MATLAB仿真環(huán)境搭建與核心代碼解析理論鋪墊完畢現(xiàn)在進(jìn)入實(shí)戰(zhàn)環(huán)節(jié)。我們將用 MATLAB 一步步構(gòu)建這個(gè)仿真系統(tǒng)。我個(gè)人的習(xí)慣是先搭建一個(gè)清晰、模塊化的框架這樣調(diào)試和擴(kuò)展都會(huì)很方便。3.1 參數(shù)初始化與需求數(shù)據(jù)生成首先我們創(chuàng)建一個(gè)腳本文件比如叫newsvendor_simulation.m。開(kāi)頭先定義所有基礎(chǔ)參數(shù)。%% 1. 參數(shù)設(shè)置 clear; clc; close all; % 清空環(huán)境好習(xí)慣 % 經(jīng)濟(jì)參數(shù) unit_cost 2; % c: 每份報(bào)紙批發(fā)成本元 unit_price 5; % p: 每份報(bào)紙零售價(jià)格元 unit_salvage 0.5; % s: 每份未售出報(bào)紙的殘值元 penalty_cost 1; % g: 每份缺貨的懲罰成本元可選設(shè)為0則為經(jīng)典模型 % 需求分布參數(shù) - 這里假設(shè)需求服從正態(tài)分布 demand_mean 100; % 平均日需求 demand_std 20; % 日需求標(biāo)準(zhǔn)差 % 仿真參數(shù) num_days 10000; % 模擬的天數(shù)天數(shù)越多結(jié)果越穩(wěn)定 order_quantity 90; % Q: 我們要測(cè)試的訂購(gòu)量可以先設(shè)一個(gè)值跑跑看接下來(lái)是生成隨機(jī)需求。MATLAB 的統(tǒng)計(jì)工具箱提供了豐富的隨機(jī)數(shù)生成器。%% 2. 生成隨機(jī)需求序列 % 使用正態(tài)分布生成需求。注意需求應(yīng)為非負(fù)整數(shù)所以需要取整和取最大值。 daily_demand max(round(normrnd(demand_mean, demand_std, num_days, 1)), 0); % normrnd生成正態(tài)分布隨機(jī)數(shù)round四舍五入取整max(...,0)確保非負(fù)。 % 可視化一下需求分布可選但強(qiáng)烈推薦 figure; subplot(2,1,1); histogram(daily_demand, Normalization, probability); xlabel(日需求量); ylabel(頻率); title(模擬日需求分布直方圖); grid on; subplot(2,1,2); cdfplot(daily_demand); % 繪制經(jīng)驗(yàn)累積分布函數(shù) xlabel(日需求量); ylabel(F(x)); title(需求的經(jīng)驗(yàn)CDF); grid on;實(shí)操心得生成需求時(shí)round和max(...,0)這兩個(gè)處理很重要。現(xiàn)實(shí)中需求是整數(shù)四舍五入更合理。雖然正態(tài)分布理論上可能產(chǎn)生負(fù)數(shù)但我們的demand_mean100,demand_std20產(chǎn)生負(fù)數(shù)的概率極低max(...,0)是一個(gè)安全的保護(hù)措施。如果你模擬的需求均值很小比如接近0則需要考慮使用嚴(yán)格非負(fù)的分布如泊松分布poissrnd(lambda, num_days, 1)。3.2 單周期利潤(rùn)計(jì)算與仿真循環(huán)核心的計(jì)算邏輯封裝成一個(gè)函數(shù)會(huì)非常清晰。我們先寫(xiě)一個(gè)計(jì)算單日利潤(rùn)的函數(shù)。function profit calculate_daily_profit(Q, D, p, c, s, g) % 計(jì)算報(bào)童模型單日利潤(rùn) % 輸入: Q - 訂購(gòu)量, D - 當(dāng)日實(shí)際需求, p,c,s,g - 經(jīng)濟(jì)參數(shù) % 輸出: profit - 當(dāng)日利潤(rùn) sales min(Q, D); % 實(shí)際銷量 leftover max(0, Q - D); % 剩余庫(kù)存 shortage max(0, D - Q); % 缺貨量 revenue p * sales; % 銷售收入 salvage_income s * leftover; % 殘值收入 procurement_cost c * Q; % 采購(gòu)成本 shortage_penalty g * shortage; % 缺貨懲罰 profit revenue salvage_income - procurement_cost - shortage_penalty; end然后在主腳本中我們進(jìn)行仿真循環(huán)計(jì)算長(zhǎng)期平均利潤(rùn)。%% 3. 仿真計(jì)算 daily_profits zeros(num_days, 1); % 預(yù)分配數(shù)組提升效率 for day 1:num_days current_demand daily_demand(day); daily_profits(day) calculate_daily_profit(order_quantity, ... current_demand, ... unit_price, ... unit_cost, ... unit_salvage, ... penalty_cost); end % 計(jì)算關(guān)鍵績(jī)效指標(biāo) (KPIs) average_daily_profit mean(daily_profits); profit_std std(daily_profits); service_level sum(daily_demand order_quantity) / num_days; % 需求滿足率庫(kù)存覆蓋概率 fprintf(仿真結(jié)果訂購(gòu)量 Q%d\n, order_quantity); fprintf( 平均日利潤(rùn): %.2f 元\n, average_daily_profit); fprintf( 利潤(rùn)標(biāo)準(zhǔn)差: %.2f 元\n, profit_std); % 衡量風(fēng)險(xiǎn) fprintf( 服務(wù)水平需求滿足率: %.2f%%\n, service_level * 100);3.3 結(jié)果可視化與分析數(shù)字有了但圖表更能說(shuō)明問(wèn)題。我們來(lái)繪制利潤(rùn)的分布和收斂情況。%% 4. 結(jié)果可視化 figure; % 子圖1日利潤(rùn)分布 subplot(2,2,1); histogram(daily_profits, 50, FaceColor, [0.2 0.6 0.8]); xlabel(日利潤(rùn)元); ylabel(頻數(shù)); title(sprintf(日利潤(rùn)分布 (Q%d), order_quantity)); grid on; hold on; % 標(biāo)記平均利潤(rùn)線 yl ylim; plot([average_daily_profit, average_daily_profit], [yl(1), yl(2)], r--, LineWidth, 2); legend(利潤(rùn)分布, 平均利潤(rùn), Location, best); hold off; % 子圖2累積平均利潤(rùn)看仿真收斂性 subplot(2,2,2); cumulative_avg_profit cumsum(daily_profits) ./ (1:num_days); plot(1:num_days, cumulative_avg_profit, b-, LineWidth, 1.5); xlabel(模擬天數(shù)); ylabel(累積平均利潤(rùn)元); title(平均利潤(rùn)隨仿真天數(shù)的收斂過(guò)程); grid on; hold on; plot([1, num_days], [average_daily_profit, average_daily_profit], r--); legend(累積平均, 最終平均, Location, southeast); hold off; % 子圖3利潤(rùn)與需求的關(guān)系散點(diǎn)圖 subplot(2,2,3); scatter(daily_demand, daily_profits, 10, filled, MarkerFaceAlpha, 0.6); xlabel(日需求量); ylabel(日利潤(rùn)); title(需求與利潤(rùn)關(guān)系散點(diǎn)圖); grid on; % 可以添加趨勢(shì)線或分界線 hold on; plot([order_quantity, order_quantity], ylim, k--, LineWidth, 1.5); % 標(biāo)記訂購(gòu)量 hold off; % 子圖4不同需求下的利潤(rùn)構(gòu)成示例取一天 subplot(2,2,4); sample_day find(daily_demand round(demand_mean), 1); % 找一個(gè)需求接近均值的天 if isempty(sample_day) sample_day 1; end sample_demand daily_demand(sample_day); [sales, leftover, shortage] deal(min(order_quantity, sample_demand), ... max(0, order_quantity - sample_demand), ... max(0, sample_demand - order_quantity)); profit_breakdown [unit_price*sales, unit_salvage*leftover, -unit_cost*order_quantity, -penalty_cost*shortage]; labels {銷售收入, 殘值收入, 采購(gòu)成本, 缺貨懲罰}; bar(profit_breakdown); set(gca, XTickLabel, labels); ylabel(金額元); title(sprintf(第%d天利潤(rùn)構(gòu)成 (需求%d), sample_day, sample_demand)); grid on;運(yùn)行這段代碼你就能得到一個(gè)完整的單點(diǎn)仿真結(jié)果。但我們的目標(biāo)是找到最優(yōu)的Q*所以下一步是進(jìn)行敏感性分析。4. 尋找最優(yōu)訂購(gòu)量仿真與理論對(duì)比現(xiàn)在我們讓Q動(dòng)起來(lái)觀察平均利潤(rùn)如何隨Q變化并嘗試找到那個(gè)最高點(diǎn)。4.1 遍歷搜索與利潤(rùn)曲線繪制我們?cè)O(shè)定一個(gè)Q的搜索范圍比如從demand_mean - 3*demand_std到demand_mean 3*demand_std覆蓋需求的絕大部分可能區(qū)間。%% 5. 尋找最優(yōu)訂購(gòu)量 Q* % 定義搜索范圍 Q_range floor(demand_mean - 3*demand_std) : ceil(demand_mean 3*demand_std); Q_range Q_range(Q_range 0); % 確保非負(fù) num_Q length(Q_range); avg_profit_list zeros(num_Q, 1); service_level_list zeros(num_Q, 1); fprintf(開(kāi)始掃描 %d 個(gè)不同的Q值...\n, num_Q); % 對(duì)每個(gè)Q進(jìn)行仿真。注意這里為了速度復(fù)用之前生成的需求序列。 % 如果追求絕對(duì)準(zhǔn)確應(yīng)對(duì)每個(gè)Q重新生成獨(dú)立的需求序列但計(jì)算量會(huì)大很多。 % 在Q值掃描中使用同一組需求序列是標(biāo)準(zhǔn)做法保證了比較的公平性。 for i 1:num_Q current_Q Q_range(i); temp_profits zeros(num_days, 1); for day 1:num_days temp_profits(day) calculate_daily_profit(current_Q, daily_demand(day), ... unit_price, unit_cost, ... unit_salvage, penalty_cost); end avg_profit_list(i) mean(temp_profits); service_level_list(i) sum(daily_demand current_Q) / num_days; end % 找到仿真下的最優(yōu)Q [sim_max_profit, sim_opt_idx] max(avg_profit_list); sim_opt_Q Q_range(sim_opt_idx); fprintf(【仿真結(jié)果】最優(yōu)訂購(gòu)量 Q* %d對(duì)應(yīng)平均日利潤(rùn) %.2f 元服務(wù)水平 %.2f%%\n, ... sim_opt_Q, sim_max_profit, service_level_list(sim_opt_idx)*100);繪制利潤(rùn)-訂購(gòu)量曲線。% 可視化利潤(rùn)曲線 figure; yyaxis left; plot(Q_range, avg_profit_list, b-o, LineWidth, 1.5, MarkerSize, 4); hold on; plot(sim_opt_Q, sim_max_profit, r*, MarkerSize, 15, LineWidth, 2); xlabel(訂購(gòu)量 Q); ylabel(平均日利潤(rùn)元); yyaxis right; plot(Q_range, service_level_list*100, g--s, LineWidth, 1.5, MarkerSize, 4); ylabel(服務(wù)水平 (%)); title(平均利潤(rùn)與服務(wù)水平隨訂購(gòu)量變化曲線); grid on; legend(平均利潤(rùn), sprintf(最優(yōu)點(diǎn) (Q%d), sim_opt_Q), 服務(wù)水平, ... Location, best);你會(huì)看到一條經(jīng)典的凹曲線利潤(rùn)先隨Q增加而上升因?yàn)槟茏プ「噤N售機(jī)會(huì)達(dá)到一個(gè)頂峰后開(kāi)始下降因?yàn)闇N損失開(kāi)始超過(guò)新增銷售的收益。那個(gè)頂峰對(duì)應(yīng)的Q就是我們的仿真最優(yōu)解。4.2 理論解計(jì)算與對(duì)比現(xiàn)在我們用前面提到的臨界分位數(shù)公式來(lái)計(jì)算理論最優(yōu)解并與仿真結(jié)果對(duì)比。%% 6. 理論解計(jì)算與對(duì)比 % 計(jì)算臨界分位數(shù) critical_ratio (unit_price - unit_cost penalty_cost) / ... (unit_price - unit_salvage penalty_cost); fprintf(臨界分位數(shù) (p - c g) / (p - s g) %.4f\n, critical_ratio); % 由于我們假設(shè)需求服從正態(tài)分布 N(mu, sigma^2) % 理論最優(yōu)Q*是滿足 F(Q*) critical_ratio 的值即逆CDF % 使用 norminv 函數(shù) theory_opt_Q norminv(critical_ratio, demand_mean, demand_std); theory_opt_Q round(theory_opt_Q); % 取整因?yàn)镼是整數(shù) fprintf(【理論解】最優(yōu)訂購(gòu)量 Q*_theory %.2f (取整后為 %d)\n, ... norminv(critical_ratio, demand_mean, demand_std), theory_opt_Q); % 計(jì)算理論解對(duì)應(yīng)的仿真利潤(rùn)用同一組需求數(shù)據(jù)評(píng)估 theory_profits zeros(num_days, 1); for day 1:num_days theory_profits(day) calculate_daily_profit(theory_opt_Q, daily_demand(day), ... unit_price, unit_cost, ... unit_salvage, penalty_cost); end theory_avg_profit mean(theory_profits); theory_service_level sum(daily_demand theory_opt_Q) / num_days; fprintf(理論解Q%d對(duì)應(yīng)的仿真評(píng)估平均利潤(rùn)%.2f服務(wù)水平%.2f%%\n, ... theory_opt_Q, theory_avg_profit, theory_service_level*100); % 對(duì)比分析 comparison_table table([sim_opt_Q; theory_opt_Q], ... [sim_max_profit; theory_avg_profit], ... [service_level_list(sim_opt_idx); theory_service_level]*100, ... VariableNames, {最優(yōu)訂購(gòu)量Q, 平均日利潤(rùn), 服務(wù)水平_百分比}, ... RowNames, {仿真搜索, 理論公式}); disp(comparison_table);正常情況下仿真搜索得到的Q*和理論公式計(jì)算的Q*應(yīng)該非常接近。如果差異較大可能的原因有1) 仿真天數(shù)num_days不夠多結(jié)果有波動(dòng)2) 需求分布不是完美的正態(tài)分布因?yàn)槲覀冏隽巳≌腿》秦?fù)處理3) 搜索的步長(zhǎng)不夠精細(xì)。增加num_days和縮小Q_range的步長(zhǎng)例如以1為步進(jìn)可以改善。注意事項(xiàng)norminv函數(shù)要求critical_ratio在 (0,1) 開(kāi)區(qū)間內(nèi)。如果您的成本參數(shù)設(shè)置導(dǎo)致critical_ratio非常接近0或1例如售價(jià)遠(yuǎn)低于成本norminv可能會(huì)返回-Inf或Inf。在實(shí)際業(yè)務(wù)中這通常意味著最優(yōu)策略是“不訂購(gòu)”或“訂購(gòu)極大數(shù)量”需要在實(shí)際代碼中加入邊界判斷。5. 深入分析與擴(kuò)展應(yīng)用場(chǎng)景基礎(chǔ)仿真跑通后我們可以玩點(diǎn)更花的讓模型更貼近復(fù)雜的現(xiàn)實(shí)情況。5.1 敏感性分析參數(shù)如何影響決策最優(yōu)訂購(gòu)量Q*對(duì)成本參數(shù)非常敏感。我們可以系統(tǒng)地改變一個(gè)參數(shù)比如單位成本c觀察Q*和最大利潤(rùn)的變化。%% 7. 敏感性分析示例單位成本c的影響 cost_range 1.5:0.1:2.5; % 單位成本從1.5元到2.5元變化 num_costs length(cost_range); opt_Q_vs_cost zeros(num_costs, 1); max_profit_vs_cost zeros(num_costs, 1); % 固定其他參數(shù)和需求序列 for i 1:num_costs current_cost cost_range(i); % 計(jì)算當(dāng)前成本下的臨界分位數(shù)和理論Q* current_cr (unit_price - current_cost penalty_cost) / ... (unit_price - unit_salvage penalty_cost); % 防止cr超出(0,1)范圍 current_cr max(min(current_cr, 0.999), 0.001); current_opt_Q round(norminv(current_cr, demand_mean, demand_std)); opt_Q_vs_cost(i) current_opt_Q; % 評(píng)估該Q下的仿真利潤(rùn) temp_profits zeros(num_days, 1); for day 1:num_days temp_profits(day) calculate_daily_profit(current_opt_Q, daily_demand(day), ... unit_price, current_cost, ... unit_salvage, penalty_cost); end max_profit_vs_cost(i) mean(temp_profits); end figure; subplot(2,1,1); plot(cost_range, opt_Q_vs_cost, b-o, LineWidth, 1.5); xlabel(單位成本 c (元)); ylabel(最優(yōu)訂購(gòu)量 Q*); title(最優(yōu)訂購(gòu)量隨單位成本變化); grid on; subplot(2,1,2); plot(cost_range, max_profit_vs_cost, r-s, LineWidth, 1.5); xlabel(單位成本 c (元)); ylabel(最大期望利潤(rùn) (元)); title(最大期望利潤(rùn)隨單位成本變化); grid on;你可以清晰地看到隨著批發(fā)成本c上升最優(yōu)訂購(gòu)量Q*會(huì)下降因?yàn)槊糠莘e壓的損失風(fēng)險(xiǎn)變大同時(shí)最大期望利潤(rùn)也會(huì)下降。類似的你可以分析售價(jià)p、殘值s或需求波動(dòng)demand_std的影響。5.2 需求分布誤判的風(fēng)險(xiǎn)現(xiàn)實(shí)中我們可能錯(cuò)誤地估計(jì)了需求分布。假設(shè)真實(shí)需求是泊松分布但我們誤以為是正態(tài)分布并據(jù)此制定了訂購(gòu)策略結(jié)果會(huì)怎樣%% 8. 需求分布誤判的風(fēng)險(xiǎn)分析 % 假設(shè)真實(shí)需求服從泊松分布均值 lambda 100 lambda_true 100; true_demand poissrnd(lambda_true, num_days, 1); % 決策者誤以為需求是正態(tài)分布并用歷史數(shù)據(jù)擬合了參數(shù)這里假設(shè)擬合出的均值和標(biāo)準(zhǔn)差恰好也是100和sqrt(100)10 demand_mean_wrong 100; demand_std_wrong sqrt(100); % 泊松分布方差等于均值 % 基于錯(cuò)誤的正態(tài)分布假設(shè)計(jì)算“理論最優(yōu)Q” critical_ratio (unit_price - unit_cost penalty_cost) / ... (unit_price - unit_salvage penalty_cost); Q_decision_wrong round(norminv(critical_ratio, demand_mean_wrong, demand_std_wrong)); % 基于真實(shí)的泊松分布計(jì)算真正的最優(yōu)Q通過(guò)仿真搜索 Q_range_poisson floor(lambda_true - 3*sqrt(lambda_true)) : ceil(lambda_true 3*sqrt(lambda_true)); Q_range_poisson Q_range_poisson(Q_range_poisson 0); profit_poisson zeros(length(Q_range_poisson), 1); for i 1:length(Q_range_poisson) temp_profits zeros(num_days, 1); for day 1:num_days temp_profits(day) calculate_daily_profit(Q_range_poisson(i), true_demand(day), ... unit_price, unit_cost, ... unit_salvage, penalty_cost); end profit_poisson(i) mean(temp_profits); end [true_max_profit, true_opt_idx] max(profit_poisson); Q_decision_true Q_range_poisson(true_opt_idx); % 評(píng)估錯(cuò)誤決策在真實(shí)世界中的表現(xiàn) profits_wrong zeros(num_days, 1); for day 1:num_days profits_wrong(day) calculate_daily_profit(Q_decision_wrong, true_demand(day), ... unit_price, unit_cost, ... unit_salvage, penalty_cost); end avg_profit_wrong mean(profits_wrong); fprintf(\n 需求分布誤判分析 \n); fprintf(真實(shí)需求分布泊松(λ%d)\n, lambda_true); fprintf(決策者誤認(rèn)為正態(tài)(μ%.1f, σ%.1f)\n, demand_mean_wrong, demand_std_wrong); fprintf(基于錯(cuò)誤模型決策的訂購(gòu)量 Q_wrong %d\n, Q_decision_wrong); fprintf(基于真實(shí)模型的最優(yōu)訂購(gòu)量 Q_true %d\n, Q_decision_true); fprintf(錯(cuò)誤決策在真實(shí)環(huán)境下的平均利潤(rùn)%.2f 元\n, avg_profit_wrong); fprintf(正確決策可達(dá)到的最大平均利潤(rùn)%.2f 元\n, true_max_profit); fprintf(因模型誤判導(dǎo)致的利潤(rùn)損失%.2f 元/天 (損失率 %.2f%%)\n, ... true_max_profit - avg_profit_wrong, ... (true_max_profit - avg_profit_wrong)/true_max_profit*100);這個(gè)分析能讓你直觀地感受到錯(cuò)誤的需求模型會(huì)帶來(lái)真金白銀的損失。這也說(shuō)明了在現(xiàn)實(shí)中使用更魯棒的預(yù)測(cè)方法或采用數(shù)據(jù)驅(qū)動(dòng)的仿真優(yōu)化而非依賴強(qiáng)分布假設(shè)的重要性。5.3 擴(kuò)展到多周期與動(dòng)態(tài)規(guī)劃思想經(jīng)典的報(bào)童問(wèn)題是單周期的。但現(xiàn)實(shí)中庫(kù)存可以跨期持有。我們可以做一個(gè)簡(jiǎn)單的兩周期擴(kuò)展思考今天沒(méi)賣(mài)完的報(bào)紙可以留到明天賣(mài)但可能貶值或完全報(bào)廢而明天的需求又是隨機(jī)的。這就變成了一個(gè)動(dòng)態(tài)規(guī)劃問(wèn)題。雖然用MATLAB實(shí)現(xiàn)完整的動(dòng)態(tài)規(guī)劃求解稍復(fù)雜但我們可以用仿真來(lái)近似評(píng)估一個(gè)簡(jiǎn)單的(s, S)策略當(dāng)庫(kù)存低于s時(shí)補(bǔ)貨到S。%% 9. 簡(jiǎn)單多周期仿真思路兩周期帶庫(kù)存結(jié)轉(zhuǎn) % 假設(shè)當(dāng)天未售出報(bào)紙可以以更低的殘值 s2 s 留到第二天銷售。 % 第二天報(bào)紙的批發(fā)價(jià)和售價(jià)不變。 num_periods 2; initial_inventory 0; % 期初庫(kù)存 holding_cost 0.1; % 每份報(bào)紙每周期持有成本如倉(cāng)儲(chǔ)費(fèi) salvage_period2 0.2; % 第二周期末的殘值比第一周期末s更低 % 策略每周期初如果庫(kù)存低于 reorder_point則訂購(gòu)到 order_up_to_level reorder_point 20; order_up_to_level 100; total_profit_multi 0; current_inv initial_inventory; for period 1:num_periods % 本期決策是否補(bǔ)貨補(bǔ)多少 if current_inv reorder_point order_qty order_up_to_level - current_inv; current_inv current_inv order_qty; procurement_cost_this_period unit_cost * order_qty; else order_qty 0; procurement_cost_this_period 0; end % 生成本期需求 period_demand max(round(normrnd(demand_mean, demand_std)), 0); % 計(jì)算本期銷售、剩余等 sales min(current_inv, period_demand); leftover max(0, current_inv - period_demand); shortage max(0, period_demand - current_inv); revenue unit_price * sales; shortage_penalty penalty_cost * shortage; % 本期利潤(rùn)不考慮期末庫(kù)存價(jià)值 period_profit revenue - procurement_cost_this_period - shortage_penalty; total_profit_multi total_profit_multi period_profit; % 庫(kù)存結(jié)轉(zhuǎn)剩余庫(kù)存進(jìn)入下一期但產(chǎn)生持有成本并可能貶值 if period num_periods holding_cost_this holding_cost * leftover; total_profit_multi total_profit_multi - holding_cost_this; current_inv leftover; % 庫(kù)存結(jié)轉(zhuǎn)到下期 else % 最后一期計(jì)算期末殘值 salvage_income salvage_period2 * leftover; total_profit_multi total_profit_multi salvage_income; end end fprintf(\n 簡(jiǎn)單兩周期(s,S)策略仿真 \n); fprintf(策略(s%d, S%d)\n, reorder_point, order_up_to_level); fprintf(兩周期總利潤(rùn)%.2f 元\n, total_profit_multi);這個(gè)簡(jiǎn)單的多周期仿真框架可以很容易地?cái)U(kuò)展到更多周期并用于評(píng)估不同的庫(kù)存策略參數(shù)(s, S)通過(guò)網(wǎng)格搜索或優(yōu)化算法來(lái)尋找長(zhǎng)期最優(yōu)策略。6. 常見(jiàn)問(wèn)題、調(diào)試技巧與性能優(yōu)化在仿真過(guò)程中你可能會(huì)遇到各種問(wèn)題。這里分享一些我踩過(guò)的坑和總結(jié)的技巧。6.1 仿真結(jié)果不穩(wěn)定或與理論值偏差大問(wèn)題每次運(yùn)行程序找到的仿真最優(yōu)Q*都不一樣或者與理論解差距較大。排查與解決增加仿真天數(shù) (num_days)這是最直接有效的方法。大數(shù)定律要求樣本足夠多才能收斂到期望值。對(duì)于報(bào)童問(wèn)題我建議至少num_days10000對(duì)于更精細(xì)的分析可以增加到100000甚至更多。檢查隨機(jī)數(shù)種子在調(diào)試階段為了結(jié)果可復(fù)現(xiàn)可以在腳本開(kāi)頭固定隨機(jī)數(shù)種子rng(12345); % 設(shè)置隨機(jī)種子。這樣每次運(yùn)行都會(huì)生成相同的隨機(jī)需求序列。驗(yàn)證需求分布繪制生成的需求數(shù)據(jù)的直方圖并與你假設(shè)的理論分布概率密度函數(shù)PDF進(jìn)行對(duì)比。使用histfit函數(shù)或ksdensity函數(shù)。figure; histfit(daily_demand, 50, normal); % 擬合正態(tài)分布 title(生成的需求數(shù)據(jù)與正態(tài)分布擬合對(duì)比);細(xì)化搜索步長(zhǎng)在尋找最優(yōu)Q*時(shí)確保Q_range的步長(zhǎng)是1整數(shù)。如果步長(zhǎng)太大可能會(huì)錯(cuò)過(guò)真正的峰值。6.2 代碼運(yùn)行速度慢當(dāng)num_days很大或者需要掃描很多Q值時(shí)循環(huán)嵌套會(huì)導(dǎo)致運(yùn)行變慢。優(yōu)化技巧向量化操作這是 MATLAB 性能提升的關(guān)鍵。避免在循環(huán)內(nèi)進(jìn)行逐元素計(jì)算。例如計(jì)算所有天數(shù)利潤(rùn)的循環(huán)可以改寫(xiě)為% 向量化計(jì)算針對(duì)固定的Q sales_vec min(order_quantity, daily_demand); % 向量與標(biāo)量的min生成向量 leftover_vec max(0, order_quantity - daily_demand); shortage_vec max(0, daily_demand - order_quantity); profit_vec unit_price * sales_vec unit_salvage * leftover_vec ... - unit_cost * order_quantity - penalty_cost * shortage_vec; average_daily_profit mean(profit_vec);這種方法比f(wàn)or循環(huán)快一個(gè)數(shù)量級(jí)。預(yù)分配數(shù)組在循環(huán)前使用zeros()預(yù)分配存儲(chǔ)結(jié)果的大數(shù)組避免數(shù)組在循環(huán)中動(dòng)態(tài)增長(zhǎng)這能顯著提升速度。我們的代碼中已經(jīng)這樣做了。并行計(jì)算如果掃描多個(gè)Q值可以使用parfor循環(huán)需要 Parallel Computing Toolbox。注意并行循環(huán)內(nèi)部的操作需要是獨(dú)立的。avg_profit_list zeros(num_Q, 1); parfor i 1:num_Q % 將 for 改為 parfor current_Q Q_range(i); % ... 計(jì)算 temp_profits ... avg_profit_list(i) mean(temp_profits); end6.3 理論公式計(jì)算報(bào)錯(cuò)NaN或Inf問(wèn)題使用norminv(critical_ratio, mu, sigma)時(shí)返回NaN或Inf。原因與解決norminv函數(shù)的第一個(gè)參數(shù)必須在 (0,1) 開(kāi)區(qū)間內(nèi)。檢查critical_ratio的計(jì)算公式是否正確。確保(p - s g)不為零分母為零意味著模型無(wú)意義。在計(jì)算前對(duì)critical_ratio進(jìn)行鉗制critical_ratio max(min(critical_ratio, 0.9999), 0.0001);這能保證數(shù)值穩(wěn)定性。如果critical_ratio被鉗制到極端值說(shuō)明你的成本參數(shù)設(shè)置導(dǎo)致最優(yōu)策略是“永不訂購(gòu)”或“無(wú)限訂購(gòu)”需要重新審視業(yè)務(wù)參數(shù)。6.4 如何將模型應(yīng)用于實(shí)際數(shù)據(jù)仿真模型的強(qiáng)大之處在于能處理實(shí)際數(shù)據(jù)。假設(shè)你有一份歷史日銷量數(shù)據(jù)historical_sales.csv。數(shù)據(jù)導(dǎo)入與處理data readtable(historical_sales.csv); demand_data data.SalesQuantity; % 假設(shè)列名為SalesQuantity % 注意歷史銷量可能受庫(kù)存限制存在缺貨并非真實(shí)需求。 % 更嚴(yán)謹(jǐn)?shù)淖龇ㄐ枰褂眯枨蠊烙?jì)技術(shù)來(lái)還原未觀測(cè)到的需求。經(jīng)驗(yàn)分布替代理論分布不再假設(shè)正態(tài)分布直接用歷史數(shù)據(jù)的經(jīng)驗(yàn)分布來(lái)生成隨機(jī)需求。% 方法1自助法 (Bootstrap) - 有放回地隨機(jī)抽取歷史數(shù)據(jù) num_days_sim 10000; bootstrap_demand datasample(demand_data, num_days_sim); % 方法2使用經(jīng)驗(yàn)累積分布函數(shù) (ecdf) 和逆變換采樣 [f, x] ecdf(demand_data); % f是累積概率x是對(duì)應(yīng)的需求值 % 生成均勻分布隨機(jī)數(shù)然后插值得到需求 u rand(num_days_sim, 1); ecdf_demand interp1(f, x, u, linear, extrap); ecdf_demand max(round(ecdf_demand), 0); % 取整并確保非負(fù)然后用bootstrap_demand或ecdf_demand替代之前代碼中normrnd生成的需求序列進(jìn)行仿真。這種方法完全由數(shù)據(jù)驅(qū)動(dòng)避免了錯(cuò)誤指定理論分布的風(fēng)險(xiǎn)。通過(guò)這個(gè)從理論到實(shí)踐、從基礎(chǔ)到擴(kuò)展的完整仿真流程你不僅掌握了用 MATLAB 解決報(bào)童問(wèn)題的方法更獲得了一套處理不確定性庫(kù)存決策的建模與分析框架。這個(gè)框架的核心——定義參數(shù)、建立利潤(rùn)模型、生成隨機(jī)場(chǎng)景、評(píng)估策略、優(yōu)化搜索——可以遷移到無(wú)數(shù)類似的運(yùn)營(yíng)決策問(wèn)題中去比如航空公司的超售決策、零售商的季節(jié)性商品采購(gòu)、甚至金融領(lǐng)域的風(fēng)險(xiǎn)管理。真正理解了這個(gè)簡(jiǎn)單的“報(bào)童”你就拿到了打開(kāi)運(yùn)籌優(yōu)化世界大門(mén)的一把鑰匙。