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