Optimizzazione dell’inventario con la programmazione dinamica in meno di 100 righe di codice Python

Ottimizzazione dell'inventario con la programmazione dinamica in meno di 100 righe di codice Python

Come padroneggiare le decisioni di inventario utilizzando la programmazione dinamica per massimizzare la redditività

Foto di Christin Hume su Unsplash

Introduzione

L’ottimizzazione dell’inventario è un problema ampio che si presenta in molti settori. La domanda centrale riguarda la ricerca di una risposta a:

“Quanti prodotti ordinare periodicamente per il tuo negozio al fine di massimizzare la redditività?”

Credo che tu sia un gestore di un negozio di biciclette. Ogni giorno, devi contattare il tuo fornitore per ordinare un certo numero di biciclette per il tuo negozio. Se ordini troppe biciclette ogni giorno, spenderai troppi soldi per il mantenimento e la conservazione delle biciclette nel negozio. D’altra parte, se ordini meno biciclette, potresti non avere abbastanza biciclette per soddisfare le richieste dei clienti. Hai bisogno di una “strategia di ordinazione” che bilanci questi due fattori.

Suddivideremo questo problema e troveremo la “strategia di ordinazione” ottimale. La “strategia di ordinazione ottimale” si basa su tre pilastri:

  • Processo di Markov,
  • Rewards Process di Markov
  • e Processi Decisionali di Markov.

In questo e in questo, ho affrontato come strutturare il problema di “ottimizzazione dell’inventario”, il Processo di Markov e il Processo di Ricompensa di Markov.

In questo blog, riassumeremo tutto e collegheremo il Processo di Markov e il Processo di Ricompensa di Markov con il Processo di Decisione di Markov per arrivare alla “strategia di ordinazione” ottimale, il tutto con codici in Python.

Processo di Decisione di Markov:

Il Processo di Decisione di Markov (MDP) è un modello matematico il cui componente chiave sono gli stati, le azioni, le probabilità di transizione e le ricompense.

Nel Processo di Ricompensa di Markov, le decisioni erano fisse e non stavamo considerando tutte le possibili azioni per ciascuno stato. Nel Processo di Decisione di Markov (MDP), dobbiamo costruire una struttura dati che, per ogni stato, consideri tutte le possibili azioni e, seguendo il $(S_t,A_t)$, dobbiamo sapere cosa sarà $(S_{t+1}, R_{t+1})$.

Esempio di Struttura Dati del Processo di Decisione di Markov

Qui ti do un esempio di un dizionario MDP, in cui puoi vedere che per ogni stato devono essere considerate tutte le possibili azioni:

MarkovDecProcessDict = {"Stato Corrente A":{"Azione 1":{("ProssimoS1daAact1", "Ricompensa1"): "PProssimoS1daAact1"                                           ,("ProssimoS2daAact1", "Ricompensa2"): "PProssimoS2daAact1"},                                           "Azione 2":{("ProssimoS1daAact2", "Ricompensa2"): "PProssimoS1daAact2"                                           ,("ProssimoS2daAact2", "Ricompensa2"): "PProssimoS2daAact2"}},                                         "Stato Corrente B":{"Azione 1":{("ProssimoS1daBact1", "Ricompensa1"): "PProssimoS1daBact1"                                           ,("ProssimoS2daBact1", "Ricompensa2"): "PProssimoS2daBact1"},                                           "Azione 2":{("ProssimoS1daBact2", "Ricompensa2"): "PProssimoS1daBact2"                                           ,("ProssimoS2daBact2", "Ricompensa2"): "PProssimoS2daBact2"}}}for stato_corrente, azioni in MarkovDecProcessDict.items():    print(f"Stato Corrente: {stato_corrente}")        for azione, transizioni in azioni.items():        print(f"  Azione: {azione}")                for (prossimo_stato, ricompensa), probabilità in transizioni.items():            print(f"  ({prossimo_stato},{ricompensa}): {probabilità}")
Output del codice Python — fonte: Autore

Processo decisionale di Markov per l’ottimizzazione dell’inventario

