Models

• 1: RLC-Elements
• 2: Branches
• 3: Transformer
• 4: Ideal Transformer Model
• 5: Switches
• 6: Loads
• 7: Sources
• 8: Network Injection and Compensation
• 9: Synchronous Generator
• 9.1: Reduced Order Machine Models
• 10: Synchronous Generator Regulators
• 10.1: Exciters and Power System Stabiliser
• 10.2: Turbines and Governors
• 11: Power Electronics
• 11.1: EMT Ph3 Averaged Voltage Source Inverter
• 11.2: DP Ph1 Averaged Voltage Source Inverter
• 11.3: DP Ph3 Averaged Voltage Source Inverter
• 11.4: EMT Ph3 Grid-Forming Inverter
• 12: Converter Control
• 13: State-Space Nodal Components
• 14: Signal Processing Blocks

1 - 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 $\boldsymbol{x}$ and the source vector $\boldsymbol{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

2 - Branches

Both line models below are composite components: they do not stamp the system matrix directly but are built from resistor, inductor and capacitor subcomponents, each of which contributes its own stamp. See subcomponents for how that composition works, and RLC elements for the stamps of the individual elements.

The transformer is documented separately under transformer.

RX-Line

The RX line represents a line by its series resistance and series inductance only, ignoring the shunt admittance. It is the appropriate choice for short lines, where the charging current is negligible, and it is what the CIM reader produces for an ACLineSegment when no shunt data is present.

The model is composed of a series resistor and a series inductor between the two terminals:

An additional resistor from the inductor terminal to ground is present to make initialisation well posed. It is not part of the physical model.

RxLine exists in DP::Ph1, EMT::Ph3, SP::Ph1 and SP::Ph3.

PI-Line

The PI line adds the shunt admittance of the line, split evenly between the two terminals, which matters once the line is long enough for the charging current to affect the result. The name comes from the shape of the equivalent circuit: a series branch with one shunt branch at each end.

The series branch carries the resistance and inductance as above. Each terminal additionally carries half of the total shunt capacitance and half of the total shunt conductance:

The shunt capacitance and conductance are specified as totals for the line, and the halving between the two ends is part of the model rather than something the user does.

Decoupling Line

The decoupling line is a distributed parameter line based on the Bergeron travelling wave method. Unlike the two models above it is not primarily a fidelity improvement: its purpose is to remove the direct coupling between the two terminals so that the network on either side can be solved as an independent system, which is what makes splitting a network across solvers or across simulators possible.

The method rests on the behaviour of a lossless line. For a line with distributed inductance and capacitance, the quantity $v + Z_c, i$ observed at one end reappears unchanged at the other end one travel time later, and likewise in the opposite direction. Nothing propagates faster than that travel time, so the two ends cannot influence each other within it. The surge impedance and the travel time follow from the line’s total inductance and capacitance,

Each terminal is then represented by a resistance to ground in parallel with a current source. The resistance is $Z_c + R/4$, and the current source carries the history term, whose value depends on the voltage and current recorded at the other terminal one travel time ago. Because that value is already known when the step begins, it enters the system as a constant injection rather than as a coupling into the admittance matrix, and the matrix separates into two blocks that can be factorised and solved independently.

The series resistance is not distributed along the line. It is lumped, with $R/4$ placed at each end and the remainder in the middle of the equivalent, which is why the terminating resistance and the history coefficients carry $R/4$ terms rather than the full $R$.

The travel time is not required to be a whole number of time steps. The recorded quantities are held in a buffer of $\lceil \tau / \Delta t \rceil$ samples and the value one travel time ago is recovered by linear interpolation between the two nearest entries. The one hard requirement is that the travel time exceed the time step; a line whose $\tau$ is shorter than $\Delta t$ cannot decouple anything, and setting one up is rejected rather than silently approximated.

In the dynamic phasor domain the history terms carry an additional rotation $e^{-j \omega_s \tau}$. This is a direct consequence of working with envelopes: a delay of $\tau$ applied to an instantaneous waveform becomes, for the envelope, the same delay together with a phase rotation of the carrier over that interval, as described under dynamic phasors. Note that this rotation is currently evaluated at a fixed 50 Hz rather than at the system frequency in use.

The decoupling is exact for the lossless travelling wave line it is derived from. The error introduced in practice comes from the lumped treatment of the series resistance and from the interpolation of the delayed quantities, and it grows as the time step approaches the travel time.

Choosing between them

Use the RX line when the shunt admittance can be neglected and you want the smaller system matrix, since the PI line introduces additional nodes for its shunt branches. Use the PI line when the line is long enough that its charging current matters, or when you are comparing against a reference tool that models the shunt branch.

Both are lumped parameter models and therefore do not reproduce travelling wave behaviour. Use the decoupling line when you need that behaviour, or when the reason for reaching for a line model is to split the network in the first place. For the domains each model is available in, see model availability.

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.

Why the ideal part needs an extra equation

The ideal transformer imposes two constraints at once: the voltages are in a fixed ratio and the powers on the two sides are equal, which makes the currents inversely proportional to the same ratio,

Neither is a current balance at a node, so neither can be written as an admittance. This is the same situation as an ideal voltage source described under sources: the system is extended with the branch current as an unknown, the constraint occupies the added row, and the added diagonal entry is zero. The asymmetry of the stamp, $-1$ against $T$ in the added column and $1$ against $-T$ in the added row, is exactly the statement that voltage scales by $T$ while current scales by $1/T$ with opposite sign.

Making $T$ complex adds a phase shift, which is how a delta-wye connection is represented without modelling the windings. The magnitude and the angle then carry the tap ratio and the vector group respectively.

Series impedance and the direction of the ratio

The winding resistance and leakage inductance are lumped into one series branch on one side of the ideal part rather than split between the two sides. Referring an impedance across an ideal transformer scales it by $T^2$, so the choice of side is a choice of reference, not an approximation, and the parameters have to be given consistently with it.

The ratio is defined greater than one, from high voltage to low. Supplying it the other way round describes the same physical device but with the two ends exchanged, so a transformer given an inverted ratio has to have its terminal assignment inverted with it to remain the same transformer.

Numerical damping

Connecting an inductive branch between two nodes that have no other path to ground leaves those nodes weakly defined, and the resulting matrix can be poorly conditioned or singular. Small shunt elements at each terminal remove that, at the cost of a negligible current that would not exist in the physical device.

Those elements are sized from the transformer’s rated power, which makes the rating a required parameter rather than documentation. Without a positive rating there is no scale to size them against, and the natural result is an infinite resistance and a zero capacitance whose admittance is not a number. A single such entry propagates through the factorisation and destroys the whole solution, not merely the transformer, so the rating cannot be treated as optional.

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

Why the resistors are necessary

The two failure cases they prevent are the same one seen twice. A voltage source directly across a capacitor over-determines that node: both impose a voltage, and the capacitor’s companion model and the source’s constraint row describe the same quantity. A current source in series with an inductor under-determines the branch in the dual way, since both impose a current. In each case the system matrix becomes singular rather than merely ill-conditioned, so the resistors are a condition for the method to work at all and not a refinement of it.

Their values are a compromise of the kind described under switches. Small enough to be electrically negligible, large enough not to dominate the condition number.

What the delay costs

The exchanged signals are always at least one step old, because each side computes from what the other produced previously. That delay is the reason the two subcircuits can be solved separately at all, and it is also the entire error of the method: the coupled system is not the original circuit but the original circuit with a transport delay inserted at the point of common coupling.

The consequence is that accuracy is governed by how much the interface signals change within one exchange interval, not by how accurately either side is solved internally. Refining the step inside a subcircuit while holding the macro-step fixed improves nothing at the interface.

A macro-step larger than the simulation step makes this explicit, which is the point of the second ring buffer: it is the co-simulation case, where the two sides may be different tools exchanging at a rate neither controls. Zero-order hold holds the last received value for the whole interval; first-order hold extrapolates linearly from the last two. The first is safe and lags; the second tracks a smoothly varying signal better and overshoots at a discontinuity, which is exactly what a fault produces.

Contrast with the alternatives

Three ways of splitting a network appear in this documentation and they differ in what they cost.

Tearing, described under alternative solution methods, is exact: the removed branches are restored within the same step, so the answer matches the intact network. It does not allow the parts to be advanced independently.

The travelling-wave line under branches is exact for the lossless line it derives from, and its delay is physical rather than introduced. It requires that a real line with a travel time longer than the step exists at the splitting point.

The ideal transformer model requires no such line and can split anywhere, and pays for that with a delay that has no physical counterpart. It is the general method and the least accurate of the three.

5 - Switches

A switch in a nodal formulation is not an ideal open or short. Both would be singular: an ideal short shorts two node equations together, and an ideal open leaves a node with no path to ground. Switches are therefore represented by a finite resistance that takes one of two values.

The two-resistance model

The switch contributes a single admittance between its two terminals,

stamped as a conductance between the two terminal nodes, with the usual reduction when one terminal is grounded. Typical values are far apart, of the order of milliohms closed and megohms open, so the switch is a near short or a near open without ever being singular.

The consequence of this choice is that the ratio $R_{open} / R_{closed}$ lands directly in the condition number of the system matrix. Making the contrast arbitrarily large to approach an ideal switch degrades the accuracy of every node voltage in the network, not only those near the switch. The values are a numerical compromise, not a physical measurement.

Because the admittance appears in the system matrix rather than in the right hand side, changing state requires the matrix to be refactorised. This is why a network that switches often costs more than one that does not, even though the model itself is trivial.

Why a step change in resistance is a problem

Opening a switch that carries inductive current asks the network to interrupt that current within one time step. The inductor opposes it, and with the trapezoidal companion model the result is a numerical oscillation across the switch: the current alternates sign at the step frequency and decays slowly, contaminating the solution for many steps after the event.

This is a property of the discretisation, not of the physical circuit. The physical arc that would form across real contacts dissipates that energy; a two-valued resistance has no equivalent mechanism.

The variable-resistance switch

The variable-resistance switch removes the oscillation by refusing to make the change in a single step. On opening, the resistance is multiplied by a fixed factor each step,

until it reaches the target open value, after which it is held there. The current therefore decays geometrically over several steps rather than being interrupted at once, which is close to what an arc does and which the trapezoidal companion model can follow without ringing.

The growth factor is tied to the step size so that the transition covers a comparable interval of time rather than a comparable number of steps. Closing is not ramped: the resistance is taken straight to its closed value, because energising a path through a small resistance does not produce the same interruption problem.

The cost is that the system matrix changes on every step of the transition rather than once, so each of those steps requires a refactorisation. The switch is worth its cost where the interruption is severe, typically a fault applied at a machine terminal or a transformer winding, and unnecessary for ordinary load switching.

Series switching

Where a switch is combined with the series resistance it energises, the two are represented as one element rather than as a switch plus a resistor. This keeps the branch to a single admittance and avoids introducing an internal node that carries no physical meaning and adds an equation to the system.

6 - Loads

A load is specified as an active and a reactive power at a nominal voltage, but a nodal solver needs either an admittance or a current. The two ways of making that conversion behave differently as the terminal voltage moves away from nominal, and the difference matters more than the model’s simplicity suggests.

Constant impedance

The powers are converted once, at the nominal voltage, into a resistance and a reactance,

and the reactance becomes an inductance or a capacitance according to its sign,

The load is then an ordinary passive branch to ground, and it is stamped exactly as the elements it is built from.

Both conversions divide by a power, so a load with zero active power has no defined resistance and one with zero reactive power has no defined reactance. Such a branch is simply absent rather than infinite, which is the correct behaviour but means a load specified with one of the two set to zero is not the load a reader might expect.

The assumption is that consumption follows the square of the voltage. At nominal voltage the load draws exactly $P$ and $Q$; at 0.9 per unit it draws 81 percent of them. For a genuinely impedance like load this is right, and for anything regulated it understates the demand during a depression.

Constant current

The alternative injects a current derived from the specified power,

held fixed as the terminal voltage varies. Consumption then falls linearly with voltage rather than quadratically, which is closer to the behaviour of many aggregated loads.

Note what this is not. Because the current is computed from the nominal voltage and not from the present terminal voltage, this is a constant current model and not a constant power one. A true constant power load would require the current to be recomputed from the solved voltage at every step, making the component nonlinear and the nodal solve iterative. The linear model is used because it keeps the system matrix constant.

Which to use

The three canonical load characteristics are constant impedance, constant current and constant power, differing in whether demand follows the square of voltage, the voltage, or neither. Only the first two are available as linear models. For a voltage excursion of a few percent the choice changes little; for a deep depression during a fault it changes the answer materially, and the constant impedance model is the optimistic one because it sheds load exactly when the network is weakest.

Shunts

A shunt is specified directly as a conductance and a susceptance rather than as a power, so no conversion is involved. It is the natural representation for a capacitor bank or a reactor, where the rating is an admittance and the consumed power is a consequence of the voltage rather than the specification.

7 - Sources

A source imposes a quantity on the network. Which quantity it imposes, and whether it does so exactly, determines how it enters the system of equations and what it costs.

Current sources are free, voltage sources are not

A current source imposes a known current into a node. Its contribution is entirely on the right hand side of the nodal equations, and the system matrix does not know it exists.

A voltage source imposes a relation between two node voltages, which is not a nodal equation at all. Nodal analysis has one equation per node expressing current balance, and there is no current variable for an ideal voltage source to appear in. The system is extended with the source current as an unknown and with the constraint that fixes the voltage difference, as described under nodal analysis. The matrix grows by one row and column per source, and the added diagonal entry is zero, so the extended matrix is no longer positive definite and cannot be factorised by methods that assume it is.

This asymmetry is the reason so many models are formulated as current injections even when what they physically represent is a voltage behind an impedance.

The Norton equivalent

A voltage source with a series resistance can avoid the extension entirely. Source transformation replaces a voltage $V$ behind a resistance $R$ with a current $V/R$ in parallel with the same resistance,

which contributes a conductance to the matrix and a current to the right hand side. No extra unknown, no zero on the diagonal, and the matrix stays the shape it would have had without the source.

The two representations are equivalent at the terminals, exactly, for any $R$ that is not zero. The choice is therefore numerical rather than physical, and the cost is that the source is no longer ideal: its terminal voltage falls with the current drawn. Where a genuinely stiff source is wanted, $R$ has to be made small, and a small $R$ means a large conductance, which is the same conditioning trade-off that appears in switches.

Sources that change over time

The simplest time-varying source takes its value from a generator, as described under signal processing blocks.

A ramp source is more specific: it holds one value, then moves to a second over a defined interval, and holds that. The subtlety is what happens when the ramp changes not only the magnitude and phase but also the frequency. Interpolating a frequency linearly and applying it as if it had always been in force produces a phase discontinuity at both ends of the ramp, because phase is the integral of frequency and not its product with time. Blending the frequency contribution in and out smoothly over the ramp interval avoids that, at the price that the frequency during the transition is not the linear interpolation it appears to be.

A profile source takes its value from a recorded sequence instead of from a formula, stepping through samples as the simulation advances. It is the right choice when the excitation comes from a measurement, and it carries the obvious constraint that the sample rate and the simulation step must be reconciled: a profile is silent about what happens between its samples, and the simulation will ask.

Controlled sources

A controlled source takes its reference from another quantity in the simulation rather than from a parameter or a clock. This is what allows a component to be built out of sources: a converter imposes a voltage its control law computed, and an interface between two solvers imposes a value the other side produced.

The distinction from a time-varying source is that the reference is not known in advance. Since the reference is read as an input rather than solved simultaneously, it is the value from the previous step, which introduces a delay of one step into whatever loop the source closes. For a control loop that is usually acceptable and always worth knowing about; for a coupling between two solvers it is the central property of the method, and it is the subject of branches where the same delay is used deliberately.

8 - Network Injection and Compensation

Every simulated network stops somewhere. What lies beyond the boundary has to be represented by something, and once it is, the question of holding voltage at a bus becomes a question about what that representation will support.

The external network

The rest of the grid is represented as an ideal voltage source behind no impedance: a bus whose voltage is imposed and whose current is whatever the network draws. This is the slack of the powerflow carried into the time domain, and it supplies unlimited power at a fixed voltage and frequency.

That idealisation is the right one when the external system is genuinely much stiffer than what is being studied, and it is misleading when it is not. A stiff boundary suppresses exactly the behaviour that a weak grid study is about: it holds the voltage the compensator is supposed to be regulating and absorbs the power swings the converters are supposed to be sharing. Placing an impedance between the source and the network is what makes the boundary finite, and the short circuit ratio it produces is a modelling decision rather than a detail.

Because the imposed voltage comes from a signal generator rather than a constant, the boundary can also be driven: a frequency ramp to study the response to rate of change of frequency, or a modulated frequency to probe a control loop. The boundary then becomes the disturbance source rather than the reference.

Static reactive compensation

A static compensator regulates bus voltage by varying a shunt susceptance. It exchanges reactive power only, so it can raise or lower voltage but supplies no energy.

The regulator measures the bus voltage through a first-order lag, forms the per-unit error against a reference, and drives the susceptance through a further first-order lag with gain $K_r$ and time constant $T_r$,

both lags integrated with the trapezoidal rule so that the controller and the network advance consistently.

Two properties follow from the physics rather than from the controller. The susceptance is bounded at both ends by the installed capacitive and inductive ratings, and the regulator saturates against those bounds rather than failing; a compensator sitting on its limit is providing everything it has and the voltage error persists. And because the device is a susceptance rather than a source, the reactive power it delivers falls with the square of the voltage. It is weakest exactly when the voltage is lowest, which is when it is most needed.

The measurement lag matters more than it appears. It sits inside the regulator loop, so it is not merely a smoothing of the reported value; making it small to track faster couples the compensator to noise, and making it large delays the response into a range where it can interact with nearby machine controls.

Discrete compensation

Where the compensation is switched rather than continuous, the control is a different kind. The regulator compares the voltage error against a deadband, and only if the error exceeds it does it move by one discrete step, in the direction that reduces the error, subject to end stops.

The deadband is not a refinement but the central element. Without it any measurement noise drives continual switching, and the switching is mechanical and finite in life. With it, the steady-state voltage is not the reference but anywhere within a band around it, which is the accepted cost.

Power flow control

A device that connects two systems through a converter pair rather than through a magnetic circuit does not transfer voltage; it transfers power. Both sides are then specified as power exchanges rather than by a turns ratio: an active power to be moved from one side to the other, and a reactive power at each terminal set independently.

The distinction from a conventional transformer is that the sides are decoupled. The reactive power on one side is not a consequence of the other, the transfer does not depend on the angle across the device, and the two systems need not share a frequency. What is not free is the active power, which is common to both terminals up to losses; specifying it independently on each side would ask the device to create or destroy energy.

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

Reduced order machines from third to sixth order are derived under reduced order machine models. How any of these are stamped into the solver is covered under reduced order generator implementation.

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

9.1 - Reduced Order Machine Models

The full dq0 machine keeps every rotor winding as a state. Reduced order models keep the rotor flux linkages that matter on the timescale of interest and represent the rest algebraically, which removes the fastest states and lets the machine be stepped at the same rate as the network.

This page derives the equations. Nothing here depends on how they are arranged in software.

Per unit system and operational parameters

All quantities are in the reciprocal per unit system referred to the direct axis mutual inductance [Kundur1994]. With rated apparent power $S_n$ and rated line to line voltage $V_n$ as the stator base, the base peak phase voltage, base current and base impedance are

The machine is described by operational parameters rather than by winding data: the synchronous inductances $L_d$ and $L_q$, the transient inductances $L_d’$ and $L_q’$, the subtransient inductances $L_d’’$ and $L_q’’$, the corresponding open circuit time constants $T_{d0}’$, $T_{q0}’$, $T_{d0}’’$, $T_{q0}’’$, and the inertia constant $H$. A sixth parameter $T_{aa}$, the armature to field coupling time constant, distinguishes the two sixth order variants.

Which states each order retains

The state variables are the voltages behind the transient and subtransient reactances, $E_d’$, $E_q’$, $E_d’’$ and $E_q’’$, together with the two mechanical states. They are defined from the terminal quantities by

and analogously for the subtransient pair with $L_d’’$ and $L_q’’$.

Every order carries the same two mechanical states, so the third order model has five states in total and the sixth order models have eight.

Voltage behind reactance form

Written directly, the stator equations couple the machine currents to the network currents, and the machine inductances appear in the axis frame while the network is solved in the phase frame. The voltage behind reactance form removes that coupling: the machine is expressed as an internal voltage in series with a reactance that is constant in the axis frame, so the only quantity that changes between steps is the internal voltage.

The internal voltage is not a free variable. Applying the trapezoidal rule to the rotor flux equations over one step $\Delta t$ gives it as a recursion in quantities already known at the start of the step. For the transient states,

where $E_f$ is the field voltage supplied by the excitation system. The coefficients follow from the trapezoidal integration and depend only on the parameters and the step size:

The reactance differences are $Z_q’ = L_d - L_d’ - Y_d$ and $Z_d’ = L_q - L_q’ - Y_q$. For the orders without subtransient states $Y_d = Y_q = 0$ and $T_f = 0$. For the sixth order variant with armature coupling,

The subtransient states obey a recursion of the same shape with $T_{d0}’’$ and $T_{q0}’’$ in place of the transient time constants, and with the transient history entering as a forcing term.

Two properties of this form matter. The coefficients are computed once for a fixed step size, since they contain no state. And because $B_d’$ and $B_q’$ are the trapezoidal amplification factors $(2T - \Delta t)/(2T + \Delta t)$, they lie strictly inside the unit interval for any positive time constant, so the flux recursion is unconditionally stable regardless of step size.

Mechanical equations

The rotor obeys the swing equation in per unit,

with the electrical torque taken from the axis frame quantities as $T_e = V_d I_d + V_q I_q$. The load angle follows from the initial operating point as the phase of $V + j L_q I$, which is the standard construction of the q axis position from terminal conditions.

Solution schemes

The recursion above evaluates the internal voltage from quantities at the previous step, so the machine and the network are solved once per step in sequence. Two refinements exist for cases where that single pass is not accurate enough.

The predictor corrector method takes the single pass result as a prediction, re-evaluates the flux recursion using the corrected terminal quantities, and repeats until the change between successive passes falls below a tolerance. It converges to the solution of the implicit trapezoidal step rather than to its explicit approximation, at the cost of repeated network solutions.

The two stage predictor method splits the step differently: it advances the machine state on a predicted terminal voltage, then applies a single correction derived from the network solution, without iterating to convergence. It costs one extra network solve per step and removes most of the one step delay error.

Both are schemes for solving the same equations. They do not change the model order or the retained states.

Validity

The reduced order models assume the stator transients are fast enough to be neglected, so the stator is treated as algebraic. This is the assumption that makes the model valid for electromechanical studies and invalid where stator dynamics matter, such as during the first cycles of a close-in fault or for subsynchronous interaction. Neglecting the q axis rotor entirely, as the third order model does, additionally removes damping that is physically present, so a third order machine oscillates more than the same machine represented at fourth order.

References

• [Kundur1994] P. Kundur, Power System Stability and Control. New York: McGraw-Hill, 1994. Source of the reciprocal per unit system and of the operational parameter definitions used throughout this page.

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

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

10.1 - Exciters and Power System Stabiliser

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

10.2 - Turbines and Governors

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). $$

