現代庫存管理：模型、算法與Python實現第1章現代庫存管理概述1.1什么是庫存庫存是指為了使生產正常且不間斷的開展或為了更加及時地滿足顧客的需求，在各個生產階段或分銷環節之間設置的物品儲備庫存是供應鏈的六大功能模塊之一，在促進企業降本增效、供需匹配中發揮關鍵作用對于企業來說，庫存是一項流動資產，會給投資帶來相應的回報庫存的四方面功能：供需匹配：持有安全庫存應對不確定性規模效應：促使企業批量補貨來降低補貨的成本快速履約：持有更多的庫存，庫存分散前置庫存分離：供應鏈網絡一定程度的解耦1.2庫存管理解決的問題通過適量的庫存，實現超儲與缺貨之間的最佳或最經濟合理的平衡，是庫存管理的核心目標企業的存貨不足，可能會造成供貨不及時，供應鏈斷裂、喪失交易機會或市場占有率商品庫存會產生一定的持貨成本，同時存在由于商品擠壓或者損壞而產生的庫存風險因此，在庫存管理中，既要保持合理的庫存數量，防止缺貨、實現一定的規模效應，同時又要避免庫存過量，發生不必要的庫存費用1.2庫存管理解決的問題以簡單兩級分銷網絡為例，中央倉從供應商補貨，然后為下游的區域倉供應庫存，每個區域倉對相應區域的顧客需求進行履約1.3現代庫存管理的挑戰在具體業務場景下，實際因素使得庫存管理目標的實現更具挑戰庫存管理場景的多樣化啤酒飲料行業：重量大、體積大、貨值低的重泡貨售賣鮮食的零售商超、餐飲門店：保質期較短消費電子行業代工廠：采用訂單生產模式危險化學品的補貨過程：特殊運輸策略，保障安全運輸需求來源的多樣化與庫存網絡的復雜性零售企業線上線下全渠道融合，企業面臨的訂單來源多樣化為了快速滿足全渠道的需求，企業庫存供應網絡的結構日漸復雜1.3現代庫存管理的挑戰需求與供應不確定性高企業促銷手段層出不窮在售的商品種類繁多，很大一部分是長尾品疫情的反復變化消費者對缺貨的容忍度降低消費者的轉移成本降低，粘性變弱，消費者缺乏耐?1.4

現代庫存管理的發展趨勢在消費、產業升級和技術變革的作用之下，庫存網絡逐步發展為依賴人工智能、機器學習、云計算等新興技術的數字化矩陣式庫存網絡，呈現出以下幾點發展趨勢：庫存網絡扁平化代理商的去中間化庫存設置點前移布局前置倉，提升用戶體驗物流配送格局逐漸轉變為電商平臺+前置倉+即時物流（或消費者）1.4

現代庫存管理的發展趨勢數字化建設互聯網、物聯網、大數據、人工智能等新一代信息技術和現代化管理理念方法作為基礎以價值創造為導向、以數據為驅動，對供應鏈活動進行整體規劃設計與運作1.5現代庫存管理優化案例近年來，國內涌現出了像杉數科技（北京）有限公司（簡稱杉數科技）這樣為數字化轉型的企業提供智慧供應鏈解決方案的公司接下來介紹杉數科技的兩個庫存管理優化案例1.5現代庫存管理優化案例背景介紹企業B是目前全球頂級的食品制造商除了線下門店，公司布局了多個線上渠道、供應中心和工廠全渠道履約網絡零售企業全渠道庫存優化案例1.5現代庫存管理優化案例面臨的問題針對需求特性不一的渠道，如何實現智能預測并提高履約能力隨著每年產品數量的增加，以及各類產品繁復的市場特征，庫存缺貨和冗余現象頻繁發生，導致庫存成本持續增加解決方案根據渠道特性定制需求預測算法模型，為企業量身打造了新一代電商智能計劃系統可視化展示系統儲存的數據信息實現全渠道需求的預測和計劃制定提供算法輸出參考以及批量編輯功能幫助用戶提高計劃效率采用插件式的設計架構，可以基于場景定制更新算法包方案落地后，測試商品的月需求預測準確率提升了約6%，助力企業供應鏈智能化轉型升級1.5現代庫存管理優化案例背景介紹H公司是全球領先的信息與通信基礎設施和智能終端提供商生產制造網絡中包含數十萬個節點，數百萬條邊，存在大量共享組件及可替代原材料過去的安全庫存策略流程煩瑣，成品交付滿足率較低，也沒有深度挖掘全網絡安全庫存共享的價值面臨的問題安全庫存是庫存管理中應對需求不確定性的有效工具，如何在規模龐大的復雜制造網絡上設置安全庫存，在提升終端產品的服務水平的同時最小化系統總安全庫存成本如何對整個制造網絡上的安全庫存進行全局優化，有效地刻畫需求不確定性、全局優化安全庫存布局、平衡成本與滿足率的情況下優化安全庫存量、同時保證求解效率復雜制造網絡安全庫存優化案例1.5現代庫存管理優化案例打造了一套自動化、智能化的網絡安全庫存優化解決方案解決方案現代庫存管理：模型、算法與Python實現第2章需求預測基礎2.1需求預測的基本結構一般認為需求的實現值由兩部分構成:可預測的需求模式：這部分由商品特征、經濟規律等客觀規律決定不可預測的隨機因素：由一些難以刻畫的顧客購買行為和隨機影響導致

