Pythonコードの行数を100行未満で使用した動的プログラミングによる在庫最適化

「在庫最適化のための動的プログラミングによるPythonコードの100行未満の活用」

ダイナミックプログラミングを使った在庫決定のマスタリング方法

Photo by Christin Hume on Unsplash

イントロ

在庫最適化は、多くのドメインにわたって生じる広範な問題です。中心的な問いは次のようになります:

「利益を最大化するために、定期的にストアのために何個の商品を注文すべきか?」

あなたはバイクショップのマネージャーだと思います。毎日、ストアのために自転車の数をサプライヤーに発注する必要があります。毎日多くの自転車を注文すると、ストアでの保管と維持に余計なお金がかかってしまいます。一方、少なく注文すると、顧客の要求に対応する自転車が足りなくなるかもしれません。2つの要素をバランスする「注文戦略」を持つ必要があります。

この問題を分解し、最適な「注文戦略」を見つけます。最適な注文戦略は以下の3つの柱に基づいています:

  • マルコフプロセス
  • マルコフ報酬プロセス
  • およびマルコフ決定過程

私はここで「在庫最適化」問題、マルコフプロセス、およびマルコフ報酬プロセスのフレーム化方法をカバーしました。

このブログでは、すべてをまとめ、マルコフプロセスとマルコフ報酬プロセスをマルコフ決定プロセスにつなげて、最適な「注文戦略」に到達する方法をPythonのコードとともにまとめます。

マルコフ決定プロセス:

マルコフ決定プロセス(MDP)は、状態、アクション、遷移確率、および報酬をキーとする数学モデルです。

マルコフ報酬プロセスでは、決定は固定され、各状態に対してすべての可能なアクションを考慮していませんでした。マルコフ決定プロセス(MDP)では、各状態に対してすべての可能なアクションを考慮するデータ構造を構築する必要があり、$(S_t,A_t)$に従って次の$(S_{t+1}, R_{t+1})$がわかるようにする必要があります。

マルコフ決定プロセスのデータ構造の例

ここでは、MDP辞書の例を示します。各状態について、すべての可能なアクションを考慮する必要があります:

MarkovDecProcessDict = {"現在の状態A":{"アクション1":{("次の状態1_from_A_act1","報酬1"):"次の状態1_from_A_act1の確率",                                                  ("次の状態2_from_A_act1","報酬2"):"次の状態2_from_A_act1の確率"},                                                "アクション2":{("次の状態1_from_A_act2","報酬2"):"次の状態1_from_A_act2の確率",                                                  ("次の状態2_from_A_act2","報酬2"):"次の状態2_from_A_act2の確率"}},                                            "現在の状態B":{"アクション1":{("次の状態1_from_B_act1","報酬1"):"次の状態1_from_B_act1の確率",                                                  ("次の状態2_from_B_act1","報酬2"):"次の状態2_from_B_act1の確率"},                                                "アクション2":{("次の状態1_from_B_act2","報酬2"):"次の状態1_from_B_act2の確率",                                                  ("次の状態2_from_B_act2","報酬2"):"次の状態2_from_B_act2の確率"}}}for current_state, actions in MarkovDecProcessDict.items():    print(f"現在の状態:{current_state}")        for action, transitions in actions.items():        print(f"  アクション:{action}")                for (next_state, reward), probability in transitions.items():            print(f"  ({next_state},{reward}): {probability}")
Pythonコードの出力——ソース:著者

在庫最適化のためのマルコフ決定過程

ここでは、各状態に対して、すべての可能なアクションを考慮し、各アクションに対して、すべての可能な状態と報酬を考慮して辞書を構築しています。この辞書の名前は**MDP_dict**であり、以下にそのコードを記述しました。

