一种快速求解含有非线性电阻元件的电路的方法 Claudiu Tufan, Alexandru Gabriel Gheorghe 和 George Marian Vasilescu

《Mathematics》：A Very Fast Method for Solving Electrical Circuits with Nonlinear Resistive Elements Claudiu Tufan, Alexandru Gabriel Gheorghe and George Marian Vasilescu

　　摘要：我们提出了一种改进的H?n?il?方法（HM），适用于包含多个非线性电阻元件（NRCEs）的电路。该方法用由线性电阻器和非线性受控源组成的发电机替换NRCEs，这些发电机的值通过迭代进行更新。在每次迭代中，都会在频域中求解得到的线性电路的每个谐波。我们评估了计算工作量，并将计算时间和结果与使用传统方法处理并联连接的单相半波整流器电路所获得的结果进行了比较。

1. 引言
电力电子技术和储能解决方案的发展增加了电气网络中非线性元件的数量。这一趋势导致了更高的谐波污染[1,2,3]。因此，目前需要一种快速高效的方法来求解NRCEs，以便准确评估谐波失真效应[4,5,6,7,8,9,10,11,12,13,14,15]。
其中一种方法是在计算要求较高的情况下表现出效率的方法，即HM。该方法被应用于各种技术领域，包括解决非线性介质中的电磁场问题和涡流问题，以及处理包含非线性电阻元件（NRCEs）的电路[16,17,18,19,20,21,22,23,24,25,26]。
该方法最初在[16]中提出，用于求解直流NRCEs，随后被扩展到NRCEs的周期性稳态分析以及三相NRCE配置[16,23,24,25,26]。HM能够有效处理其他方法难以解决的NRCE类型，包括三相电路、谐波共振以及具有不同序列电抗的电路。它还能处理涉及Gibbs现象的受控切换和突然信号转换。即使这些现象同时发生，该方法仍然有效。此外，它能够在合理的时间内分析多个谐波和每个周期内的采样点[23,24,25,26]。
该方法的优点包括数学上保证的收敛性、减少的计算内存和时间需求、能够以所需的精度计算结果（无需强制消除或忽略偶数或三倍谐波），以及大量的信号采样点[23,24,25,26]。
我们之前的工作在具有挑战性的条件下解决了NRCEs，通常需要加速技术或修改算法[24,25,26]。这可能会让人误认为该方法计算复杂。本文旨在澄清这一误解。对于其他方法也能处理的简单电路（例如具有周期性误差控制的单相半波整流器电路（HB, TR）来说，HM的速度非常快，甚至比商业软件包中实现的计算方法还要快[4]。本文对此进行了说明。我们对[4]中提出的单相半波整流器电路进行了比较分析，包括单个整流器和多个并联连接的整流器的情况。我们注意到，尽管在之前的工作中简要提到了求解NRCEs的可能性，但没有提供示例。
对于HM来说，计算工作量非常低。连接到非线性元件（NLE）端子的线性部分只需计算一次。使用非线性特性的校正是在时域中进行的，非常快速。计算要求最高的步骤是傅里叶级数展开，然后将其逆变换回时域。尽管如此，与替代方法的计算要求相比，这些操作仍然非常高效。当傅里叶级数截断到少数几个谐波（少于100个）时，计算工作量很低，时间也很短。由于该方法基于迭代过程，因此无需求解复杂的非线性方程组[23,24,25,26]。最初只需求解一组线性方程来计算迭代前的传输特性。
在[26]中，分析了一种改进的HM算法，该算法通过在开始时只进行一次傅里叶分析来减少计算时间，而不是在每次迭代时都进行。这里没有使用这种改进算法。对于分析的电路，原始算法[16,23,24]被证明速度极快。