2.1需求預測的基本結構

2.1需求預測的基本結構產品的需求模式同時受到內在規律和外生因素的影響內在規律可以進一步分解為趨勢性、季節性這兩種共通的需求變動規律和產品本身決定的需求特性，由產品本身的性質決定外生影響的來源則非常廣泛，往往需要針對具體的產品進行具體分析和刻畫2.2需求預測方法的分類時間序列方法

宏觀趨勢

2.2需求預測方法的分類時間序列方法在使用時間序列模型時，最重要的是理解模型采用了哪一種分解模型，以及刻畫的是序列中的哪一部分

宏觀因素微觀因素

趨勢性季節性一元多元簡單指數平滑√√

Holt-Winters√√

ARIMA

√

sARIMA

√√

VARMA

√sVARMA

√

√HTS

√prophet√√√

2.2需求預測方法的分類機器學習方法特點優勢劣勢線性模型簡單的線性回歸形式，可以通過正則化項演化為LASSO或Ridge模型形式簡單只能擬合線性關系CART樹基于二分樹的結構構建，是后續許多集成方法的基礎結構簡單，擬合能力強大對數據敏感，容易過擬合神經網絡由線性變換、非線性激活函數和網絡結構組成多層網絡理論上可以以任意精度逼近任意函數，具有強大的擬合能力網絡結構復雜，模型參數較多，容易調試不當或過擬合隨機森林以并行的決策樹為基礎模型，通過簡單加權集成最終結果相比決策樹更加穩健，可以并行運算，計算速度快相比于GBDT集成方法效果較差XGBoost以決策樹為基礎模型，不斷強化對殘差的學習基于GBDT思想改進，能在實現較好模型效果同時以較快的速度完成運算對于較大的數據運行較慢；參數設置復雜，容易過擬合lightGBM以決策樹為基礎模型，不斷強化對殘差的學習在海量數據的場景下，能夠以較小的內存和較快的速度完成計算參數設置復雜，容易過擬合Stacking以多種預測模型為基礎模型，通過機器學習模型集成最終結果模型形式簡單，效果強大，能夠集成不同模型的優勢，包容性強容易過擬合2.2需求預測方法的分類比較和總結機器學習方法是數據驅動的，時間序列方法是模型驅動的機器學習模型的效果受到特征選擇、模型選擇、參數選擇等諸多因素的影響，實際運用效果還要取決于業務場景、數據性質和模型調試等諸多方面由于強大的擬合能力，在一個數據規模相對較小或是內在規律簡單的數據集上，機器學習很容易過度擬合隨機擾動因素，從而導致較差的泛化能力；而時間序列則能夠集中于把握宏觀趨勢，較好地剔除隨機擾動的偶然影響機器學習得到的結果往往較難得到直觀理解，解釋力較弱，較難結合已有經驗進行輔助判斷和調整；時間序列模型具有簡單清晰的模型結構，非常便于使用、解讀和調整2.3數據預處理在進行需求預測之前，通常要對數據進行預處理及分析。盡管不同企業的業務場景與產品的需求規律不同，但數據處理的步驟與方法是類似的本節以食品飲料快消品公司W2018年1月1日～2020年7月30日77個產品在18個分銷中心（distributioncenter,DC）的銷售數據為對象一般而言，盡管原始數據的具體格式不盡相同，需求預測任務的數據集至少包含以下信息：

商品名，有時存在品類信息，構成兩級及以上的產品類別層級結構銷售地點，有時存在國家/區域—城市—終端銷售點的地理層級結構時間，一般以日為顆粒度實際銷量首先，從數據文件夾讀入銷售數據表格，并查看表格前十行，以獲得對數據結構的基本了解注意數據路徑的設置，相比于將文件地址作為一個完整字段，單獨指定文件夾路徑和表格名稱更具靈活性，更適用于具有較多數據表格需要管理以及可能存在數據遷移的場景2.3數據預處理數據導入及基礎操作importpandasaspd

data_dir=

'../../data/forecast_data/'

sales_data=pd.read_csv(data_dir+

'sales_data.csv')

print(sales_data.head(10))2.3數據預處理數據導入及基礎操作可以看到，本節所使用的數據集恰好對應了需求預測任務的四個基本元素，并通過dc_id和sku_id組合成unit_id，從而標識更細致的顆粒度

datedc_idsku_idunit_idsale

02018-05-01DC001SKU024DC001_SKU0240

12018-05-01DC001SKU070DC001_SKU0700

22018-05-01DC001SKU078DC001_SKU0789

32018-05-01DC001SKU044DC001_SKU044589

42018-05-01DC001SKU051DC001_SKU05143

52018-05-01DC001SKU063DC001_SKU0632751

62018-05-01DC001SKU099DC001_SKU09996

72018-05-01DC001SKU100DC001_SKU100200

82018-05-01DC001SKU101DC001_SKU10142

92018-05-01DC001SKU102DC001_SKU102102.3數據預處理數據導入及基礎操作