11 - Power Electronics

Every inverter model here is averaged: the switching is not represented, and the converter is treated as a controllable voltage behind its filter. Averaging removes the switching frequency from the problem, which is what allows a step size set by the control bandwidth rather than by the carrier. It also means these models say nothing about switching losses, harmonic injection or any behaviour that depends on the modulation itself.

The control that sits on top of each is derived separately under converter control, because the same cascade appears in more than one of these models.

Choosing among them

The models differ along two axes: which domain they are written in, and whether the converter follows the grid or forms it.

EMT Ph3 averaged VSI is the reference formulation. All fourteen states are real, the three filter phases are represented individually, and there is no carrier, so nothing is assumed about the bandwidth of what it carries.

DP Ph1 averaged VSI is the same converter as a single positive-sequence envelope. Its six real filter states become two complex envelopes, which is the saving the envelope description buys, at the cost of being unable to represent an unbalance.

DP Ph3 averaged VSI restores per-phase representation in the envelope domain, with one complex envelope per phase and a controller that keeps a single positive-sequence frame. Because three independent phase envelopes admit a negative-sequence component, it carries negative-sequence current control that the single-phase model has no need for.

EMT Ph3 grid-forming VSI is the one that differs in kind rather than in representation. It carries its own frequency and angle as states instead of tracking the grid’s, so it can energise a network with no other source. Its control is nonlinear enough that the model is linearized numerically at each operating point rather than written in closed form.

What they share

All four are solved simultaneously with the network rather than through a delayed injection, using the state-space nodal method described under SSN components. All four are therefore re-formed as the operating point moves, and all four make the system matrix change at every step, which is the cost of the approach.

11.1 - EMT Ph3 Averaged Voltage Source Inverter

This model represents a grid-following averaged voltage source inverter in the EMT domain. Because its state-space form is recomputed as the operating point moves, it is solved simultaneously with the network rather than through a delayed injection. 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.

How this is arranged in code, together with the source and the runnable examples, is covered under EMT Ph3 averaged VSI implementation.

11.2 - DP Ph1 Averaged Voltage Source Inverter

This 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 network equations.

How this is arranged in code, together with the source and the runnable examples, is covered under DP Ph1 averaged VSI implementation.

11.3 - DP Ph3 Averaged Voltage Source Inverter

This model extends the single-phase grid-following averaged inverter 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 of the single-phase model, 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 the single-phase case, 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 6 complex per-phase envelopes ahead of the 8 real control states, keeping the carrier-band and baseband blocks separate,

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, or 22 with the optional negative-sequence loop described below.

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 real and imaginary parts of the 6 complex per-phase envelopes and the 8 real control states packed into a single real 20-vector,

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

In this default configuration the controller operates in a single positive-sequence $dq$ frame, so 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.

