MPC 为什么只执行第一个动作？

模型有误差，越往后的预测越不可靠。只执行第一个动作、到新状态后用真实观测重新规划，相当于每一步都用真实反馈纠正一次，这种“滚动时域”方式对模型误差更鲁棒。

CEM 比随机打靶好在哪里？

随机打靶的候选序列彼此独立，在 25 步 × 连续动作的空间里很难采到好的序列；CEM 每轮把采样分布的均值和方差拟合到得分最高的一批序列上，逐步把搜索集中到好的区域。本站实验中，CEM 的平均回报比随机打靶高约 108。

PETS 为什么在本站实验中没有比确定性模型 + CEM 更好？

Pendulum 的动力学是确定性的、低维的，几百个样本就能学得很准，模型不确定性的作用不大；PETS 每条候选序列要用 10 个粒子传播，计算量约为 CEM 的 7 倍，在同样的时间预算下候选数还更少（200 对 300）。PETS 的优势主要体现在有噪声、高维、数据少时模型容易出错的任务上。

MPC 和学习一个策略网络（如 SAC）相比如何？

MPC 不需要学习策略，样本效率很高（本站实验约 400–600 个真实步数就能立起摆，而 SAC 需要约 4000 步），但每一步都要做大量模型推理，决策延迟高；策略网络决策只需一次前向传播。MBPO 等方法把两者结合：用模型生成数据来训练策略。

模型预测控制：随机打靶、CEM 与 PETS

一句话总结 · 要点速览

MPC 每一步都用模型评估许多候选动作序列，只执行最优序列的第一个动作，然后在新状态重新规划，因此能不断纠正模型误差。
随机打靶直接从均匀分布采样候选序列取最好的；CEM 迭代地把采样分布拟合到得分最高的“精英”序列上，用同样的模型调用次数找到更好的动作。
PETS 用概率神经网络集成表示模型：每个网络输出高斯分布刻画随机性，成员之间的分歧刻画认知不确定性，规划时用粒子在不同成员中传播。
本站 Pendulum 实测（2 个种子，第 3–10 回合平均）：随机打靶 −278.6、CEM −170.3、PETS −196.7；三者都只用 1 个随机回合（200 步）的数据就开始规划，第 2–3 回合即可把摆立起来。

术语不熟悉？查看强化学习术语表。

基于模型的强化学习有两种用模型的方式：一是用模型生成数据去训练策略（如 Dyna、MBPO），二是不训练策略，而是每一步都用模型在线规划。后者的代表就是模型预测控制（Model Predictive Control，MPC）。MPC 在工业控制中已经使用了几十年，那里的模型通常由工程师根据物理规律写出；在强化学习里，我们把模型换成从交互数据中学到的神经网络。本节在 Pendulum 上依次实现三种 MPC：随机打靶、交叉熵方法（CEM）和 PETS。

滚动时域规划

