CN1275318C - 可处理vlsi中树网混合供电结构的瞬态仿真方法 - Google Patents

Abstract

VLSI中树网混合供电结构的高速高精度瞬态仿真方法本发明属于VLSI制造技术领域，其特性在于它能够快速对任意比例的树网混合供电网进行瞬态仿真。该方法解决了常规仿真方法不能够处理不规则拓扑结构和多供电源的问题，并且具有很高的仿真精度。同时该方法完全考虑了电阻、电容和电感作为寄生参数的输入，能够达到很高的仿真精度。方法中还采用了基于不完全乔莱斯基分解的预处理技术和快速变元消除技术，方法的效率得到大幅度提高。实验结果证明了该算法具有极高的仿真精度和数值稳定性，能够在很短的时间内对极大规模的电源/地线网络进行仿真，同时由于与拓扑结构的无关性，该方法也适用于高频集成电路中的无源时钟线网、总线线网以及信号线网的仿真。

Description

可处理VLSI中树网混合供电结构的瞬态仿真方法

技术领域

超大规模集成电路(VLSI)供电网络的电特性仿真属于VLSI制造流程中物理版图设计领域。它是设计和优化集成电路电源地线网络必须依赖的一个核心子过程，该过程的效率和精度将直接影响到整个芯片的设计周期、成品芯片的性能和寿命，同时由于集成电路版图和印制电路版(PCB)的相似性，该方法也适用于高频PCB上的供电网络的仿真分析。

背景技术

1.芯片供电网设计的重要性

随着超大规模集成电路制造工艺的不断进步，成品芯片的特征线宽在不断缩小，同时芯片的集成密度和工作频率也在大幅度的提高。在这样的环境下，不良的电源/地线供电网络设计会极大地影响芯片的性能，如引起电路参数变化，在不同模块之间产生串扰，导致产生载流子效应和金属电迁移效应等。这些不良电磁现象会引起电路逻辑错误，降低芯片寿命，使芯片无法正常工作，所以芯片的供电问题已经成为深亚微米级集成电路设计中的一个重要问题。最近对深亚微米集成电路设计的调查表明，79％的设计具有电源网络的问题，需要在流片之前矫正，否则芯片将由于电源网络上的电压降过大而不能正常工作。为了解决这个难题，世界上顶级的电子产品设计自动化(EDA)公司纷纷提出了自己的解决方案，如Candance公司提供的VoltageStorm、SingalStorm、ElectronStorm仿真工具，Mentor公司提供的Eido、Eido RF、Eido Mach等。

2.芯片版图的布图模式

芯片的供电网设计在很大程度上受芯片布图模式的约束。所谓芯片的布图模式指的是芯片中的各个功能单元在硅片上的摆放规则，目前工业上使用的比较多的主要包括全定制和半定制两大类。全定制模式下，电路被分成一些子电路，每个子电路是一个模块，通常情况下，模块的外形和放置都没有限制。除了模块以外的芯片区域都是布线区。对于三层及以上布线工艺的版图设计，模块上面也可以走线，只是走线层数有所限制。因此，全定制布图模式除了要遵循最基本的几何设计规则(最小线宽、间距、覆盖、露头等)外，没有其他任何附加的物理限制，其中以手工设计布图和积木块自动布图BBL(Building BlockLayout)为代表(参见附图1)。半定制设计模式对功能单元的高度、电源线位置和单元引线端的引出方向都有一定的限制，并且单元也只能放置在规定的区域内或者必须按行进行放置。半定制布图模式又分为若干小类，如标准单元设计模式(参见附图2)、门阵列设计模式等(参见附图3)。

3.芯片供电网的物理形式

芯片供电网络在芯片上的物理形状主要有四种：网状、树状、树网混合状以及一般图。在网状供电网下，各个功能单元从网络的交叉节点上得到电压和电流，由于每个交叉节点至少有4条支路同时供电，因此可靠性比较高(参见附图4)；在树状供电网下，各个功能单元从供电树的分支节点上获得电压和电流，使用这种供电模式比较容易优化布线面积，因此在工业制造中用的比较多，但是它的可靠性不如网状供电网(参见附图5)；树网混合的供电网是前两者的结合，兼有两者的优点，但是在线网设计优化方面比较复杂，但是近一段时期，多层金属布线已经成为IC制造的发展方向，通过多层金属工艺，原先难以解决的树形结构的不可靠的缺点现在也能够通过多层子树供电得到一定的解决，所以可以预期，在芯片工艺向更高的集成度和更快的速度发展的将来，树网混合的供电模式一定会成为主流，这一点已经得到一些相关统计资料的证实，目前工业界很少有芯片是单纯采用一种供电模式实现的(参见附图6)；一般图的供电结构多由手工布图产生，没有一定的规律，在各种供电结构中最为灵活，但也最难以处理，往往用在对布线面积有较高要求的场合。

4.当前解决芯片供电问题的手段

要使设计出来的电源/地线网能够满足芯片中各个功能电路对电压和电流的要求，即在任何时刻都要能够给芯片上的各个功能单元提供稳定的满足一定大小要求的电压和电流，就必须对电源/地线网进行电磁性能上的定量分析，然后根据分析的结果进行相应的调整。对于低频芯片，通常只进行直流静态分析就能够满足对性能和可靠性的要求了。直流分析主要针对的是供电节点的电压和线网上的电流密度，一般通过提取线网的直流参数，如电阻率，线宽，线长等，然后根据基尔霍夫电压电流定律建立节点方程组，通过求解该方程组得到对电源/地线网电性能的评价，检查有没有电压和电流上违规点，有多少违规点，再根据这个评价修改原来布线的拓扑结构，不断重复这个过程，最后得到最优的布线结构，使之不仅能够满足没有违规点的要求，同时占用的布线面积也是最小的；而对于高频芯片，由于受到互连线分布参数的影响，线网的电磁性能不能单纯由直流电参数来描述。此时要设计出满足性能和可靠性要求的供电线网就必须对线网进行高频下的交流动态仿真计算，力求通过计算对电源/地线网的拓扑结构进行一系列的优化调整。描述高频下线网的电磁特性的模型有多种，如由传输线理论得出的电阻/电感/电容(RLC)模型、但是由不管使用什么计算模型，对高频下电源/地线网络的分析最终能够归结为求解一个大型线性方程组的数学问题。

5当前技术及其缺陷

近年来学术界提出了许多新的方法，通过求解由电源网络模型得到的大型线性方程组对电源网的供电情况作出仿真，从而进一步的进行优化。这样的方法很多，如多重网格法(Multi-grid method)、层次法(Hierarchical method)、预优共轭梯度法(Precondition ConiugateGradient-PCG)、层次模型降阶法(Hierarchical Model Order Reduction)，偏微分方程的交叉隐含法(Alternating Direction Implicate)等等，但是到目前为止，仍然没有一种算法能够稳定而快速地处理树网混合的供电结构。同时，在本领域，国际上许多著名EDA公司也在根据制造上出现的新问题，提供并改进他们的解决方案，通过他们自己开发的商业软件对电源/地线网进行分析，但是这些软件一般都存在求解速度慢、内存占用多、难以同时保证收敛性和计算精度等弊病，而且对于应用越来越普遍的树网混合的电源线网的结构，这些软件也缺乏相应的加速能力，有的工具目前甚至不能够仿真含有电感寄生参数的供电网络。

发明内容

本发明针对目前商用仿真软件不能够稳定快速而且准确地对树网混合结构的供电网进行仿真的现状，作出了改进，提出了一种能够快速对任意比例的树网混合供电网进行瞬态仿真的新方法。该方法解决了常规仿真方法不能够有效处理不规则拓扑结构和多供电源的问题，并且具有很高的仿真精度。同时该方法完全考虑了电阻、电容和电感作为寄生参数的输入，能够达到很高的仿真精度。方法中还采用了基于不完全乔莱斯基分解的预处理技术和快速变元消除技术，方法的效率得到大幅度提高。实验结果证明了该算法具有极高的仿真精度和数值稳定性，能够在很短的时间内对极大规模的电源/地线网络进行仿真，同时由于与拓扑结构的无关性，该方法也适用于高频集成电路中的无源时钟线网、总线线网以及信号线网的仿真。本发明的特征在于：它是以SunV880工作站为计算平台使用C语言在unix环境下实现的，它依赖于一下几个步骤：

1)程序扫描以Spice网表形式给出的电源网络的拓扑连接信息和电参数信息

2)计算机程序利用步骤1)得到的信息，建立起分别描述电导、电容、电感和电流源、电压源参数的数据表

3)程序通过扫描网表文件找到整个电源网络中的所有树形结构，并将其表示为一个层次式的节点列表，同时生成一个描述RL结构的列表

4)通过下述任何一种时域下的差分离散方法计算电源/地线网络中的电容元件以及电感元件的等效电导大小和等效瞬时电流源大小，即把电容元件和电感元件近似地转换为并联的电流源和电阻元件，同时设初始时刻的等效电流源是零

