CN107305535B - 一种加速电路网络状态方程迭代求解的方法 - Google Patents
一种加速电路网络状态方程迭代求解的方法 Download PDFInfo
- Publication number
- CN107305535B CN107305535B CN201610245472.2A CN201610245472A CN107305535B CN 107305535 B CN107305535 B CN 107305535B CN 201610245472 A CN201610245472 A CN 201610245472A CN 107305535 B CN107305535 B CN 107305535B
- Authority
- CN
- China
- Prior art keywords
- matrix
- equation
- formula
- state equation
- state
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 66
- 239000011159 matrix material Substances 0.000 claims abstract description 59
- 239000013598 vector Substances 0.000 claims abstract description 29
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 20
- 230000010354 integration Effects 0.000 claims abstract description 15
- 230000008569 process Effects 0.000 claims abstract description 15
- 238000002939 conjugate gradient method Methods 0.000 claims abstract description 8
- 238000004364 calculation method Methods 0.000 claims description 10
- 238000006467 substitution reaction Methods 0.000 claims description 5
- 238000009795 derivation Methods 0.000 claims description 3
- 230000001133 acceleration Effects 0.000 abstract description 3
- 230000001131 transforming effect Effects 0.000 abstract description 2
- 238000000354 decomposition reaction Methods 0.000 description 2
- 229910052731 fluorine Inorganic materials 0.000 description 2
- 125000001153 fluoro group Chemical group F* 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000008901 benefit Effects 0.000 description 1
- 239000003990 capacitor Substances 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000006870 function Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- General Physics & Mathematics (AREA)
- Pure & Applied Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Operations Research (AREA)
- Computing Systems (AREA)
- Complex Calculations (AREA)
- Networks Using Active Elements (AREA)
Abstract
本发明涉及一种加速电路网络状态方程迭代求解的方法，包括下述步骤：(1)列写电路的系统状态方程，并将其变换为一阶状态方程；(2)利用直接积分法对一阶状态方程进行迭代求解；(3)在直接积分法迭代过程中调用共轭梯度法或带预条件的共轭梯度算法实现加速求解。本发明提供的技术方案通过变换直接积分法迭代格式的求解步骤，使求解过程中的矩阵逆乘向量运算满足了共轭梯度迭代求解算法所要求的矩阵特性，使求解过程可以调用现已成熟的、高速的共轭梯度算法来完成求解加速。
Description
技术领域
本发明涉及一种电路网络状态方程的迭代求解方法，具体涉及一种加速电路网络状态方程迭代求解的方法。
背景技术
目前，在电力系统分析软件中一般采用LU分解等直接解法对电路网络方程进行求解。直接解法一般要求生成满阵，不能充分利用电路网络方程中系数矩阵的稀疏特性，会占用较大的内存空间。直接解法可以有效求解一千阶以下的电路方程，而对于拥有上万个节点的大规模电路网络则不具有任何优势。
共轭梯度法(以及带有预调件的共轭梯度法)是一种成熟的迭代求解稀疏矩阵方程的算法，但其要求待解矩阵方程中的系数矩阵具有正定性和对称性。使用状态方程法对电路网络进行分析时，其迭代求解式中的系数矩阵并不满足正定性和对称性，而无法直接使用共轭梯度算法。
发明内容
为解决上述现有技术中的不足，本发明提供一种加速电路网络状态方程迭代求解的方法，该方法通过变换矩阵方程的形式，使待解方程的系数矩阵满足正定性、稀疏性和对称性，从而可以调用共轭梯度等成熟的快速迭代算法，以达到求解加速的目的。
本发明的目的是采用下述技术方案实现的：
本发明提供一种加速电路网络状态方程迭代求解的方法，其改进之处在于，所述方法包括下述步骤：
(1)列写电路的系统状态方程，并通过变量代换，将其变换一阶形式的状态方程；
(2)利用直接积分法对一阶状态方程进行迭代求解；
(3)在直接积分法迭代过程中调用共轭梯度法或带预条件的共轭梯度算法实现加速求解。
进一步地，所述步骤(1)中，电路的系统状态方程表达式如下：
其中：
是n阶列向量，每一个元素值代表对应节点的电位；KC是电容系数矩阵；KR是电阻系数矩阵；KL是电感系数矩阵；I是电路中的输入向量，代表电路中的电源；上述三个系数矩阵均为n行n列的方阵，均满足正定性、对称性和稀疏性；
通过变量代换，将系统状态方程①变换为一阶状态方程，如下：
其中：
K1、K2均表示系数矩阵；x为状态变量；E是单位对角矩阵。
进一步地，所述步骤(2)中，对于标准形式的一阶动态方程②，用直接积分法进行迭代求解，直接积分法的迭代格式如下式④所示：
其中：Δt是第n个时刻与n+1时刻之间的时间步长；β是一个可供选择的系数，当β取值为0时，称为前差法(或前向欧拉法)，稳定性为条件稳定，数值精度为一阶精度；当β取值为0.5时，称为Crank-Nicolson法(或梯形法)，稳定性为无条件稳定，数值精度为二阶精度；当β取值为1时，称为后差法(或后向欧拉法)，稳定性为无条件稳定，数值精度为一阶精度；xn代表了第n步，即t＝nΔt时刻的系统状态变量；xn+1代表了第n+1步，即t＝(n+1)Δt时刻的系统状态变量；Rn和Rn+1分别代表了第n步，即t＝nΔt时刻和第n+1步，即t＝(n+1)Δt时刻的系统输入向量。
进一步地，所述步骤(3)包括下述步骤：
1)引入中间变量a和b，将迭代式④拆解为式⑤和⑥；
2)利用式⑤计算中间变量a和b；
3)式⑥拆分成式⑦和⑧，再将式⑧代入⑦中消去其中的变量Ψn+1，得到式⑨；
4)调用预调件的共轭梯度算法(共轭梯度法是一种常见的矩阵方程求解方法)对式⑨中的矩阵逆乘向量运算进行求解，得到状态变量
5)利用式⑧获得状态变量Ψn+1，最终得到状态变量
进一步地，所述步骤1)中，将式③代入式④中，并按照计算顺序将迭代式④在等号左右两边分解成两个部分，如式⑤和⑥，其中的a和b是推导过程引入的中间变量，均为n阶列向量；
其中：KC是电容系数矩阵；KR是电阻系数矩阵；KL是电感系数矩阵；E是单位对角矩阵；式⑤是利用前一时刻的状态量
和Ψn计算中间变量a和b，仅涉及矩阵-向量的乘法、向量的加法运算；式⑥是利用中间变量求解后一时刻的状态量
和Ψn+1，需要进行一次矩阵逆乘向量的运算；Δt是迭代求解时选择的时间步长；In和In+1分别是t＝nΔt和t＝(n+1)Δt时刻的系统输入向量。
进一步地，所述步骤3)包括：
其中：KC是电容系数矩阵；KR是电阻系数矩阵；KL是电感系数矩阵；E是单位对角矩阵；a和b均为中间变量；由于KC、KL和KR都满足正定性和对称性，且β和Δt均大于零，所以
为正定对称的，满足调用共轭梯度算法求解
的前提条件。
与最接近的现有技术相比，本发明提供的技术方案具有的优异效果是：
1.本发明通过变换直接积分法迭代格式的求解步骤，使求解过程中的矩阵逆乘向量运算满足了共轭梯度迭代求解算法所要求的矩阵特性，使求解过程可以调用现已成熟的、高速的共轭梯度算法来完成求解加速。
2.本发明在求解过程中保持了状态方程稀疏性的特点，可以使用成熟的稀疏存储和求解技术，可以在有限的内存空间内实现大规模节点数的电路方程仿真。
附图说明
图1是本发明提供的加速电路网络状态方程迭代求解的方法的流程图。
具体实施方式
下面结合附图对本发明的具体实施方式作进一步的详细说明。
以下描述和附图充分地示出本发明的具体实施方案，以使本领域的技术人员能够实践它们。其他实施方案可以包括结构的、逻辑的、电气的、过程的以及其他的改变。实施例仅代表可能的变化。除非明确要求，否则单独的组件和功能是可选的，并且操作的顺序可以变化。一些实施方案的部分和特征可以被包括在或替换其他实施方案的部分和特征。本发明的实施方案的范围包括权利要求书的整个范围，以及权利要求书的所有可获得的等同物。在本文中，本发明的这些实施方案可以被单独地或总地用术语“发明”来表示，这仅仅是为了方便，并且如果事实上公开了超过一个的发明，不是要自动地限制该应用的范围为任何单个发明或发明构思。
本发明提供一种加速电路网络状态方程迭代求解的方法，其流程图如图1所示，包括下述步骤：
(1)一般线性电路(由电阻、电容、电感以及电源等基本元件组成的，具有n+1个节点的电路，其中0号节点定义为参考地电位节点)列写电路的系统状态方程，并通过变量代换，将其变为一阶形式的状态方程；
电路的系统状态方程表达式如下：
其中：
是n阶列向量，每一个元素值代表对应节点的电位；KC是电容系数矩阵；KR是电阻系数矩阵；KL是电感系数矩阵；Is是电路中的输入向量，代表电路中的电源；上述三个系数矩阵均为n行n列的方阵，均满足正定性、对称性和稀疏性；
一阶状态方程表达式如下：
其中：
x为状态变量；E是单位对角矩阵。
(2)利用直接积分法对一阶状态方程进行迭代求解；对于标准形式的一阶动态方程②，用直接积分法进行迭代求解，直接积分法的迭代格式如下式④所示：
其中：Δt是第n个时刻与n+1时刻之间的时间步长；β是一个可供选择的系数，当β取不同的值时，相关的名称如下表所述：
表1β的选择与数值稳定性
xn和xn+1分别是t＝nΔt和t＝(n+1)Δt时刻的状态变量；Rn和Rn+1分别是t＝nΔt和t＝(n+1)Δt时刻的系统输入向量。
将式③代入式④中，并按照计算顺序将迭代式④在等号左右两边分解成两个部分，如式⑤和⑥，其中的a和b是推导过程引入的中间变量，均为n阶列向量；
其中：KC是电容系数矩阵；KR是电阻系数矩阵；KL是电感系数矩阵；E是单位对角矩阵；式⑤是利用前一时刻的状态量
和Ψn计算中间变量a和b，仅涉及矩阵-向量的乘法、向量的加法运算，使用计算机计算时并不占用太多计算时间。式⑥是利用中间变量求解后一时刻的状态量
和Ψn+1，需要进行一次矩阵逆乘向量的运算，这个运算是整个计算过程中最耗时的计算，一般采用LU分解等计算方法求解，需要占用大量的计算机内存和计算时间。Δt是第n个时刻与n+1时刻之间的时间步长、In和In+1分别是t＝nΔt和t＝(n+1)Δt时刻的系统输入向量。
如果能够加速式⑥的求解速度，也就加速了整个迭代的求解过程。本专利即给出一种加速式④迭代运算的方法。
(3)在直接积分法迭代过程中调用共轭梯度法或带预条件的共轭梯度算法实现加速求解，对于满足稀疏性的高阶正定对称系数矩阵方程，一般采用已经成熟的共轭梯度等算法来加速计算。迭代式④中的各系数矩阵满足稀疏性，却不满足正定性和对称性。为了加速迭代式④的计算，采用如下步骤：
1)引入中间变量a和b，将迭代式④拆解为式⑤和⑥；
2)利用式⑤计算中间变量a和b；
3)式⑥拆分成式⑦和⑧，再将式⑧代入⑦中消去其中的变量Ψn+1，得到式⑨；
其中：KC是电容系数矩阵；KR是电阻系数矩阵；KL是电感系数矩阵；E是单位对角矩阵；a和b均为中间变量；由于KC、KL和KR都满足正定性和对称性，且β和Δt均大于零，所以
为正定对称的，满足调用共轭梯度算法求解
的前提条件。
4)调用预调件的共轭梯度算法对式⑨中的矩阵逆乘向量运算进行求解，得到状态变量
5)利用式⑧获得状态变量Ψn+1，最终得到状态变量
本发明提供的加速电路网络状态方程迭代求解的方法，采用直接积分法迭代式分步求解的方法，即：通过方程变换获得满足正定对称特性的系数矩阵的方法；在直接积分法迭代过程中调用共轭梯度法(或带预条件的共轭梯度算法)等快速稀疏矩阵算法实现求解加速的方法。
以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员依然可以对本发明的具体实施方式进行修改或者等同替换，这些未脱离本发明精神和范围的任何修改或者等同替换，均在申请待批的本发明的权利要求保护范围之内。
Claims (1)
1.一种加速电路网络状态方程迭代求解的方法，其特征在于，所述方法包括下述步骤：
(1)列写电路的系统状态方程，并通过变量代换，将其变换一阶形式的状态方程；
(2)利用直接积分法对一阶状态方程进行迭代求解；
(3)在直接积分法迭代过程中调用共轭梯度法或带预条件的共轭梯度算法实现加速求解；
所述步骤(1)中，电路的系统状态方程表达式如下：
其中：
是n阶列向量，每一个元素值代表对应节点的电位；KC是电容系数矩阵；KR是电阻系数矩阵；KL是电感系数矩阵；I是电路中的输入向量，代表电路中的电源；上述三个系数矩阵均为n行n列的方阵，均满足正定性、对称性和稀疏性；
通过变量代换，将系统状态方程①变换为一阶状态方程，如下：
其中：
K1、K2均表示系数矩阵；x为状态变量；E是单位对角矩阵；
所述步骤(2)中，对于标准形式的一阶状态方程②，用直接积分法进行迭代求解，直接积分法的迭代格式如下式④所示：
其中：Δt是第n个时刻与n+1时刻之间的时间步长；β是一个可供选择的系数，当β取值为0时，称为前差法，稳定性为条件稳定，数值精度为一阶精度；当β取值为0.5时，称为Crank-Nicolson法，稳定性为无条件稳定，数值精度为二阶精度；当β取值为1时，称为后差法，稳定性为无条件稳定，数值精度为一阶精度；xn代表了第n步，即t＝nΔt时刻的系统状态变量；xn+1代表了第n+1步，即t＝(n+1)Δt时刻的系统状态变量；Rn和Rn+1分别代表了第n步，即t＝nΔt时刻和第n+1步，即t＝(n+1)Δt时刻的系统输入向量；
所述步骤(3)包括下述步骤：
1)引入中间变量a和b，将迭代式④拆解为式⑤和⑥；
2)利用式⑤计算中间变量a和b；
3)式⑥拆分成式⑦和⑧，再将式⑧代入⑦中消去其中的变量Ψn+1，得到式⑨；
4)调用预调件的共轭梯度算法对式⑨中的矩阵逆乘向量运算进行求解，得到状态变量
5)利用式⑧获得状态变量Ψn+1，最终得到状态变量
所述步骤(1)中，将式③代入式④中，并按照计算顺序将迭代式④在等号左右两边分解成两个部分，如式⑤和⑥所示，其中的a和b是推导过程引入的中间变量，均为n阶列向量；
其中：KC是电容系数矩阵；KR是电阻系数矩阵；KL是电感系数矩阵；E是单位对角矩阵；式⑤是利用前一时刻的状态量
和Ψn计算中间变量a和b，仅涉及矩阵-向量的乘法、向量的加法运算；式⑥是利用中间变量求解后一时刻的状态量
和Ψn+1，需要进行一次矩阵逆乘向量的运算；Δt是迭代求解时选择的时间步长；In和In+1分别是t＝nΔt和t＝(n+1)Δt时刻的系统输入向量；
所述步骤(3)包括：
其中：KC是电容系数矩阵；KR是电阻系数矩阵；KL是电感系数矩阵；E是单位对角矩阵；a和b均为中间变量；由于KC、KL和KR都满足正定性和对称性，且β和Δt均大于零，所以
为正定对称的，满足调用共轭梯度算法求解
的前提条件。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610245472.2A CN107305535B (zh) | 2016-04-19 | 2016-04-19 | 一种加速电路网络状态方程迭代求解的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610245472.2A CN107305535B (zh) | 2016-04-19 | 2016-04-19 | 一种加速电路网络状态方程迭代求解的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107305535A CN107305535A (zh) | 2017-10-31 |
| CN107305535B true CN107305535B (zh) | 2022-07-26 |
Family
ID=60152746
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610245472.2A Active CN107305535B (zh) | 2016-04-19 | 2016-04-19 | 一种加速电路网络状态方程迭代求解的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107305535B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113326671A (zh) * | 2021-06-15 | 2021-08-31 | 深圳华大九天科技有限公司 | 一种对电路中线性电阻计算电路方程节点电流的方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101763087A (zh) * | 2009-12-29 | 2010-06-30 | 浙江大学 | 一种基于非线性共轭梯度法的工业过程动态优化系统及方法 |
| CN102110079A (zh) * | 2011-03-07 | 2011-06-29 | 杭州电子科技大学 | 一种基于mpi的分布式共轭梯度法的调优计算方法 |
| CN103207941A (zh) * | 2013-04-27 | 2013-07-17 | 清华大学 | 集成电路供电网络全参数模型下瞬态分析方法及系统 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| GB2523342A (en) * | 2014-02-20 | 2015-08-26 | Ibm | Conjugate gradient solvers for linear systems |
-
2016
- 2016-04-19 CN CN201610245472.2A patent/CN107305535B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101763087A (zh) * | 2009-12-29 | 2010-06-30 | 浙江大学 | 一种基于非线性共轭梯度法的工业过程动态优化系统及方法 |
| CN102110079A (zh) * | 2011-03-07 | 2011-06-29 | 杭州电子科技大学 | 一种基于mpi的分布式共轭梯度法的调优计算方法 |
| CN103207941A (zh) * | 2013-04-27 | 2013-07-17 | 清华大学 | 集成电路供电网络全参数模型下瞬态分析方法及系统 |
Non-Patent Citations (1)
| Title |
|---|
| 基于预处理共轭梯度法的电力系统机电暂态仿真;林济铿等;《电工技术学报》;20080526(第05期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107305535A (zh) | 2017-10-31 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Gugercin et al. | Structure-preserving tangential interpolation for model reduction of port-Hamiltonian systems | |
| CN107944556B (zh) | 基于块项张量分解的深度神经网络压缩方法 | |
| Chen | The finite element method: its fundamentals and applications in engineering | |
| Zhao et al. | Power grid analysis with hierarchical support graphs | |
| Liu et al. | Error and its upper bound estimation between the solutions of GM (1, 1) grey forecasting models | |
| CN107305535B (zh) | 一种加速电路网络状态方程迭代求解的方法 | |
| Holod et al. | Enhanced preconditioner for JOREK MHD solver | |
| Balochian et al. | Stabilization of multi-input hybrid fractional-order systems with state delay | |
| Liu et al. | Legendre neural network for solving linear variable coefficients delay differential-algebraic equations with weak discontinuities | |
| WO2006132639A1 (en) | Circuit splitting in analysis of circuits at transistor level | |
| CN109417295B (zh) | 获取直流电力网潮流的无损耗全局线性偏心方法 | |
| Lantsov | A new algorithm for solving of harmonic balance equations by using the model order reduction method | |
| CN106934123A (zh) | 一种基于递归卷积的电路瞬态响应计算方法 | |
| Hossain et al. | Reduce order modelling of power system models using interpolatory projections technique | |
| Arun Govind Neelan et al. | Hyperbolic Runge–Kutta method using evolutionary algorithm | |
| CN109257945B (zh) | 获取直流电力网潮流的均衡电导补偿型全局线性偏心方法 | |
| Yan et al. | DeMOR: Decentralized model order reduction of linear networks with massive ports | |
| Pan et al. | Hierarchical resampling algorithm and architecture for distributed particle filters | |
| Fangzong et al. | Fast algorithm for structural dynamics problems using differential quadrature method and V-transformation | |
| TWI632544B (zh) | 四維度聲波傳播分析系統 | |
| Polcz et al. | Passivity analysis of rational LPV systems using Finsler’s lemma | |
| Wang et al. | Research on Two-Dimensional Electrostatic Particle Simulation Method Based on superLU | |
| Uddin | Revisionist integral deferred correction methods with application to the moving method of lines | |
| CN112949240B (zh) | 一种集中参数模型多物理场耦合仿真方法 | |
| Gourary et al. | Model order reduction by state vector selection (SVS) approach |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |