Chapter 17 — Modeling control systems

Modules 1 through 4 built the plant: resistors, inductors and capacitors as companion models, lines as travelling waves, transformers and machines as coupled magnetic circuits. Every one of them was assembled into a single matrix and solved simultaneously. This chapter is about the other half of a power system — the control domain, the block diagrams that sense, compute and command — and it is solved on a completely different principle. There is no matrix, no simultaneous solution, and no iteration. There is an order, and there is a delay, and getting both right is the whole job.

Chapter 4 introduced the control palette and traced one loop from a sensed power to an IGBT gate. Here we open the machinery: how a continuous block becomes a difference equation, how the engine decides what runs before what, and exactly what it costs — in time steps — for a signal to cross between the two domains. Chapters 18 and 19 then build real generator controllers on top of it.

Learning objectives

By the end of this chapter you should be able to:

• Explain why the control domain is solved as an ordered sequence of assignments rather than as a matrix, and what that buys.
• Write the discrete update NumaSim actually uses for a first-order lag, and say how it differs from the trapezoidal companion models of Chapter 2.
• State the accuracy and stability consequences of that choice, and compute the time step at which a given block goes unstable.
• Explain what a topological sort is doing in a simulator, which edges it deliberately ignores, and what happens when the graph has a cycle.
• Predict, in time steps, the lag a signal accumulates when it crosses between the control and electrical domains — in either direction, and around a loop.
• Choose between a wire, a label and a Go-To/From pair, and know which signals you can put on a scope.

17.1 The control domain is a second solver

An electrical network is a system of simultaneous equations. Kirchhoff's current law at one node involves the voltages at neighbouring nodes, those nodes' equations involve their neighbours, and there is no way to compute one voltage before another. That is why Chapter 3 assembled everything into $\mathbf{G}\mathbf{v}=\mathbf{i}$ and factorized it: the coupling is mutual, so the solution has to be simultaneous.

A control diagram has no such coupling. A gain block's output depends on its input and nothing else. An integrator's output depends on its input and its own stored state. Every block is a plain assignment

$$y_n=f({\text{inputs}}_n,{\text{state}}_{n-1}),$$

and the state is private to the block. Nothing in a control diagram requires two blocks to be solved together, so nothing needs a matrix. The engine walks a list of blocks, evaluates each one, and moves on.

This is not a simplification or an approximation — it is a structural property of signal-flow diagrams, and it is why control blocks are almost free compared with network nodes. A hundred extra control blocks add a hundred cheap function evaluations per step. A hundred extra electrical nodes change the size of the matrix.

Two consequences follow immediately, and they are the subject of the next three sections. First, if each block is an assignment, then somebody has to decide what order to evaluate them in — get it wrong and a block reads a value that has not been computed yet. Second, the control domain and the electrical domain are now two different solvers running side by side, and something has to define how they exchange values across the boundary.

17.2 Discretizing a control block

A control block specified in the frequency domain — a lag $K/(1+sT)$, an integrator $1/s$, a washout $sT/(1+sT)$ — is a differential equation, and the engine has to turn it into a recurrence. Chapter 2 did exactly this for the inductor and capacitor using the trapezoidal rule, which averages the derivative over the step.

NumaSim's control blocks do not use the trapezoidal rule. They use forward Euler. This is a real and deliberate difference between the two domains, and you need to know it to predict what a control block will do.

The first-order lag

Take the first-order lag $H(s)=K/(1+sT)$, i.e. the differential equation

$$T\dot{y}+y=Ku.$$

Forward Euler evaluates the derivative using values already in hand and steps forward:

$$\frac{y_n-y_{n-1}}{\mathrm{\Delta }t}=\frac{K{u}_n-y_{n-1}}{T},$$

which rearranges to the update the engine actually performs:

$$\boxed{{\displaystyle y_n=y_{n-1}+\frac{\mathrm{\Delta }t}{T}(K{u}_n-y_{n-1}).}}$$

There is one multiply and one add. Notice what is not there: no division by a step-dependent denominator, no averaging of $u_n$ with $u_{n-1}$, no history term to maintain separately. The block's entire state is the single number $y_{n-1}$.

You can verify this against the engine without running anything. NumaSim keeps golden reference outputs for its blocks; the one for the first-order transfer function drives a unit step into $K=1$, $T=0.1$ s at $\mathrm{\Delta }t=1$ ms, so $\mathrm{\Delta }t/T=0.01$, and the recurrence predicts