from typing import Dict, Tuple# poisson is used to find pdf of Poisson distribution from scipy.stats import poissonMDP_dict: Dict[tuple, Dict[tuple, tuple]] = {}user_capacity = 2user_poisson_lambda = 1holding_cost = 1missedcostumer_cost = 10for alpha in range(user_capacity+1):                                                                                               for beta in range(user_capacity + 1 - alpha):                # This is St, the current state        state = (alpha, beta)                                           # This is initial inventory, total bike you have at 8AM         init_inv = alpha + beta                                         base_reward = -alpha* holding_cost                action = {}        # Consider all possible actions        for order in range(user_capacity-init_inv +1):                        #action = {}            dict1 = {}            for i in range(init_inv +1):            # if initial demand can meet the deman                if i <= (init_inv-1):                                # probality of specifc demand can happen                    transition_prob = poisson.pmf(i,user_poisson_lambda)                    dict1[((init_inv - i, order), base_reward)] = transition_prob                                     # if initial demand can not meet the demand                else:                                    transition_prob = 1- poisson.cdf(init_inv -1, user_poisson_lambda)                                # probability of not meeting the demands                    transition_prob2 = 1- poisson.cdf(init_inv, user_poisson_lambda)                                # total reward                                    reward = base_reward - missedcostumer_cost*((user_poisson_lambda*transition_prob) - \                                                  init_inv*transition_prob2)                                    dict1[((init_inv - i, order),reward)] = transition_prob                    #if state in MDP_dict:            action[order] = dict1        MDP_dict[state]= actionMDP_dict# Constants

“MDP_dict”は、各キーが「現在の状態」を表し、値が「その状態で可能なすべてのアクション」を表すPythonの辞書です。次の状態、即時報酬、および次の状態の確率に対応します。以下に、概略の説明があります:

MDP_dictの概略図 — 出典: 著者

動的計画法

リチャード・ベルマン(1950年代)は、動的計画法(DP)と呼ばれる用語を最初に提唱しました。ダイナミックプログラミングは、複雑な問題をより単純な部分問題に分割して解決するための方法です。

DPは一般的に、マルコフ決定過程の一般的な理論と、MDPで最適ポリシーを見つけるためのアルゴリズムを指します。これはベルマン方程式に大きく依存しています。

この記事の文脈では、在庫最適化問題の最適ポリシーを見つけるために、動的計画法という用語を使用します。

一般的には、次の2つの重要なDPアルゴリズムがあります:

価値関数繰り返しアルゴリズムベルマン1957

– ポリシー繰り返しアルゴリズム(ハワード1960)

この記事では、ポリシー繰り返しアルゴリズムに焦点を当て、Pythonで実装します。

在庫最適化のためのポリシー繰り返しアルゴリズム:

ポリシー繰り返しアルゴリズムは、与えられたMDPの最適ポリシーを見つけるための方法です。アルゴリズムは以下のアイデアに基づいています:

– 1) 初期ポリシー $π_0$ から始めます。

– 2) ポリシー $π_0$ を評価し、状態価値関数 $V^{π_0}$ を計算します。

– 3) $V^{\pi_0}$ に従って利己的に行動して新しいポリシー $\pi_1$ を改善します。

アルゴリズムは、ポリシーが収束するまで(ポリシーが変わらなくなるまで)上記の手順を繰り返します。以下のセクションで、3つのステージすべてを説明します。

ポリシー反復アルゴリズムを説明するための単純な図 — 出典:著者

1) 初期ポリシーで始める

ポリシー反復アルゴリズムでは、初期ポリシー(注文戦略)が必要です。任意のポリシーを使用できます。この記事では、2番目の記事で見つけたポリシーを使用します。

初期ポリシー:

π_0=C-(α + β)

初期ポリシーは、在庫の各段階で在庫の容量と初期アイテム($α$)と明日の朝に入荷するアイテム($β$)の合計の差である$C – (α + β)$の量を注文します。

以下は初期ポリシーのコードです。

user_capacity_val = 2def policy_0_gen(user_capacity: int):                # 初期ポリシーを生成        return {(alpha, beta): user_capacity - (alpha + beta)                 for alpha in range(user_capacity + 1)                 for beta in range(user_capacity + 1 - alpha)}policy_0 = policy_0_gen(user_capacity_val)policy_0

