Models

• 1: Power Electronics
• 2: Synchronous Generator Regulators
• 3: Transformer
• 4: Branches
• 5: Induction Machine
• 6: RLC-Elements
• 7: Synchronous Generator

1 - Power Electronics

DPsim provides several averaged power-electronic inverter models for simulation using EMT, DP, and SP network modeling domains.

Three-Phase Averaged Voltage Source Inverter with State-Space Nodal Interface

The EMT::Ph3::AvVoltSourceInverterStateSpace model represents a grid-following averaged voltage source inverter in the EMT domain. It is implemented as a variable state-space nodal component and can therefore be directly stamped into the MNA system. The model includes a PLL, filtered active/reactive power measurement, outer power control, inner current control, and an LC filter with coupling resistance to the grid node.

The terminal input is the PCC voltage vector

and the state vector is

The model output is the interface current injected into the MNA system,

Model equations

The controller uses the opposite current direction, i.e. positive current denotes inverter injection into the grid,

The Park transformation with PLL angle $\theta_{\mathrm{PLL}}$ is used to obtain dq quantities,

The instantaneous active and reactive powers are calculated as

The PLL and power-filter dynamics are

The outer power-control integrators and current references are

The inner current-control integrators and voltage references are

The reference voltage is transformed back to abc coordinates,

The LC filter dynamics are

At each simulation step, the nonlinear model is locally linearized into the affine state-space form

which is then discretized and stamped into the EMT MNA system.

Source code and examples

• Source code: header, implementation
• C++ example
• Python notebook

2 - Synchronous Generator Regulators

In DPSim, synchronous generator control systems are solved separately from the electric network. The outputs of the electric network (active and reactive power, node voltages, branch currents and rotor speed of synchronous generators) at time $k- \Delta t$ are used as the input of the controllers to calculate their states at time $k$. Because of the relatively slow response of the controllers, the error in the network solution due to the time delay $\Delta t$ introduced by this approach is negligible.

Exciter

DC1 type model is the standard IEEE type DC1 exciter, whereas the other model is a simplified version of the IEEE DC1 type model. The inputs of the exciters are the magnitude of the terminal voltage of the generator connected to the exciter $v_h$ and the voltage reference $v_{ref}$, which is defined as a variable since other devices such as over-excitation limiters or power system stabilizers (PSS) modify such reference with additional signals. At the moment, no over-excitation limiters have been implemented in DPSim so that the reference voltage is given by: $$ v_{ref}(t) = v_{ref,0} + v_{pss}(t) $$ where $v_{ref,0}$ is initialized after the power flow computations and $v_{pss}(t)$ is the output of the (optional) PSS connected to the exciter. The output of the exciter systems is the induced emf by the field current at $t=k + \Delta t$: $v_{ef}(k + \Delta t)$ (sometimes the alternative notation $e_{fd}(k + \Delta t)$ is used).

IEEE Type DC1 exciter model

$$ T_{R} \frac{d}{dt} v_{R}(t) = v_{h}(t) - v_{R}(t) $$

$$ T_{b} \frac{d}{dt} v_{b}(t) = v_{ref} - v_{R}(t) - v_{f}(t) - v_{b}(t), $$ $$ T_{a} \frac{d}{dt} v_{a}(t) = K_{a} v_{in}(t) - v_{a}(t), $$ $$ T_{f} \frac{d}{dt} v_{f}(t) - K_{f} \frac{d}{dt} v_{ef}(t) = -v_{f}(t), $$ $$ T_{ef} \frac{d}{dt} v_{ef}(t) = v_{a}(t) - (K_{ef} + sat(t)) v_{ef}(t), $$ where $v_h$ is the module of the machine’s terminal voltage, and $v_{in}$ is the amplifier input signal, which for the IEEE Type DC1 is given by: $$ v_{in}(t) = T_{c} \frac{d}{dt} v_b(t) + v_b(t). $$

The ceiling function approximates the saturation of the excitation winding: $$ sat(t) = A_{ef} e^{(B_{ef} | v_{ef}(t) | )} $$

The set of differential equations are discretized using forward euler in order to solve it numerically, which leads to the following set of algebraic equations: $$ v_R(k + \Delta t) = v_R(k) + \frac{\Delta t}{T_R} ( v_h(k) - v_R(k) ), $$ $$ v_b(k + \Delta t) = v_b(k)(1 - \frac{\Delta t}{T_b}) + \frac{\Delta t}{T_b} ( v_{ref}(k) - v_R(k) - v_f(k)), $$ $$ v_{in}(k + \Delta t) = \Delta t \cdot \frac{T_c}{T_b} (v_{ref}(k) - v_R(k) - v_{f}(k) - v_b(k)) + v_b(k+1), $$ $$ v_a(k + \Delta t) = v_a(k) + \frac{\Delta t}{T_a} ( v_{in}(k) K_a - v_a(k) ), $$ $$ v_f(k + \Delta t) = (1 - \frac{\Delta t}{T_f}) v_f(k) + \frac{\Delta t K_f}{T_f T_{ef}} ( v_{a}(k) - (K_{ef} + sat(k)) v_{ef}(k) ), $$ $$ v_{ef}(k + \Delta t) = v_{ef}(k) + \frac{\Delta t}{T_{ef}} ( v_{a}(k) - (sat(k) + K_{ef}) v_{ef}(k)), $$ $$ sat(k) = A_{ef} e^{(B_{ef} | v_{ef}(k) | )} $$

Since the values of all variables for $t=k$ are known, $v_{ef}(k+1)$ can be easily calculated using the discretised equations, which is carried out in the preStep function of the generator connected to each exciter.

The initial values of all variables, which are used in the first simulation step, are calculated assuming that the simulation starts in the steady. This is equivalent to assume that all derivative are equal to zero, which leads to: $$ v_R(k=0) = v_h(k=0), $$ $$ v_f(k=0) = 0 $$ $$ v_a(k=0) = K_{ef} v_{ef}(k=0) + A_{ef} e^{B_{ef} |v_{ef} (k=0)|} v_{ef}(k=0), $$ $$ v_{in}(k=0) = \frac{v_a(k=0)}{K_a}, $$ $$ v_b(k=0) = v_{in}(k=0), $$ $$ v_{ref}(t=0) = v_{in}(t=0) + v_b(t=0), $$ where $v_h(k=0)$, $v_{ef}(k=0)$ are calculated after the power flow analysis and after the initialization of synchronous machines (see section initialization of SG).

Simplified IEEE Type DC1 exciter model (DC1Simp)

Because the time constants $T_b$ and $T_c$ of the IEEE Type DC1 exciter model are frequently small enough to be neglected, in DPSim a simplified model of this exciter which neglect these time constants is also implemented. The control diagram of this exciter is depicted in Fig. 2 and it is described by the following set of differential equations: $$ T_R \frac{d}{dt} v_R(t) = v_h(t) - v_R(t) $$ $$ T_a \frac{d}{dt} v_a(t) = - v_a(t) + K_a v_{in}(t) $$ $$ T_f \frac{d}{dt} v_f(t) - K_f \frac{d}{dt} v_{ef}(t) = -v_f(t), $$ $$ T_e \frac{d}{dt} v_{ef}(t) = v_a(t) - v_{ef}(t) (sat(t) + K_{ef}) $$ where $v_h$​ is the module of the machine’s terminal voltage, and $v_{in}$​ is the amplifier input signal, which is given by: $$ v_{in}(t) = v_{ref} (t) - v_R(t) - v_f(t) $$ The set of differential equations are discretized using forward euler in order to solve it numerically, which leads to the following set of algebraic equations: $$ v_R(k + \Delta t) = v_R(k) + \frac{\Delta t}{T_R} ( v_h(k) - v_R(k) ), $$ $$ v_{in}(k) = v_{ref}(k) - v_R(k) - v_f(k), $$ $$ v_a(k + \Delta t) = v_a(k) + \frac{\Delta t}{T_a} ( v_{in}(k) K_a - v_a(k) ), $$ $$ v_f(k + \Delta t) = (1 - \frac{\Delta t}{T_f}) v_f(k) + \frac{\Delta t K_f}{T_f T_{ef}} ( v_{a}(k) - (K_{ef} + sat(k)) v_{ef}(k) ), $$ $$ v_{ef}(k + \Delta t) = v_{ef}(k) + \frac{\Delta t}{T_{ef}} ( v_{a}(k) - (sat(k) + K_{ef}) v_{ef}(k)), $$ $$ sat(k) = A_{ef} e^{(B_{ef} | v_{ef}(k) | )} $$

Since the values of all variables for $t=k$ are known, $v_{ef}(k+1)$ can be easily calculated using the discretised equations, which is carried out in the preStep function of the generator connected to each exciter.

The initial values of all variables, which are used in the first simulation step, are calculated assuming that the simulation starts in the steady. This is equivalent to assume that all derivative are equal to zero, which leads to: $$ v_R(k=0) = v_h(k=0), $$ $$ v_f(k=0) = 0, $$ $$ v_a(k=0) = K_{ef} v_{ef}(k=0) + A_{ef} e^{B_{ef} |v_{ef} (k=0)|} v_{ef}(k=0), $$ $$ v_{in}(k=0) = \frac{v_a(k=0)}{K_a}, $$ $$ v_{ref}(t=0) = v_R(t=0) + v_{in}(t=0), $$ where $v_h(k=0)$, $v_{ef}(k=0)$ are calculated using the power flow analysis and after the initialization of synchronous machines (see section initialization of SG).

Static Exciter

$$ C_{a} = \frac{T_{a}}{T_{b}}, \quad C_{b} = \frac{T_{b}-T_{a}}{T_{b}}. $$ and it is described by the following set of differential equations: $$ T_{R} \frac{d}{dt} v_{r}(t) = v_{h}(t) - v_{r}(t) $$ $$ T_{b} \frac{d}{dt} x_{b}(t) = v_{in}(t) - x_{b}(t) $$ $$ T_{e} \frac{d}{dt} e_{fd}(t) = K_{a} v_{e}(t) - e_{fd}(t), $$

Then, the set of differential equations are discretized using forward euler in order to solve it numerically, which leads to the following set of algebraic equations:

$$ v_r(k + \Delta t) = v_r(k) + \frac{\Delta t}{T_R} ( v_h(k) - v_r(k) ), $$ $$ v_{in}(k) = v_{ref}(k) - v_{r}(k), $$ $$ X_b(k + \Delta t) = \frac{\Delta t}{T_{b}} (v_{in}(k) - x_{b}(k)) + x_{b}(k), $$ $$ v_e(k) = K_{a} (C_{b} x_{b}(k) + C_{a} v_{in} (k)) , $$ $$ e_{fd}(k + \Delta t) = \frac{\Delta t}{T_{e}}(v_{e}(k) - e_{fd}(k)) + e_{fd}(k). $$

To consider the saturation of $e_{fd}$ there are two different implementations, which is automatically selected depending of value of the parameter $K_{bc}$:

Standard ($K_{bc}=0$):

$$ e^{}{fd} = e{fd, max} \quad \quad if \quad \quad e^{}{fd} > e{fd, max} \ e^{}{fd} = e{fd, min} \quad \quad if \quad \quad e^{}{fd} < e{fd, min}, $$

where $e^{*}_{fd}$ represents the output of the exciter.

Anti-windup ($K_{bc}>0$): for controllers with an integral component, i.e. also for PID controllers, the so-called “windup effect” can occur when using the standard saturation function. A strategy for limiting the anti-windup effect is shown in Fig. 5.

which means that the input of the differential equation describing $e_{fd}$, $v_{e}$, takes now the following form:

$$ v_{e} = C_{a} v_{in} + C_{b} x_{b} - K_{bc} (e_{fd} - e_{fd}^{*}) $$

The initial values of all variables, which are used in the first simulation step, are calculated assuming that the simulation starts in the steady. This is equivalent to assume that all derivative are equal to zero, which leads to:

$$ v_{r}(t=0) = v_{h}(t=0), $$

$$ v_{e}(t=0) = \frac{e_{fd}(t=0)}{K_{a}}, $$

$$ v_{in}(t=0) = \frac{v_{e}(t=0)}{C_{a}+C_{b}}, $$

$$ x_{b}(t=0) = v_{in}(t=0), $$

$$ v_{ref}(t=0) = v_{in}(t=0) + v_{r}(t=0) $$

Power System Stabilizer (PSS)

PSS is a controller of synchronous generators used to enhance damping of electromechanical oscillations. The PSS1A implemented in DPSim accepts three optional input signals: rotor speed $\omega$, active power $P$, and terminal voltage magnitude $V_h$. The combined input signal is: $$ s(t) = K_w \omega(t) + K_p P(t) + K_v V_h(t) $$ Setting $K_p = K_v = 0$ recovers the speed-only special case. The PSS output $v_{pss}$ at time $t=k$ is a signal used as the input of the AVR to calculate the field voltage at $t=k+\Delta t$, $v_{fd}(k+\Delta t)$. At present, only one PSS is implemented in DPSim which is a simplified version of the IEEE PSS1A type model.

IEEE PSS1A type PSS

The control diagram of this PSS is depicted in Fig. 6. It includes a washout filter and two lead-lag blocks and is described by the following set of differential equations: $$ T_w \frac{d}{dt} v_1(t) = -(s(t) + v_1(t)), $$ $$ T_2 \frac{d}{dt} v_2(t) = (1 - \frac{T_1}{T_2})(s(t) + v_1(t)) - v_2(t), $$ $$ T_4 \frac{d}{dt} v_3(t) = (1 - \frac{T_3}{T_4})\left(v_2(t) + \frac{T_1}{T_2}(s(t) + v_1(t))\right) - v_3(t), $$ $$ v_{pss}(t) = v_3(t) + \frac{T_3}{T_4}\left(v_2(t) + \frac{T_1}{T_2}(s(t) + v_1(t))\right), $$

where $s(t) = K_w \omega(t) + K_p P(t) + K_v V_h(t)$ is the combined input signal and $v_{pss}(t)$ is the output signal used to modify the reference voltage of the AVR.

The set of differential equations are discretized using forward euler in order to solve it numerically, which leads to the following set of algebraic equations: $$ v_1(k + \Delta t) = v_1(k) - \frac{\Delta t}{T_w} (s(k) + v_1(k)), $$ $$ v_2(k + \Delta t) = v_2(k) + \frac{\Delta t}{T_2} \left((1-\frac{T_1}{T_2})(s(k) + v_1(k)) - v_2(k)\right), $$ $$ v_3(k + \Delta t) = v_3(k) + \frac{\Delta t}{T_4} \left((1-\frac{T_3}{T_4})\left(v_2(k) + \frac{T_1}{T_2}(s(k) + v_1(k))\right) - v_3(k)\right), $$ $$ v_{pss}(k) = v_3(k) + \frac{T_3}{T_4} \left(v_2(k) + \frac{T_1}{T_2} (s(k) + v_1(k))\right) $$

Since the values of all variables for $t=k$ are known, $v_{pss}(k)$ can be easily calculated using the discretised equations, which is carried out in the preStep function of the generator connected to each exciter. Then, $v_{pss}(k)$ is used as input of the AVR to calculate the field voltage at time $k+1$. The values $v_1(k+1)$, $v_2(k+1)$, $v_3(k+1)$ are stored and used to calculate the PSS output of the next time step.

The initial values of all variables, which are used in the first simulation step, are calculated assuming that the simulation starts in steady state. This is equivalent to assuming that all derivatives are equal to zero, which leads to: $$ v_1(k=0) = -s(k=0), $$ $$ v_2(k=0) = (1 - \frac{T_1}{T_2})(s(k=0) + v_1(k=0)), $$ $$ v_3(k=0) = (1 - \frac{T_3}{T_4})\left(v_2(k=0) + \frac{T_1}{T_2}(s(k=0) + v_1(k=0))\right), $$ $$ v_{pss}(k=0) = v_3(k=0) + \frac{T_3}{T_4}\left(v_2(k=0) + \frac{T_1}{T_2}(s(k=0) + v_1(k=0))\right), $$

where $s(k=0) = K_w \omega(k=0) + K_p P(k=0) + K_v V_h(k=0)$ is evaluated after the power flow analysis and initialization of synchronous machines (see section initialization of SG). In steady state $\omega(k=0) = 1.0$ (pu), and if $K_p = K_v = 0$ then $v_2 = v_3 = v_{pss} = 0$.

Turbine Governor Models

In DPsim there are two types of Turbine Governor implementations. The Turbine Governor Type 1 implements both the turbine and the governor in one component. In contrast, Steam Turbine and Steam Turbine Governor are implemented as two separate classes and their objects are created independently. Steam/Hydro Turbine and Steam/Hydro Turbine Governor are two blocks that must be connected in series.

