Partie 2 — REINFORCE (Policy Gradient)

Imports et configurations¶

from collections.abc import Callable

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

Paramètres du problème¶

Nous conservons les mêmes paramètres que dans la partie précédente.

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

# Paramètres de la distribution de demande (géométrique tronquée)
q = 0.1

Refactorisation de l’environnement¶

Nous encapsulons la logique du magasin de vélos dans une classe propre. Contrairement au code de base, nous vectorisons la génération de la demande pour plus d’efficacité si besoin, et nous typons explicitement les méthodes.

class BikeStoreEnv:
    """Environnement de gestion de stock de vélos.

    State space: 0..M (stock en début de semaine)
    Action space: 0..M (quantité commandée)
    """

    def __init__(
        self,
        capacity: int = M,
        q_param: float = q,
        h_cost: float = h,
        c_cost: float = c,
        c0_cost: float = c0,
        p_price: float = p,
    ):
        self.capacity = capacity
        self.h = h_cost
        self.c = c_cost
        self.c0 = c0_cost
        self.p = p_price
        # Pré-calcul de la distribution de demande
        self.p_demand = np.array(
            [q_param * (1 - q_param) ** m for m in range(capacity + 1)]
        )
        # Le reste de la probabilité est concentré sur M (demande >= M)
        self.p_demand[capacity] += 1 - np.sum(self.p_demand)

        self.state = 0
        self.rng = np.random.default_rng(42)

    def reset(self, initial_state: int = 10, seed: int | None = None) -> int:
        """Réinitialise l'environnement."""
        if seed is not None:
            self.rng = np.random.default_rng(seed)
        self.state = initial_state
        return self.state

    def step(self, action: int) -> tuple[int, float, bool, dict]:
        """Exécute une action (commande) et retourne le nouvel état et la récompense.

        Returns:
            next_state: int
            reward: float
            done: bool (toujours False car horizon infini simulé)
            info: dict

        """
        # Contrainte physique : on ne peut pas commander plus que la place disponible
        # Note: L'agent doit apprendre à ne pas commander trop, mais l'environnement
        # peut aussi clipper ou interdire l'action. Ici on clippe pour être robuste.
        # Cependant, pour REINFORCE, on masquera les actions invalides dans la politique.
        effective_action = min(action, self.capacity - self.state)

        # Génération de la demande
        demand = self.rng.choice(len(self.p_demand), p=self.p_demand)

        # Calcul du stock après livraison (immédiate)
        stock_after_order = self.state + effective_action

        # Ventes réalisées
        sales = min(stock_after_order, demand)

        # Coûts et Revenus
        # Coût de commande : proportionnel + fixe si commande > 0
        order_cost = self.c * effective_action
        if effective_action > 0:
            order_cost += self.c0

        # Coût de stockage (sur le stock restant de la semaine précédente, ici self.state)
        maintenance_cost = self.h * self.state

        revenue = self.p * sales

        reward = revenue - order_cost - maintenance_cost

        # Transition d'état
        self.state = max(0, stock_after_order - demand)

        return self.state, reward, False, {"demand": demand, "sales": sales}

Implémentation de la politique Softmax¶

Nous utilisons une politique paramétrée par une matrice $\theta$ de taille $(M+1) \times (M+1)$ . $\theta_{s,a}$ représente la préférence pour l’action $a$ dans l’état $s$ .

La probabilité de choisir l’action $a$ dans l’état $s$ est donnée par :

\pi_\theta(a|s) = \frac{e^{\theta_{s,a}}}{\sum_{a' \in \mathcal{A}(s)} e^{\theta_{s,a'}}}

(1)

où $\mathcal{A}(s)$ est l’ensemble des actions valides (telles que $s + a \le M$ ).