4.1)梯形欧拉法

对电容而言 $I_c^{t+h}=G_c\·V_c^{t+h}-(G_c\·V_c^t+I_c^t)G_c=\frac{2C}{h}$

对电感而言 $I_L^{t+h}=G_L\·V_L^{t+h}+(G_L\·V_L^t+I_L^t)G_L=\frac{h}{2L}$

其中h是用户输入的仿真步长，单位是秒

I_c ^t+h、V_c ^t+h分别是t+h时刻电容支路上的电流和端电压

I_c ^t、V_c ^t分别是t时刻电容支路上的电流和端电压；G_c是容元件离散化后的等效

电导

I_L ^t、V_L ^t分别是t时刻电感支路上的电流和端电压；G_L是电容元件离散化后的等效电导

若令 ${\mathrm{IS}}_c^t=-(G_c\·V_c^t+I_c^t)$ 则 $I_c^{t+h}=G_c\·V_c^{t+h}+{\mathrm{IS}}_c^t$

${\mathrm{IS}}_L^t=(G_L\·V_L^t+I_L^t)$ 则 $I_L^{t+h}=G_L\·V_L^{t+h}+{\mathrm{IS}}_L^t$

4.2)线性多步法

$I_c^{t+h}=G_c\·V_c^{t+h}-(G_c\·V_c^t+a_1{I}_c^t+a_2{I}_c^{t-h}+a_3{I}_c^{t-2h})$

$G_c=\frac{24C}{9h},a_1=\frac{19}{9},a_2=-\frac{5}{9},a_3=\frac{1}{9}$

$I_L^{t+h}=G_L\·V_L^{t+h}+(b_1{V}_L^t+b_2{V}_L^{t-h}+b_3{V}_L^{t-2h}+I_L^t)$

$G_L=\frac{9h}{24L},b_1=\frac{19h}{24L},b_2=-\frac{5h}{24L},b_3=\frac{h}{24L}$

若令 ${\mathrm{IS}}_c^t=-(G_c\·V_c^t+a_1{I}_c^t+a_2{I}_c^{t-h}+a_3{I}_c^{t-2h})$ 则 $I_c^{t+h}=G_c\·V_c^{t+h}+{\mathrm{IS}}_c^t$

若令 ${\mathrm{IS}}_L^t=(b_1{V}_L^t+b_2{V}_L^{t-h}+b_3{V}_L^{t-2h}+I_L^t)$ 则 $I_L^{t+h}=G_L\·V_L^{t+h}+{\mathrm{IS}}_L^t$

5)按照步骤3)中得到的森林节点层次表中自底向上的顺序，通过查找节点连接关系索引表和电容、电感元件底离散数据表，把连接在各个树节点上，需要改变的电容和电感替换成相应的电导和电流源

6)按照上述森林节点层次表中节点的顺序，从每棵树的底端开始，反复使用线性含源两端电阻网络的等效电路原理，把替换后的电流源和电导元件作层次式的合并，使供电网上所有的树都归并为根节点上的电导和电流源，它含有以下顺序进行的两个步骤：

6.1)把所有的电容和电感的等效电导合并到根节点上，这一步可以在步骤3)中进行

6.2)在每个时刻更新完所有电感和电容的等效电流源之后，再把更新后树节点上电容和电感的等效电流源再进行合并，最终得到根结点上的等效电流源

7)程序按照RL结构索引表给出的顺序用上述等效原理进行变换，消除电源网中的中间节点，生成一个以RL结构索引表的顺序为排列顺序的RL结构化参数表，其中m，n，o按下式计算：

                         m＝o·G，n＝1-m， $o=\frac{1}{G+G_L}$

其中，G是RL结构上的自带电导

      GL是RL结构中L等效后得到的电导

8)重新对网络中所有的节点进行排序，形成新的节点连接关系索引表

9)根据简化后的电源/地线拓扑结构，生成计算矩阵和相应的电流源向量，形成供电节点电压、电流的线性方程组，它依赖于以下步骤：

9.1)把接地节点标记为零节点，把接地节点从节点连接关系索引表中出去。

9.2)基于基尔霍夫定律和欧姆定律得到供电节点电压、电流的线性方程组，如下式所示 $\underset{j\∈P}{\mathrm{\Σ}}{g}_{\mathrm{jk}}(v_j-v_k)+\underset{i=1,i\≠k}{\overset{N^{\′}}{\mathrm{\Σ}}}{g}_{\mathrm{ik}}(v_i-v_k)=i_k$

形式为G·V＝I

G为一个N×M的矩阵，称为计算矩阵，其中： $\left\{\begin{array}{c}{g}_{\mathrm{ij},i\≠j}=g_{\mathrm{ij}}\\ g_{\mathrm{ii}}=-(\underset{j\∈P}{\mathrm{\Σ}}{g}_{\mathrm{ij}}+\underset{j=1,j\≠i}{\overset{N^{\′}}{\mathrm{\Σ}}}{g}_{\mathrm{ij}})\end{array}\right.$

I为一个N维列向量，其中的元素代表每一个节点到地的电流值，

其中N’是供电网中网结构的节点数，M是网结构的支路数

    v_k是矩阵中一行所代表的一个节点的电压值

    v_i是其他节点的电压值

    v_j是电压源的电压值

    P是和电压源相连的节点的下标集合

    g_ik是节点k和其他节点i相连的电导值

    g_jk是节点k和电压源相连的电导值

10)通过调用基于不完全乔莱斯基分解的预优共轭梯度算法求解引擎得到简化后节点供电电压的分布

$V_{\mathrm{RLnode}}^t=m\·V_{\mathrm{Gnode}}^t+n\·V_{\mathrm{Lnode}}^t-o\·{\mathrm{IS}}_L^t$

11)在得到所有RL结构的两端节点在某个时刻的电压之后，根据RL结构化简参数表，根据下式逐一恢复RL节点此时刻的电压值，再按中间节点链进行四则运算就可以恢复所有中间节点的电压

${V^t}_{\mathrm{middle}}=c\·V_{\mathrm{top}}^t-{\mathrm{eI}}_{\mathrm{middle}}^t$

$V_{\mathrm{bottom}}^t=a\·V_{\mathrm{middle}}^t+d\·(I_L^t-I_{\mathrm{bottom}}^t)$

a＝G_L·d，b＝1-a，c＝G_top·e， $d=\frac{1}{G_L+G_{\mathrm{bottom}}},e=\frac{1}{G_{\mathrm{top}}+G_{\mathrm{middle}}}$

其中G、G_L的定义如前述。

12)根据森林节点层次表中节点的顺序和树结构化简参数表，按下式还原所有的树节点电压值

13)根据以上得到的全部节点电压分布，重新计算电流源向量

14)如果没有达到预定的仿真时间，跳转到步骤10)继续执行。

本发明和已有软件的性能比较

1.计算精度比较

参见附图19，其中带有锯齿状的曲线是商用软件Hspice给出的一个供电网中某个供电节点电压值的仿真结国，比较光滑的为本发明程序给出的仿真结果。可以看出，不仅本发明给出的仿真结果和Hspice相比误差很小，而且结果也更加光滑。

2.计算效率比较

下表是一个本发明和商用软件Hspice以及ICCG仿真程序在仿真时间上的比较，可以看出本发明的计算时间几乎只有ICCG的十分之一，而和Hspice相比则几乎快了两个数量级。

其中TME的耗时为本发明程序的耗时。

节点总数	树节点百分比(％)	ICCG计算耗时(秒)	TME计算耗时(秒)	Hspice计算耗时(秒)	加速比
							比ICCG	比
					2,653	25			6.23	1.85	未测	3.37	未测
9,313	25	43.88	7.31	未测			6.00	未测
					41,413	25			301.40	51.80	8597.6	5.82	165.97
73,171	25	970.10	94.74	16945.75			10.24	178.86
					154,843	25			3639.29	219.39	22594.34	16.59	102.98
208,393	25	7737.32	322.55	未测			23.99	未测
					357,007	25			19622.79	544.64	未测	36.03	未测
660,157	25	50157.83	1104.36	未测			45.42	未测
					1,011,0	25			98038.98	1706.75	未测	57.44	未测

附图说明

图1：积木块自动布图模式(BBL)示意

图2：标准单元布图模式示意

图3：门阵列布图模式示意

图4：树形供电网示意

图5：网形供电网示意

图6：树状和网状供电网示意

图7：一个2×2网状供电结构的仿真模型

图8：离散化后2×2网状供电结构的仿真模型

图9：电源网中树形结构层次扫描过程示意图

图10：诺顿等效电路原理图

图11：戴维南等效电路原理图

图12：诺顿等效电路于戴维南等效电路的转换

图13：RL结构及其化简示意图

图14：电容和电感元件差分离散等效示意图

图15：树形结构中二阶元件的替换示意图

图16：树形结构中电导元件的合并等效示意图

图17：树形结构中电流源合并等效示意图