The input of the turbine governor models is the mechanical omega at time $t=k-\Delta t$ and the output is the mechanical power at time $t=k$. This variable is then used by the SG to predict the mechanical omega at time $t=k+\Delta t$.

Turbine Governor Type 1

This model includes a governor, a servo and a reheat block. The control diagram of this governor is depicted in Fig. 7 and it is described by the following set of differential equations: $$ p_{in}(t) = p_{ref} + \frac{1}{R} (\omega_{ref} - \omega(t)), $$ $$ T_s \frac{d}{dt} x_{g1}(t) = p_{in}(t) - x_{g1}(t), $$ $$ T_c \frac{d}{dt} x_{g2}(t) = \left(1 - \frac{T_3}{T_c}\right) x_{g1}(t) - x_{g2}(t), $$ $$ T_5 \frac{d}{dt} x_{g3}(t) = \left(1 - \frac{T_4}{T_5}\right) \left(x_{g2}(t) + \frac{T_3}{T_c} x_{g1}(t)\right) - x_{g3}(t), $$ $$ \tau_m(t) = x_{g3}(t) + \frac{T_4}{T_5} \left(x_{g2}(t) + \frac{T_3}{T_c} x_{g1}(t)\right), $$ where $\omega(t)$ is the input signal and $\tau_m(t)$ is the output signal of the governor.

The differential equations are discretized using the forward Euler method, which leads to the following set of algebraic equations: $$ p_{in}(k-\Delta t) = p_{ref} + \frac{1}{R} (\omega_{ref} - \omega(k-\Delta t)), $$ $$ x_{g1}(k) = x_{g1}(k-\Delta t) + \frac{\Delta t}{T_s} \left(p_{in}(k-\Delta t) - x_{g1}(k-\Delta t)\right), $$ $$ x_{g2}(k) = x_{g2}(k-\Delta t) + \frac{\Delta t}{T_c} \left(\left(1 - \frac{T_3}{T_c}\right) x_{g1}(k-\Delta t) - x_{g2}(k-\Delta t)\right), $$ $$ x_{g3}(k) = x_{g3}(k-\Delta t) + \frac{\Delta t}{T_5} \left(\left(1 - \frac{T_4}{T_5}\right) \left(x_{g2}(k-\Delta t) + \frac{T_3}{T_c} x_{g1}(k-\Delta t)\right) - x_{g3}(k-\Delta t)\right), $$ $$ \tau_m(k) = x_{g3}(k) + \frac{T_4}{T_5} \left(x_{g2}(k) + \frac{T_3}{T_c} x_{g1}(k)\right). $$ Since all variables at $t=k-\Delta t$ are known, $\tau_m(k)$ is computed in the preStep of the generator and used to approximate the mechanical equations at time $k+\Delta t$.

Steam Governor

The control diagram of this model is depicted in Fig. 8. This model receives as input the frequency deviation $\Delta\omega = \omega_{ref} - \omega$ from the nominal frequency (normally $50,\text{Hz}$ or $60,\text{Hz}$) and produces the valve opening signal $p_{gv}$ for the turbine. $p_{ref}$ is the mechanical power produced at nominal frequency. The governor implements a lead-lag controller $\frac{K(1+sT_2)}{(1+sT_1)}$ where $K=1/R$ and $R$ is the droop coefficient, followed by a PT1 integrator with embedded rate limiters and an anti-windup loop. To avoid unnecessary dead-beat behaviour, complex transfer functions with more than one pole and zero are decomposed via partial fraction expansion into parallel PT1 elements, as shown in Fig. 9.

Analogous to the static exciter model, the integrator uses an anti-windup strategy as shown in Fig. 10.

