附录D - 算法对比实验¶

import numpy as np
import matplotlib.pyplot as plt
import time
from collections import defaultdict
import random

# ==================== 环境定义 ====================

class GridWorld:
    """Grid World环境"""

    def __init__(self, size=5):
        self.size = size
        self.start_pos = (0, 0)                     # 起点位置
        self.goal_pos = (size-1, size-1)             # 终点位置（右下角）
        self.obstacles = {(1, 1), (2, 2), (3, 3)}   # 障碍物坐标集合

        # 动作映射：0=上, 1=下, 2=左, 3=右
        self.action_effects = {
            0: (-1, 0), 1: (1, 0), 2: (0, -1), 3: (0, 1)
        }

    def reset(self):
        """重置环境，智能体回到起点"""
        self.agent_pos = self.start_pos
        return self._get_state(), {}

    def _get_state(self):
        """将二维坐标转换为一维状态编号"""
        return self.agent_pos[0] * self.size + self.agent_pos[1]

    def step(self, action):
        """执行动作，返回(新状态, 奖励, 是否结束, _, _)"""
        x, y = self.agent_pos
        dx, dy = self.action_effects[action]
        # 计算新位置，clip确保不越界
        new_x = np.clip(x + dx, 0, self.size - 1)
        new_y = np.clip(y + dy, 0, self.size - 1)
        new_pos = (new_x, new_y)

        # 根据新位置计算奖励
        if new_pos in self.obstacles:
            reward = -5.0       # 碰到障碍物：大惩罚
            done = False
        elif new_pos == self.goal_pos:
            reward = 10.0       # 到达终点：大奖励
            done = True
        else:
            reward = -0.1       # 普通移动：小惩罚（鼓励尽快到达）
            done = False

        self.agent_pos = new_pos
        return self._get_state(), reward, done, False, {}

    def get_state_space(self):
        return self.size * self.size

    def get_action_space(self):
        return 4

# ==================== 算法实现 ====================

