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