The forward-Euler discretised equations are: $$ \Delta \omega (k-\Delta t) = \omega_{ref} - \omega (k-\Delta t), $$ $$ p_{1}(k) = p_{1}(k-\Delta t) + \frac{\Delta t}{T_{1}} \left(\Delta \omega (k-\Delta t) \cdot \frac{T_{1} - T_{2}}{T_{1}} - p_{1}(k-\Delta t)\right), $$ $$ p(k-\Delta t) = \frac{1}{R} \left(p_{1}(k-\Delta t) + \Delta \omega(k-\Delta t) \cdot \frac{T_{2}}{T_{1}}\right), $$ $$ \dot{p}(k-\Delta t) = \frac{1}{T_{3}}\left(p(k-\Delta t) + p_{ref} - p_{gv}(k-\Delta t)\right) - K_{bc} \left(p_{gv}^{}(k-\Delta t) - p_{gv}(k-\Delta t)\right), $$ $$ p_{gv}^{}(k) = p_{gv}^{*}(k-\Delta t) + \Delta t \cdot \dot{p}(k-\Delta t), $$

and

$$ p_{gv}(k) = p_{gv}^{}(k) \quad \text{if} \quad P_{m,\min} \leq p_{gv}^{}(k) \leq P_{m,\max}, \ p_{gv}(k) = P_{m,\max} \quad \text{if} \quad p_{gv}^{}(k) > P_{m,\max}, \ p_{gv}(k) = P_{m,\min} \quad \text{if} \quad p_{gv}^{}(k) < P_{m,\min}. $$

If $T_1 = 0$ the $p_1(k)$ equation is skipped and $p(k)$ is instead: $$ p(k-\Delta t) = \frac{1}{R} \left(\Delta \omega(k-\Delta t) + \frac{T_{2}}{\Delta t} \left(\Delta \omega(k-\Delta t) - \Delta \omega(k-2\Delta t)\right)\right). $$

Assuming the simulation starts in steady state (all derivatives zero, $\Delta\omega(0)=0$), the initial values are: $$ p_{1}(t=0) = 0, \quad p(t=0) = 0, \quad p_{ref} = p_{gv}^{*}(t=0) = p_{gv}(t=0). $$

Steam Turbine

The steam turbine receives the valve opening signal $p_{gv}$ from the Steam Governor and outputs mechanical power $p_m$ to the synchronous generator. It is divided into high-pressure (HP), intermediate-pressure (IP), and low-pressure (LP) stages, each modelled as a first-order lag with time constants $T_{CH}$, $T_{RH}$, $T_{CO}$ respectively. Setting a time constant to zero disables that lag element. The total mechanical power is a weighted sum of each stage: $F_{HP} + F_{IP} + F_{LP} = 1$ must hold. The forward-Euler discretised equations are:

$$ p_{hp}(k) = p_{hp}(k-\Delta t) + \frac{\Delta t}{T_{CH}} \left(p_{gv}(k-\Delta t) - p_{hp}(k-\Delta t)\right), $$ $$ p_{ip}(k) = p_{ip}(k-\Delta t) + \frac{\Delta t}{T_{RH}} \left(p_{hp}(k-\Delta t) - p_{ip}(k-\Delta t)\right), $$ $$ p_{lp}(k) = p_{lp}(k-\Delta t) + \frac{\Delta t}{T_{CO}} \left(p_{ip}(k-\Delta t) - p_{lp}(k-\Delta t)\right), $$ $$ p_{m}(k) = F_{HP} \cdot p_{hp}(k) + F_{IP} \cdot p_{ip}(k) + F_{LP} \cdot p_{lp}(k). $$

Assuming the simulation starts in steady state (all derivatives zero), the initial values are: $$ p_{hp}(t=0) = p_{gv}(t=0), \quad p_{ip}(t=0) = p_{hp}(t=0), \quad p_{lp}(t=0) = p_{ip}(t=0), $$ $$ p_{m}(t=0) = F_{HP} \cdot p_{hp}(t=0) + F_{IP} \cdot p_{ip}(t=0) + F_{LP} \cdot p_{lp}(t=0). $$

Hydro Turbine Governor

The Hydro Turbine Governor receives the frequency deviation $\Delta\omega = \omega_{ref} - \omega$ as input and produces the valve/gate opening signal $p_{gv}$ for the turbine. $p_{ref}$ is the mechanical power produced at nominal frequency. The controller transfer function is $K\frac{1+sT_2}{(1+sT_1)(1+sT_3)}$, where $K=\frac{1}{R}$ and $R$ is the droop coefficient. The transfer function is decomposed into two parallel PT1 blocks as shown in Fig. 13.