$$y_1=0.01,y_2=0.01+0.01(1-0.01)=0.0199,$$$$y_3=0.0199+0.01(1-0.0199)=0.029701.$$

Those are exactly the first three rows of the stored reference. Had the block been trapezoidal, the same step would have produced $0.00995$, $0.019801$, $0.029554$ — close, but not the same numbers, and the discrepancy is a useful way to tell the two schemes apart when you are checking someone's model against a reference.

Why the two domains differ

The trapezoidal rule earns its cost in the network because the network is solved as a matrix anyway: the companion model's history source folds into the right-hand side at no extra cost, and second-order accuracy comes free. A control block has no matrix to hide in, and forward Euler is the cheapest thing that works — one multiply-add, no per-block bookkeeping. Given that control time constants are typically milliseconds to seconds while EMT time steps are microseconds, the accuracy the explicit method gives up is usually far below the uncertainty in the parameters themselves.

What forward Euler costs you

Two properties matter in practice.

It is first-order accurate. The local error per step is $O(\mathrm{\Delta }{t}^2)$, against $O(\mathrm{\Delta }{t}^3)$ for trapezoidal, so halving the step roughly halves the accumulated error instead of quartering it. For a lag with $T=0.1$ s stepped at 1 ms, the response is visibly but slightly fast; at the 20 µs steps typical of an EMT run, the error is invisible.

It has a stability limit. Write the update with $u=0$:

$$y_n=(1-\frac{\mathrm{\Delta }t}{T})y_{n-1}.$$

The block decays only if $|1-\mathrm{\Delta }t/T|<1$, which requires

$$\boxed{{\displaystyle \mathrm{\Delta }t<2T.}}$$

Step at exactly $2T$ and the block oscillates forever without decaying; step beyond it and the output grows without bound. The trapezoidal rule has no such limit — it is unconditionally stable for this equation — so this really is a constraint the control domain imposes and the network does not.

In an EMT study the margin is usually enormous. A voltage transducer with $T_R=20$ ms simulated at $\mathrm{\Delta }t=20$ µs has a factor of two thousand in hand. The bound bites in exactly one situation: a time constant you set smaller than your time step. If you type $T=5$ µs into a block and run at 20 µs, that block is unstable, and the symptom is an output that alternates sign and doubles every step until the run fills with infinities. The fix is never to reduce the time step to chase it — it is to ask whether a 5 µs time constant means anything in a study whose resolution is 20 µs. Usually the honest answer is to set it to zero.

Zero means "not there"

Which brings us to a convention every model in the library relies on: a zero time constant removes its block rather than dividing by zero. A lag with $T=0$ becomes a pass-through, $y_n=K{u}_n$. A lead-lag $(1+s{T}_C)/(1+s{T}_B)$ with $T_B=0$ falls back to a discrete lead, and with both zero becomes a wire. A gain of zero disables its branch entirely.

This matters because the standard controller models in Chapters 18 and 19 are published with tables of parameters in which zeros are common and meaningful — a manufacturer who does not use transient gain reduction publishes $T_B=T_C=0$, and the model is expected to simply not have that stage. You will see this convention in the datasheet defaults throughout Module 5.

The rest of the palette

Not every block is a lag. The engine uses whichever discretization suits the block, and it is worth knowing which is which:

Block family	Discretization	Examples
Lags, integrators, washouts	Forward Euler	integrator, tf1, the internal stages of every controller in Chapters 18-19
Derivatives, degenerate leads	Backward difference $(u_n-u_{n-1})/\mathrm{\Delta }t$	the D term of pid, a lead-lag with a zero denominator
Pure transport delay	Integer-sample ring buffer	the dead time of tf1, the engine delay of a gas turbine governor
One-step shift	The memory arena itself, no arithmetic	unit_delay
Sliding window	Ring buffer over one cycle	rms
Algebraic, no state	Evaluated directly at the current sample	gain, the Park transform abc_dq0, the PWM comparator pwm

Two of these deserve a note. A transport delay is rounded to a whole number of samples, $N=\mathrm{round}(L/\mathrm{\Delta }t)$, so the realized dead time is the requested one snapped to the time-step grid — request 350 µs at a 100 µs step and you get 400 µs. And the integrator's saturation limits clamp the state, not just the output, which is a deliberate anti-windup: once the state is pinned at a limit, further input in the same direction cannot accumulate a debt that has to be paid back before the output moves.