该方法有两种变体：电压控制或电流控制的非线性源校正。它包括用由非线性源组成的发电机（电压或电流）替换非线性元件，这些发电机的值通过迭代进行更新，同时还包括内部电阻。这样，非线性就被转移到了源上，从而得到一个包含受控非线性源和被动线性电路元件的电路[16,24,25,26]。这样的电路可以通过应用标准线性电路定理在频域中轻松求解每个谐波。内部电阻值的选择范围确保了计算算法的收敛性。通过在该范围内选择最优值，可以最大化收敛速度[23,24,25,26]。非线性源的值通过迭代进行更新，直到达到指定的目标精度。一旦确定了非线性源的值，就可以在频域中轻松计算电路中的所有电流和电压。
本研究证明，与HB或TR等传统技术相比，HM具有更快的计算速度。我们将HM的应用扩展到包含多个非线性元件的电路。通过各种整流器电路配置提供了数值验证。结果表明，所提出的方法显著减少了计算工作量。此外，HM在整个分析过程中保持了高精度。这确立了HM作为复杂非线性电路仿真的一个稳健且高效的替代方案。

2. HM用于NRCEs的分析
以下部分概述了[23,24,25]中应用于包含单个NLE的电路的HM的主要特性，并强调了将其扩展到NRCEs时的差异。分析重点介绍了使用电压校正的变体。简而言之，任何非线性源的算法为[23,24,25]：
(1)
其中 是非线性源的源电压，其值需要在时域中计算。 是使用傅里叶级数展开及其截断得到的源的复数谐波值向量——该操作用 表示。 是替换NLE的发电机端子的频域电压，而 是相应的时间域电压。 是 的逆操作，用于返回时域。 表示线性电路的谐波解以及基于电路中源的值（迭代过程中保持不变的初始电源值和在前面步骤中计算的非线性源）计算出的发电机端子处的电压。要解决的电路是线性的，可以使用叠加原理按每个谐波进行分析。操作符 在时域中使用元件的非线性特性进行校正， 是迭代次数。算法流程图如图1所示。
图1. HM算法流程图。
分析包含单个NLE的电路与包含多个非线性元件的电路之间的区别体现在操作符 上。替换NLE的发电机端子处的电压不仅取决于电路中已存在的源和发电机自身的非线性源，还取决于替换电路中其他非线性元件的其他发电机的非线性源[24,25,26]。从每个非线性源的电压到参与计算的非线性发电机的电压的传输特性仅在开始时计算一次。在每次迭代中，这些预先计算出的传输特性被应用于上一步获得的非线性源的频域值，并使用叠加定理对结果进行求和。
实际上，在每个NLE的端子上，等效发电机是在频域中计算的。此时，电路由线性元件和非线性源组成，可以逐个谐波轻松分析。处理多个非线性时的关键区别在于，每个元件的计算考虑了替换其余非线性元件的所有其他非线性发电机的贡献。
由于分析将每个频率分量分解为单独的线性问题，叠加定理使得解决方案直接且高效。无论电路的复杂性或组件数量（线性或非线性）如何，最终都在每个非线性发电机的端子上确定了等效发电机。此外，传输特性仅在预处理阶段计算一次，这是HM的主要优势之一。
通过将非线性特性转移到源中，只需要在频域中求解线性电路[23,24,26]。无论拓扑结构的复杂性如何，都可以使用现有的仿真工具（例如基于基尔霍夫定律的改进节点分析）非常快速地确定等效电压或电流发电机。涉及非线性特性的校正仅在时域中进行，计算效率高且所需的工作量最小。
图2展示了用电压发电机替换非线性元件以及计算 和 操作符的过程。
为了使算法收敛，(1)中的至少一个操作符 必须是收缩的，而其他操作符必须是非扩张的[23,24,25,26]。在我们之前的工作中详细讨论了选择非线性发电机的内部电阻 以确保 是收缩的过程，以及傅里叶级数展开和非展开的非扩张性质[23,24,25,26]。对于 是收缩的一个充分条件是——NLE在时域中的特征函数是利普希茨连续且均匀单调的[23,24,25,26]。这导致了(2)。
对于通常遇到的非线性元件，满足这些条件，因为它们在有限的电压和电流范围内工作。因此，电压和电流的比率有定义的最大值和最小值，这意味着它们的电阻保持在最大值 和最小值 之间[23,24,25,26]。
NLE的特征是 。替换NLE的发电机遵循构成方程。由此我们定义[23,24,25]：
(3)
替换NLE的发电机的电阻选择为在周期函数的希尔伯特空间中是收缩的。收缩因子 必须存在，以确保不等式(4) 对于所有 都成立。确保这一不等式成立的有效区间是[16,23,24,25,26]。
尚需确定对于包含多个非线性元素的电路，是否仍然是非扩张的。证明 是非扩张的，以及截断傅里叶级数甚至可以使其成为收缩的过程与[23,24,25,26]中介绍的过程类似。对于每个谐波，在每个发电机的端子（具有源组件 和内部电阻 ），我们可以计算包含源 和等效阻抗 的等效发电机。然后我们计算等效发电机在谐波 处的电压：
(5)
由此可知，满足 是非扩张的条件：
(6)
我们定义源的误差（距离）为非线性源在两个连续频域值之间的范数[23,24,26]：
(7)
我们定义迭代中的算法误差为：
(8)
其中 是电路中非线性元素的数量。
作为迭代的停止准则，我们使用相对误差，其中 表示电路供电电压的幅度。
在将频域中的非线性源值计算到所需精度后，可以轻松确定频域中所有电路元件的电流和电压[23,24,25,26]。
由于该算法基于皮卡-巴拿赫序列，总是可以保证收敛到固定点——即电路解。与牛顿-拉夫森（NR）方法不同，后者在某些条件下容易发散，HM提供了更好的数值鲁棒性。
总之，HM根据方程(1)生成一个收敛的皮卡-巴拿赫序列。我们从非线性发电机的时域值开始（初始值可以设为零）。应用傅里叶级数展开，得到频域中的源值。级数截断到所需的谐波数量。
在所有源值已知的情况下，结合初始电路源，可以求解每个谐波的电路，并获得所有发电机端子处的电压。这些值随后被转换到时域。对于每个非线性元件（NLE），根据其值，在时域中应用基于非线性特性的校正，从而得到一组新的非线性源值。然后从开始重复该算法，直到达到所需的相对误差。3. 说明性示例提出的HM算法的计算效率和准确性（图1和图2）与ADS仿真进行了基准测试。比较包括ADS提供的频域HB和时域TR求解器。我们使用带有电压校正的HM分析了[4]中提出的电路，考虑了图3所示的单半波整流器配置、图4所示的两个相同并联整流器的配置以及十个相同整流器并联的配置。图3. 单相半波整流器。图4. 分别有2个和10个单相半波整流器并联连接。电源是正弦波，幅度为325 V，频率为50 Hz，初始相位为0度。[4]中的非线性二极管特性使用由导通电阻和阻断电阻定义的两段式分段线性模型进行近似。非线性特性非常陡峭，比率非常大。图5中显示的非线性特性被应用于图3中的二极管和图4中的二极管。如[4]中所提到的，多整流器电路包括用于并联耦合的额外互连导体，这使它们与单整流器配置区分开来。因此，图3中的单整流器电路没有整合到图4所示的统一示意图中。图5. 非线性二极管特性。未来的工作将把模型扩展到包含指数部分的多段式分段线性特性。然而，预计计算时间不会显著增加。非线性特性校正是在时域中进行的，速度非常快——大约为10 μs。傅里叶级数被截断到9个、20个、40个和100个谐波，加上直流分量（分别表示为9+1、20+1、40+1和100+1）。HM的仿真是在GNU Octave 6.2.0 [27]中使用的，该软件工具基于图1和图2中呈现的算法和流程图开发的。我们使用了一台配备2.3 GHz 8核Intel Core i9处理器（第9代）和16 GB 2667 MHz DDR4内存的笔记本电脑。表1总结了计算时间（以毫秒为单位）、所需的迭代次数以及每个周期的采样点数。当相对误差达到时，假设收敛。分析使用了原始的HM算法[16,17,18,19,20,21,22,23,24,25]，而不是[26]中的修改版本。表1. 使用HM的各种配置的计算时间和迭代次数。没有使用[23,24,25,26]中提出的任何特定加速技术，因为它们是不必要的。计算使用了最佳内部电阻，这本质上提高了方法的收敛速度[23,24,25,26]。值得注意的是，可以在NLE端点计算等效电阻[24]，这可以进一步减少计算时间。如果需要更高的计算速度，可以减少谐波的数量、采样点的数量或相对误差。相反，如果需要更高的精度，可以增加这些参数。此外，可以应用加速技术[21,22,23,24,25]或修改后的算法[26]来提高性能，展示了该方法的多功能性。分析考虑了现有电路的对称性和非线性特性的相同性。使用非线性特性进行校正以及应用运算符和运算符仅每次迭代执行一次，从而显著节省了时间。如果每个整流器的非线性特性不同，则需要针对每个元件分别执行这些操作。在最坏的情况下，计算时间将与整流器的数量成线性比例增加。这种情况将在未来的研究中讨论。表2展示了使用ADS 2026 [28]和HB方法对相同电路和二极管特性进行的计算时间。为了保持一致性，傅里叶级数的截断和每个周期的采样点数保持不变。表2. 使用HB的ADS 2026 [28]的计算时间。表1中的计算时间低于[4]和表2中的时间。图6、图7和图8展示了使用100+1个谐波和256个采样点计算的10个整流器情况下的RC负载电压及其谐波谱（峰值）。除了直流分量外，谐波幅度都很小。然而，第二和第三谐波仍然显著。在这种情况下，忽略偶数或三倍谐波的方法可能会产生错误的结果。考虑到整流器的低通滤波器配置（图3），预期不会出现更高阶谐波。图6. 10个整流器的RC电压和偏差：使用100+1个谐波的HM与ADS HB结果。图7. 10个整流器的RC电压的谐波谱（峰值和直流分量的百分比）和偏差：使用100+1个谐波的HM与ADS HB结果。图8. 1个、2个和10个整流器的RC电压波形的叠加细节：使用100+1个谐波的HM与ADS HB结果。图8显示了1个、2个和10个整流器电路的RC电压叠加细节。随着并联整流器数量的增加，观察到了预期的电压降。HM的结果得到了HB结果的验证。两种方法之间的RC电压相对偏差（图6），定义为，不超过1.8 × 10?5。TR可以在ADS中检测到达到稳态的时间，对于本文研究的电路，这一时间发生在t = 140 ms。为了获得更接近HM和HB得到的解（波形），需要将最终仿真时间增加到Tstop = 260 ms。在这种设置下，HB和TR给出的解在周期内的电压偏差（对于TR）不超过5.5 mV。在TR中，时间步长可以由模拟器自动设置或用户手动设置。为了将计算时间与使用HM和HB获得的时间进行比较，我们在与前述方法相同的时间段内施加了相同数量的样本。结果如表3所示。表3. 使用TR的ADS 2026 [28]的计算时间。如果我们选择稳态检测并让程序根据积分方法（2到6阶的梯形或齿轮）选择时间步长，则达到稳态解的时间会更晚（190 ms到210 ms）。表1中的计算时间低于[4]和表2以及表3中的时间。在当前的GNU Octave [27]实现中，HM的传递特性是解析推导出来的，而程序仅执行数值计算。正在进行的发展旨在集成一个SPICE模块，以自动分析电路并计算任何用户定义的拓扑的传递特性。对于和运算符，使用了经典的傅里叶级数展开及其逆变换。系数在预处理阶段计算一次。周期被划分为相等的间隔，假设信号在点之间是线性的，从而允许解析计算积分。为了减少浮点误差，正弦和余弦的参数被映射到第一象限[23,24,25,26]。所有预处理操作都包含在总计算时间中。GNU Octave中可用的标准FFT和IFFT例程，如ADS中使用的，并未使用。虽然该方法支持FFT/IFFT，但在某些情况下，例如由共振或受控开关元件（例如晶闸管）引起的突然信号变化[21,22,23,24,25,26]，这些例程可能不可靠。在[24]中，使用FFT/IFFT显著增加了达到目标误差所需的迭代次数，尽管减少了每次迭代的时间内。尽管所选择的解析方法较慢，但它提供了更高的精度、多功能性和鲁棒性，适用于任何类型的电路。对于具有少量显著谐波的简单电路，FFT/IFFT可能会进一步减少计算时间。值得注意的是，ADS——无论是在频域还是时域中操作——都使用NR来处理非线性。虽然NR的特点是二次收敛，并以其速度著称，但它也容易出现收敛问题。尽管在ADS的情况下迭代次数减少了，但总体计算时间仍然很长。该方法具有多功能性，可以根据应用在计算速度和高精度之间进行权衡。通过减少分析的谐波数量或采样点数量，或通过提高收敛阈值，可以进一步最小化计算时间。相反，如果需要更高的精度，可以增加这些参数。此外，可以应用加速技术[21,22,23,24,25]或修改后的算法[26]来提高性能，展示了该方法的多功能性。分析考虑了现有电路的对称性和非线性特性的相同性。使用非线性特性进行校正以及应用运算符和运算符仅每次迭代执行一次，从而显著节省了时间。如果每个整流器的非线性特性不同，则需要分别对每个元件执行这些操作。在最坏的情况下，计算时间将与整流器的数量成线性比例增加。这种情况将在未来的研究中讨论。表2展示了使用ADS 2026 [28]和HB方法对相同电路和二极管特性进行的计算时间。为了保持一致性，傅里叶级数的截断和每个周期的采样点数保持不变。表2. 使用HB的ADS 2026 [28]的计算时间。表1中的计算时间低于[4]和表2中的时间。图6、图7和图8展示了使用100+1个谐波和256个采样点计算的10个整流器情况下的RC负载电压及其谐波谱（峰值）。除了直流分量外，谐波幅度都很小。然而，第二和第三谐波仍然显著。忽略偶数或三倍谐波的方法在这种情况下可能会产生错误的结果。鉴于整流器的低通滤波器配置（图3），预期不会出现更高阶谐波。图6. 10个整流器的RC电压和偏差：使用100+1个谐波的HM与ADS HB结果。图7. 10个整流器的RC电压的谐波谱（峰值和直流分量的百分比）和偏差：使用100+1个谐波的HM与ADS HB结果。图8. 1个、2个和10个整流器的RC电压波形的叠加细节：使用100+1个谐波的HM与ADS HB结果。图8显示了1个、2个和10个整流器电路的RC电压叠加细节。随着并联整流器数量的增加，观察到了预期的电压降。HM的结果得到了HB结果的验证。两种方法之间的RC电压相对偏差（图6），定义为，不超过1.8 × 10?5。TR可以检测到达到稳态的时间，对于本文研究的电路，这一时间发生在t = 140 ms。为了获得更接近HM和HB得到的解（波形），需要将最终仿真时间增加到Tstop = 260 ms。在这些设置下，HB和TR给出的解在周期内的电压偏差（对于TR）不超过5.5 mV。在TR中，时间步长可以由模拟器自动设置或用户手动设置。为了将计算时间与使用HM和HB获得的时间进行比较，我们在之前提到的方法中施加了相同数量的样本。结果如表3所示。表3. 使用TR的ADS 2026 [28]的计算时间。如果我们选择稳态检测并让程序根据积分方法（2到6阶的梯形或齿轮）选择时间步长，则达到稳态解的时间会更晚（190 ms到210 ms）。表1中的计算时间低于[4]和表2以及表3中的时间。在当前的GNU Octave [27]实现中，HM的传递特性是解析推导出来的，而程序仅执行数值计算。正在进行的发展旨在集成一个SPICE模块，以自动分析电路并计算任何用户定义的拓扑的传递特性。对于和运算符，使用了经典的傅里叶级数展开及其逆变换。系数在预处理阶段计算一次。周期被划分为相等的间隔，假设信号在点之间是线性的，从而允许解析计算积分。为了减少浮点误差，正弦和余弦的参数被映射到第一象限[23,24,25,26]。所有预处理操作都包含在总计算时间中。GNU Octave中可用的标准FFT和IFFT例程，如ADS中使用的，并未使用。虽然该方法支持FFT/IFFT，但在某些情况下，例如由共振或受控开关元件（例如晶闸管）引起的突然信号变化[21,22,25,26]，这些例程可能不可靠。在[24]中，使用FFT/IFFT显著增加了达到目标误差所需的迭代次数，尽管减少了每次迭代的时间。尽管所选择的解析方法较慢，但它提供了更高的精度、多功能性和鲁棒性，适用于任何类型的电路。对于具有少量显著谐波的简单电路，FFT/IFFT可能会进一步减少计算时间。值得注意的是，ADS——无论是在频域还是时域中操作——都使用NR来处理非线性。虽然NR的特点是二次收敛，并以其速度著称，但它也容易出现收敛问题。尽管在ADS的情况下迭代次数减少了，但总体计算时间仍然很长。该方法具有多功能性，可以根据应用在计算速度和高精度之间进行权衡。通过减少分析的谐波数量或采样点数量，或通过提高收敛阈值，可以进一步最小化计算时间。相反，通过增加这些参数，可以获得更高的精度。对于分析的电路，显著谐波的数量非常少（最多到第4个谐波）。因此，即使使用每个周期较少的采样点，也可以实现良好的精度。即使不利用电路对称性，在最坏的情况下，计算时间也会与整流器的数量成线性比例增加，同时计算成本仍然较低。HM所需的计算工作量很小。通过连接到NLE端点的电路的线性部分确定的传递特性仅在开始时计算一次。使用非线性特性进行的校正是在时域中进行的，速度非常快。最耗时的操作是傅里叶级数展开和随后的逆变换到时域。当截断到较少数量的谐波（少于100个）和较少的采样点（少于1000个）时，计算工作量很小，计算时间非常短，即使在标准计算硬件上也是如此。由于该方法基于迭代过程，因此无需解决复杂的非线性方程组[23,24,25,26]。只需最初解决一个线性方程组来计算传递特性，然后在迭代过程中重复使用这些特性。4. 结论本文的结果表明，HM非常快。它的性能超过了商业软件中实现的传统计算方法。这种性能是在不使用加速技术或修改算法的情况下实现的。无论是否存在高度非线性（陡峭）特性，该方法都保持高效。它适用于解决NRCEs，并能够快速准确地评估谐波失真。该方法的优点包括保证的、数学上证明的收敛性，以及低计算工作和时间。该方法还可以通过使用大量谐波和采样点来实现高精度，而无需强制消除偶数或三倍谐波。值得注意的是，HM已经成功解决了其他方法未能解决的特定NRCEs。这包括三相非线性配置、谐波共振以及具有不同序列电抗的电路。此外，它还处理了受控开关元件，如晶闸管，以及容易发生Gibbs现象的突然信号变化。该方法同时处理这些现象的组合，并允许在合理的计算时间内分析大量谐波[23,24,25,26]。到目前为止，尚未观察到HM的任何显著缺陷。对于性能可能下降的特定情况，我们在之前的论文[24,25,26]中开发的加速技术已成功实施，以解决这些限制。然而，这对于本文分析的电路来说并不是必需的。未来的研究将专注于全波整流器和使用更复杂的分段线性近似的非线性元件，并补充指数部分。在未来的工作中，我们打算对HM与ADS和其他商业软件包进行比较分析，以评估更复杂的拓扑和具有更多线性和非线性元件的电路，包括使用实验室测量数据的实验验证。应当注意的是，我们之前的研究[24,25,26]中获得的结果已经通过LTspice [29]中的时域仿真得到了验证。此外，我们打算将HM与其他文献中提出的计算和数值方法进行比较，并欢迎在这个领域的研究合作。我们的目标是证明HM具备集成到商业电路仿真软件中所需的特性和优势。这可能包括为现有的开源、免费软件或专有仿真平台开发一个专门的附加模块。