电路的计算机仿真 电磁学小论文 演示文稿 邱哲儒 PB15000034.

Presentation on theme: "电路的计算机仿真 电磁学小论文 演示文稿 邱哲儒 PB15000034."— Presentation transcript:

1 电路的计算机仿真 电磁学小论文 演示文稿 邱哲儒 PB

2 目录 一.理论讨论 二.计算机实现 1.编写的程序XLWSpice 1.从最简单的情况开始 2.对恒压源的处理 3.受控电源的处理
4.非线性电路的迭代求解 5.电路的时域动态分析 二.计算机实现 1.编写的程序XLWSpice 2.应用举例 i)包含二极管的算例 ii)RC,RLC电路暂态过程的算例 iii)大型电阻网络的计算 3.进一步发展的思路

3 最简单的情况：只有电流源与线性电阻 用节点法处理 𝐼 𝑠𝑟 𝑐 𝑘 + 𝑗 𝐼 𝑗𝑘 =0
回路法处理图的问题，确定独立回路的过程较为复杂，节点法建立方程的过程很直观，便于自动化处理 对每个节点列出基尔霍夫电流方程： 对第k个节点： 𝐼 𝑠𝑟 𝑐 𝑘 指电流源注入节点k的电流，取流入方向为正， 𝐼 𝑗𝑘 是从j流向k节点的电流 对每个电阻列出欧姆定律方程： G 𝑗𝑘 𝑙 是j,k节点间的第l个线性电阻的电导，用0表示节点间无连接,取第0节点作为电势参考点 𝑢 0 ≡0 整理可得对于第k个节点的方程 𝐼 𝑠𝑟 𝑐 𝑘 + 𝑗 𝐼 𝑗𝑘 =0 𝐼 𝑗𝑘 =( 𝑢 𝑘 − 𝑢 𝑗 ) 𝑙 𝐺 𝑗𝑘 𝑙 − 𝑢 𝑘 𝑗,𝑙 𝐺 𝑗𝑘 𝑙 + 𝑗 𝑙 𝐺 𝑗𝑘 𝑙 𝑢 𝑗 = 𝐼 𝑠𝑟 𝑐 𝑘

4 写成矩阵的形式 𝐴𝑢=i 𝑢= 𝑢 1 𝑢 2 ⋮ 𝑢 𝑛 ，i= 𝐼 𝑠𝑟 𝑐 1 𝐼 𝑠𝑟 𝑐 2 ⋮ 𝐼 𝑠𝑟 𝑐 𝑛 ， 𝐴 𝑘𝑝 = − 𝑗,𝑙 𝐺 𝑗𝑘 𝑙 ,𝑘=𝑝 𝑙 𝐺 𝑝𝑘 𝑙 ,𝑘≠𝑝 可以根据各电阻元件参数构建A矩阵 此时的A矩阵是对称的 实际上可以进一步简化处理 case R://R //Building A if (pe.node0) Spmat.coeffRef(pe.node0 - 1, pe.node0 - 1) += 1.0 / pe.value; //对 𝑨 𝒏𝒐𝒅𝒆𝟎,𝒏𝒐𝒅𝒆𝟎 的贡献 if (pe.node1) Spmat.coeffRef(pe.node1 - 1, pe.node1 - 1) += 1.0 / pe.value; //对 𝑨 𝒏𝒐𝒅𝒆𝟏,𝒏𝒐𝒅𝒆𝟏 的贡献 if (pe.node0&&pe.node1){ Spmat.coeffRef(pe.node0 - 1, pe.node1 - 1) += -1.0 / pe.value; Spmat.coeffRef(pe.node1 - 1, pe.node0 - 1) += -1.0 / pe.value; //对 𝑨 𝒏𝒐𝒅𝒆𝟎,𝒏𝒐𝒅𝒆𝟏 ， 𝑨 𝒏𝒐𝒅𝒆𝟏,𝒏𝒐𝒅𝒆𝟎 的贡献}