Optional negative-sequence current control

A second, negative-sequence current-control loop can be added alongside the positive-sequence one, giving the dual-sequence structure of Yazdani and Iravani, chapter 8. The two configurations answer different questions: without the loop the model has the same 20 states and the same eigenvalue count as its EMT::Ph3 counterpart, which is what a cross-domain comparison requires, while with it the model gains 2 states and can regulate an unbalanced terminal.

The negative-sequence quantities are obtained by projecting the same three envelopes onto the conjugate sequence set,

A negative-sequence component rotates backwards relative to the PLL frame, so in envelope terms its $dq$ image follows from conjugating the projected phasor and rotating by $+\psi$ rather than $-\psi$,

The loop itself is the same PI structure as the positive-sequence inner loop,

reusing the inner-loop gains $K_{p,I}$ and $K_{i,I}$. Its output is redistributed to the per-phase bridge voltages through the sequence-orthogonal set, and adds to the positive-sequence command of the previous section,

The two references $i_{nd,\mathrm{ref}}$ and $i_{nq,\mathrm{ref}}$ default to zero, which makes the loop a negative-sequence suppressor. A non-zero pair commands a deliberate negative-sequence injection instead, as required by some unbalanced fault ride-through grid codes.

The state vector grows to 22 by appending the two integrators after the control block, so that the envelope and positive-sequence control indices are unaffected. Under a single-line-to-ground fault, enabling the loop suppresses the negative-sequence component of the injected current by about 40 percent while moving the positive-sequence component by less than 0.1 percent.

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