class SoftmaxPolicy:
    def __init__(
        self,
        num_states: int,
        num_actions: int,
        rng: np.random.Generator | None = None,
    ):
        self.num_states = num_states
        self.num_actions = num_actions
        # Initialisation aléatoire faible pour briser la symétrie, mais proche de l'uniforme
        self.theta = np.zeros((num_states, num_actions))
        self.rng = rng if rng is not None else np.random.default_rng()

    def get_valid_actions_mask(self, state: int) -> np.ndarray:
        """Retourne un masque booléen des actions valides (s + a <= M)."""
        mask = np.zeros(self.num_actions, dtype=bool)
        # On ne peut pas commander plus que la capacité restante
        max_order = self.num_states - 1 - state  # num_states = M + 1
        mask[: max_order + 1] = True
        return mask

    def probs(self, state: int) -> np.ndarray:
        """Calcule les probabilités d'action pour un état donné."""
        preferences = self.theta[state]
        mask = self.get_valid_actions_mask(state)

        # Pour éviter les problèmes numériques (overflow), on soustrait le max des préférences valides
        valid_preferences = preferences[mask]
        if len(valid_preferences) == 0:
            # Cas théoriquement impossible ici sauf si state > M
            return np.zeros(self.num_actions)

        exp_pref = np.exp(valid_preferences - np.max(valid_preferences))
        sum_exp = np.sum(exp_pref)

        probabilities = np.zeros(self.num_actions)
        probabilities[mask] = exp_pref / sum_exp

        return probabilities

    def sample_action(self, state: int) -> int:
        """Échantillonne une action selon la politique courante."""
        probabilities = self.probs(state)
        return self.rng.choice(self.num_actions, p=probabilities)

    def compute_gradient(self, state: int, action: int) -> np.ndarray:
        """Calcule le gradient du log-vraisemblance par rapport aux paramètres theta[state].

        nabla_theta log pi(a|s) = 1_a - pi(a|s) pour les actions valides.
        Le gradient est de même taille que theta[state] (vecteur de taille num_actions).
        """
        probabilities = self.probs(state)
        grad = -probabilities  # - pi(a'|s)
        grad[action] += 1.0  # 1 - pi(a|s)

        # Le gradient est nul pour les actions invalides (déjà 0 dans probabilities)
        # mais on s'assure que le masque est respecté si jamais
        mask = self.get_valid_actions_mask(state)
        grad[~mask] = 0.0

        return grad

Algorithme REINFORCE¶

L’algorithme se déroule comme suit :

Générer une ou plusieurs trajectoires avec la politique courante.
Pour chaque pas de temps $t$ , calculer le retour $G_t = \sum_{k=t}^T \gamma^{k-t} r_k$ .
Mettre à jour les paramètres $\theta$ : $\theta \leftarrow \theta + \alpha G_t \nabla_\theta \log \pi(a_t|s_t)$ .

def generate_episode(env: BikeStoreEnv, policy: SoftmaxPolicy, horizon: int) -> dict:
    """Simule un épisode complet."""
    states = []
    actions = []
    rewards = []

    state = env.reset()
    for _ in range(horizon):
        action = policy.sample_action(state)
        next_state, reward, _, _ = env.step(action)

        states.append(state)
        actions.append(action)
        rewards.append(reward)

        state = next_state

    return {
        "states": np.array(states),
        "actions": np.array(actions),
        "rewards": np.array(rewards),
    }


def reinforce(
    env: BikeStoreEnv,
    num_iterations: int = 1000,
    episodes_per_batch: int = 10,
    horizon: int = 100,
    alpha: float = 0.01,
    gamma: float = 0.97,
    seed: int = 42,
) -> tuple[SoftmaxPolicy, list[float]]:

    rng = np.random.default_rng(seed)
    policy = SoftmaxPolicy(env.capacity + 1, env.capacity + 1, rng=rng)

    avg_rewards_history = []

    for i in range(num_iterations):
        grad_accumulator = np.zeros_like(policy.theta)
        batch_rewards = 0.0

        for _ in range(episodes_per_batch):
            episode = generate_episode(env, policy, horizon)
            rewards = episode["rewards"]
            states = episode["states"]
            actions = episode["actions"]

            batch_rewards += np.sum(rewards)  # Somme non actualisée pour le monitoring

            # Calcul des retours G_t (Backward pass)
            T = len(rewards)
            returns = np.zeros(T)
            G = 0
            for t in reversed(range(T)):
                G = rewards[t] + gamma * G
                returns[t] = G

            # Accumulation du gradient
            for t in range(T):
                s_t = states[t]
                a_t = actions[t]
                G_t = returns[t]

                # Baseline optionnelle (non demandée explicitement mais aide la stabilité)
                # Ici on utilise REINFORCE vanilla comme demandé

                local_grad = policy.compute_gradient(s_t, a_t)
                grad_accumulator[s_t] += alpha * G_t * local_grad

        # Mise à jour des poids (moyenne sur le batch pour stabiliser)
        # Note: alpha est déjà appliqué ci-dessus. On divise par le nombre d'épisodes ?
        # Standard REINFORCE: theta += alpha * mean(gradients)
        policy.theta += grad_accumulator / episodes_per_batch

        avg_rewards_history.append(batch_rewards / episodes_per_batch)

        if (i + 1) % 100 == 0:
            print(
                f"Iteration {i + 1}/{num_iterations}, Avg Reward: {avg_rewards_history[-1]:.2f}"
            )

    return policy, avg_rewards_history