Qui stiamo costruendo il dizionario, dove per ogni stato, stiamo considerando tutte le azioni possibili e, per ogni azione, stiamo considerando tutti gli stati e le ricompense possibili. Il nome di questo dizionario è **MDP_dict** e ho scritto il codice di seguito:

from typing import Dict, Tuple# poisson è utilizzato per trovare la pdf della distribuzione di Poisson from scipy.stats import poissonMDP_dict: Dict[tuple, Dict[tuple, tuple]] = {}user_capacity = 2user_poisson_lambda = 1holding_cost = 1missed_customer_cost = 10for alpha in range(user_capacity+1):                                                                                               for beta in range(user_capacity + 1 - alpha):                # Questo è St, lo stato corrente        state = (alpha, beta)                                           # Questo è l'inventario iniziale, il totale di biciclette che si hanno alle 8 del mattino         init_inv = alpha + beta                                         base_reward = -alpha* holding_cost                action = {}        # Considera tutte le azioni possibili        for order in range(user_capacity-init_inv +1):                        #action = {}            dict1 = {}            for i in range(init_inv +1):            # se la domanda iniziale può soddisfare la domanda                if i <= (init_inv-1):                                # probabilità che una specifica domanda possa accadere                    transition_prob = poisson.pmf(i,user_poisson_lambda)                    dict1[((init_inv - i, order), base_reward)] = transition_prob                                     # se la domanda iniziale non può soddisfare la domanda                else:                                    transition_prob = 1- poisson.cdf(init_inv -1, user_poisson_lambda)                                # probabilità di non soddisfare i requisiti                    transition_prob2 = 1- poisson.cdf(init_inv, user_poisson_lambda)                                # ricompensa totale                                    reward = base_reward - missed_customer_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# Costanti

“MDP_dict” è un dizionario Python in cui ogni chiave rappresenta lo “stato corrente” e il valore è “tutte le azioni possibili in quello stato”, corrispondenti allo stato successivo, alla ricompensa immediata e alla probabilità dello stato successivo. Ecco spiegato in modo schematico:

Schema di MDP_dict - Fonte: Autore

Programmazione dinamica

Richard Bellman (anni ’50) ha coniato per la prima volta il termine chiamato Programmazione dinamica (DP). La programmazione dinamica è un metodo per risolvere problemi complessi scomponendoli in sottoproblemi più semplici.

DP si riferisce in generale alle teorie generali del processo decisionale di Markov e agli algoritmi per trovare una politica ottimale in MDP, affidandosi pesantemente all’Equazione di Bellman.

Nel contesto di questo articolo, utilizziamo il termine “programmazione dinamica” con l’obiettivo di trovare la politica ottimale per il problema dell’ottimizzazione dell’inventario.

In generale, esistono due importanti algoritmi di DP:

Algoritmo di iterazione della funzione di valore (Bellman 1957)

– Algoritmo di iterazione di politica (Howard 1960)

In questo articolo, ci concentreremo sull’algoritmo di iterazione di politica e lo implementeremo in Python.

Algoritmo di iterazione di politica per l’ottimizzazione dell’inventario:

L’algoritmo di iterazione di politica è un metodo per trovare la politica ottimale per un determinato MDP. L’algoritmo si basa sulla seguente idea:

– 1) Inizia con una politica iniziale $π_0$.

– 2) Valuta la politica $π_0$ calcolando la funzione di valore di stato $V^{π_0}$.

– 3) Migliora la politica agendo avidamente rispetto a $V^{\pi_0}$ per ottenere una nuova politica $\pi_1$.

L’algoritmo itera sui passaggi precedenti finché la politica converge (la politica smette di cambiare). Andremo attraverso tutte e tre le fasi nelle seguenti sezioni.

Un semplice diagramma per spiegare l'Algoritmo di Iterazione delle Politiche - Fonte: Autore

1) Iniziare con la Politica Iniziale

L’algoritmo di iterazione delle politiche ha bisogno di una politica iniziale (strategia di ordinamento), che può essere qualsiasi politica arbitraria. In questo articolo, useremo la politica che abbiamo trovato nel secondo articolo, che è la seguente:

Politica iniziale:

π_0=C-(α + β)

La politica iniziale prevede che ad ogni fase dell’inventario, ordiniamo la quantità di $C-(α + β)$, che è la differenza tra la capacità dell’inventario e la somma degli articoli iniziali nell’inventario ($α$) e quelli che arriveranno domattina ($β$) .

Ecco il codice per la politica iniziale:

user_capacity_val = 2def policy_0_gen(user_capacity: int):                # Genera una politica iniziale        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) Valutare la politica π_0 calcolando il valore dello stato V^{π_0}.

Si noti che ogni processo decisionale di Markov con una politica fissa porta a un **processo di ricompensa di Markov implicito**. Quindi, come nel precedente articolo, se abbiamo il processo di ricompensa di Markov, possiamo trovare la funzione di valore di ogni stato.

In altre parole:

Riferimento: Autore

La seguente funzione prenderà una **politica fissa** (come input) e restituirà un processo di ricompensa di Markov implicito:

def MRP_using_fixedPolicy(full_MDP, policy):        # Calcola il processo di ricompensa di Markov utilizzando una politica fissa        MRP_policy = {}        for state in full_MDP.keys():            action = policy[state]            MRP_policy[state] = full_MDP[state][action]        return MRP_policy

A titolo di esempio, possiamo fornire la politica iniziale alla seguente funzione e restituirà il processo di ricompensa di Markov implicito:

MRP_p0=MRP_using_fixedPolicy(MDP_dict, policy_0)MRP_p0
Output del codice Python - Fonte: Autore

Avendo il processo di ricompensa di Markov, è abbastanza facile trovare le Ricompense Immediate per ogni stato e anche la funzione di valore dello stato per ogni stato:

def calculate_expected_immediate_rewards(MRP_policy):        # Calcola le ricompense immediate attese dalla politica 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        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):        # Crea la matrice delle probabilità di transizione        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):        # Calcola la funzione di valore dello stato        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
Output del codice Python — Fonte: Autore

3) Migliorare la policy attuando la strategia di massimo guadagno per ottenere una nuova policy $π_1$.

La parte finale dell’ Algoritmo di Iterazione della Policy è migliorare la policy attuando la strategia di massimo guadagno rispetto a $V^{\pi_0}$ per ottenere una nuova policy $\pi_1$.

La formula dell’approccio a massimo guadagno si basa sull’Equazione di Bellman e, in effetti, consiste nel trovare l’azione con il valore di stato più alto per ogni stato.

Equazione di Ottimalità di Bellman — Fonte: Autore

Ecco il codice per l’operazione di massimo guadagno:

def operazione_massimo_guadagno(MDP_full, state_val_policy, old_policy, gamma=0.9):        # Esegui l'operazione di massimo guadagno per migliorare la 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 + 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 = operazione_massimo_guadagno(MDP_full=MDP_dict,                               state_val_policy=state_val_p0,                               old_policy=policy_0)

La nuova policy dopo la prima iterazione dell’algoritmo di iterazione della policy è la seguente:

new_policy

Nell’algoritmo di iterazione della policy, continuiamo a iterare sui tre passaggi sopra elencati fino a quando la policy si stabilizza. In altre parole, continuiamo ad iterare sui tre passaggi sopra elencati fino a quando la policy smette di cambiare. Ecco il codice per l’Algoritmo di Iterazione della Policy:

def iterazione_della_policy():        # Esegui l'iterazione della policy per trovare la policy ottimale        policy = policy_0_gen(user_capacity_val)        while True:            MRP_policy_p0 = MRP_usando_Policy_Fissa(MDP_dict, policy)            expected_immediate_rew = calcola_ricompense_immediate_attese(MRP_policy_p0)            trans_prob_mat_val = crea_matrice_probabilità_di_transizione(MRP_policy_p0)            value_function = calcola_funzione_valore_stato(trans_prob_mat=trans_prob_mat_val,                                                            expected_immediate_rew=expected_immediate_rew,                                                            gamma=0.9)            new_policy = operazione_massimo_guadagno(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

Qual è l’Ordine Ottimale?