數據類型轉換時間數據在csv等文件中常以字符串形式儲存，對這種格式的數據無法使用一些針對時間戳格式的便捷操作。因此，在讀入帶有時間的數據時，往往需要確認時間數據類型并進行數據類型的轉換sales_data['date']=pd.to_datetime(sales_data['date'])2.3數據預處理數據導入及基礎操作

數據整理當拿到的原始數據的排列并不按照商品、DC和時間順序進行排列時，我們并不知道數據從何時開始、到何時結束，也不知道究竟有多少種商品和DC。所以接下來我們進行對表格的排序、時間范圍的識別、商品種類和DC個數的統計工作表格的排序#排序

sales_data=sales_data.sort_values(

by=['dc_id','sku_id','date']).reset_index(drop=True)

print(sales_data.head())

datedc_idsku_idunit_idsale

02018-02-02DC001SKU001DC001_SKU0017

12018-02-03DC001SKU001DC001_SKU0010

22018-02-04DC001SKU001DC001_SKU0010

32018-02-05DC001SKU001DC001_SKU0012

42018-02-06DC001SKU001DC001_SKU0010時間范圍的識別商品種類及DC個數的統計2.3數據預處理數據導入及基礎操作

#商品種類和DC個數統計

all_sku_list=sales_data['sku_id'].unique().tolist()

print(f'商品個數為:{len(all_sku_list)}')

all_dc_list=sales_data['dc_id'].unique().tolist()

print(f'DC個數為:{len(all_dc_list)}')#時間范圍的識別

print(f'數據起始時間為:{sales_data["date"].min()}')

print(f'數據結束時間為:{sales_data["date"].max()}')數據起始時間為:2018-01-0100:00:00

數據結束時間為:2020-07-3000:00:00

商品個數為:77

DC個數為:182.3數據預處理數據導入及基礎操作

部分數據篩選上述數據表格包含了所擁有的全部信息，但當我們僅需要其中部分信息時，可以將這部分數據存放到一個新的數據表中，便于之后取用unit_id=

'DC001_SKU044'

unit_df=sales_data[sales_data['unit_id']==unit_id]

print(unit_df.head())

datedc_idsku_idunit_idsale

212282018-01-01DC001SKU044DC001_SKU04479

212292018-01-02DC001SKU044DC001_SKU044492

212302018-01-03DC001SKU044DC001_SKU0441261

212312018-01-04DC001SKU044DC001_SKU0441067

212322018-01-05DC001SKU044DC001_SKU0445352.3數據預處理數據導入及基礎操作

描述性統計pandas自帶的數據分析函數describe()可以自動統計數值型數據的均值方差等信息，也可以自動統計非數值型數據的個數、唯一值數量、最高頻率等信息print(unit_df['sale'].describe())count940.000000

mean1199.408511

std913.354405

min0.000000

25%628.000000

50%1106.000000

75%1661.750000

max7777.000000

Name:sale,dtype:float642.3數據預處理數據導入及基礎操作

數據聚合對于非時間類型的類別信息，我們一般可以通過groupby方法，將原始數據按照某一類進行聚合下面的代碼展示了將原始的銷售數據按照商品進行聚合，統計每一類產品的日度銷售情況，并按照日期和商品類別進行排序print(sales_data.groupby(['sku_id','date'])['sale']\

.agg(['sum','mean','median','std'])\

.sort_values(by=['date','sku_id']).reset_index()[:5])2.3數據預處理數據導入及基礎操作

數據聚合由于日度銷量往往波動性較高，有時需要將原始的日銷量數據聚合為一段時間內的總銷量數據，例如周銷量、月銷量等可以通過resample方法實現，下列代碼展示了我們如何將原始數據聚合得到每個unit的周銷量數據weekly_sales_df=sales_data.groupby('unit_id')[['date','sale']]\

.resample('W',on='date').sum().reset_index()

print(weekly_sales_df.head())2.3數據預處理數據導入及基礎操作缺失值處理由于對于真實數據，往往需要考察其完整性，并對缺失值做出一定的處理一般地，我們需要針對每個商品單元考察記錄完整性，即每個商品標識下的數據是否完整覆蓋了數據集時間對于缺失數據，不同的成因對應了不同的解決方式這可能是正常的商業現象，即某天的銷售量為零，對此我們可以用0填充缺失數據；這可能是數據遺失導致的，缺失部分存在銷售，但銷量并未被記錄，此時可以用插值的方式進行填充也可能是因為產品上市和銷售時間不一致，某些新品、促銷品、季節品本身在市時間就較短，此時可以選擇將這些產品與常規品區別開，單獨構建模型，或者通過定性分析的方式單獨進行分析和處理2.3數據預處理數據導入及基礎操作數據輸出在完成上述的基本數據整理工作后，我們可以將得到的新數據表格輸出，在后續進行深入分析和模型構建時，可以直接讀取整理好的結果，從而避免重復的數據預處理工作weekly_sales_df.to_csv(data_dir+

'weekly_sales_data.csv',index=False)2.3數據預處理時序數據可視化導入Python中常用的可視化工具包，matplotlib和seaborn，并統一繪圖格式代碼示例#導入繪圖包

importmatplotlib.pyplotasplt

importseabornassns

#設置畫圖風格