17.3 Evaluation order is a topological sort

If every block is an assignment, order is everything. Evaluate a summing junction before the gain that feeds it and the junction adds last step's value instead of this step's. The result is not a crash; it is a subtly wrong answer, which is worse.

Read the control diagram as a directed graph: one vertex per block, one edge from each block that produces a signal to each block that consumes it. "Every producer runs before all of its consumers" is precisely a topological sort of that graph. NumaSim computes one with Kahn's in-degree algorithm: count each block's unsatisfied inputs, repeatedly emit any block whose count has reached zero, and decrement its consumers. The cost is linear in blocks plus wires.

The sort runs once, at compile time. It happens when the circuit is built, before $t=0$; every step then walks the same fixed schedule. This is the control-domain twin of factorizing $\mathbf{G}$ once — do the structural work at setup so the inner loop is pure arithmetic. A run of ten million steps sorts the graph exactly one time.

Three kinds of edge the sort ignores

A naive reading would build the graph from every connection in the circuit. The engine deliberately does not, and the exclusions are the interesting part.

Electrical-to-electrical connections are not edges. Two resistors sharing a node do not have a producer/consumer relationship — they are simultaneous, and the matrix solve handles their coupling. Feeding those connections to a topological sort would be meaningless (and would usually produce a cycle, since electrical connections are undirected).

Cross-domain connections are not edges either. This is the surprising one. When a control block commands an electrical device, you might expect the engine to force the control block to run first. It does not. There is no ordering constraint between the domains at all, which is what allows them to be evaluated independently — and, on a multi-core run, genuinely in parallel. Section 17.4 explains what makes that safe.

Edges into a unit delay are not edges. A unit_delay is defined to output its input from the previous step, so it does not care whether its source has run yet. Dropping its incoming edges is what lets it break a feedback loop.

Underneath the sort there is a second, coarser mechanism: each executable carries an execution phase, with control blocks at phase 0, the electrical pre-solve (which stamps switch conductances and injections) at 10, the network solve at 40, and the electrical post-solve (which reads node voltages back out) at 50. These are a tie-break, not a dependency: when several blocks are ready to run, the lowest phase goes first. The phases keep the electrical sub-loop in a sensible internal order — you cannot read node voltages before solving for them — but they impose nothing between control and electrical.

Feedback needs a delay somewhere

A topological order exists only for an acyclic graph, and a feedback loop is a cycle by definition. Two cases, with different outcomes.

The loop closes through the network. A controller senses a voltage, computes, and commands a device; the device changes the voltage. On the graph this is not a cycle at all, because neither cross-domain connection is an edge. The loop sorts cleanly with no work from you, and the delay that makes it physically sensible is supplied by the mechanism of §17.4.

The loop stays inside the control domain. A block's output reaches its own input through control blocks only. That is a genuine cycle — the classic algebraic loop, in which a block would need its own output in order to compute its own output. Break it with an explicit unit_delay on the feedback path. Its output is its input from the previous step, so it imposes no ordering constraint, and you choose where the delay lands.

Leave the cycle unbroken and the run still proceeds. When Kahn's algorithm stalls with blocks remaining, the engine picks the one with the fewest unsatisfied inputs, forces it into the schedule, logs a warning to the log panel, and carries on. The forced block reads a stale value for its unsatisfied input, which is to say the engine has silently inserted a one-step delay — on an edge chosen by a tie-break rather than by you. The numbers you get are plausible and the run completes, which is exactly what makes it worth heeding the warning. If a loop needs a delay, put it where you want it.

17.4 The delay between control and electrical

Here is the mechanism that makes everything above hold together, and it is worth stating precisely because it is the most common source of "why is my controller a step behind?"

Every signal value in a simulation lives in a slot with two cells: a current cell and a prev cell. Producers always write current. At the start of every step the engine copies the whole current plane into the prev plane in one bulk memory copy, so for the entire duration of a step, prev holds a frozen, consistent snapshot of the previous step while current is being filled in.

When the circuit is compiled, each connection is classified once — never again at run time — into which of the two cells the consumer will read:

• Same-domain connections (control to control, or electrical to electrical) read current. The topological order guarantees the producer already ran this step, so the value is fresh and the block behaves exactly as its transfer function says.
• Cross-domain connections read prev. Both directions. Control to electrical, and electrical to control.

That last rule is the whole story. It is why the topological sort can ignore cross-domain edges: a consumer reading prev has no dependency on the producer's execution this step, because it is not reading anything the producer writes this step. The two domains genuinely cannot race, whatever order they run in.

Counting the steps

The arithmetic follows directly. A controller commanding a device sees its command take effect one step later. A controller sensing the network sees measurements that are one step old. And a loop that leaves the control domain and comes back — the normal case for any real regulator — accumulates two steps of lag.

Where the two steps come from. Each step begins by copying the current plane into prev. Same-domain reads take current and are fresh; cross-domain reads take prev. A measurement written by the network in step $k-1$ is read by the controller in step $k$, and the command the controller writes in step $k$ is not seen by the network until step $k+1$.

You can see all three numbers on a scope. The control loop delay sample is built for nothing else: a square-wave command drives a breaker, and the resulting load current is brought back into the control domain. At a 1 ms step the command falls at $t=0.050$ s, the breaker state and load current follow at $t=0.051$ s, and the control-domain copy of that current follows at $t=0.052$ s. One step out, one step back, two steps around.

Is this a defect?

No — it is the honest representation of a system in which the controller and the plant are separate machines. A real regulator samples a transducer, computes for a finite time, and drives an actuator; it is never in instantaneous communication with the flux in an iron core. A simulator that pretended otherwise would have to iterate the two domains to convergence within every step, at a cost far out of proportion to the physics.

What matters is keeping the lag small compared with the dynamics you are studying. Two steps at 20 µs is 40 µs against an exciter time constant of half a second — four orders of magnitude, and utterly negligible. Two steps at 1 ms against a switching converter running at 10 kHz is 2 % of a switching period, which is starting to matter. The rule of thumb is the same one that governs every other choice of $\mathrm{\Delta }t$: make the step small enough that the fastest loop you care about sees many steps per cycle, and the cross-domain delay disappears along with every other discretization artifact.

Three related cases

The same prev rule shows up in three other places, and they are worth recognizing so they do not surprise you:

• Portless named inputs — signals a component picks up by name rather than through a drawn wire — always read prev, regardless of domain. They are always one step delayed.
• Across a rate boundary. In a multi-rate subsystem, a signal read from a different rate group is forced to prev even within the control domain, giving a one-substep delay. See multi-rate monitoring.
• A unit delay feeding across the domain boundary costs two steps, not one: the unit_delay contributes its own step, and the domain crossing contributes another.

17.5 Naming, observables and wire-less routing

The last piece of the control domain is bookkeeping, and it is what makes a model readable and measurable.

Any control signal you give a name becomes an observable — something you can put on a scope, plot, or export. There are three ways to name one, and they are equivalent as far as the engine is concerned:

• Label a wire. Give the wire between two blocks a signal name and the value travelling along it is recorded.
• Name a component's monitored output. Most components have measure_… toggles with an accompanying name field — a resistor's current, a breaker's state, a controller's field-voltage output. Turning one on publishes that internal quantity under the name you choose.
• Publish with a Go-To. A goto block takes a signal in and publishes it under a name, with no wire leaving it.

Whatever publishes a name, a from block reads it back somewhere else on the sheet, with no wire drawn between them. This is not cosmetic. A generator, its exciter, its governor and its stabilizer exchange half a dozen signals, and drawing all of them as wires across a one-line diagram makes the diagram unreadable. Naming lets you keep the power diagram looking like a power diagram and the control diagram looking like a control diagram.

One caution follows from §17.4: a Go-To/From pair inside the control domain is a plain alias and costs nothing, but a From reading a signal published by an electrical component is a domain crossing and costs its step, exactly as a wire would. Naming changes how the sheet looks, never what the engine computes.

17.6 Lab: counting the steps

Lab 17A — measure the cross-domain delay

Open the control loop delay demo in simulator →

(Full description: control loop delay.)

The circuit is a square-wave command driving a breaker that feeds a 10 Ω load from a 10 V source, with the load current brought back into the control domain through a From tag and a unity gain. The run uses a deliberately coarse 1 ms step so that one step is one visible millisecond.

Run it and plot cmd, brk, i_load and i_ctl on the same axes. Zoom in on the edge at $t=0.050$ s and confirm the ladder:

Signal	Domain	Edge time	Steps behind cmd
cmd	control	0.050 s	0
brk, i_load	electrical	0.051 s	1
i_ctl	control	0.052 s	2

Then change the time step to 100 µs and re-run. Every gap shrinks by a factor of ten, because the delay is two steps and not a fixed span of time. That is the single most useful thing to internalize about it.

Lab 17B — find the loop in a working controller

Open the PWM power controller in simulator →

(Full description: PWM power control.)

This is the loop from Chapter 4, revisited with the machinery of this chapter. Trace it and answer three questions on the canvas. Which connections are same-domain and therefore exact within the step? Which two are the domain crossings that give the loop its lag? And why does this circuit sort without needing a unit_delay anywhere, even though it is unambiguously a feedback loop? (Answers in §17.8.)

17.7 Summary

• The control domain is a second solver with a different principle: every block is a plain assignment $y_n=f(\text{inputs},\text{state})$, so there is no matrix and no simultaneous solution — only an order.
• Continuous control blocks are discretized with forward Euler, not the trapezoidal rule used for the network's companion models. The first-order lag is $y_n=y_{n-1}+(\mathrm{\Delta }t/T)(K{u}_n-y_{n-1})$.
• That choice is first-order accurate and conditionally stable, requiring $\mathrm{\Delta }t<2T$. The bound is irrelevant at EMT time steps unless you set a time constant smaller than the step, in which case the block diverges.
• A zero time constant removes its block rather than dividing by zero — the convention that lets published controller parameter sets contain meaningful zeros.
• Blocks are evaluated in a topological order computed once at compile time. The sort ignores electrical-to-electrical connections, cross-domain connections, and edges into a unit delay.
• An algebraic loop inside the control domain needs an explicit unit_delay. Without one the engine breaks the cycle itself, warns, and puts the delay wherever the tie-break lands.
• Every cross-domain signal reads the previous step's value, in both directions. One step out, one step back, two steps around a loop that closes through the network.
• Naming a signal — by wire label, by a component's monitor, or with a Go-To — makes it an observable, and lets a From read it anywhere on the sheet without a drawn wire.

17.8 Problems

Problem 17.1. A first-order lag with $K=2$ and $T=50$ ms is simulated at $\mathrm{\Delta }t=5$ ms with $y_0=0$ and a unit step applied from the first step. Compute $y_1$, $y_2$ and $y_3$, and state the final value the block converges to.

Solution 17.1

Here $\mathrm{\Delta }t/T=0.1$, so $y_n=y_{n-1}+0.1(2-y_{n-1})$.

$y_1=0+0.1(2-0)=0.2$. $y_2=0.2+0.1(2-0.2)=0.38$. $y_3=0.38+0.1(2-0.38)=0.542$.

The steady state is where $y_n=y_{n-1}$, which requires $Ku-y=0$, so $y_{\mathrm{\infty }}=Ku=2$. The DC gain is exact regardless of step size — forward Euler gets the final value right and only the transient's shape wrong.

Problem 17.2. A colleague sets a transducer time constant to 10 µs in a model running at $\mathrm{\Delta }t=50$ µs, and the run fills with infinities within a few milliseconds. Explain what happened and give two fixes, saying which one you would use.

Solution 17.2

Forward Euler on a lag is stable only for $\mathrm{\Delta }t<2T$. Here $2T=20$ µs and $\mathrm{\Delta }t=50$ µs, so the homogeneous factor is $1-\mathrm{\Delta }t/T=1-5=-4$: the output alternates sign and grows by a factor of four every step, reaching overflow in about a dozen steps.

Fix one: reduce the time step below 20 µs. Fix two: set the time constant to zero, which by the zero-parameter convention makes the block a pass-through.

Use the second. A 10 µs time constant carries no information in a study resolved at 50 µs — it is far faster than anything else being modelled — so representing it as instantaneous is both cheaper and more honest than dragging the whole simulation down to resolve a dynamic nobody is studying.

Problem 17.3. A sheet contains control blocks $A$, $B$, $C$ and an electrical device $E$. The connections are $A\to B$, $B\to E$ (a gate command), $E\to C$ (a measured current) and $C\to A$. Draw the dependency graph the engine actually builds, give a valid evaluation order, and say how many steps of delay the loop carries. Then repeat for the same four blocks with $E$ replaced by a control-domain gain.

Solution 17.3