图18：不完全乔莱斯基分解的过程示意图

图19：本发明的仿真结果和商用仿真软件Hspice的比较

图20：仿真程序流程图

具体实施方式

本方法已经使用C语言在unix环境下开发并实现为一个高性能的仿真程序。该仿真程序以Spice网表形式给出的供电网络电参数描述文件作为输入，通过用户交互地输入仿真时间和仿真步长等参数，最终给出每一个时刻供电网络上各个供电节点的电压电流分布作为输出。

                      具体的实现环境如下：

硬件系统：Sun V880 WorkStation with 4 CPU at 800MHz，2GB Memory

操作系统：Sun Solaris 5.0

程序设计语言：ANSIC

编译器：GNU gcc版本：3.3.1

程序调试工具：GNU gdb版本：6.0

辅助生成工具：GNU make版本：3.79.1

下面介绍该仿真程序的工作流程：

1.程序扫描以Spice网表形式给出的电源网络描述文件，得到电源网络的拓扑连接信息和电参数信息。

1.1Spice网表文件简介

Spice网表文件是通过对电路中各个节点进行编号，使用分类记录各个元件的引脚编号和相应的参数值的方法来描述电路的拓扑结构和电特性的。对于电源网络等效电路模型中的两端元件：电阻、电感、电容、电压源、电流源，该文件具有以下描述：

元件名称	描述语法
		电阻	元件类型 元件编号 左端节点号 右端节点号 电阻值
电感	元件类型 元件编号 左端节点号 右端节点号 电感值
		电容	元件类型 元件编号 左端节点号 右端节点号 电容值
直流电压源	元件类型 元件编号 左端节点号 右端节点号 电压值
		分段电流源	元件类型 元件编号 左端节点号 右端节点号 端电压类型 端电压大小PWL{时刻，电流大小}

其中电阻、电感、电容、电压源、电流源的元件类型分别是：R、L、C、V和I分段电流源中端电压类型有两种：

如果是直流电流端电压类型为DC

如果是交流电流端电压类型为AC

{}表示括号内的内容重复次数可以是大于等于1的任意次，但是不可以为0次文件中各个物理量值的默认的量量纲分别是：欧姆(Ω)、亨利(H)、法拉(F)、福特(V)和安培(A)。如果所使用的物理量大小不是默认量纲的大小，则在物理量的具体数值之前显示的加以说明，下面是常用倍乘系数的描述：

系数	f	p	n	u	m	MI
							大小	10-15	10-12	10-9	10-6	10-3	25.4×106
系数	K	MEG	G	T	DB
							大小	103	106	109	1012	20logAi

还需要说明一下的是文件中的分段电流源。分段电流源在所描述的时刻上的输入/输出电流值等于文件中给出的值，而在两个给出时刻之间的电流值则使用线性插值的方法给出。网表文件中注释行使用引导前缀：^**必要的说明信息和电路的部分统计信息都可以出现在注释行中，同时在文件开始的时候必须交代进行时域分析的时间步长和总共的分析时间使用以下语法描述：.tran时间步长总分析时间

以下是一个2×2纯网格结构电源网的网表文件内容，其对应的实际供电网参见附图7：

^**supply voltage：5

^**#row：2

^**#RC/RLC section per row：2

^**#p/g strips per row：2

^**#max depth of a tree：0

^**#max branches in each tree node：0

^**This is a 2×2 power network

.tran 0.05ns 6ns

^**row：1

^**p/g strip connection

R1 1 6 0.76032

R2 3 8 0.76032

R3 5 10 0.76032

^**end of p/g strip connection

R4 1 11 0.0018128

I1 10 DC 0PWL(0 0.01ns 0.0525mA 2ns 0.5285mA 3ns 0.5145mA 4ns 0.35mA 5ns 0.0315mA)

C1 10 1.12e-10

R5 1 2 0.18656

L1 2 3 1.46387e-12

C2 3 0 1.67e-10

I2 3 0 DC 0 PWL(0 0.01ns 0.294mA 2ns 0.665mA 3ns 0.3675mA 4ns 0.2555mA 5ns 0.0315mA)

R6 3 4 0.18128

L2 4 5 1.43701e-12

C3 5 0 1.66e-10

I3 5 0 DC 0 PWL(0 0.01ns 0.231mA 2ns 0.728mA 3ns 0.5775mA 4ns 0.4235mA 5ns 0.0245mA)

R7 5 11 0.0018304

^**row：2

R8 6 11 0.0017776

I4 6 0 DC 0 PWL(0 0.01ns 0.1575mA2ns 0.637mA 3ns 0.651mA 4ns 0.329mA 5ns 0.0175mA)

C4 6 0 1.94e-10

R9 6 7 0.17952

L3 7 8 1.35643e-12

C5 8 0 1.67e-10

I5 8 0 DC 0 PWL(0 0.01ns 0.189mA 2ns 0.6405mA 3ns 0.6335mA 4ns 0.371mA 5ns 0.028mA)

R10 8 9 0.19008

L4 9 10 1.41015e-12

C6 10 0 1.43e-10

I6 10 0 DC 0 PWL(0 0.01ns 0.2345mA 2ns 0.476mA 3ns 0.5075mA 4ns 0.3255mA 5ns0.007mA)

R11 10 11 0.0018304

^**total num of nodes：11

^**total num of trees：0

^**supply voltage

V1 11 05

.end

1.2元件参数表的建立

为了能够利用网表文件给出的所有电源网络的连接信息和电参数信息，仿真程序需要建立起分别描述电导(电阻的倒数)、电容、电感和电流源、电压源参数的数据表。同时为了简化以后的计算复杂性，还会针对供电网络中的一些比较特殊的结构建立相应的数据表。

1.2.1参数表的建立过程

1)由于供电网络中的电导、电容、电感、电压源的大小不随时间改变，所以只需要扫描一次网表文件就可以得到，而电流源表中的电流源大小一列的数据由于供电网络的负载各个时刻吸取的电流不同因此是一个动态变化的量，所以在各个不同的分析时刻需要进行更新。

       电感参数表举例，其他表类似

电感编号	头节点	尾节点	电感大小(H)
				L1	2	3	1.46387e-12
L2	4	5	1.43701e-12
				L3	7	8	1.35643e-12
L4	9	10	1.41015e-12

2)通过扫描网表文件仿真程序还要生成一个节点连接关系索引表，通过该表要能够得出每个节点和那些元件相连接，根据给出的元件编号程序就能够在不同的元件参数表中查出该元件的具体电参数。

                           节点连接关系索引表

节点编号	连接电阻	连接电容	连接电感	连接电流源	连接电压源
						1	R4，R5，R1	C1	无	I1	无
2	R4	无	L1	无	无

3)对于供电网络中的树形结构，仿真程序将通过扫描网表文件找到整个电源网络中的所有树形结构，并将其表示为一个层次式的节点列表。这样作的原因是为了提高后面计算时的查找所需数据的效率，相关内容参见3．2节。整个树的生成过程严格按照自底向上的方法进行归约，由于实际上电源网络中往往有多颗树同时存在，所以程序实际上会生成一个森林的节点层次表。假设电源网络中有如下的一个森林结构，如附图9所示，程序将能够得到如下4个表

表1：第一次被剪枝掉的节点

3

4

8

9

10

11

14

15

表2：第二次被剪枝掉的节点

7

6

13

表3：第三次被剪枝掉的节点

表4：第四次被剪枝掉的节点

通过以上的剪枝，电源网中所有的树形结构都被上推到根节点，所有的剪枝表组合到一起就形成森林的节点层次表。

4)在进行森林的节点层次表的同时，可以同时进行一些预处理的操作，目的也在于提高后期计算的效率，详细处理方法参见3.2.2节

5)由于电源网络中还存在一种比较特殊的连接结构，RL结构，如附图13所示，这种结构将能够通过后面所介绍的方法简化计算量，详细内容参见7.2节，因此仿真程序在扫描网表文件的时候也同时生成一个描述所有RL结构的表。对附图7有：

         RL结构索引表举例

RL结构编号	R端的节点号	L端的节点号
			RL1	1	3
RL2	3	5

2.特殊元件的处理

得到电源网所有的元件连接信息之后，程序要根据算法的要求对一些特殊的元件进行处理，并对整个电源网络的连接结构进行一定的改动，以加快后面的计算速度。

2.1什么是特殊元件

特殊元件指供电网络描述中所有的电容元件和电感元件

2.2为什么特殊

含有电容和电感的电路是二阶电路，如果不对电感和电容进行任何处理，则对所有节点用基尔霍夫定律所建立的电路方程组将会是一个含有节点电压的二阶常微分线性方程组；即使引入电感电流作为方程组的辅助变量，将系统方程组降阶为一个一阶的常微分线性方程组，由于通常情况下电源网络的等效电路中节点的个数是很大的，要直接求解这样的一个大规模线性微分方程组的代价是很高的，同时要在时域上对等效电路进行分析，不仅是要计算一个时间点上所有节点的电压值，而且需要计算许多时间点上这些节点的电压值，就算使用数值的方法去求解这个微分方程组，效率也是很低的。