Expérimentation 1 : Recherche d’hyperparamètres¶

Nous testons l’algorithme sur 1000 itérations avec un petit batch size pour commencer. L’horizon $T$ est choisi pour que $\gamma^T$ soit négligeable. Pour $\gamma=0.97$ , $T=200$ donne $0.97^{200} \approx 0.002$ .

env = BikeStoreEnv()
horizon = 200

# Entraînement
policy_reinforce, history = reinforce(
    env,
    num_iterations=1000,
    episodes_per_batch=5,
    horizon=horizon,
    alpha=0.005,  # Learning rate prudent
    gamma=gamma,
)

Iteration 100/1000, Avg Reward: 227.34
Iteration 200/1000, Avg Reward: 178.72
Iteration 300/1000, Avg Reward: 179.14
Iteration 400/1000, Avg Reward: 197.88
Iteration 500/1000, Avg Reward: 163.72
Iteration 600/1000, Avg Reward: 227.46
Iteration 700/1000, Avg Reward: 185.36
Iteration 800/1000, Avg Reward: 207.50
Iteration 900/1000, Avg Reward: 206.28
Iteration 1000/1000, Avg Reward: 191.02

# Visualisation de l'apprentissage
df_learning = pl.DataFrame(
    {"Iteration": np.arange(len(history)), "Average Reward": history}
)

alt.Chart(df_learning).mark_line().encode(
    x="Iteration", y="Average Reward", tooltip=["Iteration", "Average Reward"]
).properties(title="Convergence de REINFORCE (Moyenne mobile)").transform_window(
    rolling_mean="mean(Average Reward)", frame=[-50, 0]
).encode(y="rolling_mean:Q")

Expérimentation 2 : Influence du nombre de trajectoires (Batch Size)¶

Nous comparons la stabilité et la performance pour différents batch sizes : 5, 20, 50.

“100, 1000, 5000” trajectoires comme suggéré dans la consigne semble excessif pour chaque pas de gradient étant donné la simplicité du problème, cela ralentirait trop l’exécution.

Nous allons tester des ordres de grandeur raisonnables ou interpréter la consigne comme “nombre total de trajectoires simulées” ou ajuster en conséquence. Nous allons tester des batch sizes plus conséquents comme 5, 20, 50 pour voir l’effet de réduction de variance.

batch_sizes = [5, 20, 50]
results_by_batch = {}

for bs in batch_sizes:
    print(f"Training with batch size {bs}...")
    # On ajuste le nombre d'itérations pour avoir un nombre comparable d'épisodes total ou juste observer la variance
    # Ici on garde le même nombre d'itérations pour voir la vitesse de convergence par update
    _, hist = reinforce(
        env,
        num_iterations=500,  # Moins d'itérations car chaque itération est plus longue
        episodes_per_batch=bs,
        horizon=horizon,
        alpha=0.01,  # On peut augmenter un peu alpha car le gradient est moins bruité
        gamma=gamma,
    )
    results_by_batch[f"Batch {bs}"] = hist

Training with batch size 5...
Iteration 100/500, Avg Reward: 257.62
Iteration 200/500, Avg Reward: 238.02
Iteration 300/500, Avg Reward: 258.66
Iteration 400/500, Avg Reward: 280.36
Iteration 500/500, Avg Reward: 298.78
Training with batch size 20...
Iteration 100/500, Avg Reward: 270.71
Iteration 200/500, Avg Reward: 278.86
Iteration 300/500, Avg Reward: 269.79
Iteration 400/500, Avg Reward: 281.19
Iteration 500/500, Avg Reward: 282.14
Training with batch size 50...
Iteration 100/500, Avg Reward: 264.25
Iteration 200/500, Avg Reward: 272.71
Iteration 300/500, Avg Reward: 287.22
Iteration 400/500, Avg Reward: 293.58
Iteration 500/500, Avg Reward: 285.78

