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

> 来源：https://dqn.cn/advanced/mpc ｜ 作者：DQN.CN 编辑部 ｜ 发布：2026-10-30 ｜ 更新：2026-10-30 ｜ 许可：CC BY-NC-SA 4.0（文字）/ MIT（代码）

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

[基于模型的强化学习](https://dqn.cn/advanced/model-based) 有两种用模型的方式：一是用模型生成数据去训练策略（如 Dyna、[MBPO](https://dqn.cn/advanced/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}), \qquad \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 |\mathcal{A}|$ 维，均匀采样在高维空间里很难碰到好的序列。

### 交叉熵方法（CEM）

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

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

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

### 用于规划的动力学模型

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

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

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

- **偶然不确定性**（aleatoric）：环境本身的随机性。让每个网络输出高斯分布 $\mathcal{N}(\mu_\theta(s, a), \Sigma_\theta(s, a))$，用负对数似然训练：
$$
\mathcal{L}(\theta) = \sum_n \left[\mu_\theta(s_n, a_n) - \Delta s_n\right]^\top \Sigma_\theta^{-1} \left[\mu_\theta(s_n, a_n) - \Delta s_n\right] + \log \det \Sigma_\theta(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 步 |

```python
"""模型预测控制（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 页](https://dqn.cn/advanced/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 自由度机械臂等更高维、接触更复杂、模型更容易出错的任务上；在那里，单个确定性模型会被规划器“钻空子”，选出在模型中回报虚高的动作。

## 实践要点

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

## 小结

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

## 练习题

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

## 参考文献

- 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](https://arxiv.org/abs/1805.12114)
- Lakshminarayanan, B., Pritzel, A., & Blundell, C. (2017). Simple and scalable predictive uncertainty estimation using deep ensembles. *NeurIPS*.

## 下一步

继续阅读：[MBPO](https://dqn.cn/advanced/mbpo)——不在线规划，而是用模型生成短视界数据来训练 SAC 策略。

## 常见问题速查（FAQ）

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

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

**问：CEM 比随机打靶好在哪里？**

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

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

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

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

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

## 如何引用本文

DQN.CN 编辑部. 模型预测控制：随机打靶、CEM 与 PETS[EB/OL]. DQN.CN, 2026-10-30（更新于 2026-10-30）. https://dqn.cn/advanced/mpc.