2) ポリシーπ_0を評価し、状態価値V^{π_0}を計算する

固定されたポリシーを持つマルコフ決定過程は、**暗黙のマルコフ報酬過程**につながります。したがって、前の記事のようにマルコフ報酬過程を持っていれば、各状態の価値関数を求めることができます。

つまり:

参考:著者

以下の関数は、**固定されたポリシー**(入力)を取得し、暗黙のマルコフ報酬過程を返します。

def MRP_using_fixedPolicy(full_MDP, policy):        # 固定されたポリシーを使用してマルコフ報酬過程を計算        MRP_policy = {}        for state in full_MDP.keys():            action = policy[state]            MRP_policy[state] = full_MDP[state][action]        return MRP_policy

例えば、以下の関数に初期ポリシーを渡すと、暗黙のマルコフ報酬過程が返されます。

MRP_p0=MRP_using_fixedPolicy(MDP_dict, policy_0)MRP_p0
Pythonコードの出力 — 出典:著者

マルコフ報酬過程があれば、各状態の即時報酬と状態価値関数を簡単に見つけることができます。

def calculate_expected_immediate_rewards(MRP_policy):        # マルコフ報酬過程から即時報酬を計算        E_immediate_R = {}        for from_state, value in MRP_policy.items():            expected_reward = sum(reward[1] * prob for (reward, prob) in value.items())            E_immediate_R[from_state] = expected_reward        return E_immediate_RR_ime_p0 = calculate_expected_immediate_rewards(MRP_p0)R_ime_p0

import numpy as npimport pandas as pddef create_transition_probability_matrix(MRP_policy):        # 遷移確率行列を作成        states = list(MRP_policy.keys())        num_states = len(states)        trans_prob = np.zeros((num_states, num_states))        df_trans_prob = pd.DataFrame(trans_prob, columns=states, index=states)        for i, from_state in enumerate(states):            for j, to_state in enumerate(states):                for (new_state, reward) in MRP_policy.get(from_state, {}):                    if new_state == to_state:                        probability = MRP_policy[from_state].get((new_state, reward), 0.0)                        df_trans_prob.iloc[i, j] = probability                                return df_trans_probdef calculate_state_value_function(trans_prob_mat, expected_immediate_rew, gamma=0.9):        # 状態価値関数を計算        states = list(expected_immediate_rew.keys())        R_exp = np.array(list(expected_immediate_rew.values()))        val_func_vec = np.linalg.solve(np.eye(len(R_exp)) - gamma * trans_prob_mat, R_exp)        MarkRevData = pd.DataFrame({'Expected Immediate Reward': R_exp, 'Value Function': val_func_vec},                                    index=states)        return MarkRevDatatrans_prob_p0 = create_transition_probability_matrix(MRP_p0)state_val_p0 = calculate_state_value_function(trans_prob_mat=trans_prob_p0,                               expected_immediate_rew=R_ime_p0)state_val_p0
Pythonコードによる出力 — ソース: 著者

3) 方針を改善して、新たな方針$π_1$を得るために、$V^{π_0}$に従って利己的に行動する。

方針改善アルゴリズムの最後の部分は、$V^{\pi_0}$に従って利己的に行動して新しい方針$\pi_1$を得ることです。

利己的方程式はベルマン方程式に基づいており、実際には各状態に対して最も高い**状態価値**関数を持つアクションを見つけることについてです。

ベルマン最適方程式 — ソース: 著者

以下は利己的な操作のコードです:

def greedy_operation(MDP_full, state_val_policy, old_policy, gamma=0.9):        # 新しい方針を得るために利己的な操作を実行        new_policy = {}        for state in old_policy.keys():            max_q_value, best_action  = float('-inf'), None            state_val_dict = state_val_policy.to_dict(orient="index")            for action in MDP_full[state].keys():                q_value = 0                for (next_state, immediate_reward), probability in MDP_full[state][action].items():                    q_value = q_value +  probability * (immediate_reward + gamma *                        (state_val_dict[next_state]["Value Function"]))                if q_value > max_q_value:                    max_q_value, best_action = q_value, action            new_policy[state] = best_action        return new_policynew_policy = greedy_operation(MDP_full=MDP_dict,                               state_val_policy=state_val_p0,                               old_policy=policy_0)

ポリシー反復アルゴリズムの最初の反復後の新しい方針は次の通りです:

new_policy

ポリシー反復アルゴリズムでは、方針が収束するまで上記の3つのステップを繰り返します。つまり、方針が変化しなくなるまで上記の3つのステップで反復を行います。以下はポリシー反復アルゴリズムのコードです:

def policy_iteration():        # 最適方針を見つけるためにポリシー反復を実行        policy = policy_0_gen(user_capacity_val)        while True:            MRP_policy_p0 = MRP_using_fixedPolicy(MDP_dict, policy)            expected_immediate_rew = calculate_expected_immediate_rewards(MRP_policy_p0)            trans_prob_mat_val = create_transition_probability_matrix(MRP_policy_p0)            value_function = calculate_state_value_function(trans_prob_mat=trans_prob_mat_val,                                                            expected_immediate_rew=expected_immediate_rew,                                                            gamma=0.9)            new_policy = greedy_operation(MDP_full=MDP_dict,                                           state_val_policy=value_function,                                           old_policy=policy)            if new_policy == policy:                break            policy = new_policy                opt_policy = new_policy        opt_value_func = value_function                return opt_policy, opt_value_func

最適な順序とは何ですか?

さて、ポリシー反復の結果を見て、各状態に対する最適な順序を確認することができます。以下はそのためのコードです:

for state, order_quantity in opt_policy.items():    print(f"状態{state}の最適な注文数量は:{order_quantity}")
最適な最終方針 — ソース: 著者

すべてをまとめる:

Pythonクラス内のすべてのコードを1つにまとめます。