2.3本发明中对特殊元件的处理方法

发明中我们采取的方法是使用有限差分的方法(即将原来微分运算符得出的精确值用一系列的简单函数的组合表示)对电容和电感元件的特性方程进行时域离散化。这样处理之后，在每个时刻，我们可以将电容和电感中只与当前时刻电压电流相关的系数处理成等效电阻，而把和前一个或者前几个时刻电压电流相关的系数处理成为随时间变化的电流源。于是在每一个时刻，我们就可以得到一个只含有电阻和电流源的网络，而求解这个大型的电阻性网络可以使用非常高效的算法。

下面我们介绍几种在我们的程序中实现的，具有不同精度和计算复杂度的离散化方法。

1)梯形欧拉法

首先考虑描述线性电感和电容伏安特性的微分方程：

$L\·\frac{{\mathrm{di}}_L}{\mathrm{dt}}=v_L---\left(2.1\right),C\·\frac{{\mathrm{dv}}_c}{\mathrm{dt}}=i_c---\left(2.2\right)$

其中i_L是流过电感元件的电流，v_L是电感元件两端的电压，L是该元件的电感值，单位为亨利；i_c是流过电容元件的电流，v_c是电容元件两端的电压，C是该元件的电容，单位为法拉。通过将(2.1)和(2.2)使用梯形欧拉方法进行离散化处理可以得到下面的近似公式：

$I_c^{t+h}=G_c\·V_c^{t+h}-(G_c\·V_c^t+I_c^t)G_c=\frac{2C}{h}---\left(2.3\right)$

$I_L^{t+h}=G_L\·V_L^{t+h}+(G_L\·V_L^t+I_L^t)G_L=\frac{h}{2L}---\left(2.4\right)$

以上两式中h是用户输入的仿真步长，单位是秒，I_c ^t+h和V_c ^t+h是t+h时刻电容支路上的电流和端电压，I_c ^t和V_c ^t是t时刻电容支路上的电流和端电压，G_c可以看成是电容元件离散化后的等效电导；而I_L ^t+h和V_L ^t+h是t+h时刻电感支路上的电流和端电压，I_L ^t和V_L ^t是t时刻电感支路上的电流和端电压，G_L可以看成是电感元件离散化后的等效电导。

2)线性多步法

如果采用四阶隐式亚当斯线性多步法对(2.1)和(2.2)进行离散化处理，则可以得到以下近似公式：

$I_c^{t+h}=G_c\·V_c^{t+h}-(G_c\·V_c^t+a_1{I}_c^t+a_2{I}_c^{t-h}+a_3{I}_c^{t-2h})G_c=\frac{24C}{9h},a_1=\frac{19}{9},a_2=-\frac{5}{9},a_3=\frac{1}{9}---\left(2.5\right)$

$I_L^{t+h}=G_L\·V_L^{t+h}+(b_1{V}_L^t+b_2{V}_L^{t-h}+b_3{V}_L^{t-2h}+I_L^t)G_L=\frac{9h}{24L},b_1=\frac{19h}{24L},b_2=-\frac{5h}{24L},b_3=\frac{h}{24L}---\left(2.6\right)$

上面两式中电压、电流变量及其上下标含意同(2.3)、(2.4)中的相同，只是增加了和t-h时刻、t-2h时刻相关的电压、电流变量以及相应的系数。

3)如果将(2.3)-(2.6)中带负号的括号部分记为IS_c ^t，带正号的括号部分记为IS_L ^t，则经过以上离散化后的电感和电容可以看作是由电阻和时变电流源组成的线性元件，因为它们都可以用标准表达式(2.7)，(2.8)表示：

$I_c^{t+h}=G_c\·V_c^{t+h}+{\mathrm{IS}}_c^t---\left(2.7\right)$

$I_L^{t+h}=G_L\·V_L^{t+h}+{\mathrm{IS}}_L^t---\left(2.8\right)$

(2.7)和(2.8)意味着我们可以将电容元件和电感元件近似地转换为并联的电流源和电阻元件，如附图14所示。

2.4对特殊元件进行近似后的仿真精度

通过对梯形欧拉方法进行泰勒展开分析可以发现其截断误差为O(h³)，但由于存在前后多步的迭代计算，其总体误差将可以大到O(h²)；而对于四阶隐式亚当斯线性多步法而言，其截断误差为O(h⁵)，由于也存在多步的迭代计算，所以其总体误差可以达到O(h⁴)，其中h为用户输入的仿真步长。一般而言要到达更高的精度，在计算每一个时刻电压、电流值的时候就需要知道更多以前时刻的电压、电流值，这不仅增加了算法的时间开销，也增加了算法的内存开销，因而在实际的计算中，应该根据不同的仿真精度要求，选择使用不同的差分格式进行离散化。从以上的差分模型的误差分析可知，仿真步长的选择对仿真精度有很大的影响。一般而言，h选择得越小计算结果的精度就越高，但是对于有的不具有绝对稳定性的差分方法，为了保证其逼进的收敛性必须对仿真步长进行约束，如果h过小则整个计算就会很快发散掉，得不到正确结果。在本发明中，我们给出了梯形欧拉离散化方法和h的选择无关同时绝对稳定的证明。该证明解释了为什么对梯形欧拉方法而言h在纳秒(ns)数量级之下还是可以随意选择的；另一方面，h也不是选得越小越好，受到分布参数C和L的影响，过小的h不仅会导致计算中舍入误差增加，而且还会导致(2.3)-(2.6)中电容的等效电导变的很大、电感的等效电导变得很小，这将给后面的化简带来计算上的麻烦，所以在对常规的分析精度和计算速度进行了一定的折衷之后，我们实现的分析软件中使用的是梯形欧拉的差分格式，同时取h为0.05ns(0.05×10-9s)。对于实际的工程应用，由于h选择的比较小，O(h²)的精度已经非常精确了。

2.5梯形欧拉方法的稳定性

对于描述电感和电容的电压电流关系的微分方程(2.1)和(2.2)而言，由于方程右端不显示地含有微分变量，所以需要写出特殊的误差方程。考虑(2.2)式经梯形欧拉格式离散所对应的误差方程：

$\left|\frac{\mathrm{\δ}{I}_c^{t+h}}{\mathrm{\δ}{I}_c^t}\right|=\left|\frac{G_c\·\mathrm{\δ}{V}_c^{t+h}-(G_c\·\mathrm{\δ}{V}_c^t+\mathrm{\δ}{I}_c^t)}{\mathrm{\δ}{I}_c^t}\right|\≤p$

其中δI_c ^t+h、δI_c ^t是t+h和t两个时刻电容电流上的数值误差，δV_c ^t+h、δV_c ^t是t+h和t两个时刻电容电压上的数值误差，p是需要满足的误差收敛速度，如果能够达到10^-2或10^-3就非常理想了，通过化简上面的不等式我们可以得到

$\frac{2C\left|x\right|}{h(p+1)}\≤\left|\mathrm{\δ}{I}_c^t\right|\≤\frac{2C\left|x\right|}{h(1-p)}$ 其中 $\left|x\right|=|\mathrm{\δ}{V}_c^{t+h}-\mathrm{\δ}{V}_c^t|$

从上式可以看出，当单步的数值误差在阀值 以下时是会出现积累的，但是一旦超过阀值就会按预设的收敛速度衰减，所以我们可以通过选取一定的阀值来能够保证算法的稳定性。由于我们是通过求解节点电压方程得到电容上的电压值的，而求解线性方程组的数值方法可以保证|x|≤M，所以要满足 $\frac{2C\left|x\right|}{h(p+1)}<\mathrm{\&eta;},$ η为累积误差上限，只需要有 $M<\frac{\mathrm{\&eta;h}(p+1)}{2C}$ 就可以了。当取η＝10^-5，分析要求h的大小在ns数量级(10^-9)时，由于分布电容C的数量级在10^-10左右，所以只要求M＜10⁴就可以了，这个大小一般的数值计算方法很容易达到。如果要求更小的误差累积上限，对于迭代法只需要增加迭代次数提高数值精度就可以了。对于(2.1)式，进行类似的分析也能得到相同的结论，所以我们说采用梯形欧拉方法进行离散计算的稳定性和时间步长h的选取无关，只要保证节点方程组的求解精度，使用梯形欧拉的离散方法在芯片电源网络的分析计算中一定不会出现迭代发散的情况。

2.6得到特殊元件的离散参数表

由于电感和电容元件离散之后得到的等效电导大小是不随时间变化而变化的，而电流源的大小却是一个随时间变化的量，因此在每个不同的时刻需要进行刷新，而对电流源的刷新需要使用到前一个时刻或是前几个时刻的电压值，所以在扫描输入网表文件的阶段，我们可以根据用户输入的仿真步长h和各个特殊元件的电参数(电容量或电感量)，一次性的计算出特殊元件的离散化参数。初始时刻的等效电流源设为是0。