With $E$ electrical, the two cross-domain connections $B\to E$ and $E\to C$ are not edges: their consumers read the previous step's value, so there is no execution dependency. The graph the engine builds is just $A\to B$ and $C\to A$ — a chain, not a cycle. A valid order is $C,A,B$, with $E$ scheduled independently in the electrical phases. The loop carries two steps of delay, one per crossing.

With $E$ replaced by a gain, every connection is same-domain and every one is an edge. The graph is $A\to B\to E\to C\to A$: a four-node cycle with no topological order at all. This is an algebraic loop, and it needs an explicit unit_delay on whichever edge you want the delay to land — after which the loop carries exactly one step, not two.

The pair is worth comparing. Replacing a physical device with a "mathematically equivalent" control block does not preserve the timing: it changes the loop from two steps of delay to one, and turns a graph that sorts silently into one that warns.

Problem 17.4. In Lab 17B, the PWM power controller is a closed feedback loop, yet it needs no unit_delay and produces no cycle warning. Why? Contrast this with a loop built entirely from gain and add2 blocks.

Solution 17.4

The PWM loop leaves the control domain (the gate signal commands an IGBT) and re-enters it (the sensed load voltage and current come back). Neither crossing is an edge in the dependency graph, because cross-domain consumers read the previous step's value and therefore have no execution dependency on the producer. With both crossings removed the remaining graph is a simple chain, which sorts trivially. The physical delay the loop needs is supplied by the prev reads — two steps around.

A loop of pure gain and add2 blocks never leaves the control domain, so every edge is a same-domain edge reading current, every edge constrains the order, and the cycle is real. That is an algebraic loop and it needs an explicit unit_delay.

Problem 17.5. A controller reads a machine's terminal voltage, computes a field voltage, and commands the machine, all at $\mathrm{\Delta }t=20$ µs. The regulator's dominant time constant is 0.4 s. Quantify the round-trip lag as a fraction of that time constant, and comment. Then repeat for a switching converter loop at $\mathrm{\Delta }t=1$ ms with a 10 kHz switching frequency.

Solution 17.5

Generator case: the round trip is two steps, 40 µs, against a 0.4 s time constant — a ratio of ${10}^{-4}$. The lag is four orders of magnitude faster than the dynamics being studied and cannot affect the result in any way you could measure.

Converter case: two steps is 2 ms... except that a 10 kHz switching period is only 100 µs, so a 1 ms step cannot resolve the switching at all — the time step is the primary problem, not the delay. Choosing a step that does resolve switching, say 1 µs, makes the round trip 2 µs, or 2 % of a switching period, which is small but no longer negligible for a tightly-tuned inner current loop.

The general lesson: the cross-domain delay is never the binding constraint on $\mathrm{\Delta }t$. Whatever step resolves your fastest dynamic already makes two steps of lag small.

Problem 17.6. A tf1 block is given a dead time of 350 µs in a run at $\mathrm{\Delta }t=100$ µs. What dead time does the model actually realize, and what would you change to get closer to 350 µs?

Solution 17.6

Transport delay is stored in a ring buffer holding a whole number of samples, $N=\mathrm{round}(L/\mathrm{\Delta }t)=\mathrm{round}(3.5)=4$, so the realized dead time is 400 µs — about 14 % long.

To get closer, reduce the time step so the requested delay lands nearer a sample boundary: at 50 µs, $N=7$ and the realized delay is exactly 350 µs. There is no sub-sample interpolation, so the only lever is the step size.

17.9 References

• P. Kundur, Power System Stability and Control, McGraw-Hill — representation of control blocks and transfer functions in stability programs.
• J. Arrillaga and N. R. Watson, Power Systems Electromagnetic Transients Simulation, IET Power and Energy Series 39 — control-system representation and the interface between control and network solutions.
• H. W. Dommel, Electromagnetic Transients Program (EMTP) Theory Book, Bonneville Power Administration — transient analysis of control systems (TACS) and its coupling to the network solution.
• A. B. Kahn, "Topological sorting of large networks," Communications of the ACM, vol. 5, no. 11, pp. 558-562, 1962 — the in-degree algorithm used to order the control blocks.
• J. C. Butcher, Numerical Methods for Ordinary Differential Equations, Wiley — accuracy and stability regions of explicit and implicit one-step methods.

---

Previous: Chapter 16 — The induction machine · Next: Chapter 18 — Excitation systems and AVRs.