电路仿真器是如何工作的

电路仿真器是如何工作的

前言

THE DESIGNER’S GUIDE TO SPICE AND SPECTRE 的作者说：

被动用户受模拟器控制，而积极主动的用户则掌控模拟器。

被动用户运行模拟器时只是希望一切顺利，一旦出现问题，就从一堆解决办法中逐一尝试，却并不真正理解其原理，只希望其中一种能解决问题。如果问题仍未解决，要么重新设计电路以避免问题，要么干脆不用模拟器。而积极主动的用户会预料到可能出现的问题，当问题出现时，他们知道问题的原因，并且清楚该如何解决。

这句话我觉得特别有道理，是否掌握了仿真技术，就得看 TA 对于仿真器的理解，知道仿真器内部是如何运作得到结果的，而不是一味地假设仿真器正常运作而不断尝试修改参数（我以前就是这样的）。

---

仿真器流程简述

这里介绍 Spectre / SPICE 实现的流程，DC / 瞬态每一步的真实算法流程

1. 写出 MNA 非线性方程

   $$F(x) = 0$$

   • x：节点电压 + 电压源电流
   • F：KCL + 器件本构关系

2. 给定初始猜测

   $$x^{(0)}$$

3. 在当前点做一阶泰勒展开

   $$F(x^{(k)} + \Delta x) \approx F(x^{(k)}) + J(x^{(k)}) \Delta x$$

4. 构造线性 MNA 方程

   $$J(x^{(k)}) \Delta x = -F(x^{(k)})$$

5. 求解线性方程组（一次 MNA）

   • 稀疏矩阵
   • LU / Krylov

6. 更新解

   $$x^{(k+1)} = x^{(k)} + \Delta x$$

7. 判断收敛

   • 若满足 → 结束
   • 否则 → 回到第 3 步

瞬态分析只是多了一步：

把 $\frac{dx}{dt}$ 用数值积分公式转换成当前时刻的代数项，然后仍然进入同一套 Newton + MNA 流程

也就是说，每一个时间点都在解一个 $F^{(n)}(x) = 0$，用的完全是同一套牛顿线性化 + MNA

---

流程详解

上面浓缩了太多信息，下面会对流程简述里提到的概念一一解析，我们就能对仿真器内部在做什么有着更加清晰的认识。

MNA 非线性方程

MNA（Modified Nodal Analysis）是 SPICE 核心使用的组网方法，用来把电路变成数学方程。
它的形式是：

$$\mathbf{A(x)} \mathbf{x} = \mathbf{b}$$

不过这只是最简单线性的版本。对于非线性器件，它实际上是：

$$\mathbf{F(x)} = 0$$

其中：

• $\mathbf{x}$：节点电压、电流未知变量向量
• $\mathbf{F}$：由 KCL + 元件模型构成

拼装 MNA 方程的流程是：

读取 netlist → 建立拓扑 → 按 MNA 规则“stamp”每个器件 → 组装全局矩阵

首先了解一下 netlist 在仿真器眼里是什么？

以 Spectre / SPICE 为例，一行 netlist：

R1 n1 n2 1k

在仿真器眼里并不是“一个电阻”，而是一个结构化对象：

Device type: resistor
Node+: n1
Node-: n2
Value: 1000

结构化对象可以理解为一个封装好的数据对象，其中包含了该器件的类型，节点位置（表明电压方向以确定正负号），数值大小（电阻/电容/电感数值）。

因此，更一般地，netlist = 器件列表 + 节点连接关系

它提供了三样 MNA 必需的信息：

1. 节点编号（拓扑）
2. 器件类型
3. 器件参数 / 模型

MNA 的“组网”到底在做什么？

节点编号（node indexing）

仿真器第一步会：

n1 → 1
n2 → 2
gnd → 0 (reference node)

得到未知量向量：

$$x = \begin{bmatrix} v_1 \\ v_2 \\ \vdots \\ i_{V1} \\ i_{V2} \end{bmatrix}$$

相当于声明了一个数组变量，包括所有的节点电压和流过电压源 V1 V2 的支路电流（也可能有 V3 V4，看具体电路），得到了数组的大小并确定了不同的节点对应数组的第几个 index。

每个节点一条 KCL（电流守恒）

这是 MNA 的“骨架”：

对每一个非参考节点写一条 KCL

形式统一为：

