Models

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

1 - Ideal Transformer Model

The Ideal Transformer Model (ITM) is a signal component that splits a circuit into two subcircuits, using a common node as a Point of Common Coupling (PCC), in such a way that a copy of this node is found in the two subcircuits, as shown in Fig. 1, where the copies of the node are denoted as $n$ and $m$. Moreover, the circuits are coupled using a controlled voltage source and a controlled current source, which exchange their interface currents and voltages, respectively, namely the interface signals. This exchange takes place using a ring buffer, on top of which a second ring buffer has been implemented to emulate a co-simualtion using a macro-step, which means that the exchange of interface signals can be made at an interval larger than the simulation’s step size. This second ring buffer is used to implement Zero- and First-Order hold extrapolation methods, while the first ring buffer allows to linearly interpolate the value of the signal at the current time step, in case the delay between both subcircuits is not an integer multiple of the step size.

To add an ITM, users must split the cirtuit and create the copies of the PCC node. An example of this process can be found in the Notebook ITM.ipynb.

To avoid connections of the controlled voltage source with a capacitor, or the controlled current source with an inductor, the resistors $R_{\mathrm{series}}$ and $R_{\mathrm{parallel}}$ are included.

2 - 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

Single-Phase Averaged Voltage Source Inverter with State-Space Nodal Interface (Dynamic Phasor)

The DP::Ph1::AvVoltSourceInverterStateSpace model ports the same grid-following averaged inverter into the dynamic-phasor (DP) domain, as a single positive-sequence complex envelope rather than three abc waveforms. The PLL, power filter, outer power control, and inner current control are baseband and stay real; only the LC filter’s two states are genuine carrier-band envelopes and carry the $-j\omega_n$ shift described in State-Space Nodal.

The terminal input is the PCC voltage envelope

and the state vector is the mixed real/complex-envelope form

where $\psi := \theta_{\mathrm{PLL}} - \omega_n t$ is the PLL angle’s deviation from the nominal carrier phase, tracked instead of the raw, unboundedly growing $\theta_{\mathrm{PLL}}$ for relinearization accuracy, and $V_c$, $I_f$ are complex envelopes replacing EMT’s six abc filter states.

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,

Because the DP envelope already demodulates the carrier, the dq quantities are obtained by rotating the envelope by $\psi$ alone, not by the full absolute angle $\theta_{\mathrm{PLL}}$,

with $v_{c,d} = \operatorname{Re}{V_{c,dq}}$, $v_{c,q} = \operatorname{Im}{V_{c,dq}}$, and likewise for $i_{rc,d}$, $i_{rc,q}$.

The instantaneous active and reactive powers are calculated as

identical in form to EMT’s; DP::Ph1’s own voltage/current scale already represents total power directly, with no three-phase multiplier.

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 a complex envelope, rotating by $\psi$,

The LC filter dynamics carry the envelope’s carrier shift explicitly,

At each simulation step, the nonlinear model is locally linearized into the affine state-space form, packing the 8 real states and the real/imaginary parts of the 2 complex states into one real 12-vector,

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

Source code and examples

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

Three-Phase Averaged Voltage Source Inverter with State-Space Nodal Interface (Dynamic Phasor)

The DP::Ph3::AvVoltSourceInverterStateSpace model extends the grid-following averaged inverter of DP::Ph1 to the three-phase dynamic-phasor (DP) domain. Each phase of the LC filter is represented by an independent complex envelope, $V_{c,a/b/c}$ and $I_{f,a/b/c}$, in contrast to the single positive-sequence envelope employed by DP::Ph1, whereas the controller retains a single positive-sequence $dq$ frame shared by the PLL, the power filter, and the outer and inner control loops. As in DP::Ph1, the control states are baseband quantities and remain real-valued; only the six per-phase filter envelopes are carrier-band quantities, and each carries the $-j\omega_n$ frequency shift introduced in State-Space Nodal.

The terminal input is the PCC voltage envelope of the three phases,

and the state vector concatenates the 8 real control states ahead of the 6 complex per-phase envelopes,

where $\psi := \theta_{\mathrm{PLL}} - \omega_n t$ again denotes the deviation of the PLL angle from the nominal carrier phase, retained as a state to preserve relinearization accuracy. Each per-phase envelope contributes its real and imaginary parts to the packed real vector, yielding 20 real states in total.

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

Model equations

The main extension relative to DP::Ph1 is the per-phase projection onto, and redistribution from, the single positive-sequence $dq$ control frame. The three capacitor-voltage envelopes are projected onto a single positive-sequence phasor,

and the PCC input $\underline{U}$ is projected identically, so that the coupling-current envelope seen by the controller is $\underline{I}_{rc} = (\underline{V}_c - \underline{U})/R_c$, with positive current again denoting injection from the inverter into the grid. The $dq$ quantities are obtained by rotating the projected envelopes by $\psi$,