Ora possiamo osservare l’esito dell’iterazione della policy e vedere qual è l’ordine ottimale per ogni stato. Ecco il codice per farlo:

for state, order_quantity in opt_policy.items():    print(f"Per lo stato {state}, la quantità d'ordine ottimale è: {order_quantity}")
La policy finale ottimale — Fonte: Autore

Mettere tutto insieme:

Sto mettendo tutti i codici insieme in una classe Python con i suoi metodi.

import numpy as np
from scipy.stats import poisson
import pandas as pd

class MarkovDecisionProcess:
    def __init__(self, capacità_utente, poisson_lambda, costo_mantenimento, costo_esaurimento, gamma):
        # Inizializza la MDP con i parametri dati
        self.capacità_utente = capacità_utente
        self.poisson_lambda = poisson_lambda
        self.costo_mantenimento, self.costo_esaurimento = costo_mantenimento, costo_esaurimento
        self.gamma = gamma
        self.MDP_completa = self.crea_MDP_completa()  # Crea la MDP completa

    def crea_MDP_completa(self):
        # Crea il dizionario della MDP completa
        dizionario_MDP = {}
        for alpha in range(self.capacità_utente + 1):
            for beta in range(self.capacità_utente + 1 - alpha):
                stato, inv_iniziale = (alpha, beta), alpha + beta
                azione = {}
                for ordine in range(self.capacità_utente - inv_iniziale + 1):
                    dict1 = {}
                    for i in range(inv_iniziale + 1):
                        if i <= (inv_iniziale - 1):
                            probabilità_transizione = poisson.pmf(i, self.poisson_lambda)
                            dict1[((inv_iniziale - i, ordine), -alpha * self.costo_mantenimento)] = probabilità_transizione
                        else:
                            probabilità_transizione = 1 - poisson.cdf(inv_iniziale - 1, self.poisson_lambda)
                            probabilità_transizione2 = 1 - poisson.cdf(inv_iniziale, self.poisson_lambda)
                            ricompensa = -alpha * self.costo_mantenimento - self.costo_esaurimento * (
                                        (self.poisson_lambda * probabilità_transizione) - inv_iniziale * probabilità_transizione2)
                            dict1[((0, ordine), ricompensa)] = probabilità_transizione
                    azione[ordine] = dict1
                dizionario_MDP[stato] = azione
        return dizionario_MDP

    def policy_0_gen(self):
        # Genera una politica iniziale
        return {(alpha, beta): self.capacità_utente - (alpha + beta)
                for alpha in range(self.capacità_utente + 1)
                for beta in range(self.capacità_utente + 1 - alpha)}

    def MRP_utilizzando_politica_fissa(self, politica):
        # Crea la MRP utilizzando una politica fissa
        return {stato: self.MDP_completa[stato][azione]
                for stato, azione in politica.items()}

    def calcola_funzione_valore_stato(self, MRP_politica):
        # Calcola le ricompense immediate attese dalla politica MRP
        E_ricompensa_immediata = {}
        for stato_partenza, valore in MRP_politica.items():
            ricompensa_attesa = sum(ricompensa[1] * probabilità for (ricompensa, probabilità) in valore.items())
            E_ricompensa_immediata[stato_partenza] = ricompensa_attesa

        # Crea la matrice delle probabilità di transizione
        stati = list(MRP_politica.keys())
        probabilità_transizione = np.zeros((len(stati), len(stati)))
        df_probabilità_transizione = pd.DataFrame(probabilità_transizione, columns=stati, index=stati)
        for i, stato_partenza in enumerate(stati):
            for j, stato_destinazione in enumerate(stati):
                for (nuovo_stato, ricompensa) in MRP_politica.get(stato_partenza, {}):
                    if nuovo_stato == stato_destinazione:
                        probabilità = MRP_politica[stato_partenza].get((nuovo_stato, ricompensa), 0.0)
                        df_probabilità_transizione.iloc[i, j] = probabilità

        # Calcola la funzione di valore dello stato
        R_exp = np.array(list(E_ricompensa_immediata.values()))
        funzione_valore_vec = np.linalg.solve(np.eye(len(R_exp)) - self.gamma * df_probabilità_transizione, R_exp)
        MarkRevData = pd.DataFrame({'Ricompensa Immediata Attesa': R_exp, 'Funzione Valore': funzione_valore_vec},
                                   index=stati)
        return MarkRevData

    def operazione_avidà(self, MDP_completa, funzione_valore_politica, vecchia_politica):
        # Esegue l'operazione avida per migliorare la politica
        nuova_politica = {}
        for stato in vecchia_politica.keys():
            max_q_valore, migliore_azione = float('-inf'), None
            dizionario_funzione_valore_stati = funzione_valore_politica.to_dict(orient="index")
            for azione in MDP_completa[stato].keys():
                q_valore = 0
                for (stato_successivo, ricompensa_immediata), probabilità in MDP_completa[stato][azione].items():
                    q_valore = q_valore + probabilità * (ricompensa_immediata + self.gamma *
                                                          (dizionario_funzione_valore_stati[stato_successivo][
                                                               "Funzione Valore"]))
                if q_valore > max_q_valore:
                    max_q_valore, migliore_azione = q_valore, azione
            nuova_politica[stato] = migliore_azione
        return nuova_politica

    def policy_iteration(self):
        # Esegue l'iterazione di politica per trovare la politica ottimale
        politica = self.policy_0_gen()
        while True:
            MRP_politica_p0 = self.MRP_utilizzando_politica_fissa(politica)
            funzione_valore = self.calcola_funzione_valore_stato(MRP_politica_p0)
            nuova_politica = self.operazione_avidà(self.MDP_completa, funzione_valore, politica)
            if nuova_politica == politica:
                break
            politica = nuova_politica
        politica_ottimale, funzione_valore_ottimale = nuova_politica, funzione_valore
        return politica_ottimale, funzione_valore_ottimale