电容元件的离散数据表(假设用户输入仿真步长为0.05ns)

元件编号	等效电导Gc	时刻	等效电流源ISc t
				C1	4.48	0ns	0mA
C2	6.68	0ns	0mA

3.树形结构的处理

3.1电容、电感元件的替换

完成离散化的计算之后，程序并不直接进行下一步的计算，而是先将电源网中所有的电容元件和电感元件根据2.6中得到的元件离散数据表进行替换，替换成相应的电导和电流源。

如附图15所示

3.2替换后对树形结构的合并

3.2.1为什么需要进行树形结构的合并

我们知道使用迭代法求解电路的节点方程的时候，算法的复杂度近似的与节点数的二次方成正比，要是能够大幅度地减少电路中的节点数，那就肯定可以极大地提高求解的速度。考虑在每一个时间点上供电网上的节点电压值构成了系统电路方程的解向量，而这个解向量只是一个1×n的矩阵，它的秩只有1，也就是说在n个解元中有n-1个解元是线性相关的。如果我们能够知道这些解元之间的关系，那么我们只要解出其中一个解的值，就可以得到所有解的值。然而要想得到这些解元之间得所有关联系数并不容易，如果使用高斯消元法将系数矩阵通过一系列的运算化为上三角或是下三角矩阵的话将会带来很大的运算开销。虽然如此，但如果我们能够比较容易的得到部分耦合比较紧的解元之间的系数关系，将这部分节点从待求方程组中去掉，只留下一个特征变元，则通过求解节点数大大化简后的方程组，我们就可以得到特征变元的值，再通过特征变元和归并节点之间的系数关系，便能够得到原方程组的完整解。这种思想正是我们实现的快速变量消除算法的核心。

基于变量消除的方法能不能行之有效关键在于能不能在很小运算开销的情况下找到紧密耦合解元的相关系数以及能不能在知道特征变元值的情况下很快地计算出被约简节点的电压值。通过分析电源地线网的等效电路模型，我们发现，使用电路理论中等效电路的方法正好能够高效地约简树网混合供电网中大量的中间节点和树节点，求出其相关系数，同时只要能够设计合适的数据结构，重建节点电压值的过程仅仅需要少量的运算。

3.2.2树形结构的合并过程

1)合并原理

一个经过适当限制的线性含源二端电阻网络R可以用一个独立电流源与一个独立电阻器的并联组合来等效，这种等效成为诺顿等效。其中，独立电流源电流是源网络的短路电流，即该网的负载电阻为0时的电流；而等效电导则是网络的内电导，即将网络R中所有独立电压源和独立电压全部置0后计算得到的网络的端口电导，参见附图10；一个经过适当限制的线性含源二端电阻网络R可以用一个独立电压源与一个独立电阻器的串联组合来等效，这种等效称为戴维南等效。其中，等效电压源的电压是原网络R的开路电压，即端口不接负载时的电压；而等效电阻则是原网络的内电阻，即将网络R中所有独立电压源和独立电压全部置0后计算得到的网络的端口电阻，参见附图11。

一个电阻网络的诺顿等效电路和戴维南等效电路可以相互转化，参见附图12。

a.戴维南等效电路中的电阻值和诺顿等效电路中的电导值相乘等于1G_ab·R_ab＝1

b.诺顿等效电路中的独立电流源的大小等于戴维南等效电路中独立电流源的值和诺顿等效电路中电导值的积。I_sc＝U_oc·G_ab

对于供电网中的所有树形结构，我们可以从每棵树的底端开始作归并运算。我们的目标是通过反复使用上面提到的等效电路原理，将供电网中所有的树都归并成为根结点上的电导和电流源。由于在1.2.1的第三个小步已经得到了整个供电网中森林结构的层次列表，所以此时没有必要重新扫描整个供电网络，只需要按照森林节点层次表中给出的节点顺序进行处理就可以了。

从效率上考虑，我们的化简将分成两步进行。因为电导的合并不受电流源大小的影响，而电流源的合并却和各个电流源相关，所以我们可以在扫描网表文件生成森林结构的层次列表时就提前合并树枝上的电导，将所有电导(包括电感和电容的等效电导)合并到根节点上，这个阶段的化简结果能在后面每一个时间点的计算中使用到，可以看成是一个预处理步骤。由于在每一个时刻需要同时得到此时根节点上的等效电导和等效电流源的值才能将整棵树的所有节点用诺顿等效电路等效掉，所以在每个时刻更新完所有电感和电容的等效电流源大小之后，我们必须将更新后树节点上电感和电容的等效电流源再次进行合并，最终得到根结点上新的等效电流源。以上提到的两个化简步骤都必须使用森林结构的层次列表，否则就需要每次重新作一次全网表的扫描。因此，一个优秀的层次列表设计可以大大的提高计算的效率。

2)电导的预合并

下面我们先来看看如何得到树结构根结点上的等效电导。本化简过程是层次式的，严格按照1.2.1中第三小步生成的森林的节点层次表进行。考虑一个树枝上节点的化简过程，参见附图16。附图中G₁、G₂…G_k是以前各次化简后得到的连接在BottomNode节点上各个树枝的对地电导，G_c、G_L、G_top分别是电容的等效电导、电感的等效电导以及树枝的金属连线电导。首先对连接在BottomNode节点上的所有电导进行合并，得到BottomNode节点的合并对地电导G_bottom；然后再对连接在MiddleNode节点和地之间的电导进行合并，消去BottomNode节点得到G_middle；最后对连接在TopNode节点和地之间的电导进行合并，消去MiddleNode节点得到G^*；其中

$G_{\mathrm{bottom}}=\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{G}_i+G_c,G_{\mathrm{middle}}=\frac{G_L\&CenterDot;G_{\mathrm{bottom}}}{G_L+G_{\mathrm{bottom}}},G^{*}=\frac{G_{\mathrm{top}}\&CenterDot;G_{\mathrm{middle}}}{G_{\mathrm{top}}+G_{\mathrm{middle}}}---\left(3.1\right)$

不断重复以上使用的合并方法就可以在树的根节点得到总的等效电导。在合并的同时，仿真程序会按照归约的顺序根据每一个树枝上的G_bottom、G_middle、G^*，作为后来进行电流源合并计算时的参数。同时程序会计算每一个树形节点上a，b，c，d，e五个参数，生成一个以森林节点层次表中节点的顺序为排列顺序的化简参数表。其中a，b，c三个参数将用于电流源的后期合并(参见接下来的第3)小节)；而d，e两个参数将用于重建被消除的节点(参见7.1节)。

                          一个典型的树结构化简参数表

节点序号	参数a	参数b	参数c	参数d	参数e	Gbottom	Gmiddle	G*
									1
2

α＝G_L·d，b＝1-a，c＝G_top·e， $d=\frac{1}{G_L+G_{\mathrm{bottom}}},e=\frac{1}{G_{\mathrm{top}}+G_{\mathrm{middle}}}---\left(3.2\right)$

3)电流源的后期合并

下面我们再来看如何在每个时间点上求得根节点的等效电流源的值，参见附图17。附图中IS^t _L、IS^t _c、I^t _node分别是电感、电容的在t时刻的等效电流以及功能模块在t时刻的吸取电流，其中只有IS^t _c的方向是取从地流入节点为正方向，其他都取流向地的方向为正方向；I^t ₁、I^t ₂…I^t _k是C节点上子树枝在以前各步化简后得到的t时刻的对地等效电流源，都取流入地的方向为其正方向。我们首先对连接在BottomNode节点上的所有电流源进行合并，得BottomNode节点的等效对地电流源I_bottom ^t

$I_{\mathrm{bottom}}^t=\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}{I^t}_j-{\mathrm{IS}}_c^t+I_{\mathrm{node}}^t---\left(3.3\right)$

然后我们将节点到地的电流源和电导从诺顿形式变换成戴维南形式，得到V1和V2

$V1=G_L\&CenterDot;{\mathrm{IS}}_L^t,V2=G_{\mathrm{bottom}}\&CenterDot;I_{\mathrm{bottom}}^t---\left(3.4\right)$

再合并的电压源V1、V2，电导G_L、G_bottom，反变幻到诺顿形式就得到MiddleNode到地的等效电导G_middle和到地的等效电流源I^t _middle，由于G_middle已经在(3.1)式中计算得出，所以此步仅需计算I^t _middle，可以证明

$I_{\mathrm{middle}}^t=a\&CenterDot;I_{\mathrm{bottom}}^t+b\&CenterDot;{\mathrm{IS}}_L^t---\left(3.5\right),$

其中a和b是(3.2)中定义的系数。

最后我们对连接在Middle节点上的对地电导和对地电流源进行从诺顿形式到戴维南形式的转换，合并G_top之后再反变换到诺顿形式，这样就可以得到TopNode节点上t时刻的等效对地电导G^*和等效对地电流源I^t _*，其中 $I_{*}^t=c\&CenterDot;I_{\mathrm{middle}}^t---\left(3.6\right)$