with $v_{c,d} = \operatorname{Re}{V_{c,dq}}$, $v_{c,q} = \operatorname{Im}{V_{c,dq}}$, and analogously for $i_{rc,d}$ and $i_{rc,q}$. Taken together, the $1\times 3$ projection, the scalar $dq$ rotation, and the $3\times 1$ redistribution defined below constitute a rank-one $3\times 3$ Park mapping on the envelope triple, which reduces to the single-envelope relation of DP::Ph1 under balanced operation.

The positive-sequence active and reactive power measurements used by the controller are

with the projection scaling chosen so that $p$ and $q$ match the total three-phase active and reactive powers under balanced operation; under unbalanced operation they are the positive-sequence components seen by the single-frame controller.

The control chain from the PLL through the inner current loop is identical in form to that of DP::Ph1 and operates on the single positive-sequence $dq$ pair. The PLL and power-filter dynamics read

The outer power-control integrators and current references are

and the inner current-control integrators and voltage references are

The single $dq$ voltage reference $V_{\mathrm{ref},dq} = v_{d,\mathrm{ref}} + j v_{q,\mathrm{ref}}$ is redistributed to the per-phase bridge-voltage envelopes through the inverse projection,

so that all three phases are driven by the same positive-sequence command.

The LC-filter dynamics are decoupled per phase within the plant and carry the carrier shift of the envelope explicitly,

the phases being coupled only through the shared control chain, that is, through $V_{\mathrm{ref},p}$.

At each simulation step the nonlinear model is linearized about the current operating point into the affine state-space form, with the 8 real control states and the real and imaginary parts of the 6 complex per-phase envelopes packed into a single real 20-vector,

which is subsequently discretized and stamped into the DP MNA system.

Because the controller operates in a single positive-sequence $dq$ frame, only the positive-sequence component of an unbalanced terminal is regulated. The negative-sequence response is present in the per-phase filter envelopes but is not itself a control state, and the $2\omega_n$ ripple it would otherwise induce in the $dq$ frame is therefore not represented. A dual-sequence controller with a dedicated negative-sequence frame remains the subject of future work.

References

• M. Mirz, S. Vogel, G. Reinke, and A. Monti, “DPsim: A dynamic phasor real-time simulator for power systems,” SoftwareX, vol. 10, art. 100253, 2019. https://doi.org/10.1016/j.softx.2019.100253
• A. Yazdani and R. Iravani, Voltage-Sourced Converters in Power Systems: Modeling, Control, and Applications. Hoboken, NJ: Wiley-IEEE Press, 2010. https://ieeexplore.ieee.org/book/5237659
• X. Gao, D. Zhou, A. Anvari-Moghaddam, and F. Blaabjerg, “Stability Analysis of Grid-Following and Grid-Forming Converters Based on State-Space Model,” in Proc. 2022 International Power Electronics Conference (IPEC-Himeji 2022, ECCE Asia), pp. 422–428. https://ieeexplore.ieee.org/document/9806927

Source code and examples

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

Three-Phase Averaged Grid-Forming Inverter with State-Space Nodal Interface

The EMT::Ph3::SSN_GFM model represents a grid-forming averaged voltage source inverter in the EMT domain. The control structure follows the state-space grid-forming converter of Gao2022 (VSG algorithm loop, voltage loop, current loop with active damping), whose grid-following counterpart in the same paper is the basis for the averaged inverter above; the inner voltage/current control and LC filter modeling follow Yazdani2010. Like the grid-following inverter above it is a variable state-space nodal component stamped directly into the MNA system, but instead of a PLL that locks to the grid it carries its own virtual synchronous machine (VSG): the internal angle and voltage magnitude are states driven by active- and reactive-power balance, so the inverter imposes a voltage and can run islanded. The model includes the VSG swing dynamics, a reactive-power/voltage excitation loop, filtered active/reactive power measurement, a cascaded voltage and current controller, a first-order converter/digital-delay approximation, and an LC filter with coupling resistance to the grid node.

The terminal input is the PCC voltage vector

and the 17-element state vector is

where $\theta$ is the VSG angle (there is no PLL), $E$ is the excitation-controlled voltage magnitude, $\xi_{v}$, $\xi_{i}$ are the voltage- and current-loop integrators, and $v_{\mathrm{del}}$ are the two delay states.

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

Control structure

The virtual synchronous machine sets the internal angle from the active-power balance and the internal magnitude from the reactive-power/voltage loop; the cascaded voltage and current controllers then track that internal reference through the LC filter. The dashed grid-connected extensions (virtual impedance, feed-forward scaling, Q-V droop) are described below.

Model equations

The physical grid current, positive for injection into the grid, is

All dq quantities use the VSG angle $\theta$ (amplitude-invariant Park transform $\mathbf{T}(\theta)$),

and the capacitor current is $\mathbf{i}{\mathrm{cap},dq} = \mathbf{i}{f,dq} - \mathbf{i}_{g,dq}$. Because the Park transform is amplitude invariant, three-phase instantaneous power carries the factor $3/2$,

and the PCC voltage magnitude is $U_{\mathrm{pcc}} = \sqrt{v_{c,d}^2 + v_{c,q}^2}$.

The measurement filters are first-order lags,

The VSG swing equation sets the angle from the active-power balance,

with virtual inertia $J$ and damping $D$. The reactive-power/voltage excitation controller sets the internal magnitude,