# Comparaison des courbes d'apprentissage
data_list = []
for name, hist in results_by_batch.items():
    for i, val in enumerate(hist):
        data_list.append({"Batch Size": name, "Iteration": i, "Reward": val})

df_batch_comp = pl.DataFrame(data_list)

alt.Chart(df_batch_comp).mark_line(opacity=0.5).encode(
    x="Iteration",
    y="Reward",
    color="Batch Size",
    tooltip=["Batch Size", "Iteration", "Reward"],
).properties(title="Effet de la taille du batch sur la convergence")

Comparaison avec la politique optimale (Bellman)¶

Pour valider nos résultats, nous ré-implémentons rapidement l’itération de la valeur (Value Iteration) pour obtenir la vraie politique optimale $\pi^*$ .

def get_mdp_matrices(env: BikeStoreEnv) -> tuple[np.ndarray, np.ndarray]:
    """Extraction des matrices de transition et récompense pour DP."""
    rewards = np.zeros((env.capacity + 1, env.capacity + 1))
    transitions = np.zeros((env.capacity + 1, env.capacity + 1, env.capacity + 1))

    for s in range(env.capacity + 1):
        for a in range(env.capacity + 1):
            if s + a > env.capacity:
                continue  # Action invalide

            # Coût immédiat déterministe
            cost = env.c * a + (env.c0 if a > 0 else 0) + env.h * s
            stock_after = s + a

            exp_sales = 0
            for d, prob in enumerate(env.p_demand):
                # Transition
                s_next = max(0, stock_after - d)
                transitions[s, a, s_next] += prob

                # Ventes
                sales = min(stock_after, d)
                exp_sales += prob * env.p * sales

            rewards[s, a] = exp_sales - cost

    return rewards, transitions


def solve_bellman(
    env: BikeStoreEnv, gamma: float = 0.97
) -> tuple[np.ndarray, np.ndarray]:
    R, P = get_mdp_matrices(env)
    num_states = env.capacity + 1
    V = np.zeros(num_states)

    for _ in range(1000):
        V_prev = V.copy()
        # Q[s, a] = R[s,a] + gamma * sum(P[s,a,s'] * V[s'])
        Q = R + gamma * np.dot(P, V)

        # Masquage actions invalides
        for s in range(num_states):
            Q[s, (env.capacity - s + 1) :] = -np.inf

        V = np.max(Q, axis=1)
        if np.max(np.abs(V - V_prev)) < 1e-6:
            break

    # Politique déterministe optimale
    Q = R + gamma * np.dot(P, V)
    for s in range(num_states):
        Q[s, (env.capacity - s + 1) :] = -np.inf
    Pi = np.argmax(Q, axis=1)

    return Pi, V


pi_star, v_star = solve_bellman(env, gamma)

Visualisation de la politique apprise vs Optimale¶

Nous comparons la probabilité d’action de REINFORCE avec l’action déterministe de Bellman.

# Récupération de la meilleure politique apprise (celle avec le grand batch par exemple)
# On ré-entraîne une bonne fois pour la figure finale
final_policy, _ = reinforce(env, num_iterations=1000, episodes_per_batch=50, alpha=0.01)

# Construction de la Heatmap de la politique REINFORCE
heatmap_data = []
for s in range(env.capacity + 1):
    probs = final_policy.probs(s)
    for a, prob in enumerate(probs):
        heatmap_data.append({"State": s, "Action": a, "Probability": prob})

df_heatmap = pl.DataFrame(heatmap_data)

chart_heatmap = (
    alt.Chart(df_heatmap)
    .mark_rect()
    .encode(
        x=alt.X("State:O", title="Stock actuel"),
        y=alt.Y("Action:O", title="Commande"),
        color=alt.Color("Probability:Q", scale=alt.Scale(scheme="viridis")),
        tooltip=["State", "Action", "Probability"],
    )
    .properties(title="Politique apprise par REINFORCE (Probabilités)")
)