$$\sum I_{\text{connected}}(x) = 0$$

对纯电阻电路，在数据处理上表现为导纳矩阵与节点电压向量的乘积。

器件“stamp”进方程（最关键一步）

每个器件都有一套固定的 MNA stamp 规则，仿真器做的事情就是：

for device in netlist:
    apply device.stamp() to global matrix

线性电阻的 stamp，R1 n1 n2 R 对节点 n1 的 KCL 为

$$+\frac{1}{R} v_1 - \frac{1}{R} v_2$$

对节点 n2 的 KCL：

$$-\frac{1}{R} v_1 + \frac{1}{R} v_2$$

仿真器的操作是像贴邮票一样在导纳矩阵中贴数据：

G[n1,n1] += 1/R
G[n1,n2] -= 1/R
G[n2,n1] -= 1/R
G[n2,n2] += 1/R

然后导纳矩阵与节点电压向量的乘积刚好就是 NA 方程。但此时的方程还不够完整，因为存在电压源，导纳矩阵法只能表达为 $I=GV$ 的形式，电压源不服从这个关系。于是未知向量变成：

$$\mathbf{x} = \begin{bmatrix} \mathbf{v} \\ \mathbf{i} \end{bmatrix} = \begin{bmatrix} v_1 \\ v_2 \\ \vdots \\ i_{V1} \end{bmatrix}$$

这一步是“Modified”的本质，也就是 MNA 中 M 的来源。这时候矩阵自然“长大”了，完整的 MNA 矩阵写成分块形式：

$$\begin{bmatrix} \mathbf{G} & \mathbf{A} \\ \mathbf{A}^T & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{v} \\ \mathbf{i} \end{bmatrix} = \begin{bmatrix} \mathbf{I} \\ \mathbf{E} \end{bmatrix}$$

其中 $G$ 是导纳矩阵，$A$ 是拓扑连接关系，是“电压源（或电感）与节点之间的关联矩阵”

MNA 矩阵可以看成两部分方程：

• 节点 KCL 方程（上半部分）

  $$\mathbf{G}\mathbf{v} + \mathbf{A}\mathbf{i} = \mathbf{I}$$

  含义是 节点电流 = 电阻电流 + 电压源支路电流