The forward-Euler discretised equations are: $$ x_{1}(k) = x_{1}(k-\Delta t) + \frac{\Delta t}{T_{1}} \left(\Delta\omega(k-\Delta t) - x_{1}(k-\Delta t)\right), $$ $$ x_{2}(k) = x_{2}(k-\Delta t) + \frac{\Delta t}{T_{3}} \left(\Delta\omega(k-\Delta t) - x_{2}(k-\Delta t)\right), $$ $$ p^{*}{gv}(k) = \frac{1}{R}\left(A \cdot x{1}(k) + B \cdot x_{2}(k)\right) + p_{ref}, $$

where $$ A = \frac{T_{1}-T_{2}}{T_{1}-T_{3}}, \qquad B = \frac{T_{2}-T_{3}}{T_{1}-T_{3}}, $$

and the output limiter is applied as: $$ p_{gv}(k) = \begin{cases} P_{m,\max} & \text{if } p^{}{gv}(k) > P{m,\max}, \ P_{m,\min} & \text{if } p^{}{gv}(k) < P{m,\min}, \ p^{*}_{gv}(k) & \text{otherwise.} \end{cases} $$

Assuming the simulation starts in steady state (all derivatives zero, $\Delta\omega(t=0)=0$), the initial values are: $$ x_{1}(t=0) = 0, \quad x_{2}(t=0) = 0, \quad p_{ref} = p_{gv}(t=0). $$

Hydro Turbine

The Hydro Turbine receives the gate opening signal $p_{gv}$ from the Hydro Turbine Governor and outputs mechanical power $p_m$ to the synchronous generator. The transfer function is specified by the water starting time $T_W$ and can be represented as the sum of two parallel blocks as shown in Fig. 15.

The forward-Euler discretised equations are: $$ x_{1}(k) = x_{1}(k-\Delta t) + \frac{\Delta t}{0.5,T_{W}} \left(p_{gv}(k-\Delta t) - x_{1}(k-\Delta t)\right), $$ $$ p_{m}(k) = 3,x_{1}(k) - 2,p_{gv}(k). $$

Assuming the simulation starts in steady state (all derivatives zero), the initial values are: $$ x_{1}(t=0) = p_{gv}(t=0), \quad p_{m}(t=0) = p_{gv}(t=0). $$

References

• [1] “IEEE Recommended Practice for Excitation System Models for Power System Stability Studies,” in IEEE Std 421.5-2016 (Revision of IEEE Std 421.5-2005) , vol., no., pp.1-207, 26 Aug. 2016, doi: 10.1109/IEEESTD.2016.7553421.
• [2] F. Milano, “Power system modelling and scripting,” in Power System Modelling and Scripting. London: Springer-Verlag, 2010, ISBN: 978-3-642-13669-6. doi: 10.1007/978-3-642-13669-6.
• [3] F. Milano, A. Manjavacas, “Frequency Variations in Power Systems: Modeling, State Estimation, and Control”. ISBN: 978-1-119-55184-3.
• [4] F. Milano, “Power System Analysis Toolbox: Documentation for PSAT”, ISBN: 979-8573500560.
• [5] M. Eremia; M. Shahidehpour, “Handbook of Electrical Power System Dynamics: Modeling, Stability, and Control”, https://ieeexplore.ieee.org/book/6480471
• [6] A. Roehder, B. Fuchs, J. Massman, M. Quester, A. Schnettler, “Transmission system stability assessment within an integrated grid development process”.

3 - Transformer

2-Winding Transformer

The transformer model is composed of an RL-segment and an ideal transformer. The single line diagram is depicted in the figure below.

If node reduction is not applied, two virtual nodes are created to stamp this model into the system matrix.

Furthermore, the ideal transformer has an additional equation, which requires an extension of the system matrix. The complete matrix stamp for the ideal transformer is

The variable $j$ denotes the high voltage node while $k$ is the low voltage node. $l$ indicates the inserted row and column to accommodate the relation between the two voltages at the ends of the transformer. The transformer ratio is defined as $T = V_{j} / V_{k}$. A phase shift can be introduced if $T$ is considered as a complex number.

4 - Branches

RX-Line

PI-Line

Transformer

5 - Induction Machine

6 - RLC-Elements

EMT Equations and Modified Nodal Analysis