How this is arranged in code, together with the source and the runnable examples, is covered under DP Ph3 averaged VSI implementation.

11.4 - EMT Ph3 Grid-Forming Inverter

This 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}$ and $\mathbf{D}$ automatically. The affine offsets fix the model to the current operating point,

giving the affine state-space form

The dq/abc transformations and the nonlinear controls make this local model time varying, so it holds only in a neighbourhood of the operating point it was formed at.

How the linearization is carried out and stamped, together with the source and the runnable examples, is covered under EMT Ph3 grid-forming VSI implementation.

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.

12 - Converter Control

A converter model needs an angle to transform between the phase frame and its control frame, and a control law that decides what to synthesise. The two questions are separable, and the answer to the first is what distinguishes a grid-following converter from a grid-forming one.

Tracking an angle: the phase-locked loop

A phase-locked loop drives the estimated angle until the measured voltage sits on the chosen axis of the control frame. The error signal is the off-axis component, which is zero exactly when the frame is aligned, and it is fed to a proportional-integral controller whose output is a frequency correction.

With the integrator state $\phi$ and the frequency error input $e$, the loop is

The nominal frequency enters as a feed-forward term rather than being learned, so the loop only has to supply the deviation from it. That keeps the integrator near zero in normal operation and is why a loop initialised at nominal frequency locks quickly.