sns.set_theme(style='darkgrid')

#設置畫布大小

plt.rcParams['figure.figsize']=(12.0,6.0)

#設置清晰度

plt.rcParams['figure.dpi']=

500#導入繪圖包

importmatplotlib.pyplotasplt

importseabornassns

2.3數據預處理時序數據可視化（接上頁）#設置畫圖風格

sns.set_theme(style='darkgrid')

importmatplotlib

matplotlib.rc('font',family='SongtiSC')

parameters={'figure.figsize':[12.0,6.0],

'figure.dpi':500,

'axes.labelsize':14,

'xtick.labelsize':14,

'ytick.labelsize':14,

'legend.fontsize':14,

'font.weight':'bold'

}

plt.rcParams.update(parameters)2.3數據預處理時序數據可視化下列代碼定義了針對某一個產品銷量的季節性繪制函數，通過將不同年份的數據疊加，可以更直觀地看到不同年份間是否存在相似的季節性變動

#季節性繪制函數

defplot_year_seasonality(data,date_col,qty_col):

"""

data:pd.DataFrame格式數據

date_col:data中標識日期的列名,需轉換為日期格式

qty_col:data中標識銷量的列名

"""

#以1-12月設置橫坐標

month_all_list=list(range(1,12

+

1))

plt.xticks(range(len(month_all_list)),month_all_list,rotation=60)

#設置曲線顏色和曲線標識

color=['#787878','#1c79d9','black']

marker=['o','^','s']

#初始化標簽列表

legend_list=[]

2.3數據預處理時序數據可視化color=['#787878','#1c79d9','black']

marker=['o','^','s']

#初始化標簽列表

legend_list=[]

#將傳入數據按照月份聚合

data=data.resample('M',on=date_col).agg({qty_col:'sum'})

#提取數據涉及年份，并逐年遍歷

plot_year=list(set(data.index.year))

color_idx=

0

foryearinplot_year:

#提取當年所有數據，以月份為橫坐標，月總銷量為縱坐標繪制折線圖

data_year=data[data.index.year==year]

data_inds=data_year.index.month

data_val=data_year[qty_col]

plt.plot(data_inds-

1,data_val,c=color[color_idx],

marker=marker[color_idx])

#為曲線添加標簽

legend_list.append(year)

#動態改變下一次繪圖曲線顏色

color_idx+=

1

2.3數據預處理時序數據可視化

#展示繪制圖像

plt.xlabel('月份')

plt.ylabel('銷量')

plt.legend(legend_list,loc='upperright')

plt.show()

2.3數據預處理時序數據可視化分別考察’SKU044’和’SKU008’在’DC001’的銷量季節性?？梢园l現，’SKU044’的銷量季節性并不明顯，相比之下’SKU008’則呈現出一定的季節規律：在3-7月呈現出銷量上升的共同趨勢，在9月～次年2月呈現出銷量下降的共同趨勢2.3數據預處理時序數據可視化下列代碼定義了針對某一個產品銷量的趨勢性繪制函數，利用產品的月銷量數據繪制折線圖，來直觀地考察其上升下降趨勢#趨勢性繪制函數

defplot_trend(data,date_col,qty_col):

data=data.resample('M',on=date_col).sum()

plt.plot(data['sale'],color='#1c79d9',marker='o')

plt.xlabel('時間')

plt.ylabel('銷量')

plt.show()

2.3數據預處理時序數據可視化選擇’DC001_SKU044’和’DC017_SKU002’兩個unit分別繪制其趨勢圖像?？梢钥吹?，前者呈現出一定的上升趨勢，同時疊加了一定的季節性波動；而后者則出現了較為明顯的下降趨勢。在針對這兩個unit進行需求預測時，為了獲得更好的預測效果，需要將產品的趨勢性納入模型考慮2.3數據預處理時序數據可視化當一個商品在多個地區進行銷售時，其在各區域之間的需求可能存在一定的相關性，這是由商品屬性決定的共同趨勢。同時，不同產品之間也可能存在相互關聯不管是在模型構建還是特征構建中，我們都可以將相關因素納入其中，以提升預測模型的準確性這一操作在經過格式轉換的數據上可以通過調用內置函數corr()實現，同時，也可以通過繪制熱力圖，更直觀地觀察相關性sku_df=sales_data[sales_data[‘sku_id’]==

‘SKU076’]

corr=sku_df[[‘date’,‘dc_id’,‘sale‘]]\

.pivot(index='date',columns='dc_id',values='sale').corr()

sns.heatmap(corr,annot=True,cmap="Blues")<AxesSubplot:xlabel='dc_id',ylabel='dc_id'>

2.3數據預處理時序數據可視化熱力圖解讀可以通過圖中標注出的相關系數考察變量相關性，相關系數的取值范圍在-1到1之間，越接近于1表明正相關程度越高，越接近于-1表明負相關程度越高，越接近于0表明相關性越弱也可以結合圖例，從色塊顏色直觀地考察變量之間的相關性2.3數據預處理異常點檢驗

2.3數據預處理異常點檢驗下面兩個函數定義了上述兩種異常值檢測方法，輸入DataFrame格式的數據并指定要考察異常值的數據列名，函數便會返回離群值的標識defsigma_outlier(data,qty_col):