• 电压约束方程（下半部分）

  $$\mathbf{A}^T\mathbf{v}=\mathbf{E}$$

  结合拓扑连接关系

  A[n1, iV1] += 1
  A[iV1, n1] += 1  ← 实际上是 A^T
  b[iV1] = 1.8

  对 V1 来说就是电压源的定义：

  $$v_{n1} - v_0 = 1.8$$

  为了更形象地说明电压约束方程中拓扑连接关系 A 是如何与电压源联系起来的，可以从一个简单的例子出发

  

  该电路有四个节点，我们假定左下角节点为参考节点（接地），那么该节点电压为0，其余三个节点电压分别为 $V_1,V_2,V_3$，定义顺时针方向为支路电流正方向，此时拓扑连接关系 A 可以写为

  $$A=\quad \begin{array}{ccc} & E_1 & E_2 \\ n_1 & -1 & 0 \\ n_2 & 0 & 1 \\ n_3 & 0 & -1 \\ \end{array}$$

  $n_x,E_x$ 表示节点的 index，也描述了指向关系，例如 $A(n_1,E_1)=-1$ 代表了节点1与电压源1相连，但方向与支路电流方向相反。

  因此电压约束方程为

  $$\begin{aligned} \mathbf{A}^T\mathbf{v}=\mathbf{E}\quad\Rightarrow\quad \left( \begin{array}{ccc} -1 & 0 & 0\\ 0 & 1 & -1 \end{array} \right)\times \left( \begin{array}{c} V_1\\ V_2\\ V_3\\ \end{array} \right)= \left( \begin{array}{c} E_1\\ E_2\\ \end{array} \right)\quad\Rightarrow\quad \left\{ \begin{aligned} -V_1&=E_1\\ V_2-V_3&=E_2 \end{aligned} \right. \end{aligned}$$

  从上面就可以很清晰地看出，拓扑连接关系描述了节点之间存在的电势差关系，这就是电压源带来的电压约束方程。该方程的本质就是能量守恒，即环路上所有器件电压的代数和为零。

• G：描述“器件的 I–V 关系”（物理特性）

• A：描述“谁和谁连在一起”（拓扑与约束）

• MNA = 用 G 写电流 + 用 A 写电压约束

• 纯电阻网络时，A 自动消失

这里额外提到一个问题，为什么导纳矩阵中是自加而不是直接赋值，拓扑连接矩阵也是自加？

首先解释为什么导纳矩阵中是自加，这是因为一个矩阵元素 ≠ 一个器件。在 MNA 中，每一个矩阵元素，通常是“多个器件贡献的叠加结果”，相当于是将所有导纳归类到一起，所有电流源归类到一起，所有电压源归类到一起，整体形成一个庞大统一的 MNA 方程。但为什么拓扑连接矩阵也是自加，这难道不会出现 -1+1=0 连接断开的问题么？A 矩阵里用 +=，不是因为“可能重复连线”，而是因为“统一 stamping 机制 + 支持复杂器件”。在合法电路中，每个 (节点, 电压源) 对只贡献一个 ±1。

前面仅提到了线性电阻，如果加上电容/电感/非线性器件，MNA 方程会如何变化呢？

电容本质上也按“导纳” stamp。在 DC 分析中，$\omega=0$，理想电容对应导纳为 $j\omega C=0$，因此相当于开路，都没必要 stamp 进 MNA 里。在瞬态分析中，用 Backward Euler（Spectre 常用）：

$$i_C^n = C \frac{v^n - v^{n-1}}{\Delta t}$$

重写成：

$$i_C^n = \underbrace{\frac{C}{\Delta t}}_{G_{eq}} v^n \;-\; \underbrace{\frac{C}{\Delta t} v^{n-1}}_{I_{hist}}$$

因此在当前时间步，电容支路的 MNA 等效为一个导纳加一个历史电流源，stamp MNA 矩阵的方式和电阻一模一样，唯一的区别是电阻只需要 stamp 导纳矩阵 $G$，电容还需要 stamp 当前时间步节点电流 $I$。导纳 $G$ 自加归类到一起，历史电流源与节点电流 $I$ 归类到一起，电容融入进了 MNA 方程之中。

电感情况就不一样了，因为电感的节点分析量是电压，在 MNA 中的地位和电压源完全一样。因此对于电感支路，需要引入支路电流的未知量，矩阵维度 +1，用 A 矩阵连结点表示拓扑关系。

瞬态分析中，电感的 Backward Euler：

$$v^n = L \frac{i^n - i^{n-1}}{\Delta t}$$

整理为：

$$(v_{n1} - v_{n2}) \;-\; \frac{L}{\Delta t} i_L^n = -\frac{L}{\Delta t} i_L^{n-1}$$

拓扑项 → A / Aᵀ

A[n1, iL] += +1
A[n2, iL] += -1

对应：

$$+v_{n1} - v_{n2}$$

在理想电压源中：

$$A^T v = E$$

但在电感中，这一行变成：

$$A^T v \;\;+\;\; \left(-\frac{L}{\Delta t}\right) i_L = \text{hist}$$

这意味着：

• MNA 的右下角块不再是 0

• 而是：

  $$D = -\frac{L}{\Delta t}$$

所以更一般的 MNA 结构其实是：

$$\begin{bmatrix} G & A \\ A^T & D \end{bmatrix} \begin{bmatrix} v \\ i \end{bmatrix} = \begin{bmatrix} b \\ h \end{bmatrix}$$

行	方程类型	右端向量
上半部分	节点 KCL 方程（每个节点一行）	b（节点注入电流）
下半部分	支路约束方程（每个电压源/电感一行）	h（电压/历史项）

这正是电感的历史电压 stamp 进右端向量的体现

在代码中的表现为

rhs[iL] += - L/dt * i_prev

这里值得注意的是自加，这里的自加不是时间上的自加，而是空间上的。MNA 的矩阵和右端项，是“每一个时间步重新组装的”，典型流程是：

for each time step n:
    zero G, A, b, h
    for each device:
        device.stamp(current_state, previous_state)
    solve MNA

也就是说：

• += 只是在“当前时间步的组装过程”中叠加
• 不会跨时间步保留旧的 h

这里之所以自加，是因为一个时间步内，可能有多个器件同时往 h 里贡献，例如：

• 一条支路上可能有多个电感
• 或一个“复合器件”内部拆成多个等效支路

同样的道理对于电容也是一样，每个时间步 stamp 都要先清零 b/h 向量，自加的是同一个节点上多个器件的历史项。

通常每个电感一条支路方程，不会和别的电感混在同一个 h，因此这里自加主要还是实现上的统一写法或复合/受控/等效模型。

对于非线性器件，其伏安特性曲线由支路特征方程所描述，只依赖于已有的节点电压

$$I = f(V)$$

对其进行 Newton 线性化：

$$I(V) \approx I(V_0) + g(V_0) \cdot (V - V_0)$$

其中：

• $g = \frac{dI}{dV}$
• $I_{eq} = I(V_0) - g V_0$

在每一次 Newton 迭代都会得到一个：

• 等效导纳 $g$

• 等效电流源 $I_{eq}$

• 然后像电阻一样 stamp

MNA 方程结构不变，数值在变

1. netlist 决定矩阵结构（谁连谁、维度多大）
2. 器件模型决定 stamp 的数值
3. 瞬态用数值积分把时间吃掉
4. 非线性用 Newton 把非线性吃掉
5. 最终永远在解线性 MNA 矩阵

---

Newton 迭代法

Newton 法在 MNA 里不是“数值技巧”，而是“非线性器件被线性化并 stamp 进 MNA”的统一框架。，本质就是：

把非线性电路在当前工作点线性化成一个等效 MNA，反复修正，直到 KCL 不再被违反

为什么 MNA 一定会变成 Newton 迭代？

对纯线性电路：

$$A x = b$$

• $x$：节点电压 + 电压源电流
• 一次求解，完事

非线性器件一加入，灾难就来了，比如一个二极管：

$$i_D = I_S \left( e^{\frac{v_D}{V_T}} - 1 \right)$$

它让 MNA 变成：

$$f(x) = 0$$

其中

• $f(x)$：KCL + 器件本构关系
• 非线性、指数型

Newton–Raphson 的本质就是在当前工作点，把非线性器件“拉直”，解一个线性电路。

给定：

$$f(x) = 0$$

Newton 迭代：

$$J(x_k)\,\Delta x_k = -f(x_k)$$

其中：

• $x_k$：第 k 次猜测的解
• $J = \frac{\partial f}{\partial x}$：Jacobian

于是 Newton 步就是：

$$\underbrace{A(x_k)}_{\text{等效线性 MNA 矩阵}} \cdot \Delta x = \underbrace{b(x_k)}_{\text{等效激励}}$$

非线性器件在 Newton 中是如何处理的？

以二极管为例，原始关系：

$$i(v) = I_S e^{v/V_T}$$

在 $v_k$ 处一阶泰勒展开：

$$i(v) \approx i(v_k) + \left.\frac{di}{dv}\right|_{v_k}(v - v_k)$$

整理成：

$$i(v) = g_d \, v + I_{eq}$$

其中：

• $g_d = \frac{I_S}{V_T} e^{v_k/V_T}$（等效电导）
• $I_{eq} = i(v_k) - g_d v_k$（等效电流源）

每一轮 Newton：

1. 根据 x_k 计算所有非线性器件的：
   - g_eq
   - I_eq

2. 清空 MNA 矩阵 A 和 RHS（Right-Hand Side）b

3. stamp：
   - 线性器件（R, L, C）
   - 非线性器件的等效 g_eq, I_eq
   - 独立源

4. 解线性方程：
   A Δx = b

5. 更新：
   x_{k+1} = x_k + Δx

6. 判断收敛

初始值是如何猜测的？

在 MNA 求解里，常见有三类“初始值”，

层级	初始的是什么	用在什么时候
① 拓扑/线性 MNA	❌不需要猜	纯线性 DC / AC
② 非线性器件 Newton 迭代	结点电压、电流的初值	DC operating point
③ 瞬态分析	上一个时间步的解	Transient

只有 ② 是真的“要猜”的，③ 根本不是猜，是“继承”

如果是线性 MNA，形式为

$$A x = b$$

而且 A 不依赖 x（只有 R、独立源、线性 L/C 的离散形式），直接一次 LU / Cholesky（矩阵分解方法）解完

如果是非线性 MNA（找 DC 工作点），形式为

$$f(x) = 0$$

则求解下一个迭代位置

$$x_{k+1} = x_k - J^{-1}(x_k) f(x_k)$$

从图中就可以明白牛顿迭代法到底在干什么，就是根据合理的猜测值不断一阶近似找到其 $x$ 轴上的截距 $x_{k+1}$，进而不断逼近真实解。

初始值怎么给？SPICE / Spectre 的策略：

1. 所有结点电压 = 0 V
2. 独立电压源 → 强制已知
3. MOS，$V_{gs}$ = 0，全部当关断
4. BJT，$V_{be}$ = 0，全部当截止
5. 二极管初始电流 = 0

如果初始值乱给，结果会无法收敛。

Spectre 和 HSPICE 中直接进行牛顿迭代往往因为初始猜想值偏离实际解太远而导致不收敛。为了解决这个问题，常用的三种主要策略：

• 连续法 / 步进法 (Continuation / Stepping Methods)

  当标准的牛顿法失败时，仿真器会自动切换到步进法。其核心思想是将一个“难解”的问题通过一个参数变量转化成一系列“易解”的子问题。

  • 源步进 (Source Stepping)： 将电路中所有的独立源（电压源、电流源）从 0 开始逐渐增加到设定值。在每个步进点，使用前一个点的解作为牛顿迭代的初始值。
  • Gmin 步进 (Gmin Stepping)： 在电路的每个非线性节点（如二极管、晶体管的 PN 结）并联一个非常大的电导（Gmin）。初始时 Gmin 较大，使矩阵变得“容易收敛”，随后逐渐减小 Gmin 直至其对电路影响消失，每一步都利用上一步的结果作为初值。

• 用户提供的初始值 (User-Provided Initial Values)

  用户可以通过以下命令手动干预，为牛顿迭代提供更接近真实解的起点：

  • .NODESET： 给定节点的初始估计值。仿真器会在初始迭代时在这些节点接上电压源以强迫收敛，之后再移除这些源进行最终求解。
  • .IC (Initial Conditions)： 主要用于瞬态分析的起点，但在某些情况下也可辅助 DC 分析。

• 伪瞬态分析 (Pseudo-Transient Analysis)

  当步进法也失败时，Spectre 等仿真器会尝试伪瞬态分析。

  • 它给电路中的节点添加虚拟电容和电感，将求解静态 DC 方程转化为求解动态演变过程。
  • 随着时间推进，电路最终会达到稳定状态（即导数为零的点），这个稳定点就是 DC 解。

在 Spectre 或 HSPICE 执行 dc 分析时，通常遵循以下优先级逻辑：

1. 直接牛顿迭代： 使用默认初值（通常是全 0 或基于用户 .NODESET）直接求解。
2. 自动步进： 若失败，尝试 Gmin Stepping，再失败则尝试 Source Stepping。
3. 备选算法： 最终可能会调用 Pseudo-Transient 或是 Homotopy（同伦法） 等高级算法。

瞬态分析

原始的瞬态 MNA 是这样的（考虑电容）：

$$\mathbf{C}(x)\,\dot x(t) + \mathbf{G}(x)\,x(t) = \mathbf{u}(t)$$

问题在于：

• Newton 只能解：

  $$f(x)=0$$

• 但现在 MNA 里有：

  $$\dot x(t)$$

原始瞬态方程的一项：

$$C\,\dot x(t_n)$$

代入 BE：

$$C \frac{x_n - x_{n-1}}{\Delta t}$$

整理：

$$\underbrace{\frac{C}{\Delta t}}_{\text{等效导纳}} x_n \;-\; \underbrace{\frac{C}{\Delta t}x_{n-1}}_{\text{历史电流源}}$$

这和前面提到的电容分析模式是一样的，Newton 现在要解的当前时间步的 MNA 方程为

$$f(x_n) = \Big( \mathbf{G} + \frac{\mathbf{C}}{\Delta t} \Big) x_n - \mathbf{b}_{\text{hist}} - \mathbf{u}(t_n) = 0$$

本质就变为了一个线性方程，即便内部包含了非线性器件，每次迭代求解的依然是一个线性方程。因为你只要找到了非线性器件的 IV 关系 $I=f(V)$，你就能将之转化为某个工作点 $x_k$ 下的等效导纳和等效电流源，然后整合进 MNA 线性方程之中。这就是前面提到的：

1. 瞬态用数值积分把时间导数吃掉
2. 非线性用 Newton 把非线性吃掉