从3.3、3.5、3.6中我们可以看到，求出根节点上的等效电流源大小的过程实际上就是对电路网中所有更新后的电流源进行线形组合运算的过程，只不过该线形组合必须按树的层次来进行，如果我们在预处理阶段计算并存储下每一个树枝的组合系数，则在以后的每一个时间点求各个树结构根结点的等效对地电流源大小只需对所有树结构从最底层到根按层次遍历一次就可以了。

4.RL结构体的消除

消除所有树形结构之后，原来的电源网络转换为一个单纯的网状结构，但是这个网中还有大量的RL结构，可以通过消除进一步提高分析效率。在这一步，仿真程序将按照RL结构索引表给出的顺序按算法进行变换，消除电源网中的中间节点。化简过程中仿真程序将生成一个以RL结构索引表的顺序为排列顺序的RL结构化简参数表，这些参数将在后面的恢复过程中用到，参见7.1节。

                     一个典型的RL结构化简参数表

其中参数计算如下：m＝o·G，n＝1-m， $o=\frac{1}{G+G_L}$

G是RL结构上自带的电导，G_L是RL结构中L等效得到的电导，参见附图13通过参数的计算，程序可以得到RL结构的等效电导G_RL和等效电流大小I^t _RL

G_RL＝m·G_L， $I_{\mathrm{RL}}^t=\frac{{\mathrm{IS}}_L^t}{m}$

5.计算矩阵和电流源向量的建立

5.1节点方程组的矩阵形式

使用梯形欧拉法将电源地线网的RLC模型进行时域简化之后，我们就要对每个时间点上的直流网络进行分析，而要求解这样的直流网络，我们必须对其中每一个节点使用基尔霍夫电压定律和电流定律列出节点方程组，同时考虑支路上电压和电流必须符合的欧姆定律从而得到关于供电节点电压、电流的线性方程组。根据电网络图理论可知，任何由电路元件构成的电路网的节点方程都能写成矩阵方程BΛB^T·V＝I(5.1)的形式，其中A矩阵是拓扑结构和该电路网相同的图的关联矩阵，如果电路网是具有N个节点，M条支路的电路，A矩阵就是一个N×M的矩阵，且有

(5.1)式中的Λ为一M×M的对角矩阵，对角线上的元素是相应编号支路上的电导值；V是N×1的向量，其中的元素代表相关节点的电压值；I也是N×1的向量，其中的元素代表每一个节点到地的电流值。如果记BΛB^T＝A，V＝x，I＝b，则(5.1)式可以化成线形方程组的标准型Ax＝b。

对于这个标准的线性方程组，由于它的对称正定特性，我们可以进一步的进行优化，提高求解速度。首先由于方程组右边的电流源向量一般数量值很小，在计算的时候受数值误差的影响比较大，所以我们将原来存在于方程组左边的和电压源相连的电导项移动到方程组右边，以增加右边电流向量的绝对值，提高方程组抗数值误差干扰的特性。同时由于在方程组左边的对角线元素的系数中出现了相应的和电压源相连的电导项，大大改善了该矩阵的条件数，而线性方程组的求解速度和矩阵的条件数关系密切，同时该改善方法将随着电网络中电压源的增加而变得更好。下面两个公式中的(5.3)为实际生成系数矩阵的公式。其中v_k是矩阵中一行所代表的一个节点的电压值，g_ik是该节点和电路中其他节点相连的电导值，v_i是其他节点的电压值，g_jk是该节点和电压源相连接的电导值，v_j是电压源的电压值。P是和电压源相连的节点下标集合，N’是为原网络中网状结构上的节点总数。

$\underset{j\&Element;P}{\mathrm{\&Sigma;}}{g}_{\mathrm{jk}}(v_j-v_k)+\underset{i=1,i\&NotEqual;k}{\overset{N^{\&prime;}}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ik}}(v_i-v_k)=i_k---\left(5.2\right)$

$-v_k(\underset{j\&Element;P}{\mathrm{\&Sigma;}}{g}_{\mathrm{jk}}+\underset{i=1,i\&NotEqual;k}{\overset{n}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ik}})+\underset{i=1,i\&NotEqual;k}{\overset{N^{\&prime;}}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ik}}{v}_i=i_k-\underset{j\&Element;P}{\mathrm{\&Sigma;}}{g}_{\mathrm{jk}}{v}_j---\left(5.3\right)$

可以将(5.3)写成矩阵的形式：

G·V＝I

G为一个N’×N’的矩阵，其中 $\left\{\begin{array}{c}{g}_{\mathrm{ij},i\&NotEqual;j}=g_{\mathrm{ij}}\\ g_{\mathrm{ii}}=-(\underset{j\&Element;P}{\mathrm{\&Sigma;}}{g}_{\mathrm{ij}}+\underset{j=1,j\&NotEqual;i}{\overset{N^{\&prime;}}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ij}})\end{array}\right.$

I为一个N’维列向量，其中 $I_k=i_k-\underset{j\&Element;P}{\mathrm{\&Sigma;}}{g}_{\mathrm{jk}}{v}_k$

其中G矩阵就是我们需要在计算初始时得到的计算矩阵

5.2计算矩阵的生成步骤

1)新的节点连接关系索引表

完成所有特殊元件的处理工作之后，电源网中只剩下电导、电流源和电压源构成的网结构，由于已经消除了许多节点，在这一步，仿真程序重新对网络中所有的节点进行编号，并形成新的节点连接关系索引表，此表中仅含有电导、电流源和电压源，同时仿真程序将根据这个表和算法规则，建立描述电源网络电特性的特殊矩阵和相应的电流源向量。

             新建的节点连接关系索引表

节点编号	连接电导	连接电流源	连接电压源
				1	G1，G3	I1	无
2	G4	无	V1

2)对于化简后的电源网络而言，程序将生成一个N’×N’的对称正定的计算矩阵，N’为原网络中网状结构上的节点总数。生成该矩阵的方法如下：

a)将接地节点标记为0节点，将接地节点从节点连接关系索引表中除去

b)按照节点连接关系索引表中节点排列的顺序，以该顺序作为生成矩阵的行顺序。根据(5.3)，生成矩阵的一行时，按节点索引表中连接电导一栏中列出的电导，通过查找电导参数表，找到和该节点有连接关系的其他节点在新索引表中的编号，在矩阵该行的对应列上填入相关电导的大小。

c)如果一个节点不和电压源相连，该行主元的大小等于和该行主元节点相连的电导大小之和的相反数，如果该节点和电压源相连，该行主元的大小等于和该行主元节点相连的电导大小之和加上连接在电压源和该节点间电导的和的相反数。

3)在生成矩阵的同时，仿真程序也同时生成对应仿真时刻的电流源向量。如果一个节点不和电压源相连，则该矩阵行所对应的电流源等于节点连接关系索引表中电流源一栏中所有电流源的代数和(流入为正，流出为负)；如果一个节点和电压源相连，那么该矩阵行所对应的电流源除了要取节点连接关系索引表中电流源一栏中所有电流源的代数和之外，还需要加上所有与之相连的电压源和连接电压源电导的乘积之和。

6.求解仿真矩阵

上面得到的G矩阵在仿真过程中不随时间变化而变化，而电流源向量I则需要根据每个仿真时刻进行更新。仿真程序得到该矩阵和向量之后，就通过调用一个求解引擎，计算出该时刻各个节点上电压的分布情况。本程序使用的求解引擎是基于不完全乔莱斯基分解的预优共轭梯度算法的，该引擎具有速度快，稳定性高的特点，是本仿真程序的一个核心部件。

6.1基于不完全乔莱斯基分解的预优共轭梯度求解算法

对于(5.1)式给出的代数方程组，如果其系数矩阵A是一个对称正定的矩阵，就可以借助于泛函的变分原理将代数方程组求解问题转化为一个求解泛函极值的问题。

设A是阶对称正定矩阵，b是已知向量，x^*是方程组解的充分必要条件为满足

(x^*)＝min((x))x∈Rⁿ

其中 其中(，)符号为内积符号

该原理称为Ritz原理，是共轭梯度法及其演化算法的基础。

由变分原理可知，计算(x)极小的算法也就给出了求解原方程组的算法，而求得(x)极小的关键当然是如何使(x)的值下降。

首先考虑经典的最速下降法，设x₀是Rⁿ中任意一点，则(x)在该点下降最快的方向是该点的负梯度方向

经计算可知，

其中r₀称为残差向量，令x₁＝x₀+αr₀(6.1)

其中α待定，由于(x)沿r₀方向是下降的，于是有(x₁)＜(x₀)

为了使下降的效果最佳，自然要使

这样求得的α可以使x₁沿r₀方向下降最多，而(6.2)式的解为 ${\mathrm{\&alpha;}}_1=\frac{(r_0,r_0)}{({\mathrm{Ar}}_0,r_0)}$