class DynamicProgramming:
    """动态规划（值迭代）"""

    def __init__(self, env, gamma=0.9):
        self.env = env
        self.gamma = gamma
        self.n_states = env.get_state_space()
        self.n_actions = env.get_action_space()
        self.V = np.zeros(self.n_states)
        self.policy = np.zeros(self.n_states, dtype=int)
        self.history = []

    def solve(self, theta=1e-8, max_iter=1000):
        """值迭代求解：反复更新状态值直到收敛"""
        start_time = time.time()

        for iteration in range(max_iter):
            delta = 0  # 记录本轮最大值变化

            # 遍历所有状态，更新价值函数
            for s in range(self.n_states):
                v = self.V[s]  # 保存旧值用于计算变化量

                # 计算所有动作的Q值，选最大值作为V(s)
                q_values = []
                for a in range(self.n_actions):
                    # 模拟执行动作：设置位置后调用step
                    self.env.agent_pos = (s // self.env.size, s % self.env.size)
                    next_s, r, done, _, _ = self.env.step(a)
                    q = r + (0 if done else self.gamma * self.V[next_s])
                    q_values.append(q)

                self.V[s] = max(q_values)  # 贝尔曼最优方程更新
                delta = max(delta, abs(v - self.V[s]))

            self.history.append({
                'iteration': iteration,
                'max_delta': delta,
                'time': time.time() - start_time
            })

            # 值变化小于阈值，认为已收敛
            if delta < theta:
                break

        # 从最优值函数中提取最优策略
        for s in range(self.n_states):
            q_values = []
            for a in range(self.n_actions):
                self.env.agent_pos = (s // self.env.size, s % self.env.size)
                next_s, r, done, _, _ = self.env.step(a)
                q = r + (0 if done else self.gamma * self.V[next_s])
                q_values.append(q)
            self.policy[s] = np.argmax(q_values)  # 选Q值最大的动作

        return self.V, self.policy

class MonteCarlo:
    """蒙特卡洛方法"""

    def __init__(self, env, gamma=0.9, epsilon=0.1):
        self.env = env
        self.gamma = gamma
        self.epsilon = epsilon
        self.n_states = env.get_state_space()
        self.n_actions = env.get_action_space()
        self.Q = np.zeros((self.n_states, self.n_actions))
        self.returns = defaultdict(list)  # defaultdict访问不存在的键时返回默认值
        self.history = []

    def select_action(self, state):
        """ε-贪婪策略"""
        if random.random() < self.epsilon:
            return random.randint(0, self.n_actions - 1)
        return np.argmax(self.Q[state])

    def generate_episode(self):
        """生成一个episode"""
        episode = []
        state, _ = self.env.reset()
        done = False

        while not done:
            action = self.select_action(state)
            next_state, reward, done, _, _ = self.env.step(action)
            episode.append((state, action, reward))
            state = next_state

        return episode

    def train(self, n_episodes=10000):
        """训练：通过大量episode采样来估计Q值"""
        start_time = time.time()

        for episode_num in range(n_episodes):
            episode = self.generate_episode()  # 按当前策略生成完整轨迹
            G = 0       # 累计回报
            visited = set()  # 记录已访问的(状态,动作)对

            # 从后往前计算每个时刻的累计回报 G_t
            for t in range(len(episode) - 1, -1, -1):
                state, action, reward = episode[t]
                G = self.gamma * G + reward  # G_t = r_t + γ * G_{t+1}

                # 首次访问MC：只使用每个(s,a)对在episode中第一次出现的回报
                if (state, action) not in visited:
                    visited.add((state, action))
                    self.returns[(state, action)].append(G)
                    self.Q[state, action] = np.mean(self.returns[(state, action)])  # 取历史回报的平均值

            if (episode_num + 1) % 100 == 0:
                avg_reward = self.evaluate()
                self.history.append({
                    'episode': episode_num + 1,
                    'avg_reward': avg_reward,
                    'time': time.time() - start_time
                })

        return self.Q

    def evaluate(self, n_episodes=10):
        """评估当前策略"""
        total_rewards = []
        for _ in range(n_episodes):
            state, _ = self.env.reset()
            done = False
            episode_reward = 0

            while not done:
                action = np.argmax(self.Q[state])
                next_state, reward, done, _, _ = self.env.step(action)
                episode_reward += reward
                state = next_state

            total_rewards.append(episode_reward)

        return np.mean(total_rewards)

class TDLearning:
    """时序差分学习（Q-Learning）"""

    def __init__(self, env, alpha=0.1, gamma=0.9, epsilon=0.1):
        self.env = env
        self.alpha = alpha
        self.gamma = gamma
        self.epsilon = epsilon
        self.n_states = env.get_state_space()
        self.n_actions = env.get_action_space()
        self.Q = np.zeros((self.n_states, self.n_actions))
        self.history = []

    def select_action(self, state):
        """ε-贪婪策略"""
        if random.random() < self.epsilon:
            return random.randint(0, self.n_actions - 1)
        return np.argmax(self.Q[state])

    def train(self, n_episodes=10000):
        """训练：每步即时更新Q值（无需等待episode结束）"""
        start_time = time.time()

        for episode_num in range(n_episodes):
            state, _ = self.env.reset()
            done = False
            episode_reward = 0

            while not done:
                action = self.select_action(state)  # ε-贪婪选择动作
                next_state, reward, done, _, _ = self.env.step(action)
                episode_reward += reward

                # Q-Learning更新（off-policy：用max Q选择目标值）
                if done:
                    target = reward  # 终止状态没有未来奖励
                else:
                    target = reward + self.gamma * np.max(self.Q[next_state])  # TD目标

                # 沿TD误差方向更新Q值
                self.Q[state, action] += self.alpha * (target - self.Q[state, action])
                state = next_state

            if (episode_num + 1) % 100 == 0:
                avg_reward = self.evaluate()
                self.history.append({
                    'episode': episode_num + 1,
                    'avg_reward': avg_reward,
                    'time': time.time() - start_time
                })

            # 衰减探索率
            self.epsilon = max(0.01, self.epsilon * 0.995)

        return self.Q

    def evaluate(self, n_episodes=10):
        """评估当前策略"""
        total_rewards = []
        for _ in range(n_episodes):
            state, _ = self.env.reset()
            done = False
            episode_reward = 0

            while not done:
                action = np.argmax(self.Q[state])
                next_state, reward, done, _, _ = self.env.step(action)
                episode_reward += reward
                state = next_state

            total_rewards.append(episode_reward)

        return np.mean(total_rewards)

# ==================== 实验运行 ====================

def run_comparison():
    """运行对比实验"""
    env = GridWorld(size=5)

    print("=" * 60)
    print("算法对比实验")
    print("=" * 60)

    # 1. 动态规划
    print("\n1. 运行动态规划（值迭代）...")
    dp = DynamicProgramming(env, gamma=0.9)
    dp.solve()
    print(f"   收敛迭代次数: {len(dp.history)}")
    print(f"   最终max_delta: {dp.history[-1]['max_delta']:.6f}")  # [-1]负索引取最后元素
    print(f"   总时间: {dp.history[-1]['time']:.3f}秒")

    # 2. 蒙特卡洛
    print("\n2. 运行蒙特卡洛方法...")
    mc = MonteCarlo(env, gamma=0.9, epsilon=0.1)
    mc.train(n_episodes=5000)
    print(f"   训练episode数: {len(mc.history) * 100}")
    print(f"   最终平均奖励: {mc.history[-1]['avg_reward']:.2f}")
    print(f"   总时间: {mc.history[-1]['time']:.3f}秒")

    # 3. 时序差分
    print("\n3. 运行Q-Learning...")
    td = TDLearning(env, alpha=0.1, gamma=0.9, epsilon=0.1)
    td.train(n_episodes=5000)
    print(f"   训练episode数: {len(td.history) * 100}")
    print(f"   最终平均奖励: {td.history[-1]['avg_reward']:.2f}")
    print(f"   总时间: {td.history[-1]['time']:.3f}秒")

    # 绘制对比图
    plot_comparison(dp, mc, td)

    return dp, mc, td

def plot_comparison(dp, mc, td):
    """绘制对比图：四个子图分别展示学习曲线、收敛速度、值函数和策略"""
    fig, axes = plt.subplots(2, 2, figsize=(14, 10))

    # 子图1：MC与TD的学习曲线对比
    ax1 = axes[0, 0]
    mc_episodes = [h['episode'] for h in mc.history]
    mc_rewards = [h['avg_reward'] for h in mc.history]
    td_episodes = [h['episode'] for h in td.history]
    td_rewards = [h['avg_reward'] for h in td.history]

    ax1.plot(mc_episodes, mc_rewards, label='Monte Carlo', alpha=0.7)
    ax1.plot(td_episodes, td_rewards, label='Q-Learning', alpha=0.7)
    ax1.axhline(y=8.0, color='r', linestyle='--', label='Optimal (~8.0)')
    ax1.set_xlabel('Episode')
    ax1.set_ylabel('Average Reward')
    ax1.set_title('Learning Curve Comparison')
    ax1.legend()
    ax1.grid(True, alpha=0.3)

    # 2. 收敛速度对比
    ax2 = axes[0, 1]
    dp_iters = [h['iteration'] for h in dp.history]
    dp_deltas = [h['max_delta'] for h in dp.history]

    ax2.semilogy(dp_iters, dp_deltas, label='Value Iteration')
    ax2.set_xlabel('Iteration')
    ax2.set_ylabel('Max Delta (log scale)')
    ax2.set_title('DP Convergence')
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    # 子图3：DP求解的值函数热力图（每个格子的期望回报）
    ax3 = axes[1, 0]
    V_matrix = dp.V.reshape(5, 5)  # 重塑张量形状
    im = ax3.imshow(V_matrix, cmap='viridis')  # 用颜色深浅表示值大小
    ax3.set_title('DP Value Function')
    for i in range(5):
        for j in range(5):
            text = ax3.text(j, i, f'{V_matrix[i, j]:.1f}',
                           ha="center", va="center", color="w")
    plt.colorbar(im, ax=ax3)

    # 子图4：TD学习得到的最优策略（箭头表示每个位置应走的方向）
    ax4 = axes[1, 1]
    action_arrows = {0: '↑', 1: '↓', 2: '←', 3: '→'}
    policy_matrix = np.argmax(td.Q, axis=1).reshape(5, 5)  # 取每个状态Q值最大的动作

    for i in range(5):
        for j in range(5):
            action = policy_matrix[i, j]
            ax4.text(j, i, action_arrows[action],
                    ha='center', va='center', fontsize=20)

    ax4.set_xlim(-0.5, 4.5)
    ax4.set_ylim(4.5, -0.5)
    ax4.set_aspect('equal')
    ax4.set_title('TD Policy')
    ax4.grid(True)

    plt.tight_layout()
    plt.savefig('algorithm_comparison.png', dpi=150)
    plt.show()

# 运行实验
if __name__ == "__main__":
    dp, mc, td = run_comparison()

gammas = [0.5, 0.9, 0.99]
for gamma in gammas:
    dp = DynamicProgramming(env, gamma=gamma)
    V, policy = dp.solve()
    print(f"γ={gamma}: 迭代次数={len(dp.history)}")

epsilons = [0.01, 0.1, 0.3]
for eps in epsilons:
    td = TDLearning(env, epsilon=eps)
    td.train(n_episodes=5000)
    print(f"ε={eps}: 最终奖励={td.history[-1]['avg_reward']:.2f}")

sizes = [3, 5, 7, 10]
for size in sizes:
    env = GridWorld(size=size)
    # 运行实验...