设已有动力学模型 $\hat{f} (s, a) \approx s^{'}$ 和奖励函数 $r (s, a)$ 。在当前状态 $s_{t}$ ，MPC 求解一个长度为 $H$ 的开环优化问题：

max_{a_{t}, \dots, a_{t + H - 1}} \sum_{k = 0}^{H - 1} r ({\hat{s}}_{t + k}, a_{t + k}), {\hat{s}}_{t} = s_{t}, {\hat{s}}_{t + k + 1} = \hat{f} ({\hat{s}}_{t + k}, a_{t + k})

然后只执行第一个动作 $a_{t}$ ，观察真实的下一个状态 $s_{t + 1}$ ，再从 $s_{t + 1}$ 重新求解。这种做法叫滚动时域（receding horizon）：

每一步都用真实观测重新出发，模型的多步预测误差不会一直累积；
规划视界 $H$ 只需要覆盖“足够看清后果”的长度，不必覆盖整个回合；
不需要价值函数，也不需要策略网络——代价是每一步都要做一次优化。

在强化学习的循环中，模型还要不断更新：

text

用随机动作收集少量数据 D
重复（每个回合）：
    用 D 训练动力学模型 f̂
    对回合中的每一步 t：
        用 f̂ 规划，得到动作序列，执行其第一个动作
        把 (s_t, a_t, s_{t+1}) 加入 D

动作序列优化

神经网络模型高度非线性，上面的优化问题一般没有解析解。常用的是无梯度的采样优化：只需要能批量调用模型，实现简单，也很适合 GPU 或向量化计算。

随机打靶（Random Shooting）

从动作范围内均匀采样 $N$ 条长度为 $H$ 的动作序列，用模型分别预测回报，选回报最高的一条。它的问题是维度灾难：动作序列有 $H \times | A |$ 维，均匀采样在高维空间里很难碰到好的序列。

交叉熵方法（CEM）

CEM 把“采样—挑选”做成迭代：维护一个高斯分布 $N (μ, diag (σ^{2}))$ 作为动作序列的采样分布，每轮：

从分布中采样 $N$ 条序列（并裁剪到动作范围）；
用模型评估回报，取最好的 $K$ 条作为精英；
用精英的均值和标准差更新 $μ, σ$ （可以与旧值做平滑）。

迭代几轮后， $μ$ 就是规划结果。另一个常用技巧是热启动：执行完第一个动作后，把 $μ$ 向前平移一步作为下一时刻的初始均值，因为相邻时刻的最优序列高度相似。

用于规划的动力学模型

本节的模型输入 $(s, a)$ 、预测状态增量 $Δ s = s^{'} - s$ （比直接预测 $s^{'}$ 更容易学），输入和输出都做标准化。Pendulum 的奖励可以从状态直接算出，我们把它当作已知函数，只学习转移。

PETS：概率集成与粒子传播

Chua 等（2018）提出的 PETS（Probabilistic Ensembles with Trajectory Sampling）针对两类不确定性：

偶然不确定性（aleatoric）：环境本身的随机性。让每个网络输出高斯分布 $N (μ_{θ} (s, a), Σ_{θ} (s, a))$ ，用负对数似然训练：

L (θ) = \sum_{n} {[μ_{θ} (s_{n}, a_{n}) - Δ s_{n}]}^{⊤} Σ_{θ}^{- 1} [μ_{θ} (s_{n}, a_{n}) - Δ s_{n}] + \log det Σ_{θ} (s_{n}, a_{n})

实现时要把预测的对数方差限制在一个可学习的上下界之间，避免在数据外的区域方差爆炸或塌缩。

认知不确定性（epistemic）：数据不足导致的模型不确定。训练 $B$ 个网络组成集成，每个网络用自助采样（bootstrap）的数据训练；在有数据的地方它们的预测一致，在数据稀少的地方会出现分歧。

规划时用粒子传播不确定性：每条候选动作序列配 $P$ 个粒子，每个粒子从当前状态出发、在每一步从所选网络的输出分布中采样下一个状态。“TS∞”方式让每个粒子在整个视界内固定使用同一个集成成员。候选序列的得分取所有粒子回报的平均，因此“只在某个模型里看起来好”的序列会被惩罚。

实验：Pendulum-v1

设置（三种方法相同的部分）：第 1 个回合用随机动作收集 200 步数据；之后每个回合结束时用全部数据重训模型，共 10 个回合；规划视界 $H = 25$ ，只执行第一个动作。各方法的差别：

方法	模型	优化器	每步模型调用量
随机打靶	1 个确定性网络	300 条均匀随机序列	300 条 × 25 步
CEM	1 个确定性网络	5 轮 × 300 条，精英 30 条，热启动	1500 条 × 25 步
PETS	5 个概率网络集成	同 CEM，但每轮 200 条 × 10 个粒子	10000 条 × 25 步

"""模型预测控制（MPC）：在 Pendulum-v1 上比较三种基于“学到的模型”的规划方法
  1. 随机打靶（Random Shooting）：确定性神经网络模型 + 随机采样动作序列取最好的一条
  2. CEM：确定性模型 + 交叉熵方法迭代优化动作序列
  3. PETS：概率神经网络集成（5 个）+ 粒子轨迹采样（TS∞）+ CEM
奖励函数视为已知（由预测的状态直接计算），只学习状态转移。
用法：python mpc_pendulum.py；DQN_QUICK=1 python mpc_pendulum.py 只跑很短的版本用于检查代码。"""
import os
import sys
import time

import gymnasium as gym
import numpy as np
import torch
import torch.nn as nn

QUICK = os.environ.get("DQN_QUICK") == "1"
torch.set_num_threads(1)
H, POP, ITERS, ELITE = 25, 300, 5, 30          # 规划视界、候选序列数、CEM 迭代次数、精英数
N_ENS, PARTICLES, POP_PETS = 5, 10, 200         # PETS：集成大小、每条候选序列的粒子数、候选数（计算量较大，略减）
EPISODES = 3 if QUICK else 10                   # 第 1 个回合用随机动作收集初始数据
A_MAX = 2.0


def reward_fn(obs, a):
    """Pendulum 的奖励：-(θ² + 0.1·θ̇² + 0.001·u²)，θ 由 cos、sin 恢复。"""
    th = torch.atan2(obs[..., 1], obs[..., 0])
    return -(th ** 2 + 0.1 * obs[..., 2] ** 2 + 0.001 * a[..., 0] ** 2)


class Model(nn.Module):
    """输入 (s, a)，输出 Δs 的均值（概率模型还输出对数方差）。多个成员用一个批量维度并行计算。"""

    def __init__(self, n_members, probabilistic):
        super().__init__()
        self.n, self.prob = n_members, probabilistic
        out = 6 if probabilistic else 3
        def P(*shape):
            return nn.Parameter(torch.randn(n_members, *shape) * (1 / np.sqrt(shape[0])))
        self.W = nn.ParameterList([P(4, 200), P(200, 200), P(200, out)])
        self.b = nn.ParameterList([nn.Parameter(torch.zeros(n_members, 1, d)) for d in (200, 200, out)])
        self.max_logvar = nn.Parameter(torch.full((1, 1, 3), 0.5)); self.min_logvar = nn.Parameter(torch.full((1, 1, 3), -10.0))
        self.mu = self.sd = self.dmu = self.dsd = None

    def set_norm(self, X, D):
        self.mu, self.sd = X.mean(0), X.std(0) + 1e-6
        self.dmu, self.dsd = D.mean(0), D.std(0) + 1e-6

    def forward(self, x):                        # x: [n, B, 4]（已归一化）
        for i, (W, b) in enumerate(zip(self.W, self.b)):
            x = torch.baddbmm(b, x, W)
            if i < 2:
                x = x * torch.sigmoid(x)         # Swish 激活
        if not self.prob:
            return x, None
        mean, logvar = x[..., :3], x[..., 3:]
        logvar = self.max_logvar - nn.functional.softplus(self.max_logvar - logvar)
        logvar = self.min_logvar + nn.functional.softplus(logvar - self.min_logvar)
        return mean, logvar

    def fit(self, X, D, rng, epochs=50, batch=256):
        X, D = torch.as_tensor(X, dtype=torch.float32), torch.as_tensor(D, dtype=torch.float32)
        self.set_norm(X, D); Xn, Dn = (X - self.mu) / self.sd, (D - self.dmu) / self.dsd
        opt = torch.optim.Adam(self.parameters(), lr=1e-3, weight_decay=1e-5)
        N = len(X)
        boot = torch.as_tensor(rng.integers(0, N, (self.n, N)))          # 每个成员各自的自助采样
        for _ in range(epochs):
            perm = boot[:, torch.randperm(N)]
            for i in range(0, N, batch):
                j = perm[:, i:i + batch]
                mean, logvar = self(Xn[j])
                if self.prob:                    # 高斯负对数似然
                    loss = (((mean - Dn[j]) ** 2) * torch.exp(-logvar) + logvar).mean()
                    loss = loss + 0.01 * (self.max_logvar.sum() - self.min_logvar.sum())
                else:
                    loss = ((mean - Dn[j]) ** 2).mean()
                opt.zero_grad(); loss.backward(); opt.step()

    @torch.no_grad()
    def predict(self, s, a, gen):                # s: [n, B, 3], a: [n, B, 1]
        x = (torch.cat([s, a], -1) - self.mu) / self.sd
        mean, logvar = self(x)
        if self.prob:
            mean = mean + torch.exp(0.5 * logvar) * torch.randn(mean.shape, generator=gen)
        return s + mean * self.dsd + self.dmu


@torch.no_grad()
def plan(model, s0, method, init_mean, gen):
    """返回规划出的动作序列（只执行第一个动作），以及下一步用于热启动的序列。"""
    def rollout_return(seqs):                    # seqs: [POP, H]
        P = PARTICLES if model.prob else 1
        n_pop = seqs.shape[0]
        # 粒子 p 固定使用集成成员 p % N_ENS（TS∞）；确定性模型只有 1 个成员
        s = torch.as_tensor(s0, dtype=torch.float32).expand(model.n, n_pop * P // model.n if model.prob else n_pop, 3).clone()
        a_all = seqs.repeat_interleave(P, 0) if model.prob else seqs
        a_all = a_all.view(model.n, -1, H) if model.prob else a_all.view(1, -1, H)
        ret = torch.zeros(a_all.shape[:2])
        for t in range(H):
            a = a_all[..., t:t + 1]
            ret += reward_fn(s, a)
            s = model.predict(s, a, gen)
        return ret.reshape(n_pop, P).mean(1) if model.prob else ret.reshape(n_pop)

    pop = POP_PETS if model.prob else POP
    if method == "rs":
        seqs = (torch.rand(POP, H, generator=gen) * 2 - 1) * A_MAX
        return seqs[rollout_return(seqs).argmax()], None
    mean, std = init_mean.clone(), torch.full((H,), A_MAX)
    for _ in range(ITERS):
        seqs = (mean + std * torch.randn(pop, H, generator=gen)).clamp(-A_MAX, A_MAX)
        elite = seqs[rollout_return(seqs).topk(ELITE).indices]
        mean, std = 0.9 * elite.mean(0) + 0.1 * mean, 0.9 * elite.std(0) + 0.1 * std
    return mean, torch.cat([mean[1:], torch.zeros(1)])


def run(method, seed):
    rng = np.random.default_rng(seed); torch.manual_seed(seed); gen = torch.Generator().manual_seed(seed)
    env = gym.make("Pendulum-v1"); env.action_space.seed(seed)
    model = Model(N_ENS, True) if method == "pets" else Model(1, False)
    X, D, returns = [], [], []
    for ep in range(EPISODES):
        s, _ = env.reset(seed=seed * 100 + ep); G = 0.0; warm = torch.zeros(H); t0 = time.time()
        for t in range(200):
            if ep == 0:
                a = env.action_space.sample()
            else:
                seq, nxt = plan(model, s, method, warm, gen)
                a = np.array([float(seq[0])], np.float32)
                warm = nxt if nxt is not None else warm
            s2, r, term, trunc, _ = env.step(a)
            X.append(np.concatenate([s, a])); D.append(s2 - s); G += r; s = s2
        returns.append(G)
        model.fit(np.array(X), np.array(D), rng, epochs=20 if QUICK else 50)
        print(f"  {method:4s} 种子{seed} 回合{ep}  回报 {G:8.1f}  （{time.time() - t0:.0f}s）", flush=True)
    return returns


if __name__ == "__main__":
    seeds = [0] if QUICK else [0, 1]
    methods = ("rs", "cem", "pets")
    if len(sys.argv) == 3:                        # 也可以单独运行：python mpc_pendulum.py pets 0
        methods, seeds = (sys.argv[1],), [int(sys.argv[2])]
    summary = {}
    for method in methods:
        summary[method] = [run(method, s) for s in seeds]
    print("\n最后 3 个回合的平均回报（每个种子）：")
    for m, rs in summary.items():
        print(f"{m:4s}: " + "  ".join(f"{np.mean(r[-3:]):.1f}" for r in rs))

每个回合的真实回报（2 个种子；两种子的初始状态序列对三种方法相同，因此同一列可以直接比较）：

回合	1（随机）	2	3	4	5	6	7	8	9	10
随机打靶种子0	−1071.9	−127.7	−239.1	−488.2	−512.1	−237.7	−1.4	−123.1	−123.0	−483.6
CEM 种子0	−1071.9	−1.5	−117.2	−227.0	−238.1	−115.3	−0.7	−123.9	−122.4	−242.0
PETS 种子0	−1071.9	−2.6	−121.8	−364.6	−527.7	−116.6	−0.6	−123.4	−123.8	−242.4
随机打靶种子1	−1360.5	−1233.4	−589.5	−119.9	−237.3	−122.6	−566.7	−122.8	−363.8	−126.0
CEM 种子1	−1360.5	−1103.1	−254.1	−120.2	−246.1	−122.1	−304.7	−122.2	−240.2	−127.9
PETS 种子1	−1360.5	−821.3	−240.4	−121.4	−241.4	−122.0	−309.0	−123.3	−242.0	−127.0

汇总（第 3–10 回合、2 个种子的平均）：随机打靶 −278.6，CEM −170.3，PETS −196.7。单个回合的耗时：随机打靶约 1–2 秒，CEM 约 7 秒，PETS 约 50–170 秒（单线程 CPU；与其他实验同时运行时偏慢）。

结果解读：

样本效率极高：只用 1 个随机回合（200 步）训练出的模型，种子 0 的 CEM 和 PETS 在第 2 回合就几乎完美地立起了摆（−1.5、−2.6）；种子 1 到第 3 回合也稳定下来。作为对照，MBPO 页中的无模型 SAC 需要约 4000 个真实步数才达到 −130 左右。
CEM 明显好于随机打靶：同一初始状态下，随机打靶时常需要多摆几次才立起来（例如种子 1 第 7 回合 −566.7，CEM 为 −304.7），偶尔失败（种子 0 第 4、5、10 回合约 −490）。
PETS 在这里没有优于 CEM：两者的大多数回合几乎相同，PETS 在种子 0 的第 4、5 回合反而更差（−364.6、−527.7），平均略低于 CEM。原因见下面的讨论；这一结论只针对 Pendulum 这样确定性、低维、模型容易学的任务。
回报差异主要来自初始状态：Pendulum 的初始角度是随机的，从底部附近出发时，即使是最优控制也要先来回摆动积累能量，回报约 −120 到 −250；所以比较时要看同一列（相同初始状态），而不要和其他页面上不同初始状态的数字直接对比。

为什么 PETS 没有带来提升

Pendulum 的真实动力学是确定性的，偶然不确定性为零；状态只有 3 维，几百个样本就能把模型学得很准，认知不确定性也很快消失。
PETS 每条候选序列要用 10 个粒子评估，计算量约为 CEM 的 7 倍；为了控制时间，我们把每轮候选数从 300 降到了 200，优化本身变得更粗糙。
原论文中 PETS 的优势出现在 HalfCheetah、7 自由度机械臂等更高维、接触更复杂、模型更容易出错的任务上；在那里，单个确定性模型会被规划器“钻空子”，选出在模型中回报虚高的动作。

实践要点

预测增量并标准化：预测 $Δ s$ 并对输入输出做标准化，比直接预测 $s^{'}$ 稳定得多。
视界不宜过长： $H$ 太短看不到后果，太长则模型误差累积、优化维度也更高。Pendulum 用 25 步（约 1.25 秒）就足够。
热启动和分布平滑：CEM 中把上一时刻的解平移作为初始均值，并对均值、方差做指数平滑，能显著减少所需的迭代轮数。
计算预算：MPC 的决策延迟取决于“候选数 × 视界 × 粒子数 × 迭代轮数”次模型调用。实际部署时常用 GPU 批量推理，或者把规划结果蒸馏成一个策略网络。
奖励函数是否已知：本例假设奖励可以从状态算出；若奖励未知，需要同时学习一个奖励模型，误差会进一步累积。

小结

方法	优化方式	模型	优点	缺点
随机打靶	独立采样取最优	任意	实现最简单、可完全并行	高维动作序列空间中效率低
CEM	迭代拟合精英分布	任意	同样预算下找到更好的序列	多轮串行迭代
PETS	CEM + 粒子传播	概率网络集成	显式处理两类不确定性，数据少时更稳健	计算量大

练习题

把规划视界 $H$ 分别设为 5、15、25、50，观察 CEM 的回报和耗时变化，并解释视界过短和过长各自的问题。
关掉 CEM 的热启动（每一步都从均值 0 开始），需要多少轮迭代才能达到相近的效果？
在环境的转移中加入观测噪声或动作噪声（例如对执行的动作加 $N (0, {0.3}^{2})$ 噪声），重新比较 CEM 与 PETS。
在 PETS 中把“TS∞”（粒子固定使用一个成员）改为“TS1”（每一步随机换成员），比较规划效果。
把奖励函数也改成用神经网络学习，观察对样本效率的影响。

参考文献

García, C. E., Prett, D. M., & Morari, M. (1989). Model predictive control: Theory and practice—A survey. Automatica, 25(3), 335–348.
Rubinstein, R. Y. (1999). The cross-entropy method for combinatorial and continuous optimization. Methodology and Computing in Applied Probability, 1, 127–190.
Nagabandi, A., Kahn, G., Fearing, R. S., & Levine, S. (2018). Neural network dynamics for model-based deep reinforcement learning with model-free fine-tuning. ICRA.
Chua, K., Calandra, R., McAllister, R., & Levine, S. (2018). Deep reinforcement learning in a handful of trials using probabilistic dynamics models. NeurIPS. arXiv:1805.12114
Lakshminarayanan, B., Pritzel, A., & Blundell, C. (2017). Simple and scalable predictive uncertainty estimation using deep ensembles. NeurIPS.

下一步

继续阅读：MBPO——不在线规划，而是用模型生成短视界数据来训练 SAC 策略。

模型预测控制：随机打靶、CEM 与 PETS ​

滚动时域规划 ​

动作序列优化 ​

随机打靶（Random Shooting） ​

交叉熵方法（CEM） ​

用于规划的动力学模型 ​

PETS：概率集成与粒子传播 ​

实验：Pendulum-v1 ​

为什么 PETS 没有带来提升 ​

实践要点 ​

小结 ​

练习题 ​

参考文献 ​

下一步 ​