Esempio di utilizzo della classe Processo Decisionale Markoviano

# Esempio di utilizzo:user_capacity = 2poisson_lambda = 1.0holding_cost = 1stockout_cost = 10gamma = 0.9MDP_Esempio = ProcessoDecisionaleMarkoviano(user_capacity, poisson_lambda, holding_cost, stockout_cost, gamma)opt_policy, opt_val = MDP_Esempio.policy_iteration()# Stampa la politica ottimaleprint("Politica Ottimale:")for state, order_quantity in opt_policy.items():    print(f"Per lo stato {state}, la quantità di ordine ottimale è: {order_quantity}")# Stampa la funzione di valore ottimaleprint("\nFunzione di Valore Ottimale:")print(opt_val)
Output del codice Python — Fonte: Autore

Riassunto:

  • Ottimizzazione dell’inventario non è un problema di ottimizzazione statico; piuttosto, è un problema di decisione sequenziale. Significa che abbiamo bisogno di una politica “adattiva” in cui una decisione (numero di ordini) dipende dallo stato dell’inventario.
  • In questo articolo, abbiamo trovato la politica di ordinazione “ottimale”, basata sull’Algoritmo di Programmazione Dinamica.
  • Lo scopo era quello di costruire una base su come pensare ai problemi di ottimizzazione dell’inventario, come formularli e come risolverli utilizzando adeguati modelli matematici.

Riferimenti:

  • [1] È possibile leggere questo esempio in modo più approfondito in “Foundation of Reinforcement Learning with Application in Finance”. Tuttavia, ho riscritto i codici Python in questo articolo per renderli più comprensibili.
  • Bellman, Richard., Un Processo Decisionale Markoviano (1957), Journal of Mathematics and Mechanics.
  • Howard, R. A. 1960. Programmazione Dinamica e Processi Markoviani, (1960.), Cambridge, MA: MIT Press.

Grazie per aver letto finora!

Spero che questo articolo abbia fornito un tutorial facilmente comprensibile su come ottimizzare l’inventario con Python.

Se pensi che questo articolo ti abbia aiutato a imparare di più sull’ottimizzazione dell’inventario e sul Processo Markoviano, ti prego di lasciare un 👏 e di seguirmi!