Partie 1 — Programmation dynamique - Reinforcement Learning

Ce notebook implémente la résolution du problème de gestion de stock stochastique par programmation dynamique. On utilise l’itération sur la valeur (Value Iteration) pour trouver la politique optimale.

import altair as alt
import numpy as np
import polars as pl

Paramètres du problème¶

Nous reprenons les constantes définies dans le modèle de base pour assurer la cohérence des comparaisons.

M = 15  # capacité maximale du stock
h = 0.3  # coût de maintenance par unité
c = 0.5  # coût de commande par unité
c0 = 0.3  # coût fixe par commande
p = 1.0  # prix de vente par unité
gamma = 0.97  # facteur d'actualisation

# Distribution de la demande (géométrique tronquée)
q = 0.1
p_dem = np.array([q * (1 - q) ** m for m in range(M + 1)])
p_dem[M] += 1 - np.sum(p_dem)

Implémentation de l’environnement MDP¶

Nous pré-calculons la matrice de transition $P(s'|s, a)$ et les récompenses attendues $R(s, a)$ .

def get_mdp_matrices(m: int, p_demand: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Pre-compute transition and reward matrices."""
    # rewards[s, a]
    rewards = np.zeros((m + 1, m + 1))
    # transitions[s, a, s_next]
    transitions = np.zeros((m + 1, m + 1, m + 1))

    for s in range(m + 1):
        for a in range(m + 1 - s):  # Action space is limited by capacity
            # Immediate costs
            cost = c * a + (c0 if a > 0 else 0) + h * s

            stock_after_order = s + a

            exp_sales = 0.0
            for d, prob in enumerate(p_demand):
                # Next state logic
                s_next = max(0, stock_after_order - d)
                transitions[s, a, s_next] += prob

                # Expected sales profit
                sales = min(stock_after_order, d)
                exp_sales += prob * p * sales

            rewards[s, a] = exp_sales - cost

    return rewards, transitions


r_matrix, p_matrix = get_mdp_matrices(M, p_dem)

Itérations de Bellman (Value Iteration)¶

L’objectif est de trouver le point fixe de l’opérateur de Bellman :

V_{k+1}(s) = \max_{a} \left( R(s, a) + \gamma \sum_{s'} P(s'|s, a) V_k(s') \right)

(1)

def value_iteration(
    rewards: np.ndarray, transitions: np.ndarray, gamma_val: float, tol: float = 1e-8
) -> tuple[np.ndarray, np.ndarray]:
    """Find the optimal value function and policy."""
    num_states = rewards.shape[0]
    v = np.zeros(num_states)
    pi = np.zeros(num_states, dtype=int)

    while True:
        v_old = v.copy()

        # Q-values: rewards[s, a] + gamma * sum_s' P[s, a, s'] * V[s']
        q_values = rewards + gamma_val * np.dot(transitions, v)

        # Mask invalid actions (a > M - s) with -infinity
        for s in range(num_states):
            q_values[s, (M - s + 1) :] = -np.inf

        v = np.max(q_values, axis=1)
        pi = np.argmax(q_values, axis=1)

        if np.max(np.abs(v - v_old)) < tol:
            break

    return v, pi


v_star, pi_star = value_iteration(r_matrix, p_matrix, gamma)

Analyse de la politique optimale¶

Nous transformons les résultats en DataFrame Polars pour la visualisation.

df_results = pl.DataFrame(
    {"Stock": np.arange(M + 1), "Value": v_star, "Policy": pi_star}
)

Visualisation de la fonction de valeur et de la politique¶

La politique optimale semble suivre une structure de type $(s, S)$ : si le stock descend en dessous d’un seuil $s$ , on commande pour atteindre le niveau $S$ .

chart_v = (
    alt.Chart(df_results)
    .mark_line(point=True)
    .encode(
        x=alt.X("Stock:Q", title="Niveau de stock initial"),
        y=alt.Y("Value:Q", title="Valeur optimale V*(s)"),
        tooltip=["Stock", "Value"],
    )
    .properties(title="Fonction de valeur optimale", width=300)
)

chart_pi = (
    alt.Chart(df_results)
    .mark_bar()
    .encode(
        x=alt.X("Stock:O", title="Niveau de stock"),
        y=alt.Y("Policy:Q", title="Quantité à commander"),
        tooltip=["Stock", "Policy"],
    )
    .properties(title="Politique optimale (Quantité de commande)", width=300)
)

chart_v | chart_pi

Comparaison avec les politiques de base¶

Nous simulons des trajectoires pour comparer la politique optimale avec les stratégies naïves.

def simulate_policy(
    policy_array: np.ndarray, steps: int = 200, s_init: int = 10, seed: int = 42
) -> list[float]:
    """Simulate a trajectory given a fixed policy array."""
    rng = np.random.default_rng(seed)
    s = s_init
    cum_reward = 0.0
    rewards_history = []

    for t in range(steps):
        a = policy_array[s]
        # Draw demand from distribution
        d = rng.choice(M + 1, p=p_dem)

        # Calculate reward (same logic as environment)
        stock_after_order = s + a
        reward = -c * a - (c0 if a > 0 else 0) - h * s + p * min(stock_after_order, d)

        cum_reward += (gamma**t) * reward
        rewards_history.append(float(cum_reward))

        # Transition
        s = max(0, stock_after_order - d)

    return rewards_history


# Définition des politiques baselines
pi_const = np.array([min(3, M - s) for s in range(M + 1)])
pi_threshold = np.array([10 - s if s <= 4 else 0 for s in range(M + 1)])

# Simulation
res_opt = simulate_policy(pi_star)
res_const = simulate_policy(pi_const)
res_thresh = simulate_policy(pi_threshold)

df_comp = pl.DataFrame(
    {
        "Semaine": np.arange(len(res_opt)),
        "Optimale": res_opt,
        "Constante (3)": res_const,
        "Seuil (4, 10)": res_thresh,
    }
).unpivot(
    index="Semaine",
    on=["Optimale", "Constante (3)", "Seuil (4, 10)"],
    variable_name="Politique",
    value_name="Gain cumulé",
)

Comparaison des performances¶

Le graphique suivant montre l’évolution du gain cumulé actualisé. La politique optimale surpasse les heuristiques simples.

alt.Chart(df_comp).mark_line().encode(
    x="Semaine:Q",
    y=alt.Y("Gain cumulé:Q", title="Gain cumulé actualisé"),
    color="Politique:N",
    tooltip=["Semaine", "Gain cumulé", "Politique"],
).properties(title="Comparaison des performances des politiques", width=600)