Оптимизация инвентаризации с применением динамического программирования в менее чем 100 строк кода на Python

Динамическое программирование для оптимизации инвентаризации код на Python в 100 строк

Как овладеть принятием решений о запасах с использованием динамического программирования для максимальной прибыльности

Введение

Оптимизация запасов – это широкая проблема, которая возникает во многих областях. Основной вопрос заключается в том, как определить:

«Сколько продуктов следует заказывать для вашего магазина в периодическом порядке для максимизации прибыльности?»

Я думаю, вы управляете веломагазином. Каждый день вам нужно связываться с поставщиком, чтобы заказать определенное количество велосипедов для вашего магазина. Если вы заказываете слишком много велосипедов ежедневно, вы потратите слишком много денег на обслуживание и хранение велосипедов в магазине. С другой стороны, если вы заказываете меньше, у вас может не быть достаточно велосипедов, чтобы удовлетворить потребности покупателей. Вам нужно иметь “стратегию заказа”, которая учитывает два фактора.

Мы разобьем эту проблему и найдем оптимальную “стратегию заказа”. “Оптимальная стратегия заказа” основана на трех принципах:

• Открывайте мир с LangChain общайтесь с ваших документами, чат-бот для перевода, общайтесь с Википедией, генератор синтетических данных
• Топ-15 баз данных векторной графики, которые вы должны попробовать в 2024 году
• Топ-15 векторных баз данных для науки о данных в 2024 году Подробное руководство

• Процесс Маркова,
• Процесс Маркова с вознаграждением
• и Процессы принятия решений Маркова.

Здесь и здесь я рассказывал о том, как сформулировать проблему «оптимизации запасов», о процессе Маркова и процессе Маркова с вознаграждением.

В этом блоге мы собираемся подытожить все и связать процесс Маркова и процесс Маркова с вознаграждением с процессом принятия решений Маркова, чтобы прийти к оптимальной “стратегии заказа”, все это с помощью кодов на Python.

Процесс принятия решений Маркова:

Процесс принятия решений Маркова (MDP) – это математическая модель, основными компонентами которой являются состояния, действия, вероятности перехода и вознаграждения.

В процессе Маркова с вознаграждением решения были фиксированными, и мы не учитывали все возможные действия для каждого состояния. В процессе принятия решений Маркова (MDP) нам необходимо построить структуру данных, в которой для каждого состояния мы должны учесть все возможные действия, и зная $(S_t,A_t)$, нам нужно знать, что будет $(S_{t+1}, R_{t+1})$.

Пример структуры данных процесса принятия решений Маркова

Вот пример словаря MDP, где вы видите, что для каждого состояния необходимо учитывать все возможные действия:

MarkovDecProcessDict = {"Текущее состояние А":{"Действие 1":{("СледующееS1изAдейств1", "Награда1"): "PСледующееS1изAдейств1"                                           ,("СледующееS2изAдейств1", "Награда2"): "PСледующееS2изAдейств1"},                                           "Действие 2":{("СледующееS1изAдейств2", "Награда2"): "PСледующееS1изAдейств2"                                           ,("СледующееS2изAдейств2", "Награда2"): "PСледующееS2изAдейств2"}},                                         "Текущее состояние B":{"Действие 1":{("СледующееS1изBдейств1", "Награда1"): "PСледующееS1изBдейств1"                                           ,("СледующееS2изBдейств1", "Награда2"): "PСледующееS2изBдейств1"},                                           "Действие 2":{("СледующееS1изBдейств2", "Награда2"): "PСледующееS1изBдейств2"                                           ,("СледующееS2изBдейств2", "Награда2"): "PСледующееS2изBдейств2"}}}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}")

Марковский процесс принятия решений для оптимизации инвентаризации

Здесь мы создаем словарь, где для каждого состояния мы рассматриваем все возможные действия, и для каждого действия мы рассматриваем все возможные состояния и вознаграждения. Это название словаря **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, где каждый ключ представляет «текущее состояние», а значение – «все возможные действия в этом состоянии», соответствующие следующему состоянию, немедленной награде и вероятности следующего состояния. Вот объяснение на схеме:

Динамическое программирование

Уже в 1950-х годах Ричард Беллман впервые ввел термин под названием Динамическое программирование (DP). Динамическое программирование – это метод решения сложных проблем путем разбиения их на более простые подпроблемы.