如果在(6.1)式中取α＝α₁，就得到一步最速下降法，一般地有

x_k＝x_k-1+α_kr_k-1k＝1，2，…

其中r_k-1＝b-Ax_k-1， ${\mathrm{\&alpha;}}_k=\frac{(r_{k-1},r_{k-1})}{({\mathrm{Ar}}_{k-1},r_{k-1})}$

由于A是对称正定的，所以它的所有特征值都是正数，设为λ₁≥λ₂≥…≥λ_n＞0定义 ${\left|\right|x\left|\right|}_A={(x,\mathrm{Ax})}^{\frac{1}{2}}$ 为一种向量范数，则最速下降法的收敛性有如下的结果

${\left|\right|x_k-x^{*}\left|\right|}_A\&le;{\left(\frac{{\mathrm{\&lambda;}}_1-{\mathrm{\&lambda;}}_n}{{\mathrm{\&lambda;}}_1+{\mathrm{\&lambda;}}_n}\right)}^k{\left|\right|x_0-x^{*}\left|\right|}_A$

其中x^*是方程组的精确解，x₀为初始点

应该注意的是虽然r_k是(x)在x_k处的最速下降方向，但这是局部的，要尽快达到(x)的极小值，这样的选择不一定是最优的，另一方面，由于计算中的舍入误差会影响实际下降方向偏离最速下降方向，所以计算的时候具有一定的不稳定性。

共轭梯度的方法则主要针对最速下降法的这几个缺点进行了改进，其主要思想是在(x)的极小化过程中，下降应该具有整体最优的特点。

如果在Rⁿ中能够取到p₁，p₂…p_k是线性无关的，并且能够使x_k满足

                     (x_k)＝min((x))x∈S_k(6.3)

其中S_k＝span{p₁，p₂…p_k}，即由p₁，p₂…p_k张成的线性空间，则k＝n时就给出了原方程组的解。当k＝1时，只要取p₁＝r₀＝b-Ax₀就可以满足(6.3)式的要求

假设有p₁，p₂…p_k-1(k≥2)使得 $x_{k-1}=\underset{i=1}{\overset{k-1}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i{p}_i$ 且满足(x_k)＝min((x))x∈S_k共轭梯度法取x_k＝x_k-1+α_kp_k，只要能选择p_k，使得(p₁，Ap_k)＝0 i＝1，2…k-1(6.4)成立(称p_k和p₁，p₂…p_k-1关于A共轭)，就能满足(6.3)。而满足(6.4)式的非0向量可以被构造出来：

p_k＝r_k-1+β_kp_k-1 ${\mathrm{\&beta;}}_k=-\frac{(p_{k-1},{\mathrm{Ar}}_{k-1})}{(p_{k-1},{\mathrm{Ap}}_{k-1})}$

这样通过反复迭代，序列{x_k}将会收敛到原方程组的解。

从理论上讲，共轭梯度法是一种直接求解线性方程组的方法，由方法的构造过程可见，若A是n阶矩阵，则x_n应当是方程组的精确解，但是由于实际计算中舍入误差的影响，p₁，p₂…p_k不可能严格保持A共轭的特性，甚至当k大到一定程度后，p₁，p₂…p_k可能变得几乎线性相关，所以共轭梯度法更多的作为迭代法使用。

由于A是对称正定的，设它的特征值为λ₁≥λ₂≥…≥λ_n＞0，用共轭梯度法给出的向量序列将满足以下估计：

${\left|\right|x_k-x^{*}\left|\right|}_A\&le;{\left(\frac{\sqrt{{\mathrm{\&lambda;}}_1}-\sqrt{{\mathrm{\&lambda;}}_n}}{\sqrt{{\mathrm{\&lambda;}}_1}+\sqrt{{\mathrm{\&lambda;}}_n}}\right)}^k{\left|\right|x_0-x^{*}\left|\right|}_A$

共轭梯度法作为一种快速收敛的迭代法其优点是可以充分利用系数矩阵的系数性，计算步骤简明，易于并行计算，但是在系数矩阵的特征值出现非均匀分布的情况时，收敛的速度将收到影响。对于经典的共轭梯度法，其收敛的快慢依赖于系数矩阵的谱分布情况，当系数矩阵的条件数很小或特征值分布很集中时它可以用很少几步得到高精度的近似解。而如果系数矩阵的特征值非常不均匀地分布在一个很长的区间上时，共轭梯度法的收敛速度就会变得很慢。为了克服这一个缺点，使之在分析差异很大电路方程的系数矩阵时能够给出比较平稳的求解速度，必须要考虑对电路方程的系数矩阵进行预处理。如果A是对称正定的，那它就一定存在乔莱斯基分解，但如果A还是大型稀疏的，则分解很可能会完全破坏A的稀疏性，而不完全乔莱斯基分解就是将A分解成

A＝LL^T+R

其中L是一个上三角矩阵，R称为剩余矩阵，由于这里R是可以变化的，所以L的稀疏性就可以的到保证，从而克服乔莱斯基分解破坏矩阵稀疏性的问题。可以令预优矩阵M＝LU，U和L有关系U＝DLT，其中D是由A的对角元素组成的矩阵。当剩余矩阵R选得合适时，预优矩阵可以具有和系数矩阵A相同的稀疏性并且和A相近。理论上可以证明，M和A越相近，预优共轭梯度法收敛得越快。我们使用类似于高斯消元的方法来获得系数矩阵的不完全乔莱斯基分解和预优矩阵。不完全乔莱斯基分解的过程如附图18所示，设系数矩阵A已经被分解到第i步，则此时元素a_ii成为主元。采用“高斯消元法”使a_ii所在列上的其它元素消为零。然而，这个“高斯消元法”是不完全的，这是因为：既然我们的目的是求得一个上(下)三角矩阵U(L)，那么下(上)三角矩阵的元素可以不加入消元过程。考虑第k列的情况，只有当a_ik和a_jk(i＜j＜k)都不为零时，高斯消元才会发生，在其它情况下a_ik和a_jk都不变。系数矩阵A的稀疏性保证了不完全的乔莱斯基分解的运算量很小。

7.恢复被消除的变元

上一步计算仅仅得到化简后节点电压得分布情况，仍然需要经过两个步骤的计算才能得出所有的节点电压分布情况

7.1约简节点的重建计算

在求解出紧密网上所有节点的电压值之后，我们首先需要重建所有中间节点的电压然后才能够重建树节点的电压，因为肯定会有树的根在中间节点上的情况。根据上面等效时进行的推导，在得到所有RL节点两端节点在某个

时刻的电压后，根据下面的公式可以逐一恢复RL节点此时刻的电压值。

$V_{\mathrm{RLnode}}^t=m\&CenterDot;V_{\mathrm{Gnode}}^t+n\&CenterDot;V_{\mathrm{Lnode}}^t-o\&CenterDot;{\mathrm{IS}}_L^t---\left(7.1\right)$

由于计算中需要的所有系数在开始合并中间的结点时都进行了计算和存储，所以按中间节点链进行加减和乘除运算就能恢复所有中间节点的电压。

恢复树节点电压的过程是从树根到树的底部的层次遍历过程。对于每一个层次上的树枝，如果我们知道了其亲节点(TopNode节点)的电压就能够求出附图15中Middle Node和Bottom Node的电压，而Bottom Node又将是下一层某个树枝的根节点，如此重复进行计算就可以重建所有树节点的电压。通过进行化简过程的逆向递推，可以证明恢复化简节点电压的迭代计算式如下：

$V_{\mathrm{middle}}^t=c\&CenterDot;V_{\mathrm{top}}^t-{\mathrm{eI}}_{\mathrm{middle}}^t---\left(7.2\right)$

$V_{\mathrm{bottom}}^t=a\&CenterDot;V_{\mathrm{middle}}^t+d\&CenterDot;(I_L^t-I_{\mathrm{bottom}}^t)---\left(7.3\right)$

其中V_top ^t+h是t+h时刻TopNode节点的节点电压，V_middle ^t+h和V_bottom ^t+h是要计算的t+h时刻的被消除掉的中间节点的电压，I_middle ^t是t时刻算出的MiddleNode节点的等效对地电流，I_bottom ^t是t时刻算出的BottonNode节点的等效对地电流，I^t _L是t时刻电感的等效电流源。a，b，c，d，e几个系数的定义和(3.2)式相同。

a＝G_L·d，b＝1-a，c＝G_top·e， $d=\frac{1}{G_L+G_{\mathrm{bottom}}},e=\frac{1}{G_{\mathrm{top}}+G_{\mathrm{middle}}}---\left(3.2\right)$

7.2化简节点的恢复步骤

1)根据RL结构索引表给出的顺序和RL结构化简参数表，按(7.1)还原所有中间节点的电压值

2)根据森林节点层次表中节点的顺序和树结构化简参数表，按(7.2)和(7.3)还原所有树节点的电压值

8.仿真的迭代