upper=data[qty_col].mean()+

3

*data[qty_col].std()

lower=data[qty_col].mean()-

3

*data[qty_col].std()

returndata[

(data[qty_col]>upper)|(data[qty_col]<lower)].index.tolist()defquantile_outlier(data,qty_col):

q1=data[qty_col].quantile(q=0.25)

q3=data[qty_col].quantile(q=0.75)

upper=q3+

1.5

*(q3-q1)

lower=q1-

1.5

*(q3-q1)

returndata[(data[qty_col]>upper)|(data[qty_col]<lower)].index

2.4評價預測模型的方法在預測任務中，兩組不同的預測結果何者更好，往往并不直觀。在某些數據集上一種方法可能比另一種方法的預測結果更接近真實值，而在其他數據集上則可能反之因此，我們需要一套規范化的指標對預測結果進行評估，以便選擇最合適的預測方法2.4評價預測模型的方法樣本內外誤差數據泄露（dataleakage）：其實質是進行模型評估的數據同時被用于模型訓練，為了更好地預測未來值，而不是最優擬合歷史數據，所以需要避免數據泄露模型在已知數據集上表現優于在未知數據集上表現的現象，被稱為過擬合，這樣的模型被稱為泛化能力差。這種現象是預測中常常會遇到的，需要盡力避免進行模型訓練和模型評估所用的數據集應該沒有交集，以避免因數據泄露導致的對模型的不公允評估與選擇2.4評價預測模型的方法樣本內外誤差在實踐中，將數據集按照一定的比例切分為訓練集（即用于模型訓練的數據）和測試集（即用于模型評估的數據）在和時間順序無關的數據集中，一般隨機選擇數據劃分為這兩個集合，以保證兩個集合具有類似的數據結構在時序場景下，需要遵循數據集的時間屬性，將靠前的數據劃分為訓練集，靠后的數據劃分為測試集。在進行切分的時候，也不再按照比例的方式決定兩個數據集合的大小，而是根據一個指定的時間，對數據進行前后切分模型在訓練集上給出的預測結果與真實數據的偏差程度稱為樣本內誤差，在測試集上的偏差則稱為樣本外誤差一般而言，以最小化樣本內誤差為目標進行各預測模型的訓練，以最小化樣本外誤差為標準進行模型間的選擇2.4評價預測模型的方法預測評估指標

2.4評價預測模型的方法預測評估指標

importnumpyasnp

defcal_mse(y_true,y_pred):

returnnp.mean([(y_pred[i]-y_true[i])**

2

foriinrange(len(y_pred))])

defcal_rmse(y_true,y_pred):

returnnp.sqrt(cal_mse(y_true,y_pred))

defcal_mae(y_true,y_pred):

returnnp.mean([abs(y_pred[i]-y_true[i])foriinrange(len(y_pred))])

2.4評價預測模型的方法預測評估指標defcal_mape(y_true,y_pred):

returnnp.mean([abs(y_pred[i]-y_true[i])/max(1.0,y_true[i])foriin

range(len(y_pred))])

defcal_wmape(y_true,y_pred):

sum_error=sum([abs(y_pred[i]-y_true[i])foriinrange(len(y_pred))])

sum_true=max(1.0,sum(y_true))

returnsum_error/sum_true

在需求預測業務中，我們常將wMAPE視為模型的錯誤率，而將1-wMAPE作為模型準確性的一個代表defpred_evaluate(y_true,y_pred):

returnround(max(0,1

-cal_wmape(y_true,y_pred)),3)2.4評價預測模型的方法誤差分析常見的誤差來源數據性質數據不完整，所收集的數據中遺漏了某些對結果具有影響力的重要因素數據的聚合程度較低時，數據中的隨機性影響較大用于數據評估的測試集性質也會影響評估效果需求性質產品需求本身具有特殊性質或可預測性較低產品正處于需求拐點期模型方法模型的錯誤選擇或是模型的不當調試現代庫存管理：模型、算法與Python實現第3章時間序列方法數據集本節采用’DC001_SKU044’的周度銷售數據作為示例導入上一章定義的預測相關的基本函數讀入并篩選數據繪制該unit的周銷量趨勢數據集選擇2020年4月1日為分割時間點，在此之前的數據為訓練集，共117條數據，在此之后的數據為測試集，共18條數據，測試集比例為13.3%#數據切分

unit_data=unit_data.set_index('date')['sale']

test_start_date=pd.Timestamp(2020,4,1)

unit_train_data=unit_data[unit_data.index<test_start_date]

unit_test_data=unit_data[unit_data.index>=test_start_date]3.1指數平滑方法移動平均法基本思想：利用過去一段時間數據的均值，作為對未來的預測需要決定的參數：選擇過去多長時間的數據，即預測窗口期代碼示例forwindow_lengthin[2,8]:

unit_fit=unit_data.rolling(window=window_length).mean()

plot_fitness(unit_data,unit_fit)3.1指數平滑方法移動平均法窗口期越短，對數據的變動跟蹤地越緊密細致；窗口期越長，則越能反應數據的平均變動趨勢3.1指數平滑方法簡單指數平滑法

3.1指數平滑方法簡單指數平滑法