# Superposition de la politique optimale (points rouges)
df_star = pl.DataFrame({"State": np.arange(len(pi_star)), "Action": pi_star})
chart_star = (
    alt.Chart(df_star)
    .mark_point(color="red", filled=True)
    .encode(x="State:O", y="Action:O")
)

chart_heatmap + chart_star

Iteration 100/1000, Avg Reward: 124.65
Iteration 200/1000, Avg Reward: 133.09
Iteration 300/1000, Avg Reward: 141.45
Iteration 400/1000, Avg Reward: 132.02
Iteration 500/1000, Avg Reward: 138.59
Iteration 600/1000, Avg Reward: 147.03
Iteration 700/1000, Avg Reward: 135.76
Iteration 800/1000, Avg Reward: 148.85
Iteration 900/1000, Avg Reward: 151.48
Iteration 1000/1000, Avg Reward: 146.73

Comparaison finale des performances¶

Nous simulons les politiques sur de nouvelles trajectoires pour comparer les gains cumulés moyens.

def evaluate_policy(
    env: BikeStoreEnv,
    policy_fn: Callable[[int], int],
    num_episodes: int = 100,
    horizon: int = 200,
    seed: int = 123,
) -> np.ndarray:
    rng = np.random.default_rng(seed)
    total_rewards = np.zeros(horizon)

    for _ in range(num_episodes):
        # On peut passer le seed au reset si on veut contrôler chaque épisode,
        # mais ici on veut juste de la variabilité entre épisodes
        # donc on laisse env utiliser son RNG interne ou on pourrait le reseed
        # Mais env.reset() reseed SI on passe un seed.
        state = env.reset(seed=int(rng.integers(0, 10000)))
        cum_reward = 0
        traj_rewards = []
        for t in range(horizon):
            action = policy_fn(state)
            next_state, reward, _, _ = env.step(action)
            cum_reward += (gamma**t) * reward
            traj_rewards.append(cum_reward)
            state = next_state
        total_rewards += np.array(traj_rewards)

    return total_rewards / num_episodes


# Wrappers pour les politiques
def policy_reinforce_fn(s: int) -> int:
    return final_policy.sample_action(s)


def policy_optimal_fn(s: int) -> int:
    return pi_star[s]


def policy_random_fn(s: int) -> int:
    valid_actions = env.capacity - s + 1
    # Utilisation d'un RNG local serait mieux mais pour une simple fonction lambda-like:
    return int(np.random.default_rng().integers(0, valid_actions))


# Évaluation
res_reinforce = evaluate_policy(env, policy_reinforce_fn)
res_optimal = evaluate_policy(env, policy_optimal_fn)
res_random = evaluate_policy(env, policy_random_fn)

df_final = pl.DataFrame(
    {
        "Semaine": np.arange(horizon),
        "REINFORCE": res_reinforce,
        "Optimale (Bellman)": res_optimal,
        "Aléatoire": res_random,
    }
).unpivot(index="Semaine", variable_name="Politique", value_name="Gain Cumulé")
df_final

alt.Chart(df_final).mark_line().encode(
    x="Semaine",
    y="Gain Cumulé",
    color="Politique",
    tooltip=["Semaine", "Gain Cumulé", "Politique"],
).properties(title="Comparaison des performances : REINFORCE vs Bellman")

Analyse Critique¶

Convergence : La méthode REINFORCE parvient à améliorer la politique au fil des itérations, mais présente une variance significative. L’augmentation du nombre de trajectoires par batch (batch size) aide à stabiliser l’estimation du gradient et fluidifie la courbe d’apprentissage.
Qualité de la politique : La politique apprise se rapproche de la structure de la politique optimale (visible sur la heatmap avec les points rouges), bien qu’elle reste stochastique (probabiliste) contrairement à la solution déterministe de la programmation dynamique.
Comparaison Modèle vs Sans-Modèle :
- L’approche Bellman (Modèle connu) est exacte et converge très rapidement vers la solution optimale précise. Elle nécessite cependant de connaître parfaitement la dynamique du système ( $P(s'|s,a)$ ).
- L’approche REINFORCE (Sans modèle) est plus générale et ne nécessite que des interactions avec l’environnement. Elle est cependant beaucoup plus lente à converger (sample inefficiency) et la solution finale est souvent une approximation bruitée de l’optimum.