DP обычно относится к общим теориям Марковского процесса принятия решений и алгоритмам для нахождения оптимальной политики в MDP, с большим использованием уравнения Беллмана.

В контексте этой статьи мы используем термин Динамическое программирование с целью нахождения оптимальной политики для проблемы оптимизации запасов.

Обычно существует два основных алгоритма DP:

– Алгоритм итерации функции значения (Беллман 1957)

– Алгоритм итерации политики (Говард 1960)

В этой статье мы сосредоточимся на алгоритме итерации политики, и мы его реализуем на Python.

Алгоритм итерации политики для оптимизации запасов:

Алгоритм итерации политики – это метод нахождения оптимальной политики для данного MDP. Алгоритм основан на следующей идее:

– 1) Начните с начальной политики $π_0$.

– 2) Оцените политику $π_0$, вычислив функцию значения состояния $V^{π_0}$.

– 3) Улучшите политику, действуя жадно относительно $V^{\pi_0}$, чтобы получить новую политику $π_1$.

Алгоритм повторяет вышеописанные шаги до достижения сходимости политики (политика перестает изменяться). Мы будем проходить через все три этапа в следующих разделах.

1) Начните с начальной политики

Алгоритм итерации политики требует начальной политики (стратегии упорядочивания), которая может быть любой произвольной политикой. В этой статье мы будем использовать политику, которую мы нашли во второй статье, которая выглядит следующим образом:

Начальная политика:

π_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

Имея Марковский процесс вознаграждения, довольно легко найти немедленные вознаграждения для каждого состояния и также функцию значения состояния для каждого состояния:

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

3) Улучшаем политику, действуя жадно по отношению к $V^{π_0}$, чтобы получить новую политику $π_1$.

Последний этап алгоритма итерационной политики заключается в улучшении политики, действуя жадно по отношению к $V^{\pi_0}$, чтобы получить новую политику $π_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]["Значение функции"]))                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

В алгоритме итерационной политики мы продолжаем итерации по вышеприведенным трём этапам до сходимости политики. Другими словами, мы продолжаем итерацию по этим трём шагам, пока политика не перестанет меняться. Вот код для алгоритма итерационной политики:

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"For state {state}, the optimal order quantity is: {order_quantity}")

Сборка всего вместе:

Я собираю все коды в одном классе Python со своими методами.

import numpy as npfrom scipy.stats import poissonimport pandas as pdclass MarkovDecisionProcess:    def __init__(self, user_capacity, poisson_lambda, holding_cost, stockout_cost, gamma):        # Инициализируем МДП с заданными параметрами        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()  # Создаем полный МДП    def create_full_MDP(self):        # Создаем полный словарь МДП        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):        # Создаем МРП с использованием фиксированной политики        return {state: self.full_MDP[state][action]                 for state, action in policy.items()}        def calculate_state_value_function(self, 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        # Создаем матрицу вероятностей перехода        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

Пример использования класса Markov Decision Process

# Пример использования:
user_capacity = 2
poisson_lambda = 1.0
holding_cost = 1
stockout_cost = 10
gamma = 0.9
MDP_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)

Резюме:

• Оптимизация запасов — это не статическая задача оптимизации, а последовательная проблема принятия решений. Это означает, что нам нужна «адаптивная» политика, где решение (количество заказов) зависит от состояния запасов.
• В этой статье мы нашли «оптимальную» политику заказов, основанную на алгоритме динамического программирования.
• Цель состояла в создании основы для мышления об оптимизации запасов, способе постановки задач и их решении с использованием соответствующего математического моделирования.

Литература:

• [1] Вы можете прочитать более подробный пример в «Foundation of Reinforcement Learning with Application in Finance». Однако, в этой статье я переписал Python-коды, чтобы их было легче понять.
• Bellman, Richard., Марковский процесс принятия решений (1957), Журнал математики и механики.
• Howard, R. A. 1960. Динамическое программирование и марковские процессы, (1960.), Cambridge, MA: MIT Press.

Спасибо за внимание!

Я надеюсь, что этот материал предоставил простой и понятный учебник по оптимизации запасов с помощью Python.

Если вы считаете, что этот материал помог вам изучить больше о оптимизации запасов и марковских процессах, пожалуйста, поставьте 👏 и подпишитесь на меня!