The proportional gain sets how fast the loop follows a phase step and the integral gain how fast it removes a standing frequency error. Making them large tracks disturbances the converter should arguably ignore: a phase-locked loop that follows a fault as fast as it can is not obviously desirable, since the converter then propagates the disturbance into its own control frame.

The important structural point is that a converter with a phase-locked loop takes its angle from the network. It cannot operate without a voltage to lock to, which is what “grid following” means.

Generating an angle: the oscillator

The alternative is to carry the angle as a state and advance it at a commanded frequency,

with no measurement involved. The converter then imposes a phase rather than following one, which is what “grid forming” means, and it continues to operate into a network with no other voltage source.

The difference between the two is one equation, but it determines whether the converter can start a de-energised network or support frequency, and whether it has any defined behaviour when the grid voltage collapses.

Cascaded control

Above the angle sits a cascade, ordered from slowest to fastest.

The outer loop compares measured active and reactive power against their references. The measurement is low-pass filtered first, because the instantaneous power computed from the terminal quantities carries components at twice the fundamental under any unbalance, and feeding those into a controller produces a modulation the converter should not emit. The filter cut-off therefore bounds how fast this loop can be, independently of its gains.

The inner loop regulates the filter current to the reference the outer loop produced. It must be substantially faster than the outer loop for the cascade to behave as intended: the outer loop is designed assuming its commanded current is achieved essentially immediately, and that assumption fails if the two bandwidths approach each other. The usual consequence is not instability but an interaction that appears as a poorly damped oscillation at neither loop’s design frequency.