Inductance

An inductance is described by

Integration results in an equation to compute the current at time $t$ from a previous state at $t - \Delta t$.

There are various methods to discretize this equation in order to solve it numerically. The trapezoidal rule, an implicit second-order method, is commonly applied for circuit simulation:

Applying the trapezoidal rule to leads to

This can be rewritten in terms of an equivalent conductance and current source and the number of time steps $k$ with size $\Delta t$.

Hence, components described by differential equations are transformed into a DC equivalent circuit as depicted in the figure below.

Capacitance

The same procedure can be applied to a capacitance. Integration on both side yields

Finally, the equivalent circuit is described by a current source and a conductance.

This equation set is visualized in the figure below.

Hence, the vector of unknowns $\bm{x}$ and the source vector $\bm{b}$ become time dependent and this leads to the system description:

To simulate the transient behavior of circuits, this linear equation has to be solved repeatedly. As long as the system topology and the time step is fixed, the system matrix is constant.

Extension with Dynamic Phasors

The dynamic phasor concept can be integrated with nodal analysis. The overall procedure does not change but the system equations are rewritten using complex numbers and all variables need to be expressed in terms of dynamic phasors. Therefore, the resistive companion representations of inductances and capacitances have to be adapted as well.

Inductance

In dynamic phasors the integration of the inductance equation yields

Applying the trapezoidal method leads to the finite difference equation:

Solving this for $\langle i_L \rangle(k)$ results in the \ac{DP} equivalent circuit model:

with

Capacitance

Similarly, a capacitance is described by as follows

Applying the trapezoidal rule for the capacitance equation leads to the finite difference equation:

The DP model for the capacitance is defined by

with

RL-series element

In dynamic phasors the integration of the inductance equation yields

Applying the trapezoidal method leads to the finite difference equation:

Solving this for $\langle i \rangle(k)$ results in the \ac{DP} equivalent circuit model:

with

7 - Synchronous Generator

Two different synchronous machine models are currently available:

• the full order dq0 reference frame model (EMT, DP) [Kundur, Power system stability and control, 1994]
• and the much simpler transient stability model (DP) [Eremia, Handbook of Electrical Power System Dynamics, 2003]

The machine model is interfaced to the nodal analysis network solver through a current source, which only affects the source vector and not the system matrix Wang2010.

Basic Equations

The equations of the stator and rotor voltages are

where

The flux linkage equations are defined as

The inductance matrices are varying with the rotor position $\theta_r$ which varies with time.

The mechanical equations are:

$\theta_r$ is the rotor position, $\omega_r$ is the angular electrical speed, $P$ is the number of poles, $J$ is the moment of inertia, $T_m$ and $T_e$ are the mechanical and electrical torque, respectively. Motor convention is used for all models.

dq0 Reference Frame 9th Order Model

For stator referred variables, the base quantities for per unit are chosen as follows:

• $v_{s base}$ peak value of rated line-to-neutral voltage in V
• $i_{s base}$ peak value of rated line current in A
• $f_{base}$ rated frequency in Hz

The synchronous generator equations in terms of per unit values in the rotor reference frame become:

where

The flux linkages are:

where

with

The mechanical equations in per unit become:

For the simulation, fluxes are chosen as state variables. To avoid the calculation of currents from fluxes using the inverse of the inductance matrix, the equation set needs to be solved for the fluxes analytically. To simplify the calculations, dq axis magnetizing flux linkages are defined [Krause, Analysis of electric machinery and drive systems, 2002]:

Using the flux linkages results in a simpler equation set for the fluxes:

Dynamic Phasor Model

The fundamental dynamic phasors are similar to the dq0 quantities for symmetrical conditions since both yield DC quantities in a rotating reference frame. The network abc dynamic phasor quantities can be converted to dq0 dynamic phasors by applying the symmetrical components transformation and a rotation.

The angle $\delta$ is the orientation of the dq0 reference frame relative to the abc frame.

The winding currents for positive and zero sequence components can be expressed as

In the dynamic phasor case, the equation for $\frac{d}{dt} \langle \lambda_{0s} \rangle_1$ has a frequency shift. To complete the state model, the magnetizing flux linkages are expressed as:

where

The mechanical equations in dynamic phasors are:

Transient Stability Model