WO2012113096A1 - 时域变换法 - Google Patents

时域变换法 Download PDF

Info

Publication number
WO2012113096A1
WO2012113096A1 PCT/CN2011/000273 CN2011000273W WO2012113096A1 WO 2012113096 A1 WO2012113096 A1 WO 2012113096A1 CN 2011000273 W CN2011000273 W CN 2011000273W WO 2012113096 A1 WO2012113096 A1 WO 2012113096A1
Authority
WO
WIPO (PCT)
Prior art keywords
signal
original
domain
network
coordinate system
Prior art date
Application number
PCT/CN2011/000273
Other languages
English (en)
French (fr)
Inventor
范圣韬
丁辉
施希
Original Assignee
Fan Shengtao
Ding Hui
Shi Xi
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Fan Shengtao, Ding Hui, Shi Xi filed Critical Fan Shengtao
Priority to PCT/CN2011/000273 priority Critical patent/WO2012113096A1/zh
Publication of WO2012113096A1 publication Critical patent/WO2012113096A1/zh

Links

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/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms

Definitions

  • the invention mainly relates to dynamic process simulation of power systems and circuit systems, in particular electromagnetic transient simulation of power systems, and can also be used for calculation and simulation of communication and control systems.
  • the electromagnetic transient simulation program was mainly used in the electromagnetic transient process calculation of power systems. Now it has a wider range of applications in the field of dynamic and static analysis of power systems. For example: For dynamic process simulation of large-scale AC/DC power systems including FACTS and HVDC devices, the simulator needs to have the same rapidity as the transient stability calculation program, and at the same time has the same accuracy as the electromagnetic transient calculation program. However, when the electromagnetic transient calculation program is applied to large-scale power system calculations, the calculation step size is too small (typically 50 microseconds), which results in a very slow calculation speed, which cannot meet the needs of actual power system calculation. In order to expand the simulation scale of the electromagnetic transient calculation program, researchers at home and abroad have proposed many methods, which can be mainly divided into four categories:
  • the equivalence method reduces the size of the calculation system and speeds up the calculation by performing equivalent simplification on the parts of the system that do not need to be carefully studied.
  • FDNE Frequency Dependent Network Equivalent
  • electromagnetic transient-electromechanical transient hybrid simulation For the part of the system that needs detailed modeling and simulation, the electromagnetic transient program is used for simulation, and for the remaining ⁇ points, the electromechanical transient program is used for calculation, and then the interface technology is used. Data exchange between the two programs is performed at some point in time.
  • the Chinese Patent Office published the invention patent application CN101382969A on March 11, 2009, "Grid Electromagnetic Transient and Electromechanical Transient Hybrid Simulation System and Its Simulation Method” and published on January 26, 2011.
  • the invention patent application with the publication number CN101957872A "a hybrid real-time simulation method for AC/DC power systems" adopts an electromagnetic transient-electromechanical transient hybrid simulation method.
  • the fundamental reason is that the electromagnetic transient program uses three-phase instantaneous values.
  • the representation mode, and the electromechanical transient program is based on the representation of the fundamental frequency phasor.
  • Parallel processing technology is an important means to improve the simulation speed of electromagnetic transient programs. By dividing the grid in some way, the amount of calculation can be allocated to different calculation units, thereby increasing the calculation speed.
  • the most basic grid division method is based on the natural decoupling method caused by the transmission line traveling wave transmission. The delay caused by the traveling wave transmission on the long transmission line is enough to decouple the systems at both ends of the line.
  • the famous electromagnetic transient real-time simulation program RTDS uses this method; domestic and foreign researchers have also proposed some more flexible methods to make the grid division on any bus. Whether in offline non-real-time simulation or real-time simulation of electromagnetic transients, parallel processing technology has been widely studied and applied.
  • an electromagnetic transient offline non-real-time parallel simulation system and simulation method uses parallel processing technology to accelerate electromagnetic transients.
  • Offline non-real-time simulation; parallel technology was used to improve the invention patent application of CN101719182A published by the Intellectual Property Office of China on June 2, 2010, "A Parallel Electromagnetic Transient Digital Simulation Method for AC/DC Power System”
  • the speed of electromagnetic transient calculations used for real-time simulation can not meet the needs of practical scale power system calculation.
  • Power systems typically operate at a fundamental frequency (50 Hz or 60 Hz), and disturbances like electromechanical oscillations typically only cause lower frequency offsets in the voltage and current in the system, making the voltage and current in the power system a fundamental frequency.
  • a narrow-band signal centered; this is similar to a signal in a communication system that has a very high carrier frequency, while a modulated signal has a very low frequency (the information we really care about is often. These are lower-frequency signals) .
  • the electromagnetic transient program uses a three-phase instantaneous value representation. In addition to these modulated signals with a frequency of 3 ⁇ 4, it also requires a carrier signal with a very high frequency.
  • the traditional electromagnetic transient calculation program can only use a smaller calculation step.
  • some scholars transform the signal in the original system into a narrowband signal centered on DC (0 Hz) by Hilbert transform, and then in the transform domain.
  • the electromagnetic transient calculation because the frequency of the signal in the variation domain is very low, it can be calculated using a large calculation step without loss of precision (usually the calculation step can be expanded several hundred times), and finally The transformation can get the calculation result of the original system. For example, (Kai Strunz, Rachel Shintaku, and Feng Gao.
  • the aliasing effect may occur; the forward transform and inverse transform of the frequency domain transform It is more complicated and inefficient; and when the signal in the original system contains DC offset component, the frequency domain transform method based on Hilbert transform is difficult to transform the signal in the original system into a narrowband signal centered on DC. :
  • the present invention proposes a new transform, which is different from the existing frequency adaptive transient simulation method based on frequency domain transform, and the new transform is a time domain transform method.
  • the new transform has clear physical concepts, high computational efficiency, no confusion in the frequency domain transform, and the ability to process signals with DC offset components.
  • a set of methods suitable for electromagnetic transient simulation programs is established, which can significantly improve the calculation speed without losing the calculation accuracy. Since the electromagnetic transient simulation program for the power system is essentially based on circuit simulation, this method applicable to the power system electromagnetic transient simulation program can also be used to improve the circuit simulation program (eg, SPICE), and Can be used for simulation and calculations in other fields (such as communication, control).
  • the invention provides a time domain transform method for improving system dynamic process simulation, and the idea thereof is shown in FIG. 1 .
  • the differential or integral of the original signal is introduced as its orthogonal signal, and a new variable that changes slowly with time is constructed by time domain transformation; the calculation accuracy can be performed without loss of calculation accuracy.
  • the transform domain a larger calculation step is used to solve the numerical solution, and then the inverse signal can be used to obtain the solution of the original signal.
  • the discretization in the original domain can be derived by the calculation formula in the transformation domain during the specific implementation. Formula, and implemented directly in the original domain. Detailed below:
  • a signal in a power system is typically a sine or cosine function modulated by a slowly varying signal, having the following form:
  • the formula can be solved directly using some numerical method (for example, trapezoidal integral method). However, because it is limited by the fundamental frequency of the signal, it is used at this time.
  • the step size can only be small.
  • the step size used to numerically solve a differential equation depends on the speed at which the solution x(t) is solved over time, ⁇ (0 the faster the rate of change, the smaller the required calculation step size; therefore the invention
  • the idea is to transform the signal x(t) to be a more slowly changing signal by time domain transformation.
  • Dt dx where, is a new simulation parameter, its function is equivalent to the frequency adaptation to the offset frequency in the transient simulation FAST, usually can be set as the fundamental frequency of the system; solving (6) will get 0 and (0) at the same time Solution.
  • the transformation matrix ⁇ ( ⁇ ) is defined as follows:
  • M (t) and v (t) are signals that change more slowly than ⁇ ( ⁇ ), so the numerical solution of equation (15) can take a larger step size than equation (3).
  • a numerical solution such as the trapezoidal integral method
  • the equation (14) to solve the solution from the transform domain (UV coordinates).
  • the system is inversely transformed into the original domain ( ⁇ - ⁇ coordinate system), so that the solution of the original differential equation can be obtained.
  • Equation (25) gives a new discretization method in the original domain (X-Y coordinate system).
  • the signal in addition to the sine or cosine function component modulated by the slowly varying signal, often contains a DC bias component, such as the stator current at the exit short circuit in the steady state operation of the generator.
  • a DC bias component such as the stator current at the exit short circuit in the steady state operation of the generator.
  • transformation matrix T(t) is defined as follows:
  • equation (36) 0 a cos(i3 ⁇ 4 $1 ⁇ ( ⁇ 8 ⁇ )
  • equation (36) 0 a cos(i3 ⁇ 4 $1 ⁇ ( ⁇ 8 ⁇ )
  • the numerical solution of equation (38) can take a larger step size than equations (3) and (32).
  • a numerical solution such as the trapezoidal integral method
  • the (38) formula to solve the solution from the transform domain (PUV coordinates).
  • Figure 1 shows the schematic of the idea of the invention.
  • Figure 2 shows the waveform of the example signal x(t) as a function of time.
  • Figure 3 shows the waveforms of w(t) and v(t) corresponding to the sample signal x(t) over time.
  • Figure 4 shows the original network of an example.
  • Figure 5 shows the first-order differential extension network corresponding to the original network.
  • Figure 6 shows the time domain simulation equivalent circuit of the example system.
  • the component model of various network components is usually established on the branch level based on the formula (25) or (51).
  • the component model provides an algebra of the voltage, current and historical values at any time of the branch. equation.
  • n-order differential network and “n-order differential extension network” of the original network are first introduced.
  • the so-called “n-order differential network” of the original network refers to the network that gives the n-th order differential solution of the original network, and the "0-order differential network” is the original network itself.
  • the "n-order differential extension network” of the original network is obtained by combining the 0 ⁇ n-order differential networks of the original network.
  • the original network is shown in Figure 4.
  • the parameters of the components in the original network are shown in Table 1.
  • the "first-order differential extension network" of the original network is shown in Figure 5.
  • the left half of the figure is exactly the same as the original network (that is, "0-order differential network"), and the right half is the "first-order" of the original network.
  • Differential network The resistance of R d , the inductance of 1 ⁇ and the capacitance of C d are respectively equal to the resistance of R in the original network, the inductance of L and the capacitance of C.
  • the other parameters of "first-order differential network” are summarized in Table 2.
  • G 1 — ⁇ ⁇ K (54) 2L 2 2 ⁇ 's,
  • the matrix R in (54) ⁇ (57) is defined by the formula (16), and E is a second-order identity matrix.
  • (54) ⁇ (58) give the equivalent component model of inductance and capacitance based on equation (25).
  • the equivalent models of other network components can also be obtained similarly.
  • the corresponding equivalent component model can also be established (the "2nd order differential expansion network" needs to be established first).
  • the parameters of the equivalent component model based on (51) are as follows:
  • the matrix R in (59 ⁇ 62) is defined by the formula (40), E is a 3rd order identity matrix, and 1 ⁇ are respectively three inductances in the "2nd order differential extension network" (inductance L in the original network, 1st order The corresponding inductance in the differential network! ⁇ and the corresponding inductance L dd in the second-order differential network, t beau ⁇ , the vector of the voltage and current at the moment, and the three capacitors in the "2nd order differential extension network” Capacitance in the network,
  • the initial current of the corresponding inductor and the initial voltage of the capacitor in the "n-order differential network” are calculated as follows: Replace the inductor in the "n-order differential network” with a constant current source (current is “n-1 order differential network” "Initial current corresponding to the inductor", and the capacitor is replaced by a constant voltage source (the voltage is the initial voltage of the corresponding capacitor in the "n-1 order differential network”). Solving the replaced network will get the voltage corresponding to the constant current source.
  • the initial current of the corresponding inductor and the initial voltage of the capacitor in the "n-order differential network” can be obtained by the following equation (the inductance value of the inductor is L) , the capacitance value of the capacitor is C):
  • Hi ) U L IL (64)
  • uc' ⁇ to) i c IC (65)
  • the two calculation points can be first calculated by (7) or (33)
  • p(t), W(t) between the two calculated points can be obtained by simple linear interpolation.
  • v curtan arctan [v(i) I (67)
  • phase angle ⁇ ⁇ arctan [ ⁇ ⁇ IA n ] (69)
  • the current phase angle ⁇ ⁇ can be calculated using the phase angle ⁇ ⁇ — ⁇ of the previous time step: Since the phase angle change between the two time points is not too large, the continuity of the phase angle calculated by this method is very good and is not limited to the interval of [ ⁇ , ⁇ ].
  • the resistance of the switch is switched from a closed resistance to an on-state resistance or vice versa in a short period of time as follows.
  • the switch is at t.
  • the resistance of the switch changes according to the following formula:
  • r 0N and r 0FF are the closing resistance and the on-state resistance of the switch, respectively.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Description

说 明 书 时域变换法 技术领域
本发明主要涉及电力系统和电路系统的动态过程仿真， 特别是电力系统的电磁暂态仿真， 亦可用于通 讯及控制系统的计算和仿真。
背景技术
电磁暂态仿真程序早期主要用于电力系统的电磁暂态过程计算， 现在它在电力系统动静态分析等领域 有着更广泛的 用。 例如： 对于包含 FACTS和 HVDC等设备的大规模交直流电力系统进行动态过程仿真时， 需要仿真器具有暂态稳定计算程序一样的快速性， 同时又要具有电磁暂态计算程序一样的精确性。 但是当 电磁暂态计算程序应用于大规模电力系统计算时， 由于其采用的计算步长太小 （典型步长为 50微秒） 导 致其计算速度极慢， 无法满足实际电力系统计算的需要。 为了扩大电磁暂态计算程序的仿真规模， 国内外 的研究者提出了很多方法， 主要可以分为四大类：
1 )等值法
等值法通过对系统中不需要仔细研究的部分进行等值化简，从而缩小计算的系统规模，提高计算速度。 为了计及被简化部分在不同频率点的动态特性，一些学者采用了频率相关网络等值方法（FDNE : Frequency Dependent Network Equivalent ), 如 (文献 "Xi Lin, A. M. Gole, and Ming Yu. A wide-band multi-port system equivalent for real-time digital power system simulators. IEEE Trans. Power Syst. , vol. 24, no. 1, pp. 237-249, Feb. 2009" )。 然而当被简化的系统中包含 非线性元件时， 等值法在理论存 在着难以克服的困难。
2) 电磁暂态-机电暂态混合仿真
电磁暂态-机电暂态混合仿真的思路是： 对于系统中需要详细建模仿真的部分采用电磁暂态程序进行 仿真， 而对于剩余的 ^分则采用机电暂态程序计算， 然后通过接口技术在某些时间点进行两种程序间的数 据交换。如： 中国知识产权局分别于 2009年 3月 11日公开的公开号为 CN101382969A的发明专利申请《电 网电磁暂态与机电暂态混合仿真系统及其仿真方法》和 2011年 1月 26日公开的公开号为 CN101957872A 的发明专利申请 《一种交直流电力系统的混合实时仿真方法》就是采用了电磁暂态-机电暂态混合仿真的 方法。 虽然这种方法在实际应用中取得了一定的效果， 但是由于缺乏可靠的理论基础， 限制了其进一步发 展成为一个完美的解决方案， 根本的原因在于电磁暂态程序使用的是三相瞬时值的表示方式， 而机电暂态 程序则是基于基频相量表示的表示方法。
3 ) 并行技术
并行处理技术是提高电磁暂态程序仿真速度的一个重要手段。 通过采用某种方式对电网进行划分， 可 以将计算量分配到不同的计算单元上， 从而提高计算速度。 最基本的电网划分方法是基于传输线行波传输 导致的自然解耦方法， 由长传输线上的行波传输所导致的时延足以将线路两端的系统解耦开， 著名的电磁 暂态实时仿真程序 RTDS就使用了这种方法； 国内外学者也提出了一些更灵活的方法使得电网的划分可以 在任意的母线上进行。 无论是在电磁暂态的离线非实时仿真还是实时仿真中， 并行处理技术都获得了广泛 的研究和应用。 如： 中国知识产权局于 2009年 9月 16日公开的公开号为 CN101533428A的发明专利申请 《一种电磁暂态离线非实时并行仿真系统及仿真方法》采用了并行处理技术来加速电磁暂态的离线非实时 仿真； 在中国知识产权局于 2010年 6月 2日公开的公开号为 CN101719182A的发明专利申请《一种交直流 电力系统分割并行电磁暂态数字仿真方法》 中使用了并行技术来提高用于实时仿真的电磁暂态计算的速 度。 但是单纯地依靠电磁暂态并行处理技术还不能够满足实际规模的电力系统计算的需要。
4)频率适应暂态仿真 (FAST : Frequency-Adaptive Simulation of Transients) 说 明 书
电力系统通常运行在一个基频上 （50Hz或 60Hz )， 类似机电震荡的扰动通常只会使系统中的电压和电 流发生较低的频率偏移， 使得电力系统中的电压和电流成为以基频为中心的窄带信号； 这类似于通讯系统 中的信号具有很高的载频， 而被调制的信号却只有很低的频率 （我们真正关心的信息却常.是.这些频率较低 的信号）。 电磁暂态程序使用的是三相瞬时值的表示方式， 除了这些频率较 ¾的被调制信号， 它还需要计 算频率很高的载波信号。 由于载波频率很高， 变化速度快， 这就使得传统的电磁暂态计算程序只能使用较 小的计算步长。 为了增大电磁暂态计算程序使用的计算步长， 提高计算速度， 一些学者通过希尔伯特变换 将原系统中的信号变换为以直流 （0Hz ) 为中心的窄带信号后， 再在变换域中进行电磁暂态计算； 由于变 化域中的信号频率很低， 可以在不损失精度的情况下使用很大的计算步长来计算 （通常可以使计算步长扩 大几百倍）， 最后通过反变换就可以得到原系统的计算结果。 如（文献 "Kai Strunz, Rachel Shintaku, and Feng Gao. Frequency - Adaptive Network Modeling for Integrative Simulation of Natural and Envelope Waveforms in Power Systems and Circuits. IEEE Trans. Circuits Syst. I, Fun dam. Theory Appl. , vol. 53, no. 12, pp. 2788 - 2803, Dec. 2006. " )讨论了这种方法在电力系统和电路仿真中的应用。 然 而这类方法的基础使用的都是基于希尔伯特变换的频域变换法， 当计算系统的快动态过程时， 可能发生混 淆现象 (aliasing effect)； 频域变换的正变换和反变换计算较复杂， 效率低； 而且当原系统中的信号中 含有直流偏置分量时， 基于希尔伯特变换的频域变换方法很难将原系统中的信号变换为以直流为中心的窄 带信号。：
本发明提出了一种新的变换， 不同于已有的频率适应暂态仿真方法以频域变换为基础， 新变换是一种 时域变换方法。 新变换具有物理概念明晰， 计算效率高， 不存在频域变换常存在的混淆现象， 能够处理含 有直流偏置分量的信号等优点。 基于此变换， 建立了一套适用于电磁暂态仿真类程序的方法， 可以在不损 失计算精度的情况下使计算速度有显著提高。 由于用于电力系统的电磁暂态仿真程序本质上是基于电路仿 真建立的， 因此适用于电力系统电磁暂态仿真程序的此方法， 也可以用于改进电路仿真程序（如： SPICE) , 并亦可用于其它领域的仿真和计算 （如通讯、 控制）。
发明内容
本发明提供了一种用于提高系统动态过程仿真的时域变换法， 其思路如图 1所示。 根据被仿真系统中 的信号是窄带信号的特点， 引入原信号的微分或积分作为其正交信号， 通过时域变换构造随时间变化更缓 慢的新变量； 可以在不损失计算精度的情况下在变换域内使用更大的计算步长进行数值求解， 再通过反变 换就可以得到原信号的解； 对于电磁暂态仿真类程序， 在具体实现时可以通过变换域内的计算公式推导原 域内的离散化公式， 并在原域内直接实现。 下面详述之：
根据被仿真系统的信号是否含有直流偏置分量， 需要采用不同的时域变换， 因此将分无直流偏置分量 和含直流偏置分量两种情况分别进行叙述。
1.无直流偏置分量
电力系统中的信号通常为被慢变信号调制的正弦或余弦函数， 其具有以下的形式：
Figure imgf000004_0001
其中：
'( | « o， (0i « 0 (2)
< 为系统的基频， 在机电暂态所使用的相量表示法中， 由幅值 ^(t)和相角 所刻画。 假定 x(0为如下微分方程的解：
Figure imgf000004_0002
(3) 式可以使用某种数值方法直接求解 （例如梯形积分法）。 然而因为受限于信号的基频， 此时所使用的 步长只能很小。 事实上， 数值求解一个微分方程所使用的步长取决于被求解的解 x(t)随时间变化的速度， χ(0变化速度越快， 所需要的计算步长就越小； 所以本发明的思路是通过时域变换， 将被求解的信号 x(t) 变换为变化更为缓慢的信号。
在给出变换之前，需要先引入一个正交于 X(t)的信号，可以通过对 X(t)进行微分或者积分的方式得到, 这里我们使用微分的方法， 通过对 (3)式两边进行微分可得:
„ df df ,
χ = 1 X (4) dt dx
定义新变量: 说 (5) 由此我们得到一个新的微分方程组:
x' = f(t,x书)
(6) y y
dt dx 其中， 为一个新的仿真参数， 它的作用等同于频率适应暂态仿真 FAST中的偏移频率， 通常可以设为系 统的基频； 求解 (6)式将会同时得到 0和 (0的解。
下一步， 我们将通过如下的时域变换构造两个新的变量：
Figure imgf000005_0001
其中， 变换矩阵 Τ(ί)的定义如下:
Figure imgf000005_0002
如果 x(t)具有 (1)式的形式， 则有: y(t) =— A'(t) cos [ω0ί + φ(ί)] - Α ΐ)
Figure imgf000005_0003
(9) 如果 ω5 =ω0， 则有： u(t) (10)
Figure imgf000005_0004
由 (2)式可得：
w( ~^( cos[^( ] (11) 同理可得：
ν{ΐ)^Α{ί) η[φ{ΐ)} (12) 说 明 书 由此可知， 新得到的信号 "(0 和 v(t) 比原信号 χ(0和^: 变化更缓慢，变换 T(t)的作用相当于使 X-Y 参考坐标系中的坐标轴以角速度 os反方向旋转， 从而抵消了原信号 和 )中的基频 ωϋ 。 例如若信号 0如下式所示： x(t) = e 2 cos(&)0i + 0.5ί) (13) 其中， ω0 = 50 χ 2π , x(t)的波形如图 2所示， 对应的 和 v(t)则如图 3所示； 从中可以看出， 相比 于 jc(t)， u(t) 和 v 的变化速度慢很多。
由 (7)式我们可以得到相应的反变换为:
x(t) u{t)
： τ一1 (o (14) yi ) At) 将 (14)式代入 (6)式可以得到一个关于 w(t)和 v(t)的微分方程组:
Figure imgf000006_0001
其中， M(t) 禾卩 v(t)是比 χ(ί)变化更加缓慢的信号， 因此相比于 (3)式， （15)式的数值求解可以采用更大的 步长。为了得到原微分方程 (3)式的解， 可以先在变换域中使用某种数值解法（如梯形积分法）求解 (15)式， 然后再使用 (14)式将解从变换域（U-V坐标系）反变换到原域（Χ-Υ坐标系）， 从而可以得到原微分方程的 解。 ，
然而， 对于电磁暂态仿真这类程序， 有更好的实现方式， 为此引入矩阵 R(t)如下：
Figure imgf000006_0002
矩阵 T(t)和 R(t)具有如下的一些性质- τ_1( = τ( (17)
R(a + b) = R(a)R(b) (18)
R"l( - R(- (19) T(t + h) = R(h)T(t) = T(t)R(-h) (20)
"(0 = ω5Ί{ί +— ) = 0sR(--)T(t) = ^T(t)R(-— -) (21)
2cos 2ω3 2ω8 假定使用梯形积分法求解关于 (t) 和 v(t)的微分方程组 (15)， 则有: 说 明 书
Figure imgf000007_0001
Figure imgf000007_0005
其中 h为计算所使用的步长。 通过对 (7)式两边进行微分有:
f(t,x(t))
"'(0
¾Τ(/ +— ) (23) ν'( 2ί¾ y{t)
cos dt dx
将 (7)和 (23)式代入 (22)式可得
Figure imgf000007_0002
+
对 (24)式应用 T(t)和 R(t)的性质 (17)~ (21)可得:
2
+
Figure imgf000007_0006
Figure imgf000007_0003
(25)式在原域 （X-Y坐标系） 内给出了一种新的离散化方法。
2.含直流偏置分量
在电力系统和电路系统中， 其信号除了含有被慢变信号调制的正弦或余弦函数分量， 还常含有直流偏 置分量， 如发电机稳态运行情况下出口短路时的定子电流。 此时， 其信号具有如下的形式：
x(t) = D(t) + A(t) cos [ω0ί + φ(ί)] (26) 其中 Z)(t)表示信号中含有的直流偏置分量， 而且满足：
| '(t)| « ω0 (27) 对于此类信号，，其处理方式与无直流偏置分量的思路是一样的， 为此除 7 t)， 还需要引入另外一个正交 的信号。 通过对 (4)式两边进行微分可得：
„, d2f d2f ,
Χ'" = ~ + ~— + (28) dt dtdx 若 /(t， )的混合二阶偏导数连续， 则有： d2f d'
(29) dtdx dx
并有:
Figure imgf000007_0004
dt dtdx 说 明 书 对于实际的系统， /(t,x)的混合二阶偏导数常是连续的， 因此下面的叙述将采用 (30)式， 如果不满足， 使 用 (28)式亦可使用同样的方式得到相似的结果。 定义新变量： z = y (31) 由此得到如下的一个新微分方程组:
Figure imgf000008_0001
和无直流偏置分量时一样， 求解 (32)式将会同时得到 x(0、 y(0和 z(t)的解， 但是为了使用更大的步长， 并不会直接求解 (32)式， 而是引入变换， 在变换域中进行数值求解。
我们将通过如下的时域变换构造三个新的变量：
Figure imgf000008_0002
其中， 变换矩阵 T(t)的定义如下:
1 0 1
τ( = 0 一 sin(i¾t) - COS(i¾t) (34) 0 -cos(a>st) sin(<¾t) 如果 x(0具有 (26)的形式， 且满足 (2)和 (27)式， 则有: p{t) » D{t)
u{t) » A{t) cos [¾0 (35) v(t) « A t) sin [φ< 由此可知， 新得到的信号 (t)、 w(t)和 v(t) 比原信号 χ( )、 ： 和 z(t)变化更为缓慢。
由 (33)式可以得到相应的反变换为：
(36)
Figure imgf000008_0003
其中:
1 cos(<¾ ) - \η(ω8ί)
τ-'( = 0 - sin(¾ ) - cos( 0 (37)
0 一 cos(i¾ $1η(ω8ί) 将 (36)式代入 (32)式可以得到一个关于 p(t)、 u{t)和 ν(ί)的微分方程组:
Figure imgf000009_0001
类似地， 相比于 (3)和 (32)式， （38)式的数值求解可以采用更大的步长。 为了得到原微分方程 (3)式的解， 可 以先在变换域中使用某种数值解法（如梯形积分法）求解 (38)式， 然后再使用 (38)式将解从变换域（P-U-V 坐标系） 反变换到原域 （X-Y-Z坐标系）， 从而可以得到原微分方程的解。
然而， 对于电磁暂态仿真这类程序， 同样有更好的实现方式， 为此分别引入矩阵 L(t)和 R(t)如下. - L(t) (39) t)
R( (40)
Figure imgf000009_0002
T(0、 L(t)和 R(t)具有如下的一些性质： R(a + b) = R(a)R(b) (41)
R-'( = R(- (42) L(a + b) = L(a)L(b) (43)
L-'( = L(-i) (44) L(/z)T(r) = T(t + h) = T(t)R(-h) (45)
"( = ¾FT(t +—-)F (46)
L(h)T(t) = T(t + h) = T(t)R(-h) (47) DTK (48) 其中:
Figure imgf000009_0003
说 明 书
0 0 0
F 0 1 0 (50) 0 0 1
类似于无直流偏置分量时的情况， 如果使用梯形积分法在变换域（P-U-V坐标系） 内求解微分方程组 (38), 并利用 (41)~ (48)的性质， 亦可以得到原域 （X-Y-Z坐标系） 内的离散化公式：
X„ = R(h)Xl ¾R(A -— )X„.1 + R(A)X:_ ¾R(-— )X„ + x; (51)
2¾ 2»s
其中:
Figure imgf000010_0001
附图说明
图 1 给出了本发明思路的原理图。
图 2 给出了示例信号 x(t)随时间变化的波形。 图 3 为示例信号 x(t)对应的 w(t)和 v(t)随时间变化的波形。
图 4给出了一个实例的原网络。
图 5 为示例原网络对应的 1阶微分扩展网络。
图 6 为示例系统的时域仿真等效电路。
具体实施方式 下面以 (25)和 (51)式的实现为例来介绍具体的实现方式， 其它的情况下 （如： ： 0采用对 0积分而 不是微分的方式构造， （15)式和 (38)式采用其它数值解法而不是梯形法求解） 可以使用类似的方式实现。
在具体实施时，通常都会以 (25)或 (51)式为基础在支路水平上建立各种网络元件的元件模型，元件模型 提供了一个联系支路任意时刻电压、 电流和历史值的代数方程。 为了便于叙述， 首先引入原网络的 "n阶 微分网络"和 "n阶微分扩展网络"的概念。 所谓原网络的 "n阶微分网络"就是指给出原网络 n阶微分 解的网络， "0阶微分网络"就是原网络本身。 将原网络的 0~n阶微分网络组合起来就得到了原网络的 "n 阶微分扩展网络"。
例如原网络如图 4所示， 原网络中元件的参数如表 1所示。 原网络的 " 1阶微分扩展网络"如图 5所 示， 图中的左半部分跟原网络是完全一样的 （即 "0阶微分网络 ")， 右半部分则是原网络的 " 1 阶微分网 络"。 其中 Rd的电阻、 1^的电感和 Cd的电容分别等于原网络中 R的电阻、 L的电感和 C的电容， " 1阶 微分网络"的其它参数被总结在表 2中。 说 明 书
表 1 原网络参数
Figure imgf000011_0005
表 2 1阶微分网络的参数
将 (25)式应用于图 5中的 " 1阶微分扩展网络"就可以得到如图 6所示的等效电路， 其中的参数如下: π
G1 =— Κ Ε K (54) 2L 2 2ω 's,
L和： 的历史电流为:
(55)
Figure imgf000011_0007
Figure imgf000011_0008
Figure imgf000011_0001
C和 (^的历史电流为:
(57)
Figure imgf000011_0002
Figure imgf000011_0009
Figure imgf000011_0010
1 0
K (58)
0 Ι / ω,
(54)~(57)式中的矩阵 R如（16)式所定义， E为 2阶单位矩阵。 （54)~(58)给出了电感和电容基于 (25)式的 等效元件模型， 其它网络元件的等效模型也可以类似地得到。
类似地， 基于 (51)式也可以建立相应的等效元件模型 （需要先建立 "2 阶微分扩展网络")， 电感和电 容基于 (51)式的等效元件模型的参数如下：
K (59)
21 I 2
Figure imgf000011_0003
GC =— K-1 E K (61) h
Figure imgf000011_0004
说 明 书
Figure imgf000012_0001
(59Η62)式中的矩阵 R如 (40)式所定义， E为 3阶单位矩阵， 和 1^分别为 " 2阶微分扩展网络" 中的三个电感（原网络中的电感 L， 1阶微分网络中对应的电感！^和 2阶微分网络中对应的电感 Ldd ) t„―、 时刻的电压和电流组成的向量， 禾 分别为 "2阶微分扩展网络"中的三个电容（原网络中的电容，
1阶微分网络中对应的电容 Cd和 2阶微分网络中对应的电容 Cdd ) t^时刻的电压和电流组成的向量。
" n阶微分网络"中对应的电感的初始电流和电容的初始电压， 采用如下的方法计算： 将 "n阶微分 网络"中的电感替换为恒流源（电流为 " n-1阶微分网络"中对应电感的初始电流）， 而电容替换为恒压源 (电压为 "n-1阶微分网络"中对应电容的初始电压）， 求解这个替换后的网络将会得到对应恒流源的电压
(记为 )和对应恒压源的电流 (记为 ie )， 则 "n阶微分网络"中相应电感的初始电流 和电容的初 始电压 可以由下式求得 （设电感的电感值为 L， 电容的电容值为 C):
h i ) = UL I L (64) uc' {to) = ic I C (65) 为了得到任意两个计算点之间的信号， 可以通过 (7)或 (33)式先将两个计算点的原信号变换到 U-V或 P-U-V坐标系下， 由于变换域中的这些量随时间的变化速度较慢， 通过简单的线性插值就可以得到两个计 算点之间的 p(t)、 W(t)和 v(t)， 并且具有很高的精度， 最后通过 (14)或 (36)式的反变换公式就可以得到两 个计算点之间的信号值。
有时常需要计算信号中正弦或余弦分量的幅值和相角， 这可以先通过 (7)或 (33)式计算 U-V或 P-U-V 坐标系下相应的；? (0、 和 v(0。 幅值^: t)可以通过如下的公式得到：
Figure imgf000012_0002
而直接通过下式求得的相角 (t)会出现间断点： φ(ί) = arctan [v(i) I
Figure imgf000012_0003
(67) 为了得到随时间连续变化的相角， 采用如下的方法： 设 ， v„为当前时刻的 U-V V„— i为上一个时 步的 U-V量， 首先使用下式求 一时步到当前时刻相角的变化量：
Figure imgf000012_0004
说 明 书
Αφη = arctan [Δνπ I A n ] (69) 当前的相角 φη就可以使用上一个时步的相角 φη—\计算得到:
Figure imgf000013_0001
由于两个时间点之间相角的变化不会太大，这种方法计算的相角的连续性很好，而且不会局限于 [→τ, π]的 区间内。
对于开关这类元件， 为了使其微分也连续变化， 开关的电阻使用如下方式在一个较短的时间内从闭合 电阻转变到开态电阻或者反之。 当开关在 t。时刻开始经过 At的时间由闭合状态变为打开状态时， 开关的 电阻按照如下的公式变化：
rON， '— ^0
(f—― r τΓ
rON + rOFF . Ί
rit) = ！ (rON - rOFF、 ) CQ^厂 π
—{t t0 < t < t0 + At (71)
2 2 At
['OFF ' t≥t。+
当开关在^时刻开始经过 At的时间由打开状态变为闭合状态时， 开关的电阻按照如下的公式变化:
Figure imgf000013_0002
其中 r0N和 r0FF分别为开关的闭合电阻和开态电阻。

Claims

权 利 要 求 书
1. 一种用于提高系统动态过程仿真的时域变换法， 其特征包括： 根据被仿真系统中的信号是窄带信号的 特点， 引入原信号的微分或积分作为其正交信号， 通过时域变换构造新的变量； 相比于原信号及原信号的 微分或积分量, 新变量随时间变化更缓慢， 从而可以在不损失计算精度的情况下在变换域内使用更大的计 算步长进行数值求解， 再通过反变换就可以得到原信号的解：
a) 当系统中的信号为被慢变信号调制的正弦或余弦函数， 而不包含直流偏置分量时， 新变量构造的特征 是采用如下的公式构造新变量：
Figure imgf000014_0003
其中， x(t)为原信号， 0为原信号的一次微分， 和 v(t)为新构造的变量， ^为偏移的频率， 通常 设定为系统的基频； 通过对变换矩阵的一些简单变形构造新变量的方式本质是一样的， 如：
一 sin(i¾t) - cos(<¾t) x(t)
— cos(i¾t) sin((¾t)
Figure imgf000014_0004
y(t)
b) 当系统中的信号包含直流偏置分量时， 新变量构造的特征是采用如下的公式构造新变量:
'pit) 'ι 0 1 一 "
u(t) - 0 一 sin(6¾t) -COS(iDst)
_ v( . 0 一 COS((¾t) sin(iy5t) 其中， x(t)为原信号， y(t)和 分别为原信号的一次微分和二次微分， ρ(ή、 w(t)和 v(t)为新构造的 变量， 同样通过对变换矩阵的一些简单变形构造新变量的方式本质是一样的。
2. 一种适用于电磁暂态计算类程序实现的方法， 其特征在于： 首先在变换域内 （U-V坐标系或 P-U-V坐 标系） 使用某种微分方程数值解法离散化 （如梯形积分法）， 得到变换域内的计算公式， 然后利用反变换 以及引入的 R(t)和 L(t)矩阵的性质得到原域 （原坐标系： X-Y坐标系或 X- Y-Z坐标系） 内的离散化公式； 设原信号由微分方程 x' = /(t, x)所描述， 贝 iJ : a) 当系统中的信号无直流偏置分量时
Figure imgf000014_0001
在变换域内使用梯形积分法求解， 得到的原域 （χ-γ坐标系） 内的离散化公式为:
Figure imgf000014_0005
Figure imgf000014_0002
b) 当系统中的信号含有直流偏置分量时 权 利 要 求 书
1 sin((»st) l— cos(i¾t) 1 0 0
R( = 0 cos(ft>st) sin(<¾t) L( = 0 cos(<¾t) sin((«st)
0 —sin(<¾t) cos(<¾t) 0 一 sin(i¾t) cos(<¾t)
在变换域内 ， 得到的原域 （X- Y-Z坐标系） 内的离散化公式为:
X
其中， Χ„
Figure imgf000015_0001
3. 适用于电磁暂态计算类程序实现的电感和电容的等效元件模型， 其特征在于： 首先得到原网络的 "n 阶微分扩展网络"， 然后利用在变换域内使用梯形积分法得到的原域（X-Y坐标系或 X- Y-Z坐标系）内的离 散化公式 （权项 2中的离散化公式）， 建立电感和电容的等效支路模型：
a) 当系统中的信号无直 电感对应的导纳为： GL
Figure imgf000015_0002
电感对应的历史电流为：
π
Figure imgf000015_0003
2 2ω5
Figure imgf000015_0005
2C .
电容对应的导纳为： Κ
2 2ων
电容对应的历史电流为:
2 2ω8
Figure imgf000015_0006
Figure imgf000015_0007
1 0
其中 K R(t)如权项 2. a中所定义， E为 2阶单位矩阵;
0 X I
b) 当系统中的信号含有直流偏置分量时 电感对应的导纳为： GL K
Figure imgf000015_0004
电感对应的历史电流为： ι^ = Κ-〔Ε - ：、 -R ( E+^C R(h)KlL n 权 利 要 求 书 电容对应的导纳为： Gc 电容对应的历史电流为：
Figure imgf000016_0001
1 0 0 0 -1 0
其中 0 ΙΙω, 0 C = 0 0 -1 R(t)如权项 2.b中所定义， E为 3阶单位矩阵, 0 0 ΙΙωί 0 1 0
4. "n阶微分网络"中电感的初始电流和电容的初始电压的计算方法， 其特征为： 将 "n阶微分网络" 中的电感替换为恒流源（电流为 "n- 1阶微分网络"中对应电感的初始电流）， 而电容替换为恒压源（电压 为 "n-1阶微分网络"中对应电容的初始电压）， 求解这个替换后的网络将会得到对应恒流源的电压（记为 uL )和对应恒压源的电流 (记为 ic , 设电感的电感值为 L， 电容的电容值为 C， 则 "n阶微分网络"中相 应电感的初始电流 iL' (t0) = uL /L , 相应电容的初始电压 uc' (t0) = ic IC。
5. 任意两个相邻计算点之间的信号计算方法，其特征为：先将两个计算点的原信号变换到 U-V或 P-U-V 坐标系下， 再通过线性插值得到两个计算点之间的； t；)、 w(t)和 ν(ί)， 最后再通过反变换公式得到两个计 算点之间的信号的值。
6. 信号中正弦或佘弦分量的幅值和相角的计算方法， 其特征为： 先通过正变换计算 U-V或 P-U-V坐标 系下相应的 p()、 w(t)和 v(t)，则幅值 (t) = "2(i) + v i)；设¾„, 为当前时刻的 U- V量， 为 上一个时步的 U-V量， 当前时刻的相角 则通过如下的公式计算：
V U
Figure imgf000016_0002
Αφη = arctan [Δν„ I Aun ] 其中 (^^为上一个时步的相角。
7. 开关元件的实现方法， 其特征在于： 当开关在 时刻开始经过 的时间由闭合状态变为打开状态时, 开关的电阻按
Figure imgf000016_0003
权 利 要 求 书 当开关在 时刻开始经过 At的时间由打开状态变为闭合状态时， 开关的电阻按照如下的公式变化:
rON + VOFF ！ irOFF rON ) CQg π .
<t<t0+At
~~ 2 ~~ 2 At
r0N, t≥t0 + ^t 其中 r0N和 分别为开关 S勺闭合电阻和开态电且。
PCT/CN2011/000273 2011-02-22 2011-02-22 时域变换法 WO2012113096A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
PCT/CN2011/000273 WO2012113096A1 (zh) 2011-02-22 2011-02-22 时域变换法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/CN2011/000273 WO2012113096A1 (zh) 2011-02-22 2011-02-22 时域变换法

Publications (1)

Publication Number Publication Date
WO2012113096A1 true WO2012113096A1 (zh) 2012-08-30

Family

ID=46720064

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2011/000273 WO2012113096A1 (zh) 2011-02-22 2011-02-22 时域变换法

Country Status (1)

Country Link
WO (1) WO2012113096A1 (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH11243644A (ja) * 1998-02-23 1999-09-07 Mitsubishi Electric Corp 電力系統の想定事故安定度評価方法
CN101382969A (zh) * 2008-10-31 2009-03-11 中国电力科学研究院 一种多步变步长电磁暂态仿真方法
CN101446991A (zh) * 2008-08-15 2009-06-03 中国电力科学研究院 一种电力系统全过程动态仿真的数值积分方法
CN101719182A (zh) * 2009-12-11 2010-06-02 中国电力科学研究院 一种交直流电力系统分割并行电磁暂态数字仿真方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH11243644A (ja) * 1998-02-23 1999-09-07 Mitsubishi Electric Corp 電力系統の想定事故安定度評価方法
CN101446991A (zh) * 2008-08-15 2009-06-03 中国电力科学研究院 一种电力系统全过程动态仿真的数值积分方法
CN101382969A (zh) * 2008-10-31 2009-03-11 中国电力科学研究院 一种多步变步长电磁暂态仿真方法
CN101719182A (zh) * 2009-12-11 2010-06-02 中国电力科学研究院 一种交直流电力系统分割并行电磁暂态数字仿真方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
WANQ CHENGSHAN ET AL., PROGRESSES ON ALGORITHM OF ELECTROMAGNETIC TRANSIENT SIMULATION FOR ELECTRIC POWER SYSTEM AUTOMATION OF ELECTRIC POWER SYSTEMS, vol. 33, no. 7, 10 April 2009 (2009-04-10), pages 97 - 103 *

Similar Documents

Publication Publication Date Title
Rao et al. Real-time electrical load emulator using optimal feedback control technique
Long et al. Gradient descent optimization based parameter identification for FCS-MPC control of LCL-type grid connected converter
Xu et al. Modeling of generators and their controls in power system simulations using singular perturbations
CN104866665A (zh) 基于接口等值与交互的含电力电子设备的混合仿真方法
CN105302956B (zh) 基于fpga的仿真系统及方法
Mudunkotuwa et al. Development of a hybrid simulator by interfacing dynamic phasors with electromagnetic transient simulation
CN109274284B (zh) 一种不平衡电网下并网逆变器的柔性功率控制方法
CN105429484A (zh) 基于任意周期延时的pwm整流器预测功率控制方法及系统
CN104809308B (zh) 一种适用于非对称运行状态的换流器开关函数建模方法
CN1943092A (zh) 转换开关装置和方法
CN110362937A (zh) 一种模块化多电平换流器电磁暂态仿真方法及系统
CN111239491A (zh) 采用实物控制器扰动注入的广义阻抗实时实验测量方法
CN102708225A (zh) 一种交直流大电网电磁暂态仿真的分片调试方法
CN108573094B (zh) 同步发电机的vbr电磁暂态仿真模型的建立方法及系统
CN108959671B (zh) 半桥和全桥型模块化多电平换流器的实时仿真建模方法
Semlyen S-domain methodology for assessing the small signal stability of complex systems in nonsinusoidal steady state
CN104375876A (zh) 一种输入量突变情况下的0+误差免疫电磁暂态仿真算法
CN105449705B (zh) 电力系统中分接头变压器实时负载调压的仿真方法及系统
CN111969639B (zh) 电力电子化电网级联型换流器多时间尺度暂态建模方法
CN103887798A (zh) 有源电力滤波器的反演全局快速终端滑模控制方法
Shujun et al. Modeling for VSC-HVDC electromechanical transient based on dynamic phasor method
WO2012113096A1 (zh) 时域变换法
CN103683874A (zh) 一种基于重复控制的双馈变流器的控制方法
kwon Park et al. Hardware in the loop (HILS) testing of a power electronics controller with RTDS
CN102540907B (zh) 并联型数模综合仿真系统接口和物理仿真子系统接口

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 11859259

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 11859259

Country of ref document: EP

Kind code of ref document: A1