Both loops are proportional-integral in the control frame, where a balanced fundamental quantity is constant, so an integrator can drive the steady-state error to zero. This is the reason for working in a rotating frame at all: the same controller applied to a sinusoid in the phase frame would leave a standing error, because an integrator cannot track a moving target.

Grid-forming voltage control

A grid-forming converter replaces the outer power loop with a voltage magnitude and frequency law. Droop characteristics relate active power to frequency and reactive power to voltage, which lets several converters share load without communicating: each responds to the same measured deviation, and the split follows from the droop gains.

Below that, a voltage loop regulates the filter capacitor voltage and hands a current reference to the same inner current loop as before. The inner loop is therefore common to both control philosophies; only what sits above it changes.

13 - State-Space Nodal Components

The state-space nodal method solves a component simultaneously with the network instead of coupling it through a delayed injection. This page covers the other half: how a single component is written so that the method applies to it, and what the resulting stamp is.

What a component supplies

Each component provides a continuous-time model in a state, an input and an output of its own choosing,

The choice of what the vectors mean is the entire modelling step. For a component that behaves as an admittance the input is the terminal voltage and the output is the terminal current; for one that behaves as an impedance the roles are exchanged. Everything after this is mechanical.

Discretisation

The state equation is integrated with the trapezoidal rule, which gives the discrete pair

Substituting the discrete state into the output equation separates the output into a part that depends on the present input and a part that does not,