3.1指數平滑方法簡單指數平滑法

defsimple_exp_smoothing(data,initialization=None,alpha=0.4):

"""

data:歷史需求序列,索引為時間,值為歷史需求值

initialization:初始化值，默認為原始數據的首個觀察值，也可設定為一段時間均值

alpha:模型參數,默認設置為0.4

"""

ifinitializationis

None:

fitting=[data.values[0]]

else:

fitting=[initialization]

foriinrange(len(data)):

fitting.append(alpha*data.values[i]+(1

-alpha)*fitting[-1])

returnpd.Series(fitting[:-1],index=data.index)3.1指數平滑方法簡單指數平滑法

2.2需求預測方法的分類簡單指數平滑法3.1指數平滑方法簡單指數平滑法

fromstatsmodels.tsa.apiimportSimpleExpSmoothing

defalpha_fitting(data):

model=SimpleExpSmoothing(data).fit()

returnmodel.model.params["smoothing_level"]3.1指數平滑方法簡單指數平滑法利用該函數，得到最小化訓練集上均方誤差的最優參數約為0.9。將最優參數下的模型帶入測試集進行預測，整體預測準確率為71.1%best_alpha=alpha_fitting(unit_train_data)

pred=simple_exp_smoothing(unit_test_data,unit_train_data[-1],best_alpha)

print(f'最優模型參數為$\\alpha$={alpha}')

plot_fitness(unit_test_data,pred,

title='預測準確率為'

+str(pred_evaluate(unit_test_data,pred)))

最優模型參數為$\alpha$=0.93.1指數平滑方法簡單指數平滑法

3.1指數平滑方法Holt-Winters方法針對簡單指數平滑法的不足，Holt和Winters分別對模型做出了改進，形成了Holt-Winters方法主要思路是將數據分解為趨勢項、季節項和水平項（可以對應于微觀因素），對于每一個構成項采用簡單指數平滑法進行預測，再逆向組合得到最終預測值根據分解方式的不同，這一方法有加法形式和乘法形式，選擇哪一種形式主要取決于季節性因素的影響方式。在實際中，我們可以同時嘗試兩種模型，并選擇評估結果更好的模型用于預測3.1指數平滑方法Holt-Winters方法

3.1指數平滑方法Holt-Winters方法通過Python對簡單指數平滑法進行額外的參數設置，便可以實現Holt-Winters模型的構造和使用fromstatsmodels.tsa.holtwintersimportExponentialSmoothing

defholt_winters(data,pred_length,trend='additive',damped_trend=True,

seasonal='additive',seasonal_periods=52):

model=ExponentialSmoothing(data,trend=trend,seasonal=seasonal,

damped_trend=damped_trend,

seasonal_periods=seasonal_periods).fit()

pred=model.forecast(steps=pred_length)

returnpred3.1指數平滑方法Holt-Winters方法加性模型預測效果#加性模型

pred=holt_winters(data=unit_train_data,pred_length=len(unit_test_data),

trend=’additive’,seasonal=’additive’)

plot_fitness(unit_test_data,pred,

title=’預測準確率為’

+str(pred_evaluate(unit_test_data,pred)))3.1指數平滑方法Holt-Winters方法乘性模型預測效果#乘性模型

pred=holt_winters(data=unit_train_data,pred_length=len(unit_test_data),

trend='additive',seasonal='multiplicative')

plot_fitness(unit_test_data,pred,

title='預測準確率為'

+str(pred_evaluate(unit_test_data,pred)))3.1指數平滑方法Holt-Winters方法在這組數據上，更復雜的Holt-Winters方法并沒有提供比簡單指數平滑法更好的預測這主要是因為，該方法是為趨勢性和季節性較強的數據集設計的，而本組數據覆蓋時間較短，沒有明顯的宏觀趨勢，也無法呈現出季節趨勢，整體而言波動較為劇烈，微觀趨勢更為明顯因此，并非越復雜的模型就一定有更好的效果，不同的模型有不同的假設和適用場景，只有根據數據選擇了合適的方法，才能最好地利用數據和模型，得到合理的預測3.2

ARIMA模型

3.2

ARIMA模型模型假設和數據檢驗某一隨機過程是（弱）平穩序列的標準定義為：隨機過程的均值和方差是與時間無關的有限常數，且兩個時間點的協方差只與時間間隔有關ARMA模型要求用于建模的數據是（弱）平穩的，直觀來看表現為序列圍繞某一確定水平上下波動，不呈現明顯的上升或下降趨勢當原始序列非平穩時，可以通過差分進行處理，直觀來看，對于一個帶有升降趨勢或季節周期的時間序列，兩期之間的變動可能是平穩的，即變動序列的均值和方差不隨時間改變，從而可以首先對需求變動進行預測，再反推得到需求預測，這種包含差分過程的ARMA模型被稱為ARIMA模型考慮到模型的實際意義，差分的次數一般不會超過兩次，即差分階數一般小于23.2

ARIMA模型模型假設和數據檢驗

3.2

ARIMA模型模型假設和數據檢驗因此，在面對一個時間序列時，可以按照如下流程判斷ARIMA模型的適用性：3.2

ARIMA模型參數選擇

3.2