an integral law on the reactive error with a voltage-droop term. The excitation defines the dq voltage reference; in the islanded model it is aligned with the d-axis,

The voltage controller integrates the voltage error and forms the current reference with the capacitor-current feed-forward and dq decoupling,

The current controller integrates the current error and forms the converter voltage reference, with inductor decoupling and optional active damping on the capacitor current,

A first-order lag approximates the converter/digital delay,

and its output, transformed back to abc as $\mathbf{v}{\mathrm{inv}} = \mathbf{T}^{-1}(\theta),[v{\mathrm{del},d}\ v_{\mathrm{del},q}]^\top$, drives the LC filter,

Grid-connected control extensions

The equations above describe the islanded inverter. Three opt-in extensions adapt it to a stiff grid; each defaults to the value that recovers the islanded model exactly, so the eigenstructure is unchanged unless a setter is called.

Virtual output impedance. A virtual impedance $Z_v = R_v + jX_v$ is subtracted from the excitation to form the voltage reference, using the filter current $\mathbf{i}_{f,dq}$,

i.e.

$Z_v = 0$ recovers $v_{d,\mathrm{ref}} = E,\ v_{q,\mathrm{ref}} = 0$. A finite $R_v$ adds a current-proportional term opposing motion, damping the power-synchronization loop on a stiff grid at the electrical timescale, an alternative to raising the mechanical damping $D$. The drop is taken off the filter-current state $\mathbf{i}_f$ rather than the algebraically reconstructed grid current $\mathbf{i}_g = (\mathbf{v}_c-\mathbf{u})/R_c$; the latter would multiply the reference by a factor $\propto 1/R_c$, amplifying state and linearization error.

Grid-current feed-forward scale. A scalar $\kappa$ scales the grid-current feed-forward in the current reference,

with $\kappa = 1$ the default full feed-forward.

Proportional reactive-power droop. When a cutoff $\omega_q > 0$ is set, the integral excitation is replaced by a proportional Q-V droop,

a first-order lag with a stable fixed point $E^* = E_{\mathrm{set}} + D_q(Q_{\mathrm{ref}} - Q)$ and pole at $-\omega_q$. On a stiff grid the network fixes $U_{\mathrm{pcc}}$, so the reactive error $Q_{\mathrm{ref}} - Q$ cannot be driven to zero and the integral law $\dot E = K_q(Q_{\mathrm{ref}} - Q) + K_u(U_n - U_{\mathrm{pcc}})$ has no reachable equilibrium (reactive windup); the proportional droop always has one. The setpoint $E_{\mathrm{set}}$ is captured at initialization as the operating magnitude, so $\dot{E} = 0$ when $Q = Q_{\mathrm{ref}}$ at $t = 0$.

Linearization and stamping

The model is nonlinear (Park transforms with the moving angle $\theta$, the $1/\omega$ swing term, the power products). It is not linearized by hand; at each simulation step the state and output Jacobians are computed by central finite differences of the nonlinear functions $\mathbf{f}(\mathbf{x},\mathbf{u}) = \dot{\mathbf{x}}$ and $\mathbf{g}(\mathbf{x},\mathbf{u}) = \mathbf{y}$,

each column $j$ evaluated as $[\mathbf{f}(\mathbf{x}+\delta_j\mathbf{e}_j,\mathbf{u}) - \mathbf{f}(\mathbf{x}-\delta_j\mathbf{e}_j,\mathbf{u})]/(2\delta_j)$ with a mixed relative/absolute step $\delta_j$. Because the grid-connected extensions above all enter through $\mathbf{f}$, they are captured in $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, $\mathbf{D}$ automatically, with no separate matrix code. The affine offsets fix the model to the current operating point,

giving the affine state-space form

which is discretized and stamped into the EMT MNA system. The dq/abc transformations and the nonlinear controls make the local model time varying, so the SSN equivalent is recomputed every simulation step.

Source code and examples

• Source code: header, implementation
• C++ example (IEEE 9-bus mixed SG/GFM/GFL)
• Python notebook

References

• [Gao2022] X. Gao, D. Zhou, A. Anvari-Moghaddam, and F. Blaabjerg, “Stability Analysis of Grid-Following and Grid-Forming Converters Based on State-Space Model,” in 2022 International Power Electronics Conference (IPEC-Himeji 2022 - ECCE Asia), 2022, pp. 422-428. Source of both the grid-following and grid-forming state-space control structures. Its eigenvalue analysis finds grid-following control better suited to a stiff grid and grid-forming control to a weak grid; the grid-connected extensions above (virtual impedance, Q-V droop) are what let the grid-forming model stay stable when connected to a stiff grid.
• [Yazdani2010] A. Yazdani and R. Iravani, Voltage-Sourced Converters in Power Systems: Modeling, Control, and Applications. Hoboken, NJ: Wiley-IEEE Press, 2010. Basis for the inner voltage/current control and LC-filter modeling of both inverters.

3 - 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”.

4 - 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.

5 - Branches

RX-Line

PI-Line

Transformer

6 - Induction Machine

7 - 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

8 - 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