This is exactly a companion model. $\boldsymbol{W}$ is an equivalent admittance that goes into the system matrix and $\boldsymbol{y}_{hist}$ is an equivalent source that goes into the right hand side. The difference from element-by-element companion models is only that $\boldsymbol{W}$ is derived from the component’s own state-space description rather than written by hand, so a component with internal states and cross-coupling between phases needs no special treatment.

The simplest case reproduces the classical result

Take a three-phase inductor. The natural choice is the current as state, the voltage as input and the current as output, giving

With $\boldsymbol{A} = \boldsymbol{0}$ the discretisation collapses to $\boldsymbol{A}_d = \boldsymbol{I}$ and $\boldsymbol{B}_d = \tfrac{\Delta t}{2}\boldsymbol{L}^{-1}$, so the equivalent admittance is $\tfrac{\Delta t}{2}\boldsymbol{L}^{-1}$ and the history term is the previous current plus the previous voltage contribution. For a single phase that is $\Delta t / 2L$, the familiar trapezoidal companion model of an inductor.

This is worth doing once because it shows the machinery adds no approximation of its own. A component whose model is a plain inductor gets exactly the stamp it would have had.

Where it earns its cost

The method is worth using when the component cannot be decomposed into independent elements. A series RLC branch written as three separate companion models introduces two internal nodes; written as one state-space model it introduces none, and the resulting stamp is a full matrix that captures the coupling directly. The same applies to any component whose phases are coupled through a non-diagonal inductance or through a control law.

The cost is that $\boldsymbol{W}$ is dense over the component’s terminals, where element models produce sparse contributions, and that a matrix inverse of the size of the state vector is required whenever the model changes.

The derivation above is the instantaneous case. The same component model is discretised differently in an envelope domain; see SSN across domains.

Fixed and varying models

If $\boldsymbol{A}$, $\boldsymbol{B}$, $\boldsymbol{C}$ and $\boldsymbol{D}$ are constant, the discrete matrices are computed once and the system matrix never changes on account of the component.

If the model depends on the operating point, it must be re-formed and re-discretised as the operating point moves, and the system matrix refactorised with it. A saturating inductor whose inductance is a piecewise linear function of flux is the simple case; a converter whose control law is nonlinear is the general one.

Initialization

The steady state at a given frequency follows from the continuous model directly,

which is the state-space equivalent of evaluating a phasor impedance. A component initialized this way starts in steady state rather than settling into it, provided the model is linear at the operating point.

14 - Signal Processing Blocks

Alongside the network components, a simulation contains blocks that carry no current and connect to no node. They compute a value from another value, and they exist because controllers and sources are built out of them.

Integration

An integrator advances a state from its input using the trapezoidal rule,

which is the same rule the network solver applies to reactive elements, so a control loop and the circuit it acts on are integrated consistently. Using a cruder rule for the controller would put an error into the loop that no amount of tuning removes.

Not every block needs that accuracy. An angle accumulator advancing at a commanded frequency is often stepped with the rectangular rule instead,

which is a step behind but adds no dependence on the previous input. The distinction is worth knowing when comparing an angle against one produced elsewhere, because the two rules differ by half a step of phase.

Finite impulse response filtering

A finite impulse response filter forms its output as a weighted sum of the most recent inputs,

holding those inputs in a circular buffer of length $N$. Because the output depends only on past inputs and never on past outputs, the filter cannot become unstable whatever the coefficients are, and its phase response can be made exactly linear. The price is that a given sharpness needs a long filter, which costs both memory and delay.

The delay is the part that matters in a control loop: a filter of length $N$ contributes roughly $N/2$ steps of it. Inside a feedback path that delay is a phase lag, and it erodes stability margin just as surely as raising a gain would.

Signal generators

A source that varies over time takes its value from a generator. Four behaviours cover most uses: a constant, a sinusoid at a fixed frequency and amplitude, a sinusoid whose frequency ramps between two values, and one whose frequency is modulated continuously.

The frequency ramp is the one with a subtlety. A ramp is described by a start frequency, an end frequency and a rate of change, and it is tempting to generate it by evaluating $\sin(\omega(t),t)$ with a time-varying $\omega$. That is wrong: the argument of the sine is the accumulated phase, not the product of the present frequency and the elapsed time, and the two differ whenever the frequency is not constant. The phase must be accumulated,

so that the instantaneous frequency is the derivative of the phase by construction. Generating a ramp the naive way produces a signal whose actual frequency sweeps at twice the intended rate.

Accumulating phase makes the result depend on the step size and on the history of the run. Where an exactly reproducible waveform is wanted, independent of when the simulation started or what steps it took, the phase can instead be computed in closed form from the ramp parameters, which for a linear ramp is a quadratic in time.