ARIMA模型參數選擇實踐中，可以安裝調用Pmdarima庫實現ARIMA模型的自動定階并完成預測，該庫同時支持對差分階數的自動擬合importpmdarimaaspm

defnonseasonal_arima_training(train,test,p,q,d=None):

"""

train:訓練數據

test:測試數據

p:自回歸項最大階數

q:移動平均項最大階數

"""

#自動選擇指定階數內最優模型

model=pm.auto_arima(train,

start_p=0,#p最小值

start_q=0,#q最小值

d=d,

test='adf',#如果未輸入差分值，ADF檢驗確認差分階數d

max_p=p,#p最大值

max_q=q,#q最大值

stepwise=True

#stepwise為False則不進行完全組合遍歷

)3.2

ARIMA模型參數選擇應用該方法時，需要提前給定最大滯后階數，在給定的范圍內選擇AIC最小的模型一般而言，模型滯后階數不應過大，這一方面是由ARIMA模型的性質決定的（例如高階AR模型可以由低階MA模型等價表示），另一方面也是考慮到算力限制（過高階數有時會造成運行時間過長或是求解失?。┰趯嶋H運用中，可以考慮實際意義限制最大滯后階數pred=[]

fornew_obintest.values.tolist():

pred.append(model.predict(n_periods=1)[0])

model.update(new_ob)

returnpred,model

（接上頁）3.2

ARIMA模型參數選擇由于本節數據是周度數據，選擇4階（對應一個月左右）作為最大滯后階數得到的ARIMA模型在測試集上的準確率約為72.5%，相較Holt-Winters方法有所提升，這主要因為ARIMA模型更適合捕捉數據的微觀特征在訓練集上表現最好的是ARIMA(1,1,1)模型，即使用滯后一階的真實數據和預測誤差，對原始數據進行一階差分后的平穩序列進行預測#Non-SeasonalARIMA

pred,model=nonseasonal_arima_training(unit_train_data,unit_test_data,

p=4,q=4,d=1)

pred=pd.Series(pred,index=unit_test_data.index)

print(model.summary())

plot_fitness(unit_test_data,pred,

title='ARIMA,Accuracy='

+str(pred_evaluate(unit_test_data,pred)))

3.2

ARIMA模型參數選擇ARIMA模型的預測效果plot_fitness(unit_test_data,pred,

title='預測準確率為'

+str(pred_evaluate(unit_test_data,pred)))

3.2

ARIMA模型考慮季節性因素的sARIMA模型sARIMA模型是對ARIMA模型無法考慮季節性的不足進行的改進。

Pmdarima庫可以通過修改參數設置實現sARIMA模型sARIMA模型的預測效果plot_fitness(unit_test_data,pred,

title='預測準確率為'

+str(pred_evaluate(unit_test_data,pred)))

3.3

Prophet模型Prophet模型

3.3

Prophet模型Prophet模型

3.3

Prophet模型Prophet模型

3.3

Prophet模型Prophet模型通過AIC信息準則實現自動定階的VARMA模型python代碼如下fromfbprophetimportProphet

unit_train_prophet=unit_train_data.reset_index()\

.rename(columns={'date':'ds','sale':'y'})

unit_test_prophet=unit_test_data.reset_index()\

.rename(columns={'date':'ds','sale':'y'})importos

os.environ["KMP_DUPLICATE_LIB_OK"]="TRUE"

model=Prophet(daily_seasonality=True,weekly_seasonality=True)\

.fit(unit_train_prophet)future=model.make_future_dataframe(periods=len(unit_test_prophet),

freq='W',include_history=False)

pred=model.predict(future)

pred_prophet=pd.Series(pred['yhat'].values,index=pred['ds'])

plot_fitness(unit_test_data,pred_prophet,

title='預測準確率為'

+str(pred_evaluate(unit_test_data,pred_prophet)))3.3

Prophet模型Prophet模型Prophet模型的預測效果3.4

考慮相互作用的VARMA模型

3.4

考慮相互作用的VARMA模型VARMA模型需要輸入一個包含了相同層級全部信息的矩陣形式的歷史數據,準備和模型效果評估過程代碼如下fromstatsmodels.tsa.statespace.varmaximportVARMAX

defvarmax(train,test):

model=VARMAX(train).fit()

pred=pd.DataFrame()

foriinrange(len(test)):

iflen(pred)==

0:

pred=model.forecast()

else:

pred=pd.concat([pred,model.forecast()])

model=model.append(pd.DataFrame(test.iloc[i]).T)

returnpred,modelpred,model=varmax(sku_train_data,sku_test_data)

unit_accuracy_dict={unit:pred_evaluate(sku_test_data[unit],pred[unit])

forunitinunit_list}

print("VARMA模型的平均預測準確率為："

+str(np.mean(list(unit_accuracy_dict.values()))))3.5

考慮層級結構的時間序列預測通過多個產品的需求序列往往呈現出一定的層級結構，即不同產品之間可以按照一定的指標進行聚合，由于這樣的結構存在，我們將同時面對單個產品的需求序列和聚合層級的需求序列因此，我們希望得到的預測需要滿足一定的實際約束，即對單產品的預測經過聚合后應該與對應聚合層級的預測相似，而一般的分步法（即分別預測再聚合）很難保證這一點本節介紹的預測模型HTS正是為了解決這一問題而提出的，通過調用scikit-hts庫可以完成HTS模型的構建和使用3.5