import numpy as npfrom scipy.stats import poissonimport pandas as pdclass MarkovDecisionProcess:    def __init__(self, user_capacity, poisson_lambda, holding_cost, stockout_cost, gamma):        # パラメータを指定してMDP(マルコフ決定過程)を初期化        self.user_capacity = user_capacity        self.poisson_lambda = poisson_lambda        self.holding_cost, self.stockout_cost = holding_cost, stockout_cost        self.gamma = gamma        self.full_MDP = self.create_full_MDP()  # フルMDPを作成    def create_full_MDP(self):        # フルMDP辞書を作成        MDP_dict = {}        for alpha in range(self.user_capacity + 1):            for beta in range(self.user_capacity + 1 - alpha):                state, init_inv = (alpha, beta), alpha + beta                 action = {}                for order in range(self.user_capacity - init_inv + 1):                    dict1 = {}                    for i in range(init_inv + 1):                        if i <= (init_inv - 1):                            transition_prob = poisson.pmf(i, self.poisson_lambda)                            dict1[((init_inv - i, order), -alpha * self.holding_cost)] = transition_prob                        else:                            transition_prob = 1 - poisson.cdf(init_inv - 1, self.poisson_lambda)                            transition_prob2 = 1 - poisson.cdf(init_inv, self.poisson_lambda)                            reward = -alpha * self.holding_cost - self.stockout_cost * (                                (self.poisson_lambda * transition_prob) - init_inv * transition_prob2)                            dict1[((0, order), reward)] = transition_prob                    action[order] = dict1                MDP_dict[state] = action        return MDP_dict    def policy_0_gen(self):        # 初期ポリシーを生成        return {(alpha, beta): self.user_capacity - (alpha + beta)                 for alpha in range(self.user_capacity + 1)                 for beta in range(self.user_capacity + 1 - alpha)}    def MRP_using_fixedPolicy(self, policy):        # 固定ポリシーを使用してMRPを作成        return {state: self.full_MDP[state][action]                 for state, action in policy.items()}        def calculate_state_value_function(self, MRP_policy):        # MRPポリシーから期待即時報酬を計算        E_immediate_R = {}        for from_state, value in MRP_policy.items():            expected_reward = sum(reward[1] * prob for (reward, prob) in value.items())            E_immediate_R[from_state] = expected_reward        # 遷移確率行列を作成        states = list(MRP_policy.keys())        trans_prob = np.zeros((len(states), len(states)))        df_trans_prob = pd.DataFrame(trans_prob, columns=states, index=states)        for i, from_state in enumerate(states):            for j, to_state in enumerate(states):                for (new_state, reward) in MRP_policy.get(from_state, {}):                    if new_state == to_state:                        probability = MRP_policy[from_state].get((new_state, reward), 0.0)                        df_trans_prob.iloc[i, j] = probability        # 状態価値関数を計算        R_exp = np.array(list(E_immediate_R.values()))        val_func_vec = np.linalg.solve(np.eye(len(R_exp)) - self.gamma * df_trans_prob, R_exp)        MarkRevData = pd.DataFrame({'期待即時報酬': R_exp, '価値関数': val_func_vec}, index=states)        return MarkRevData    def greedy_operation(self, MDP_full, state_val_policy, old_policy):        # ポリシーを改善するために貪欲操作を実行        new_policy = {}        for state in old_policy.keys():            max_q_value, best_action  = float('-inf'), None            state_val_dict = state_val_policy.to_dict(orient="index")            for action in MDP_full[state].keys():                q_value = 0                for (next_state, immediate_reward), probability in MDP_full[state][action].items():                    q_value = q_value +  probability * (immediate_reward + self.gamma *                        (state_val_dict[next_state]["価値関数"]))                if q_value > max_q_value:                    max_q_value, best_action = q_value, action            new_policy[state] = best_action        return new_policy    def policy_iteration(self):        # 最適ポリシーを見つけるためにポリシー繰り返しを実行        policy = self.policy_0_gen()        while True:            MRP_policy_p0 = self.MRP_using_fixedPolicy(policy)            value_function = self.calculate_state_value_function(MRP_policy_p0)            new_policy = self.greedy_operation(self.full_MDP, value_function, policy)            if new_policy == policy:                break            policy = new_policy        opt_policy, opt_value_func = new_policy, value_function        return opt_policy, opt_value_func

マルコフ決定過程クラスの使用例

# 使用例:ユーザー容量 = 2ポアソンランバダ = 1.0保持コスト = 1在庫切れコスト = 10ガンマ = 0.9MDP_Example = MarkovDecisionProcess(user_capacity, poisson_lambda, holding_cost, stockout_cost, gamma)opt_policy, opt_val = MDP_Example.policy_iteration()# 最適ポリシーを表示print("最適ポリシー:")for state, order_quantity in opt_policy.items():    print(f"状態 {state} の最適注文数量は:{order_quantity}")# 最適な価値関数を表示print("\n最適な価値関数:")print(opt_val)
Pythonコードの出力-出典:著者

要約:

  • 在庫最適化は静的な最適化問題ではありません。代わりに、順次の意思決定問題です。つまり、注文数の決定は在庫の状態に依存する「適応型」のポリシーが必要です。
  • この記事では、「最適な」注文ポリシーを、ダイナミックプログラミングアルゴリズムに基づいて見つけました。
  • 目標は、在庫最適化の問題をどのように考え、適切な数学モデリングを使用して解決するかについての基礎を築くことでした。

参考文献:

ここまで読んでいただきありがとうございます!

この記事がPythonを使った在庫最適化の方法をわかりやすく紹介できたと思います。

もし、この記事が在庫最適化やマルコフプロセスについて学ぶ上で役立ったと思ったら、👏を押してフォローしてください!

We will continue to update VoAGI; if you have any questions or suggestions, please contact us!

Share:

Was this article helpful?

93 out of 132 found this helpful

Discover more