得到一个时刻中所有节点的电压分布之后，为了继续进行下一个时刻的仿真，程序需要根据各个节点电压分布，重新计算电流源向量。具体分为以下几步：

1)重新计算第二步中元件离散数据表中等效电流源一栏的数据

2)根据更新后的元件离散数据表和森林节点层次表，再次按自底向上的顺序得到合并之后树结构的等效电流源(此时无需再对电导结构进行重复计算)

3)重新消除掉RL结构。

4)按照生成矩阵时的方法，重新计算电流源向量，此时维持原矩阵不变。

不断重复以上操作，仿真程序就不断分析出新时刻的仿真结果，直到到达规定的仿真时刻为止。

Claims (1)

1.VLSI中树网混合供电结构的瞬态仿真方法，其特征在于它是以SunV880工作站为计算平台使用C语言在unix环境下实现的，它依赖于以下几个步骤：

1)程序扫描以Spice网表形式给出的电源网络的拓扑连接信息和电参数信息；

2)计算机程序利用步骤1)得到的信息，建立起分别描述电导、电容、电感和电流源、电压源参数的数据表；

3)程序通过扫描网表文件找到整个电源网络中的所有树形结构，并将其表示为一个层次式的节点列表，同时生成一个描述RL结构的列表；

4)通过下述任何一种时域下的差分离散方法计算电源/地线网络中的电容元件以及电感元件的等效电导大小和等效瞬时电流源大小，即把电容元件和电感元件近似地转换为并联的电流源和电阻元件，同时设初始时刻的等效电流源是零；

4.1)梯形欧拉法

对电容而言 $I_c^{t+h}=G_c\&CenterDot;V_c^{t+h}-(G_c\&CenterDot;V_c^t+I_c^t)$ $G_c=\frac{2C}{h},$

对电感而言 $I_L^{t+h}=G_L\&CenterDot;V_L^{t+h}+(G_L\&CenterDot;V_L^t+I_L^t)$ $G_L=\frac{h}{2L},$

其中h是用户输入的仿真步长，单位是秒，

I_c ^t+h、V_c ^t+h分别是t+h时刻电容支路上的电流和端电压，

I_c ^t、V_c ^t分别是t时刻电容支路上的电流和端电压，

G_c是电容元件离散化后的等效电导，

I_L ^t+h、V_L ^t+h分别是t+h时刻电感支路上的电流和端电压，

I_L ^t、V_L ^t分别是t时刻电感支路上的电流和端电压，

G_L是电感元件离散化后的等效电导，

若令 ${\mathrm{IS}}_c^t=-(G_c\&CenterDot;V_c^t+I_c^t)$ 则 $I_c^{t+h}=G_c\&CenterDot;V_c^{t+h}+{\mathrm{IS}}_c^t,$

${\mathrm{IS}}_L^t=(G_L\&CenterDot;V_L^t+I_L^t)$ 则 $I_L^{t+h}=G_L\&CenterDot;V_L^{t+h}+{\mathrm{IS}}_L^t,$

4.2)线性多步法

$I_c^{t+h}=G_c\&CenterDot;V_c^{t+h}-(G_c\&CenterDot;V_c^t+a_1{I}_c^t+a_2{I}_c^{t-h}+a_3{I}_c^{t-2h})$

$G_c=\frac{24C}{9h},$ $a_1=\frac{19}{9},$ $a_2=-\frac{5}{9},$ $a_3=\frac{1}{9}$

$I_L^{t+h}=G_L\&CenterDot;V_L^{t+h}+(b_1{V}_L^t+b_2{V}_L^{t-h}+b_3{V}_L^{t-2h}+I_L^t)$

$G_L=\frac{9h}{24L},$ $b_1=\frac{19h}{24L},$ $b_2=-\frac{5h}{24L},$ $b_3=\frac{h}{24L}$

若令 ${\mathrm{IS}}_c^t=-(G_c\&CenterDot;V_c^t+a_1{I}_c^t+a_2{I}_c^{t-h}+a_3{I}_c^{t-2h})$ 则 $I_c^{t+h}=G_c\&CenterDot;V_c^{t+h}+{\mathrm{IS}}_c^t,$

若令 ${\mathrm{IS}}_L^t=(b_1{V}_L^t+b_2{V}_L^{t-h}+b_3{V}_L^{t-2h}+I_L^t)$ 则 $I_L^{t+h}=G_L\&CenterDot;V_L^{t+h}+{\mathrm{IS}}_L^t;$

5)按照步骤3)中得到的森林节点层次表中自底向上的顺序，通过查找节点连接关系索引表和电容、电感元件的离散数据表，把连接在各个树节点上，需要改变的电容和电感替换成相应的电导和电流源；

6)按照上述森林节点层次表中节点的顺序，从每棵树的底端开始，反复使用线性含源两端电阻网络的等效电路原理，把替换后的电流源和电导元件作层次式的合并，使供电网上所有的树都归并为根节点上的电导和电流源，它含有以下顺序进行的两个步骤：

6.1)把所有的电容和电感的等效电导合并到根节点上，

6.2)在每个时刻更新完所有电感和电容的等效电流源之后，再把更新后树节点上电容和电感的等效电流源再进行合并，最终得到根结点上的等效电流源；

7)程序按照RL结构索引表给出的顺序用上述等效原理进行变换，消除电源网中的中间节点，生成一个以RL结构索引表的顺序为排列顺序的RL结构化参数表，该参数表中含有的m，n，o为中间计算系数，且任意通过下列表达式得到的中间系数都应当被认为等同于系数m，n，o：

m＝o·G，n＝1-m， $o=\frac{1}{G+G_L}$

其中，G是RL结构上的自带电导，

      G_L是RL结构中L等效后得到的电导；

8)重新对网络中所有的节点进行排序，形成新的节点连接关系索引表；

9)根据简化后的电源/地线拓扑结构，生成计算矩阵和相应的电流源向量，形成供电节点电压、电流的线性方程组，它依赖于以下步骤：

9.1)把接地节点标记为零节点，把接地节点从节点连接关系索引表中除去；

9.2)基于基尔霍夫定律和欧姆定律得到供电节点电压、电流的线性方程组，如下式所示：

$\underset{j\&Element;P}{\mathrm{\&Sigma;}}{g}_{\mathrm{jk}}(v_j-v_k)+\underset{i=1,i\&NotEqual;k}{\overset{N^{\&prime;}}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ik}}(v_i-v_k)=i_k,$

该线性方程组的矩阵形式为G·V＝I，

其中，G为一个N’×N’的矩阵，称为计算矩阵，且 $\left\{\begin{array}{c}{g}_{\mathrm{ij},i\&NotEqual;j}=g_{\mathrm{ij}}\\ g_{\mathrm{ij}}=-(\underset{j\&Element;P}{\mathrm{\&Sigma;}}{g}_{\mathrm{ij}}+\underset{j=1,j\&NotEqual;i}{\overset{N^{\&prime;}}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ij}})\end{array}\right.,$

V为一个N’维列向量，其中每个元素代表网结构上各节点的电压值，

I为一个N’维列向量，其中每个元素代表每个节点到地的电流值，

同时，N’是供电网中网结构的节点数，

      v_k是矩阵中第k行所代表的一个节点的电压值，

      v_i是v_k节点之外的其他节点的电压值，

    v_j是电压源的电压值，

    P是和电压源相连的节点的下标集合，

    g_ik是节点k和其他节点i相连的电导值，

    g_ik是节点k和电压源相连的电导值；

10)通过调用基于不完全乔莱斯基分解的预优共轭梯度算法求解引擎得到简化后节点供电电压的分布；

11)在得到所有RL结构的两端节点在某个时刻的电压之后，根据RL结构化简参数表，通过扫描中间节点链，并按下式进行四则运算，就可以逐一恢复所有中间节点此时刻的电压值：

$V_{\mathrm{RLnode}}^t=m\&CenterDot;V_{\mathrm{Gnode}}^t+n\&CenterDot;V_{\mathrm{Lnode}}^t-o\&CenterDot;{\mathrm{IS}}_L^t;$

12)根据森林节点层次表中节点的顺序和树结构化简参数表，按下式还原所有该时刻树节点电压值：

$V_{\mathrm{middle}}^t=c\&CenterDot;V_{\mathrm{top}}^t-{\mathrm{eI}}_{\mathrm{middle}}^t$

$V_{\mathrm{bottom}}^t=a\&CenterDot;V_{\mathrm{middle}}^t+d\&CenterDot;(I_L^t-I_{\mathrm{bottom}}^t),$

其中a＝G_L·d，b＝1-a，c＝G_top·e， $d=\frac{1}{G_L+G_{\mathrm{bottom}}},$ $e=\frac{1}{G_{\mathrm{top}}+G_{\mathrm{middle}}},$

G、G_L的定义同7)所述；

13)根据以上得到的全部节点电压分布，重新计算电流源向量；

14)如果没有达到预定的仿真时间，跳转到步骤10)继续执行。

Images