考慮層級結構的時間序列預測HTS模型python代碼實現如下fromhtsimportHTSRegressor

defhts(data,pred_period,model_choice='auto_arima',revision_choice='OLS'):

"""

data:索引是時間序列，列名為DC，值周總銷量

model_choice：'auto_arima'、'holt_winters'、'sarimax'、'prophet'

revision_choice：

'OLS','WLSS','WLSV'

'OLS'-optimalcombinationusingordinaryleastsquares(Default),

'WLSS'-optimalcombinationusingstructurallyweightedleastsquares,

'WLSV'-optimalcombinationusingvarianceweightedleastsquares,

"""

node=dict()

node['total']=data.columns.to_list()

3.5

考慮層級結構的時間序列預測（接上頁）bottom=data

total=data.sum(axis=1).to_frame().rename(columns={0:'total'})

hierarchy_df=bottom.join(total).fillna(0)

model=HTSRegressor(model=model_choice,revision_method=revision_choice,

n_jobs=50).fit(hierarchy_df,node)

pred=model.predict(steps_ahead=len(pred_period))[-len(pred_period):]

pred.index=pred_period

returnpredpred=hts(sku_train_data,sku_test_data.index.to_list(),

model_choice='sarimax',revision_choice='OLS')

unit_accuracy_dict={unit:pred_evaluate(sku_test_data[unit],

pred[unit])forunitin

unit_list}

print(unit_accuracy_dict)

print("平均準確率為："

+str(np.mean(list(unit_accuracy_dict.values()))))現代庫存管理：模型、算法與Python實現第4章機器學習方法數據集本節采用’DC001_SKU044’的周度銷售數據作為示例運行代碼，展示模型效果導入本章使用的庫和數據，并提取想要考察的子數據fromchapter2_forecast_basicimport*

importpandasaspd

importnumpyasnp

importwarnings

warnings.filterwarnings("ignore")data_dir=

'../../data/forecast_data/'

weekly_sales_df=pd.read_csv(data_dir+

'weekly_sales_data.csv')

weekly_sales_df['date']=pd.to_datetime(weekly_sales_df['date'])

weekly_sales_df['dc_id']=weekly_sales_df['unit_id'].str[:5]

weekly_sales_df['sku_id']=weekly_sales_df['unit_id'].str[6:]

sku_id=

'SKU044'

sku_df=weekly_sales_df[weekly_sales_df['sku_id']==sku_id]\

.reset_index(drop=True)

print(sku_df.head())4.1特征工程時間特征機器學習是一種基于特征學習的技術，它無法像時間序列一樣直接對長期的趨勢性、季節性進行建模，因此需要在數據集上人為地構造和加入時間特征數值型和非數值型的數據并不由數據的表現形式決定，而應該由數據之間是否具有順序關系決定對于非數值型的數據，常用的處理方法是生成啞變量或獨熱編碼月份是與類別相關的非數值型的數據，2月的數據和1月的數據之間并沒有嚴格的順序關系，通過這一變量可以提取同一月份數據之間的相同特征，即季節性用月份特征生成啞變量，并將其與原始特征數據拼接。這里僅保留了前11個啞變量（即標記數據屬于1～11月中的哪個月），這是因為12個啞變量是多重共線的，即由前11個特征的結果即可以推出最后一個特征的取值4.1特征工程時間特征#時間特征：年份啞變量

defadd_year_dummy_feature(df):

df_=df.copy()

df_['year']=df_['date'].dt.year

drop_year=

'FEA_year_'

+str(df_['year'].unique().tolist()[-1])

df_=df_.join(pd.get_dummies(df_['year'],prefix='FEA_year')\

.drop(drop_year,axis=

1)).drop('year',axis=

1)

returndf_#時間特征：月份啞變量

defadd_month_dummy_feature(df):

df_=df.copy()

df_['month']=df_['date'].dt.month

drop_month=

'FEA_month_'

+str(df_['month'].unique().tolist()[-1])

df_=df_.join(pd.get_dummies(df_['month'],prefix='FEA_month')\

.drop(drop_month,axis=

1)).drop('month',axis=

1)

returndf_

具體代碼4.1特征工程時間特征

具體代碼結果sku_feature_df=sku_df.copy()

sku_feature_df=add_year_dummy_feature(sku_feature_df)

sku_feature_df=add_month_dummy_feature(sku_feature_df)

print(sku_feature_df.head())unit_iddatesaledc_idsku_idFEA_year_2018FEA_year_2019\

0DC001_SKU0442018-01-073782DC001SKU04410

1DC001_SKU0442018-01-143917DC001SKU04410

2DC001_SKU0442018-01-214008DC001SKU04410

3DC001_SKU0442018-01-284289DC001SKU04410

4DC001_SKU0442018-02-045347DC001SKU044104.1特征工程

FEA_month_1FEA_month_2FEA_month_3FEA_month_4FEA_month_5\

010000

110000

210000

310000

401000

FEA_month_6FEA_month_7FEA_month_8FEA_month_9FEA_month_10\

000000

100000