5 电压源的处理 实际上理想的电压源是不存在，理论上实际线性电源又可以用戴维南定理等效，也可以诺顿等效，看 作电流源并联电阻，可以只使用电流源进行一切计算。 这常常给建模带来麻烦，很多情况下，我们使用的电源可能等效内阻很小（比如市电，大型动力蓄电 池），更有可能是包含负反馈机制维持输出电压基本与负载无关的所谓稳压电源，其在工作电流范围 内几乎相当于理想电压源。 处理电压源的问题常用的方法包括“列表法”(2003年某篇获奖小论文所讨论）与改进节点电压法。 本质是另外添加关于电压源的约束方程，将电压源提供的电流作为独立的未知量。

6 以下采用改进节点电压法(Modified Nodal Analysis)
电压源独立的对两个节点电压提供了一个约束，设电压源连接在node0与node1之间，自然有方 程: 𝑢 𝑛𝑜𝑑𝑒1 − 𝑢 𝑛𝑜𝑑𝑒0 = 𝑈 𝑠𝑟𝑐 电压源提供的电流应该计入各节点电流方程，恒压源的电流是可以任意取值的，只能将流入电 压源正端的电流作为一个独立变量 𝐼 𝑠𝑟 𝑐 𝑙 ，可得对每个节点有： 对矩阵方程𝐴𝑢=𝑖进行扩充,以包含新的未知量与方程。 − 𝑢 𝑘 𝑗,𝑙 𝐺 𝑗𝑘 𝑙 + 𝑗 𝑙 𝐺 𝑗𝑘 𝑙 𝑢 𝑗 + 𝑙 𝐼 𝑈𝑠𝑟 𝑐 𝑙 ( 𝛿 𝑛𝑜𝑑𝑒 1 𝑙 ,𝑘 − 𝛿 𝑛𝑜𝑑𝑒 0 𝑙 ,𝑘 = 𝐼 𝑠𝑟 𝑐 𝑘 ……………式(*) 𝑀 𝑈 𝑘 𝐼 𝑈𝑠𝑟𝑐 = 𝐼 𝑠𝑟𝑐 𝑈 𝑠𝑟𝑐 𝑀= 𝐴 𝐵 𝐶 𝐷 𝐵 𝑘𝑙 = 𝑙 𝛿 𝑛𝑜𝑑𝑒 1 𝑙 ,𝑘 − 𝛿 𝑛𝑜𝑑𝑒 0 𝑙 ,𝑘 𝐶= 𝐵 𝑇 𝐷=𝐼

7 电压源的处理 理想运算放大器 理想运算放大器的处理非常简单，由“虚短”的概念有 𝑢 + − 𝑢 − =0,输出端的电压电流都没有限制。
case V: if (e.node0){//给B矩阵添加 𝒍 𝜹 𝒏𝒐𝒅𝒆 𝟏 𝒍 ,𝒌 − 𝜹 𝒏𝒐𝒅𝒆 𝟎 𝒍 ,𝒌 项 Spmat.insert(pp.nodecnt + e.no_NLDC, e.node0 - 1) += 1; //设置C矩阵, 𝒖 𝒏𝒐𝒅𝒆𝟎 项系数设为1， 𝒖 𝒏𝒐𝒅𝒆𝟏 项系数设为-1 Spmat.insert(e.node0 - 1, pp.nodecnt + e.no_NLDC) += 1;} if (e.node1){Spmat.insert(pp.nodecnt + e.no_NLDC, e.node1 - 1) += -1; Spmat.insert(e.node1 - 1, pp.nodecnt + e.no_NLDC) += -1;} break; 理想运算放大器 跳过理想运放 理想运算放大器的处理非常简单，由“虚短”的概念有 𝑢 + − 𝑢 − =0,输出端的电压电流都没有限制。 当然，理想运放的模型只在有合适负反馈机制的情况下适用，否则方程欠定，无法求解。

8 受控电源的处理 在模拟三极管，MOSFET等器件时都需要受控电源，受控源有四种，压控电压源（VCVS），流 控电压源（CCVS），压控电流源（VCCS），流控电流源（CCCS）。 在引入受控电源以后，方程要进行修改。 1.压控电流源（VCCS） 压控电流源的定义式是： 处理时相当于在式(*)中左侧加入 一项 需要在A矩阵中加入项 𝐼 𝑉𝐶𝐶 𝑆 𝑙 =𝛼( 𝑢 𝑠𝑒𝑛𝑠𝑒+ − 𝑢 𝑠𝑒𝑛𝑠𝑒− 𝑙 𝛼( 𝛿 𝑛𝑜𝑑𝑒 1 𝑙 ,𝑘 − 𝛿 𝑛𝑜𝑑𝑒 0 𝑙 ,𝑘 )( 𝑢 𝑠𝑒𝑛𝑠𝑒+ − 𝑢 𝑠𝑒𝑛𝑠𝑒− case G://VCCS if (e.node0)//output positive{ if (e.node2)//sense positive Spmat.coeffRef(e.node0 - 1, e.node2 - 1) += e.value; if (e.node3)//sense negative Spmat.coeffRef(e.node0 - 1, e.node3 - 1) += -e.value;} if (e.node1)//output negative{ Spmat.coeffRef(e.node1 - 1, e.node2 - 1) += -e.value; Spmat.coeffRef(e.node1 - 1, e.node3 - 1) += e.value;} 以VCCS为例子，略去其他受控源的讨论

9 2.流控电流源(CCCS) 流控电流源是双极型三极管BJT放大区的简单模型，理想变压器的模型，应用广泛。
定义式为 𝐼 𝐶𝐶𝐶𝑆 𝑙 =𝑐 𝐼 𝑠𝑒𝑛𝑠𝑒 控制端node2,node3相当于连接了一个0V的电压源而被短路，有 被控端node0,node1相当于增加了一个电流源，(*)式左端添加 一项， 𝑢 𝑛𝑜𝑑𝑒2 − 𝑢 𝑛𝑜𝑑𝑒3 =0 𝑙 𝑏( 𝛿 𝑛𝑜𝑑𝑒 1 𝑙 ,𝑘 − 𝛿 𝑛𝑜𝑑𝑒 0 𝑙 ,𝑘 ) 𝐼 𝑠𝑒𝑛𝑠𝑒 case F://CCCS if (e.node0) //更改B矩阵 Spmat.coeffRef(e.node0 - 1, pp.nodecnt + e.no_NLDC) += e.value;//out positive if (e.node1) Spmat.coeffRef(e.node1 - 1, pp.nodecnt + e.no_NLDC) += -e.value;//out negative if (e.node2){ Spmat.coeffRef(e.node2 - 1, pp.nodecnt + e.no_NLDC) += -1;//sense flow in Spmat.insert(pp.nodecnt + e.no_NLDC, e.node2 - 1) = 1; //更改C矩阵} if (e.node3){ Spmat.coeffRef(e.node3 - 1, pp.nodecnt + e.no_NLDC) += 1;//sense flow out Spmat.insert(pp.nodecnt + e.no_NLDC, e.node3 - 1) = -1; //更改C矩阵} 跳

10 3.压控电压源(VCVS) 压控电压源是MOSFET管工作在放大区的最简单的模型，在实际使用中价值很大。 定义是
作为电压源需要增加一个未知电流 𝐼 𝑉𝐶𝑉𝑆 𝑙 ，(*)式的左侧应该加上 一项，此时B,C矩阵都需要更新。 𝑢 𝑛𝑜𝑑𝑒0 − 𝑢 𝑛𝑜𝑑𝑒1 =𝑐( 𝑢 𝑛𝑜𝑑𝑒3 − 𝑢 𝑛𝑜𝑑𝑒4 𝑙 𝛿 𝑛𝑜𝑑𝑒 1 𝑙 ,𝑘 − 𝛿 𝑛𝑜𝑑𝑒 0 𝑙 ,𝑘 ) 𝐼 𝑉𝐶𝑉 𝑆 𝑙 case E://VCVS //Building B and C if (e.node0){//修改B矩阵 Spmat.coeffRef(pp.nodecnt + e.no_NLDC, e.node0 - 1) += 1;//out positive //修改C矩阵 Spmat.coeffRef(e.node0 - 1, pp.nodecnt + e.no_NLDC) += 1;} if (e.node1){ Spmat.coeffRef(pp.nodecnt + e.no_NLDC, e.node1 - 1) += -1;//out negative Spmat.coeffRef(e.node1 - 1, pp.nodecnt + e.no_NLDC) += -1;} if (e.node2)//修改B矩阵 Spmat.coeffRef(pp.nodecnt + e.no_NLDC, e.node2 - 1) += e.value;//sense positive if (e.node3) Spmat.coeffRef(pp.nodecnt + e.no_NLDC, e.node3 - 1) += -e.value;//sense negative 跳

11 4.流控电压源(CCVS) 其定义为 控制端node2,node3相当于连接了一个0V的电压源而被短路，
这时必须添加两个未知电流 𝐼 𝑠𝑒𝑛𝑠𝑒 , 𝐼 𝑜𝑢𝑡 来表征。(*)式的左端应当添加 此时不仅B,C矩阵需要修改，D矩阵也要插入新的元素 𝑢 𝑛𝑜𝑑𝑒0 − 𝑢 𝑛𝑜𝑑𝑒1 =𝑑 𝐼 𝑠𝑒𝑛𝑠𝑒 𝑢 𝑛𝑜𝑑𝑒2 − 𝑢 𝑛𝑜𝑑𝑒3 =0 𝑙 𝛿 𝑛𝑜𝑑𝑒 1 𝑙 ,𝑘 − 𝛿 𝑛𝑜𝑑𝑒 0 𝑙 ,𝑘 𝐼 𝐶𝐶𝑉𝑆 𝑙 + 𝑙 ( 𝛿 𝑛𝑜𝑑𝑒 2 𝑙 ,𝑘 − 𝛿 𝑛𝑜𝑑𝑒 3 𝑙 ,𝑘 )𝐼 𝑠𝑒𝑛𝑠𝑒 𝑙 case H://CCVS if (e.node0){//修改C矩阵 Spmat.insert(pp.nodecnt + e.no_NLDC, e.node0 - 1) = -1;//out positive //修改B矩阵 Spmat.insert(e.node0 - 1, pp.nodecnt + e.no_NLDC + 1) = 1;} if (e.node1) { Spmat.insert(pp.nodecnt + e.no_NLDC, e.node1 - 1) = 1;//out negative Spmat.insert(e.node1 - 1, pp.nodecnt + e.no_NLDC + 1) = -1; } if (e.node2){ Spmat.insert(e.node2 - 1, pp.nodecnt + e.no_NLDC) = 1;//sense flow in Spmat.insert(pp.nodecnt + e.no_NLDC + 1, e.node2 - 1) = 1;} if (e.node3) { Spmat.insert(e.node3 - 1, pp.nodecnt + e.no_NLDC) = -1;//sense flow out Spmat.insert(pp.nodecnt + e.no_NLDC + 1, e.node3 - 1) = -1; } //修改D矩阵 Spmat.insert(pp.nodecnt + e.no_NLDC, pp.nodecnt + e.no_NLDC) = e.value; 跳

12 举例 对于如图所示的电路 包含多种线性器件 按前述方法处理 得到的M矩阵如下 不再严格对称（不能“投机取巧”） 方程右侧的列向量为
%输出结果 %%%%%%NON LINEAR DC ANALYSIS%%%%%% %ABSOLUTE COVERGENCE LIMIT:0.001 %RELATIVE COVERGENCE LIMIT:0.0001 %HAVING NODES:6 %NON LINEAR DC ANALYSIS ITERATION 0 %NODE VOLTAGE 1,10 2, 3, 4, 5, 6,12 %VOLTAGE SOURCES CURRENT 0, ,V 1, ,V %STOP ITERATION CRITERIA MET 举例 对于如图所示的电路 包含多种线性器件 按前述方法处理 得到的M矩阵如下 不再严格对称（不能“投机取巧”） 方程右侧的列向量为 [ ] 𝑇 1 2 5 4 6 3 *输入文件 .DC 1e-3 1e-4 V R 1 2 5 R 2 3 6 I 4 0 1 R R V F

13 𝐽 𝑋 𝑛 𝑋 𝑛+1 =−𝐹 𝑋 𝑛 +𝐽 𝑋 𝑛 𝑋 𝑛 ………(**)
非线性元件的处理 对于含有非线性元件的电路，电路方程组不是线性的，不能直接通过线性代数的方法简便得解， 可以迭代近似求解。 用牛顿法近似求解非线性方程组 若有方程组 𝑓 𝑖 𝑥 1 , 𝑥 2 ,…, 𝑥 𝑛 =0(𝑖=1,2,…,𝑛) 泰勒级数展开至线性项 𝑓 𝑖 𝑥 1 + ℎ 1 , 𝑥 2 + ℎ 2 ,…, 𝑥 𝑛 + ℎ 𝑛 = 𝑓 𝑖 𝑥 1 , 𝑥 2 ,…, 𝑥 𝑛 + 𝐽 𝑥 1 , 𝑥 2 ,…, 𝑥 ℎ 1 , ℎ 2 ,…, ℎ 3 𝑇 𝑖 要求 𝑋 𝑛+1 是比 𝑋 𝑛 更好的近似，𝐹 X 𝑛+1 应该更接近0 0=𝐹 X 𝑛+1 ≈𝐹 𝑋 𝑛 +𝐽 𝑋 𝑛 ( 𝑋 𝑛−1 − 𝑋 𝑛 ) 𝐽( 𝑋 𝑛 )是方程组的Jacobian 𝐽 𝑋 𝑛 𝑋 𝑛+1 =−𝐹 𝑋 𝑛 +𝐽 𝑋 𝑛 𝑋 𝑛 ………(**)

14 𝐹 𝑋 =𝑀𝑋+ 𝐼 𝑁𝐿 (𝑢) 0 − 𝐼 𝑠𝑟𝑐 𝑈 𝑠𝑟𝑐 =0,𝑋= 𝑢 𝑖 𝑢𝑠𝑟𝑐
𝐹 𝑋 =𝑀𝑋+ 𝐼 𝑁𝐿 (𝑢) 0 − 𝐼 𝑠𝑟𝑐 𝑈 𝑠𝑟𝑐 =0,𝑋= 𝑢 𝑖 𝑢𝑠𝑟𝑐 只包含两端非线性器件的电路的 𝐼 𝑁𝐿 (𝑢)是 𝐼 𝑁𝐿 𝑗 (𝑢)= 𝑙 𝑖 𝑙 ( 𝑢 𝑛𝑜𝑑𝑒 + 𝑙 − 𝑢 𝑛𝑜𝑑𝑒 − 𝑙 )( 𝛿 𝑛𝑜𝑑𝑒 + 𝑙 ,𝑗 − 𝛿 𝑛𝑜𝑑𝑒 − 𝑙 ,𝑗 ) (器件下标为𝑙，node+,node-分别为正端，负端节点号) 其Jacobian为𝑀+𝑁(𝑋) 𝑁 𝑝𝑞 (𝑋)= 𝑙 𝛿 𝑛𝑜𝑑𝑒 + 𝑙 ,𝑝 − 𝛿 𝑛𝑜𝑑𝑒 − 𝑙 ,𝑝 𝛿 𝑛𝑜𝑑𝑒 + 𝑙 ,𝑞 − 𝛿 𝑛𝑜𝑑𝑒 − 𝑙 ,𝑞 𝑑 𝑖 𝑙 𝑑𝑢 ​ 𝑢=( 𝑢 𝑛𝑜𝑑𝑒 + 𝑙 − 𝑢 𝑛𝑜𝑑𝑒 − 𝑙 ) 利用（**）式可得 𝑀+𝑁 𝑋 𝑛 𝑋 𝑛+1 = 𝐼 𝑠𝑟𝑐 𝑈 𝑠𝑟𝑐 +𝑁 𝑋 𝑛 𝑋 𝑛 − 𝐼 𝑁𝐿 𝑢 0 ………(***) 多端非线性器件的处理类似(XLWSpice中暂不支持) 等效电导 等效电流源

15 经过线性近似后，在每次迭代中非线性器件可以等效为一个电流源并联一个电阻 由(. )式 知如果取 对电路方程进行求解与求解(
经过线性近似后，在每次迭代中非线性器件可以等效为一个电流源并联一个电阻 由(***)式 知如果取 对电路方程进行求解与求解(***)式完全相同，此时问题化为已经处理过的情况 实际上当两次迭代间的电压差足够小就可以认为已经找到解，停止迭代 默认要求 𝑢 𝑛 − 𝑢 𝑛−1 ∞ <1𝜇𝑉 𝐺 𝑒𝑞 (𝑢)= 𝑑𝑖 𝑑𝑢 𝐼 𝑒𝑞 (𝑢)=𝑖(𝑢)−𝑢 𝐺 𝑒𝑞 (𝑢 略,简单讲述图

16 包含非线性元件的电路的求解流程 组装只包括静态线性元件的M,在每次迭代中加入N

17 𝑥 𝑛+1 = 𝑥 𝑛 +𝑡𝑖𝑚𝑒𝑠𝑡𝑒𝑝× 𝑥 𝑛 + 𝑥 𝑛+1 2 ………(****)
电路的时域动态分析 对于伏安特性包含电流或电压导数的动态元件，比如电容，电感，可以简单地通过引入复数参数通过相量法 处理问题。 但是这样的方式很有局限性，一是一般只适用于处理稳态情况下的问题 二是若电路中含有非线性器件，叠加原理不再适用，电路中会产生不同频的谐波分量，（且对于具有指数伏 安特性的二极管等，频域展开收敛比较慢）处理较为困难。 对电路微分代数方程组进行数值积分是更为实用的方法，数值求解微分方程的算法种类众多。在电路微分代 数方程的求解中，隐式方法一般较显式方法收敛情况更好，对于非线性方程组也不可避免也不会加大计算量。 梯形法虽较原始，但容易实现，效果也还不错。取积分的步长为timestep，有公式 𝑥 𝑛+1 = 𝑥 𝑛 +𝑡𝑖𝑚𝑒𝑠𝑡𝑒𝑝× 𝑥 𝑛 + 𝑥 𝑛 ………(****) 以电容器为例：其伏安特性为𝐼=𝐶 𝑑𝑢 𝑑𝑡 ，代入(****)式中整理得 𝑈 𝑛+1 = 𝑈 𝑛 + 𝑡𝑖𝑚𝑒𝑠𝑡𝑒𝑝 2𝐶 ( 𝐼 𝑛 + 𝐼 𝑛+1 ) 在每一个时间步长中，其前一项可以看作一个阻值为 𝑡𝑖𝑚𝑒𝑠𝑡𝑒𝑝 2𝐶 的电阻，后一项可以看作一个电流为 − 2𝐶 𝑡𝑖𝑚𝑒𝑠𝑡𝑒𝑝 𝑈 𝑛 − 𝐼 𝑛 的电流源，两个元件并联。这样在每一次迭代中问题都可以化归为已经处理过的问题。 电感器的情况完全类似 𝐼 𝑛+1 = 2𝐶 𝑡𝑖𝑚𝑒𝑠𝑡𝑒𝑝 𝑈 𝑛+1 +(− 2𝐶 𝑡𝑖𝑚𝑒𝑠𝑡𝑒𝑝 𝑈 𝑛 − 𝐼 𝑛 ) 多些时间

18 //在BuildInitLHSSpMat_TRAN中
case C://修改A矩阵，添加视为电阻的部分 if (de.node0) Spmat.coeffRef(de.node0 - 1, de.node0 - 1) += 2.0*de.value0/pp.TimeStep; if (de.node1) Spmat.coeffRef(de.node1 - 1, de.node1 - 1) += 2.0*de.value0 / pp.TimeStep; if (de.node0&&de.node1){ Spmat.coeffRef(de.node0 - 1, de.node1 - 1) += -2.0*de.value0 / pp.TimeStep; Spmat.coeffRef(de.node1 - 1, de.node0 - 1) += -2.0*de.value0 / pp.TimeStep;} break; //在AddEquivalentCurrSrc_DYElems_NLTRAN中 case C://添加视为电流源的部分 dynellist[i].tmp = dynellist[i].lasti * dynellist[i].value0* dynellist[i].currvoltage / pp.TimeStep;//等效电流 if (dynellist[i].node0) zCurr(dynellist[i].node0 - 1) += dynellist[i].tmp; if (dynellist[i].node1) zCurr(dynellist[i].node1 - 1) += -dynellist[i].tmp; break;

19 电路的时域动态分析流程

20 实际实现：编写程序XLWSpice 特点： 使用C++语言，处理复杂问题时性能较高。 代码精简,不含外部库只有一千余行
使用高性能线性代数库Eigen进行稀疏矩阵的运算 通过PardisoSupport模块调用Intel MKL Pardiso LU进行方程求解（兼顾Eigen的友好接口，MKL的优秀性能） 将非线性元件的描述(𝑖 𝑢 , 𝐼 𝑒𝑞 , 𝐺 𝑒𝑞 )和交流信号的描述置入外部动态链接库中，方便 更改，不需重新编译主程序。 #include <Eigen/Sparse> 快速！带过

21 采用稀疏矩阵方法处理大型的𝑀矩阵 如果电路中有n个节点，m个需要求解输出电流的电压源，用于求解的M矩阵尺寸在 O( 𝑛+𝑚 2 )量级，在存储空间上的开销较大。然而实际上常见的电路中每一个节点一般 只有2或3个支路，一般不会连接超过4个支路，相当于矩阵的每行只有几个非零元素，采 用稀疏矩阵的技术处理可以大大降低存储空间上的开销，为𝑂(𝑛+𝑚)量级，远优于密集 矩阵处理方式。 实际上XLWSpice预分配空间时 //Reserving elements for each column Spmat.reserve(VectorXi::Constant(dim, 4 * (pp.passiveelemcnt + pp.unknownicnt_TRAN) / pp.nodecnt - 1)); 只为矩阵的每行分配了一般不超过7个元素 方程的分解与求解也明显更高效,直接对稀疏矩阵进行LU分解时可以通过一些方法大大减 少需要操作的元素，减少内存开销。迭代算法也比较有效。 使得XLWSpice可以毫不费力的在普通硬件上处理百万节点级别的问题 （使用Intel PARDISO 的Out-Of-Core功能可以将中间结果暂存在硬盘，在牺牲一定速度的前提下实现千万节点 级别问题的求解） solver.pardisoParameterArray()[59] = 1;

22 大型矩阵方程的求解 直接法：如LU分解M=𝐿𝑈 迭代法：如BiCGSTAB，以解𝑀𝑥=𝑏为例 1.选择初值 𝑥 0
2.取 𝑟 0 =𝑏−𝑀 𝑥 0 3.任取 𝑟 𝑟 0 ∙ 𝑟 0 ≠0 ,取 𝜌 0 =𝛼= 𝜔 0 =1, 𝒗 𝟎 = 𝒑 𝟎 =𝟎 4.对𝑖=1,2,3,…循环： 𝜌 𝑖 = 𝑟 0 ∙ 𝑟 𝑖−1 𝛽 = ( 𝜌 𝑖 / 𝜌 𝑖−1 )(𝛼/ 𝜔 𝑖−1 ) 搜索方向 𝒑 𝑖 = 𝒓 𝑖 −1 + 𝛽( 𝒑 𝑖−1 − 𝜔 𝑖−1 𝒗 𝑖−1 ) 𝑣 𝑖 = A 𝑝 𝑖 α = ρi/ 𝑟 0 ∙ 𝑟 𝑖−1 h = 𝒙 𝒊−𝟏 + α 𝒑 𝒊 如果h作为近似解已经足够精确，则退出 另外可以预先分解𝑀≈ 𝐾 1 𝐾 2 =𝐾 使用预处理的(Preconditioned)BiCGSTAB 比如使用不完全LU分解(ILU)算法 可以极大程度地加快收敛过程 𝒔= 𝒓 𝒊−𝟏 −𝛼 𝒗 𝒊 𝒕=𝑨𝒔 𝜔 𝑖 =(𝒕 ⋅𝒔)/(𝒕⋅𝒕) 𝒙 𝒊 =𝒉+ 𝜔 𝑖 𝒔 如果 𝑥 𝑖 作为近似解已经足够精确,则退出 残差 𝒓 𝒊 =𝒔− 𝜔 𝑖 𝒕 跳过

23 输入为网表数据，详情见源码注释 输出示例： .DC 1e-4 1e-3 VSRC1 1 0 10 R 1 2 5 R 2 3 6
I 4 0 1 R R VSRC F 输入为网表数据，详情见源码注释 输出示例： %%%%%%NON LINEAR DC ANALYSIS%%%%%% %ABSOLUTE COVERGENCE LIMIT:0.0001 %RELATIVE COVERGENCE LIMIT:0.001 %HAVING NODES:14 %NON LINEAR DC ANALYSIS ITERATION 0 %NODE VOLTAGE 1,1 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, %VOLTAGE SOURCES CURRENT 0, ,VSRC1 %STOP ITERATION CRITERIA MET %%%%%%NON LINEAR DC ANALYSIS%%%%%% %ABSOLUTE COVERGENCE LIMIT:0.0001 %RELATIVE COVERGENCE LIMIT:0.001 %HAVING NODES:6 %NON LINEAR DC ANALYSIS ITERATION 0 %NODE VOLTAGE 1,10 2, 3, 4, 5, 6,12 %VOLTAGE SOURCES CURRENT 0, ,VSRC1 1, ,VSRC2 %STOP ITERATION CRITERIA MET 给直观感觉即可

24 例子:包含晶体二极管的简单电路 𝑖 𝑢 = 𝐼 𝑠 𝑒 𝑢 𝑉 𝑡 −1 , I s , V t 为与器件特性有关的常数
由固体物理知识： 𝑖 𝑢 = 𝐼 𝑠 𝑒 𝑢 𝑉 𝑡 −1 , I s , V t 为与器件特性有关的常数 例子:包含晶体二极管的简单电路 𝐺 𝑒𝑞 (𝑢)= 𝑑𝑖 𝑑𝑢 = 𝐼 𝑠 𝑉 𝑡 𝑒 𝑢 𝑉 𝑡 𝐼 𝑒𝑞 (𝑢)= 𝐼 𝑠 [(1− 𝑢 𝑉 𝑡 ) 𝑒 𝑢 𝑉 𝑡 −1] *SAMPLE INPUT FILE .DC 1e-3 1e-4 I R R 1 2 1e4 *D node0 node1 Vguess Is Vt D E %SAMPLE OUTPUT FILE %%%%%%NON LINEAR DC ANALYSIS%%%%%% %ABSOLUTE COVERGENCE LIMIT:0.001 %RELATIVE COVERGENCE LIMIT:0.0001 %HAVING NODES:2 %NON LINEAR DC ANALYSIS ITERATION 0 %NODE VOLTAGE 1, 2, %VOLTAGE SOURCES CURRENT %WARNING:NEXT VOLTAGE IS LARGER THAN THE PREV BY 0.3V %LIMITING CHANGE TO 0.3V %NON LINEAR DC ANALYSIS ITERATION 1 1, 2, ……… %NON LINEAR DC ANALYSIS ITERATION 10 1, 2, %NON LINEAR DC ANALYSIS ITERATION 11 2, %STOP ITERATION CRITERIA MET 𝑉 𝑔𝑢𝑒𝑠𝑠 =0.9𝑉 绝对误差/V 简单带过

25 遇到的问题 对于有迅速增长的伏安特性的器件，若 𝑉 𝑔𝑢𝑒𝑠𝑠 < 𝑉 𝑤𝑜𝑟𝑘𝑝𝑡 第一次迭代会获得非常大的终值。
取与前文同样的电路，但初始尝试点取0.6𝑉 从图中可见，不加限制后尝试点会达到近1.55𝑉， 𝐼 𝑒𝑞 达 𝐴量级 对于方程组的数值求解是一个麻烦。 也会使收敛变得极为缓慢 目前的简单缓解方法：限制猜测值增长，单步变化不超过0.3𝑉 if (dynelements[i].currvoltage - dynelements[i].lastitervoltage>0.3) dynelements[i].currvoltage = dynelements[i].lastitervoltage + 0.3; %输出中的提示 %WARNING:NEXT VOLTAGE IS LARGER THAN THE PREV BY 0.3V %LIMITING CHANGE TO 0.3V 𝑉 𝑔𝑢𝑒𝑠𝑠 =0.6𝑉 绝对误差/V 不加电压限制 需要经过65次迭代才可以达到1𝑚𝑉的精度

26 例子：求解RC电路的暂态过程 𝑡𝑖𝑚𝑒𝑠𝑡𝑒𝑝=0.1𝑠,𝑅=1Ω,𝐶=1𝐹,𝑈=1𝑉 共求解40步 与理论解：1− e −𝑡 符合很好
M矩阵： 解得节点2的电压 绝对误差/V 时间/s .DC 1e-4 1e-3 .TRAN 1e-1 40 R 2 1 1 C V 1 0 1 简单带过

27 RLC电路的暂态过程 *RLC CIR .DC 1e-3 1e-4 .TRAN 1e-1 100 V 1 0 1 R 1 2 1
*L port0 port1 L I0 L C 0.1s步长，求解100步 电路处在欠阻尼状态 %SAMPLE OUTPUT %%%%%%LINEARIZED TRANSIT ANALYSIS%%%%%% %TOTAL STEP:200 %STEP LENGTH:0.2 %HAVING NODES:3 %LT TIMESTEP,TIME:0,0 %NODE VOLTAGE 1,1 2, 3, %VOLTAGE SOURCES CURRENT 0, ,V ……… 节点2电压/V

28 示例：运用于大型电阻网格问题 对于无限大的网格，两节点间电阻理论值应为 R 2 (不妨设R=1Ω)
在这里用有限个电阻去逼近，设正方形网格一条边上有N个节点 电压源放在(N+1) N 2 点与 N+1 N 2 +1之间 #生成输入文件的Python3小脚本 n=int(input(‘N?’)) of=open('input.cir','w') of.write('.DC 1e-3 1e-4\n*N={}\n'.format(n)) of.write('*ROWS\n') for i in range(0,n): for j in range(i*n,(i+1)*n-1): of.write('R {} {} 1\n'.format(j,j+1)) of.write('*COLS\n') for i in range(0,n-1): for j in range(i*n,(i+1)*n): of.write('R {} {} 1\n'.format(j,n+j))

29 用大型电阻网格问题的求解作为性能评价 测试环境： Intel(R) 64 Compiler 16.0.1.146
关键选项/O3 /Qparallel Eigen 3.2.8 Intel MKL 11.3 Windows Server 2012 R2 Intel(R) Core(TM) i5-5200U 2.20GHz 8G RAM 对于大型规律性的电阻网络，LU算法直接求解效果最佳。 Intel MKL支持多核并行计算，大大提高求解性能， PardisoLU较Eigen 内置SparseLU至少快一倍。 BiCGSTAB等迭代求解器，因Eigen库的限制，矢量矩阵乘积无法并行化 实验证明对于大型规律性的电阻网络效果很不好，N=750时需要超过50s 耗时较多步骤： 1.方程组求解 2.输入文件的读取（容易优化） 3.初始稀疏矩阵的组装 点数N 6 8 10 12 20 40 60 电阻/Ω 解算时间(PardisoLU)/s 0.11 0.12 0.13 0.16 强调性能令人满意，144万节点，280万+个电阻，20s以内 点数N 80 120 240 320 500 750 1000 1200 解算时间(PardisoLU)/s 0.18 0.26 0.72 1.23 3.12 7.10 13.0 19.0 解算时间(SparseLU)/s 0.08 0.17 0.93 1.69 5.31 14.9 33.7 54.1 总体性能令人满意

30 直观分析计算结果 计算相距两个电阻处的等效电阻 在N=1000时加电压源：V 500499 500501 1 共一百万节点
电流为 A，相当于 Ω 与理论值 2− 4 𝜋 Ω≈ Ω(见文献[5])一致 直观分析计算结果 电势-节点位置图 （横断面） 电势-节点位置图（截取中心部分）

31 N=12时的结果与文献[6]的结果完全相同 稍作改变，求解对角节点间电阻 N=120,将输入文件的最后一行替换为
无限电阻网络问题收敛其实不快，数值计算并不容易, 对于包含六千个电阻的网络求解仍有0.18%的误差。 稍作改变，求解对角节点间电阻 N=120,将输入文件的最后一行替换为 V 得到 %VOLTAGE SOURCES CURRENT 0, ,V 相当于阻值为 Ω,与理论值 𝜋 2 Ω(见文献[5])相近

32 三角形无限网格 在输入文件中加入对角线项 of.write('*DIAGS\n') for i in range(0,n-1): for j in range(i*n,(i+1)*n-1): of.write('R {} {} 1\n'.format(j,n+j+1)) 计算得n=20时R= Ω,n=120时R= Ω,n=200时R= Ω 与理论值为( 8 3 − 𝜋 )Ω(见文献[5]),约 Ω相吻合 求解耗时不长，性能满足要求

33 直观分析计算结果 局部放大 N=200时各节点电势 有四万节点，求解仍可以在0.5s内完成 “等势线” 可以明显看出电流扩散的各向异性
沿00,11方向电阻较小

34 两边电阻不等的正方形网格 与精确值相符 N 60 80 120 240 精确值 R(1,0)/Ω 0.60831 0.60825
R(1,1)/Ω 与精确值相符 精确值由公式 R p; l 1 , l 2 = 0 𝜋 𝑑𝑥 𝜋 1− 𝑒 − 𝑙 2 𝑠 cos⁡( 𝑙 1 𝑥) sinh⁡(𝑠) cosh 𝑠 =1+ 1 𝑝 − 1 𝑝 cos⁡(𝑥) 给出(见文献[7])

35 进一步发展的思路 牛顿法求解非线性电路方程时时有不收敛的情况 目前时域动态分析使用的梯形法相对粗糙，可以改用高阶的，更稳定的其他算法
尤其是电路中包含二级管等指数伏安特性的器件时，若初值选取不当，迭代过程中会产生非 常大的电流和 𝐼 𝑒𝑞 , 𝐺 𝑒𝑞 ，导致求解失败。 目前处理方法：限制迭代时尝试工作点的变化 目前时域动态分析使用的梯形法相对粗糙，可以改用高阶的，更稳定的其他算法 三阶，四阶隐式Adams公式 𝑦 𝑛+1 = 𝑦 𝑛 + ℎ 12 (5𝑓 𝑥 𝑛+1 , 𝑦 𝑛+1 +8𝑓 𝑥 𝑛 , 𝑦 𝑛 −𝑓 𝑥 𝑛−1 , 𝑦 𝑛−1 )(三阶) 实现时域动态分析的时间步长(𝑡𝑖𝑚𝑒𝑠𝑡𝑒𝑝)的自动选取与自适应调整 根据各元件参数量级预估时间步长 用较大的步长尝试求解，若失败或与初值相差过大，再采用更小步长 根据本次迭代中计算结果与上次迭代结果的差决定 增加对三端非线性器件的支持，BJT,MOSFET… 其他细节上的性能优化 改进用户界面 概括为：解决不收敛问题，改进时域求解算法

36 https://software.intel.com/en-us/intel-mkl
获取完整源码请访问 关于Eigen请访问 关于Intel® MKL请访问

37 Q&A

38 玩笑：（仅限班内传播！） 往年竞赛中的类似作品： 本作品还是有一点竞争力的。。。 2014二等奖：魏家一，曾若兰《自制电路分析软件》
“金玉其外，败絮其中” 界面引人注目，实际并无干货 2013一等奖：刘峻宇，赵文君《用节点电压方法求算电阻网格的等效电阻》 “经过艰难的编写，我们给出的matlab代码如下” 思路是如何从简陋的Matlab仿真跳到复述参考文献[5,7] 的。。。 本作品还是有一点竞争力的。。。

39 参考文献 [1] Kincaid D R, Cheney E W. Numerical analysis: mathematics of scientific computing[M]. American Mathematical Soc., 2002. [2] 邱关源 , 电路(第五版) [M], 高等教育出版社 [3] Vladimirescu A. The SPICE Book[M], John Wiley&Sons, New York, 1994. [4] 张运华等 , 数值计算方法与算法(第二版)[M] , 科学出版社 [5] Atkinson D, Van Steenwijk F J. Infinite resistive lattices[J]. American Journal of Physics, 1999, 67(6): [6] Denardo B, Earwood J, Sazonova V. Experiments with electrical resistive networks[J]. American Journal of Physics, 1999, 67(11): [7] Cserti J. Application of the lattice Green’s function for calculating the resistance of an infinite network of resistors[J]. American Journal of Physics, 2000, 68(10):