这里是"够用就好"的核心推导；完整版（含每一步代数细节、频域限制、Riccati 变分推导、摆起的 Lyapunov 证明）见
docs/理论推导.md，配套实验见 docs/教学实验手册.md。
1. 建模：拉格朗日方程
广义坐标 q = (x, θ)，θ 自竖直向上量起。摆杆按均匀细杆处理：
质心距转轴 l = L/2，绕质心转动惯量 I = mL²/12。摆杆质心位置为
(x + lsinθ, lcosθ)。
T = ½(M+m)ẋ² + m l cosθ · ẋθ̇ + ½(I + m l²)θ̇²
V = m g l cosθ
代入 d/dt(∂L/∂q̇) − ∂L/∂q = Q耗散/外力，得到两条耦合方程：
(M+m)ẍ + m l cosθ · θ̈ − m l sinθ · θ̇² = u − b ẋ
(I + m l²)θ̈ + m l cosθ · ẍ − m g l sinθ = − c θ̇
写成矩阵形式并求解加速度（model.js 里就是这么算的）：
⎡ M+m m l cosθ ⎤ ⎡ ẍ ⎤ = ⎡ u − bẋ + m l sinθ·θ̇² ⎤
⎣ m l cosθ I+m l² ⎦ ⎣ θ̈ ⎦ ⎣ m g l sinθ − cθ̇ ⎦
行列式 Δ = (M+m)(I+ml²) − (ml cosθ)² ≥ (M+m)I + m²l²sin²θ > 0，永不奇异 —— 这是"质量矩阵正定"的具体体现。
2. 线性化与不稳定性的量化
在 θ=0、θ̇=0、u=0 处取一阶近似（sinθ≈θ, cosθ≈1, θ̇²≈0），
记 J = I + ml² = mL²/3、D₀ = (M+m)J − (ml)²：
ṡ = A s + B u, s = [x, ẋ, θ, θ̇]ᵀ
A = ⎡ 0 1 0 0 ⎤ B = ⎡ 0 ⎤
⎢ 0 −Jb/D₀ −m²l²g/D₀ mlc/D₀ ⎥ ⎢ J/D₀ ⎥
⎢ 0 0 0 1 ⎥ ⎢ 0 ⎥
⎣ 0 mlb/D₀ (M+m)mgl/D₀ −(M+m)c/D₀ ⎦ ⎣ −ml/D₀ ⎦
无摩擦时特征值可解析求出：λ = 0, 0, ±p，其中
p = √( m g l (M+m) / D₀ ) （默认参数 ⇒ p = 5.59 rad/s）
于是"给控制器多少反应时间"有了硬指标：误差倍增时间 ln2 / p ≈ 0.124 s。
这直接决定了采样率下限——分析面板算出的临界采样周期 Ts_crit ≈ 37 ms 就来自这里
（经验关系 p·Ts_crit ≈ 0.2）。
B 的第四项 −ml/D₀ < 0：向右推小车，摆杆向左转。整个倒立摆的"反直觉"都源自这个负号。
3. 两个绕不开的结构性障碍
(a) 能控但只有一个执行器。能控性矩阵 [B AB A²B A³B] 满秩（秩 4），所以四个状态理论上都能控制；
但只有一个输入 u，任何时刻你都在"用同一只手同时扶摆和拉车"。
(b) u→x 通道是非最小相位的。忽略摩擦时可算出
Θ(s)/U(s) = −ml / (D₀s² − mgl(M+m)) （分子是常数：无有限零点）
X(s)/U(s) = (J s² − mgl) / [ s²(D₀s² − mgl(M+m)) ]
后者的零点位于 ±√(mgl/J) = ±√(3g/2L)，其中一个在右半平面。
计入转轴摩擦 c 后分子变为 J s² + c s − mgl，零点略微不对称：
z = ( −c + √(c² + 4 J m g l) ) / (2J) （默认参数：z = 4.849 rad/s，另一零点 −5.057）
分子只含 J、c、mgl，与小车摩擦 b 无关，也与 M 无关——把 b 从 0 改到 100、
把 M 改动 500 倍，零点变化都在 1e−12 以内（Python 参考实现已跨 3 个数量级数值验证）。
物理含义：这个零点是摆杆自身绕转轴的动力学特征，与小车受多少阻力、有多重都无关。
还有一个只在含摩擦模型里才看得见的零点。保留 b 重做消元可得
Θ(s)/U(s) = m l s / { m²l²s³ − [(M+m)s + b](Js² + cs − mgl) }
分子 ∝ s，即 b ≠ 0 时 u→θ 通道有一个原点零点（DC 增益为 0：常值力产生不了稳态倾角）。
后果之一：小车的漂移在角度信号里完全"看不见"——这就是 §5 单环 PID 失败的频域解释；
后果之二更强，见 §5 的"纯 PD 恒不可稳"。
右半平面零点带来的是数学上的硬限制：
- 位置阶跃响应必然先反向——想往右走，必须先往左推让摆杆向右倒（"欲进先退"）；
- 位置环带宽存在上限，经验判据 ωo ≲ z/2；超过它，"先反向"的那段就会被高增益放大成正反馈。
关键澄清：不稳定极点 p 约束的是内环（u→θ，须 ωn ≳ 2p），
右半平面零点 z 约束的是外环（θ_ref→x，须 ωo ≲ z/2）。两者作用在不同回路上，
所以并不构成"夹逼到空集"的矛盾。但可以证明
p / z = 1 / √( 1 − 3m / (4(M+m)) ) > 1 （无摩擦，均匀细杆）
即 z < p 恒成立，且这个比值只取决于质量比 m/(M+m)，与摆长 L、重力 g 完全无关
（默认参数 p/z = 1.1282，数值验证吻合到 6 位）。由此得到串级结构的理论依据：
内环下限 2p 与外环上限 z/2 之间必然相差约 4.5 倍，两个环被迫时标分离——这正是串级 PID 能成立的原因。
摆越短为什么越难？因为 p 与 z 都按 ∝ √(g/L) 同比例变大
（比值不变）：L 从 0.6 m 减到 0.2 m，p 从 5.59 涨到 9.68 rad/s，内环带宽要求、采样率要求、执行器带宽要求全部跟着抬高，
而这些都是有物理上限的。所以"短摆难控"的本质不是稳定性，而是硬件跟不上。
4. LQR：从代价函数到 Riccati 方程
目标：min J = ∫₀∞(sᵀQs + uᵀRu)dt，约束 ṡ = As + Bu，
Q ≥ 0、R > 0。用"猜一个二次型值函数"的办法最快：设 V(s) = sᵀPs
为从 s 出发的最优代价，则 Hamilton–Jacobi–Bellman 方程给出
minu [ sᵀQs + uᵀRu + 2sᵀP(As + Bu) ] = 0
对 u 求偏导得 2Ru + 2BᵀPs = 0 ⇒ u* = −R⁻¹BᵀP s ≡ −Ks，代回即得
代数 Riccati 方程 (ARE)：
AᵀP + PA − PBR⁻¹BᵀP + Q = 0
只要 (A,B) 可稳、Q ≥ 0、R > 0，ARE 存在唯一的对称正定稳定解 P，且闭环 A−BK
必定稳定——稳定性是免费赠送的，不需要再去试凑。这正是 LQR 与 PID 最本质的区别。
本平台怎么解 ARE（就在你的浏览器里）：
- Bass 算法给出保证稳定的初始增益：取 β > ‖A‖，解 MS + SMᵀ = 2BBᵀ
（M = A+βI），则 K₀ = BᵀS⁻¹ 必稳定。
证明：令 P=S⁻¹ 可得 (A−BBᵀP)ᵀP + P(A−BBᵀP) = −2βP < 0。
- Kleinman–Newton 迭代（对 ARE 做牛顿法，二次收敛）：解 Lyapunov 方程
(A−BKᵢ)ᵀP + P(A−BKᵢ) = −(Q + KᵢᵀRKᵢ)，再取 Kᵢ₊₁ = R⁻¹BᵀP。
Lyapunov 方程按元素展开成 16×16 线性方程组直接消元求解。
分析面板实时显示 ARE 残差（绝对与归一化两种口径）。这不是装饰：它是"我真的解出了方程"的证据。
谈残差必须归一化——绝对残差的量级由 ‖Q‖ 决定（此处 Q_θθ=328），理论下限约 1e−13，
连 scipy 自己的绝对残差都是 6.9e−12。本实现绝对残差 5.5e−13、归一化 1.7e−14。
另有一条完全独立的求解路径（自适应步长的 Riccati 微分方程反向积分）在自检脚本里做交叉验证，两者相对差 3.9e−11；
Python 参考实现 + scipy 构成第三、第四方复核，41 项对照全部通过。
一个可以手算的精确不变量：
K_x = −√(q_x / R)
它精确成立，与 m、L、g、摩擦、其他权重全都无关（由 Kalman 谱分解恒等式在 A 的零特征值方向取 s→0 极限可得）。
默认 q_x=4、R=0.01 ⇒ K_x = −20 精确命中；q_x/R=80000 ⇒ −282.842712475。
同理 LQI 的积分增益 k_i = −√(q_I/R)。这比"两个实现互相对照"更强——它是一把能用笔算校验任何 LQR 求解器的标尺。
Q、R 怎么选（Bryson 定则）：令 Qᵢᵢ = qᵢ / xᵢ,max²、
R = ρ / u_max²。这样 qᵢ、ρ 变成无量纲的"相对在意程度"，
量纲问题交给"允许的最大偏差"处理，调参立刻变得有物理直觉。
LQR 的白送裕度：单输入 LQR 的回路满足 |1 + L(jω)| ≥ 1，由此可证
增益裕度区间必含 [½, ∞)、相位裕度 ≥ 60°。分析面板会把实际的 α 稳定区间算给你看
（默认参数下达到 [0.06, ∞)，远好于理论下界）。
5. PID：能做到什么，做不到什么
(a) 单环 PID 有两个可证明的死穴。
死穴一：位置永远只是临界稳定。只反馈 θ 时，A 的第一列与 B 的第一元素都是零，
而反馈向量前两项也是零 ⇒ 闭环矩阵第一列仍全为零 ⇒ det = 0，必有原点极点（分析面板可以看到）。
所以角度能稳，位置临界稳定 → 任何微小扰动都让小车匀速漂走。这是结构性缺陷，调任何增益都无效。
死穴二：b > 0 时纯 PD 连角度都稳不住。把 §3 的含摩擦传函代入 PID，约简掉公共的 s 后得三次特征方程
D₀s³ + [(M+m)c + bJ + mlK_d]s² + [mlK_p − (M+m)mgl + bc]s + (mlK_i − b·mgl) = 0
三次多项式稳定要求全部系数同号，故常数项须为正：
K_i > b·g （默认参数：K_i > 0.981）
实测：b=0.1 时 K_i=0 ⇒ 闭环最大实部 +0.0308（不稳）；K_i=0.5 ⇒ +0.0151（仍不稳）；
K_i=1.079（1.1× 门限）⇒ +3e−6（稳）。而 b=0 时 K_i=0 就能稳。
教科书里"PD 就能稳住倒立摆角度"只在无摩擦模型下成立——有粘性摩擦时积分项是必需品，不是精修。
这个结论是与 Python 参考实现交叉验证时发现的，两个独立实现给出同样的门限。
(b) 串级结构。外环把"位置误差"翻译成"期望倾角"，内环去跟踪这个倾角：
θ_ref = Kp,o(x_ref − x) − Kd,oẋ + Ki,o∫(x_ref − x)dt
u = Kp(θ − θ_ref) + Ki∫(θ − θ_ref)dt + Kdθ̇
为什么符号是这样：因为 ∂θ̈/∂u < 0，所以 u 要和 θ 同号；
又因为稳态时 ẍ = g·tanθ，想往右加速就必须先向右倾斜，故外环增益取正。
(c) 解析调参（本平台默认值的来历）。忽略摩擦，内环闭环特征方程为
s² − p² + g₀(Kp + Kds) = 0（g₀ = ml/D₀），
于是
Kp = (ωn² + p²)/g₀, Kd = 2ζωn/g₀
外环把内环视为理想的 ẍ ≈ g·θ_ref，即双积分器：
Kp,o = ωo²/g, Kd,o = 2ζoωo/g
默认取 ωn=12 rad/s（> p=5.59）、ζ=0.9、ωo=1.2 rad/s（< z/2=2.5）、ζo=0.7。
这套参数是在 576 组候选里按"三种工况的总稳定时间 + 控制能量"筛出来的。
(d) 工程细节不是可选项。微分项必须一阶滤波（否则噪声直接乘以 Kd 进入执行器）；
积分必须抗饱和。本平台用"条件积分 + 积分限幅"，而没有用反计算(back-calculation)——
因为当比例项本身就远超限幅时，反计算会把积分器推向巨大的反向值，误差回落后积分项接管，造成反向飞车。
这个陷阱在开发本平台时真实发生过（小车飞到 8.9 m）。
6. 静差：为什么 LQR 也需要积分
加一个常值风扰 d（作用在小车上）。平衡态要求 θ=0 且合力为零，即 u控制 = −d，
而 u控制 = −K_x(x − x_ref)，于是
xss − x_ref = d / K_x
默认 K_x = −20 N/m、d = 1 N ⇒ 静差 −0.05 m。仿真实测值与此完全一致（自检脚本里有这一条断言）。
消除办法不是"手工外挂积分项"（符号极易搞错），而是把积分器并入状态做增广 LQR (LQI)：
ż = x − x_ref, 增广状态 [z, x, ẋ, θ, θ̇]
对增广系统解 ARE，Riccati 方程会自动给出符号正确的 ki；当积分权重 qI=0 时
解出的 ki 恰为 0，自动退化成普通 LQR。
7. 摆起：能量法（非线性控制）
直立点附近的线性控制器无法把摆从下垂位置甩起来（那里离平衡点 180°，线性模型完全失效）。
改用能量视角：取摆杆能量 E = ½Jθ̇² + mgl cosθ，目标 E_d = mgl。
对时间求导，并用运动方程消去 θ̈（a = ẍ 为小车加速度）：
Ė = −c θ̇² − m l · a · θ̇ cosθ
于是只要令 a = −k·θ̇cosθ（k>0），就有 Ė = mlk(θ̇cosθ)² > 0，
能量单调上升。取 k = k_E(E_d − E) 即得本平台使用的泵浦律：
a = k_E (E − E_d) θ̇ cosθ
它自带"能量过多时反向抽能"的性质。所需驱动力由完整非线性模型精确反解（forceForCartAccel）。
切换到 LQR 的判据（实测标定，不是拍脑袋）：本平台扫描过 LQR 在力饱和下的真实吸引域，
结论是"小车速度"才是关键约束——静止居中时 |θ|≤0.6 rad 都能接住，
但 |ẋ| ≥ 1 m/s 时即使 θ 很小也会失败（速度项 K_v·ẋ 把 10 N 的力矩全部吃光）。
所以切换条件是四条同时满足：|θ|<0.45、|θ̇|<2.5、|ẋ|<0.5、|x−x_ref|<0.7。
捕获后位置给定从"当前位置"以 0.25 m/s 斜坡滑回目标，实现无扰接管。
还有一个实践细节：摩擦会让能量停在 E_d 稍下方，摆差一点过不了顶（实测停在 0.355 rad）。
解决办法是把泵浦目标设为 E_d(1+ε)，ε=4%。
8. 一页总结：PID 与 LQR 的分工
| 串级 PID | LQR |
|---|
| 需要的先验 | 只需知道"谁快谁慢"，可现场试凑 | 需要完整状态空间模型 A、B |
| 需要的测量 | θ、x（速度可由差分近似） | 四个状态全部（实机需配观测器/Kalman） |
| 稳定性 | 自己保证，靠裕度检查 | 只要 Q≥0、R>0，数学上自动保证 |
| 调参对象 | 6 个增益，交互耦合 | 5 个权重，物理意义清晰、单调 |
| 多目标权衡 | 靠串级结构人为分离 | 写进代价函数，一次求解 |
| 鲁棒裕度 | 有限区间，需逐一验证 | 先天含 [½,∞) 增益裕度、≥60° 相位裕度 |
| 采样率要求 | Ts_crit = 111 ms（9 Hz） | Ts_crit = 37 ms（27 Hz），更苛刻 |
| 扩展性 | 加一个自由度就要重新搭结构 | 加状态只是把矩阵变大（LQI、LQG、MPC 同一路线） |
| 大范围机动 | 不适用 | 不适用（需能量法等非线性方法） |
关于采样率那一行别误会：LQR 更苛刻不是缺陷。实测两者的
T_s,crit · |λ|max（|λ|max = 最快闭环极点模长）几乎相同（1.075 vs 1.151），
差别完全来自"LQR 的最快极点快了 2.8 倍"，而那是我们自己在 Q/R 里要求的。
把 ρ 调大到 30，LQR 的 Ts_crit 就升到 159 ms。由此得到比教科书更好用的判据：
采样率要盯你设计出来的最快闭环极点，而不是被控对象的不稳定极点 p
—— 实测 p·Ts_crit 跨 12 倍（0.073→0.887）完全不可用，而 Ts_crit·|λ|max 稳定在 1.05~1.97。
工程取值 Ts ≤ 1/(10|λ|max)。详见 docs/理论推导.md §6.2 与实验 7。
教学结论：PID 是"局部经验的胜利"，LQR 是"模型知识的胜利"。
倒立摆之所以是经典教具，正因为它同时暴露了两者的边界——单环 PID 的结构性失败、
右半平面零点的带宽夹逼、LQR 对模型与全状态的依赖、以及线性方法在大范围机动上的彻底失效。