Controlling Power Systems with Renewables — A Modular Approach (Part II)

In this post we will introduce a structured model of an electrical power system that is well suited for modular control system design. This is the second post in the series on the modular control of electrical power systems — please see Part I for an overview and introduction to the problem if you didn’t already.

The focus today is simply on deriving the model and understanding its structure. Future posts will deal with how to exploit this structure to conduct modular design. Unlike before, to fully understand everything, some starting knowledge will be required:

• The use of ODEs to describe dynamical systems.
• Kirchhoff’s laws and some basic knowledge of circuits and AC signals.
• Some familiarity with matrices.

Nothing beyond the content of a basic course in control theory though, and hopefully most will be understandable with considerably less.

Model overview

Conceptually a power system can be seen as the interconnection of components through a transmission network.

In this post we will develop a structured mathematical model of this system. In particular we will construct a state-space model of the form

\[\dot{x}(t)=\mathbf{A}x(t)+\mathbf{B}P_d(t),\;x(0)\in\mathbb{R}^m.\]

This model describes how the state of the power system changes over time as a function of the operating point, the system parameters, the initial state, and an external power disturbance. These features correspond to the following parts of the model:

• The state of the model at time \(t\) is described by the \(m\)-dimensional vector \(x(t)\).
• The matrices \(\mathbf{A}\) and \(\mathbf{B}\) are determined from models of the components and the transmission network, and depend on the system parameters and operating point.
• The \(m\) dimensional vector \(x(0)\) describes the state of the power system at time \(t=0\) (the initial state).
• The \(n\) dimensional vector \(P_d(t)\) describes an external power disturbance at time \(t\).

The key feature of the model the we will obtain is that the matrices \(\mathbf{A}\) and \(\mathbf{B}\) are highly structured, and can be constructed in a modular manner. Ultimately it is this feature that will allow us to develop modular control design methods for power systems in future posts. A few comments are in order here:

• The important structural property lies in how the matrices are formed. This is distinct from concepts such as sparsity. In fact, in almost every practical case the \(\mathbf{A}\) matrix will not be sparse (nor should it be).
• The actual entries of the matrices are not important. As we will see in future posts, modular design can actually be conducted without ever explicitly finding even a single entry of these matrices; all the required modelling details can, for example, be inferred from open loop experiments on the components instead (though models should of course be used if they are available).

Now, on with the construction of \(\mathbf{A}\) and \(\mathbf{B}\).

A structured block diagram

In order to turn this conceptual picture into a structured dynamical model it is helpful to introduce the following block diagram:

Each arrow in the block diagram is associated with a signal, and each block the dynamical model of either a component or the transmission network.

The signals describe the voltages and power flows at every point in the network where a component can be connected. They consist of the stacking up of the signals at the individual points, for example:

\[P_{N}(t)=\begin{bmatrix}P_{N,1}(t)\\\vdots{}\\P_{N,n}(t)\end{bmatrix}.\]

The individual entries of these signals have the following meanings:

• \(P_{N,i}(t)\) — the amount of electrical power being transmitted to the \(i\)-th component by the transmission network at time \(t\).
• \(P_{P,i}(t)\) — the amount of electrical power being consumed by the \(i\)-th component at time \(t\). This quantity is negative if power is being produced.
• \(P_{d,i}(t)\) — an electrical power disturbance that can be used to describe unmodelled or uncertain power consumption by the component at time \(t\).
• \(V_i(t)\) — the AC voltage at the terminal of the \(i\)-th component at time \(t\).

There is a little ambiguity to here. For example, we haven’t said exactly what is meant by an AC voltage. We will specify the things more precisely later when we start introducing the dynamical models.

The blocks in the diagram denote dynamical systems. We have \(n+1\) blocks, one for each component, and one for the network. The presence of a block means that there is a dynamical relationship between the signals that flow in and out that depends on the the thing in the block. For the case of the network block, this means that the dynamics of the transmission lines enforces a relationship between the power flows between the components and the three-phase voltages at the terminals of the components. Similarly, the dynamics of the \(i\)-th component enforces relationships between \(P_{P,i}(t)\) and \(V_i(t)\).

The topology: it is clear that the dynamical model of the power system should include dynamical models of the components and the transmission network. But why should they be interconnected as shown in the particular block diagram? The reason is that this interconnection is enforcing the correct relationship between signals as dictated by the underlying conservation laws.

To understand this, first consider the interconnection point given by the circle with the plus sign. This just means that the signals that flow in sum to the signals that flow out. In our case this means that for each \(i\in\{1,\ldots{},n\}\):

\[P_{d,i}(t)+P_{N,i}(t)=P_{P,i}(t).\]

In words this equation says that at every point in time, the amount of power consumed or produced by the component is equal to the amount that is flowing in from the network plus a local disturbance. This is enforcing conservation of power.

The other point of interconnection, in which the component’s terminal voltages \(V\) flow into the transmission network, is also enforcing a conservation law, albeit more subtly. In words this says that voltage at the terminals of the components and the transmission network are the same, which is equivalent to Kirchhoff’s voltage law.

Dynamical models for frequency control

All that remains to complete our power system model is to mathematically describe the dynamical models in the blocks, and how they affect the signals in the block diagram. In this post we will describe the models suitable for designing systems to regulate the system frequency. But first, a salutary warning:

Power system research is all about the art of making the right assumptions.

Marija Ilić

Developing the standard models for frequency control will involve us making assumptions and simplifications that some readers may find unsatisfactory. However, the block diagram framework we have described is fully general. This means that a model suitable for modular design without these assumptions can be obtained by simply substituting in more general models. I therefore encourage the reader to view this as a first pass at modular design problems. After all if we cannot conduct modular design with a whole load of simplifying assumptions, what hope do we have without them!

Disclaimer finished, lets get on with the modelling.

The dynamics of a transmission line

The fundamental object in the modelling of the transmission network is the three-phase transmission line. The simplest mathematical model of such a line is given by the linear model

\[P_l(t)=Y_l\,(\theta_x(t)-\theta_y(t)).\]

In this equation \(P_l(t)\) is the power flow along the line relative to some nominal flow. The actual power flow along the line equals

\[P_{l,\textrm{eq}}+P_l(t),\]

where \(P_{l,\textrm{eq}}\) is the constant nominal flow determined by the operating point. From a dynamical perspective, \(P_l(t)\) is the important part of the power flow. In this model it depends on three things, the constant \(Y_l,\) and the voltage angles \(\theta_x(t),\theta_y(t)\) at the terminals at either end of the line (terminal is a fancy word for `point where transmission lines and components can be connected together’). We will elaborate on the meaning of voltage angle later, but in short, for this model, specifying the voltage angle is equivalent to specifying the AC voltage.

This model is obtained from the linearisation of a first principles model describing the power flow along the line:

\[P_{l,\textrm{n}}(t)=\frac{|V|_x(t)|V|_y(t)}{X_l}\sin{(\theta_{x,\textrm{eq}}+\theta_x(t)-\theta_y(t)-\theta_{y,\textrm{eq}})}\]

The derivation of this equation is for another time. Perhaps you remember it from an undergraduate course, though you certainly won’t have seen the time dependence included there (rigorously demonstrating the validity of this is not for the faint of heart, which is why we omit the derivation here). However we will derive the linearisation, because there is something fishy going on.

First though lets understand the quantities appearing in the equation. The constant \(X_l\) is the impedance of the transmission line, and depends on the physical properties of the material used in its construction, and the length of the line. The constants \(\theta_{x,\textrm{eq}},\theta_{y,\textrm{eq}}\) are the nominal values of the voltage angles, and are related to the nominal power flow in the line. Finally \(|V|_x(t),|V|_y(t)\) denote the voltage magnitudes at either end of the line.

But what are these voltage angles and magnitudes? Simply put, they are a convenient set of coordinates for describing a three phase voltage. The specifics of the coordinate transform are rather involved, but are closely related to the idea of a phasor. The details are not required for this post, but are briefly explained in the following figure:

The linear model is based on the Taylor series expansion of the \(\sin\) function in the first principles nonlinear model:

\[P_{l,\textrm{n}}(t)=\underbrace{\frac{|V|_x(t)|V|_y(t)}{X_l}\sin{(\theta_{x,\textrm{eq}}-\theta_{y,\textrm{eq}})}}_{P_{l,\textrm{eq}}}+\underbrace{\frac{|V|_x(t)|V|_y(t)}{X_l}\cos{(\theta_{x,\textrm{eq}}-\theta_{y,\textrm{eq}})}}_{Y_l}(\theta_x(t)-\theta_y(t))+\ldots{}\]

There is nothing wrong with this Taylor series, and the idea is that for small deviations in the voltage angles the linear and nonlinear descriptions will agree. That is

\[P_{l,\textrm{n}}(t)-P_{l,\textrm{eq}}\approx{}P_l(t).\]

This is a rather good approximation provided \(\theta_{x,\textrm{eq}}-\theta_{y,\textrm{eq}}\) is not too large (this angle is rarely larger than \(15^\circ\) in practice).

However, closer inspection of the equations shows that things don’t quite match. This is because the two time-invariant constants \(P_{l,\textrm{eq}}\) and \(Y_l\) in the linear model are associated with time-varying things from the Taylor series. The justification for using time-invariant constants in the linear model is that the systems for regulating the network voltages are so fast that we may essentially treat \(|V|_x(t)\) and \(|V|_y(t)\) as constants.

This is worth stating as an explicit assumption, since it directly affects how to mathematically describe the voltages in the block diagram.

Assumption I: The time varying nature of the AC voltage at the \(i\)-th terminal in the network is fully described by the voltage angle \(\theta_i(t)\) at that point.

However is this really true? A more rigorous approach would be to include the voltage magnitudes \(|V|_i(t)\) as variables in our transmission line model — i.e. use a two dimensional signal to describe the time varying AC voltage. However the experience of power system engineers (at least the ones I have spoken to) is that this is not required, at least for frequency regulation problems. And remember that existing power systems are extremely reliable, so they do probably know what they are talking about.

It is worth noting at this point that, strictly speaking, we really need to make an analogous assumption pertaining to the power flows. This is because we have messed with the notion of a voltage, which means that the notion of power in our model is not guaranteed to have a physical basis.

Assumption II: The time varying nature of the power transmitted by a transmission line is fully described by \(P_l(t)\).

While we’re on the topic of assumptions, there are a few other things that may strike you as odd in the above development:

1. How is it possible to describe the voltages in three wires (the three-phase voltage) in a coordinate system with two variables \(|V|(t)\) and \(\theta(t)\)?
2. Why is power a scalar quantity? There are three phases in each line, so surely any power balance would require equations in three quantities, not one.
3. Where are the dynamics? This looks like a phasor analysis, which is only valid in steady state. Transmission lines are dynamically described by a set of PDEs called the telegraphers equations – where did these go?

These are all valid criticisms, and we are required to assume that our model is still a good approximation of the `true’ behaviour, despite all these deficiencies. At the time of writing my Ph.D thesis the need to make all these odd looking assumptions, while standard in the power systems literature, did bother me considerably. There is a lot more to be said here, and I really hope to share further insights in future posts. For now, I will repeat the quote from earlier:

Power system research is all about the art of making the right assumptions.

Hopefully the above discussions have brought this statement into clearer perspective!

The dynamics of a transmission network

The dynamical model of the transmission network can be assembled from the models for the individual lines. To illustrate the idea, lets do an example with two lines, and two points where components can be connected to the grid.

First we do a power balance at every node:

\[
\begin{bmatrix}
P_1(t)\\P_2(t)\\0\end{bmatrix}
=\begin{bmatrix}
1&0\\0&-1\\-1&1
\end{bmatrix}
\begin{bmatrix}
P_a(t)\\P_b(t)
\end{bmatrix}.\]

Next we introduce the voltages by substituting in the transmission line models.

\[\begin{aligned}
\begin{bmatrix}
P_1(t)\\P_2(t)\\0\end{bmatrix}
&=\begin{bmatrix}
1&0\\0&-1\\-1&1
\end{bmatrix}
\begin{bmatrix}Y_a&0\\0&Y_b\end{bmatrix}
\begin{bmatrix}
\theta_1(t)-\theta_f(t)\\\theta_f(t)-\theta_2(t)
\end{bmatrix},\\
&=
\underbrace{\begin{bmatrix}
1&0\\0&-1\\-1&1
\end{bmatrix}}_{B}
\begin{bmatrix}Y_a&0\\0&Y_b\end{bmatrix}
\underbrace{\begin{bmatrix}
1&0&-1\\0&-1&1
\end{bmatrix}}_{B^T}
\begin{bmatrix}\theta_1(t)\\\theta_2(t)\\\theta_f(t)\end{bmatrix}.
\end{aligned}\]

Notice the nice symmetrical structure that is emerging. In particular, the matrix \(B\) is an oriented incidence matrix. This means that when we multiply the matrices out we get

\[
\begin{bmatrix}
P_1(t)\\P_2(t)\\0\end{bmatrix}
=\underbrace{\begin{bmatrix}
Y_a&0&-Y_a\\0&Y_b&-Y_b\\-Y_a&-Y_b&Y_a+Y_b
\end{bmatrix}}_{\mathbf{Y}}
\begin{bmatrix}
\theta_1(t)\\\theta_2(t)\\\theta_f(t)
\end{bmatrix},\]

where \(\mathbf{Y}\) is a weighted Laplacian matrix. Notice in particular that the \(Y\)-parameter for a line is only appearing in the \(i,k\)-position in this matrix if the \(i\)-th and \(k\)-th terminals are connected to the line. This pattern holds in general, allowing the \(\mathbf{Y}\) matrix for any network to be systematically filled out based on the \(Y\)-parameters and the network topology.

Finally we eliminate the node where no component is connected. Such nodes are called floating buses. This is done by algebraically eliminating \(\theta_f\) just as we would when solving simultaneous equations. We then arrive at our transmission network model:

\[
\begin{bmatrix}
P_1(t)\\P_2(t)\end{bmatrix}
=\left(\begin{bmatrix}
Y_a&0\\0&Y_b
\end{bmatrix}-\begin{bmatrix}
-Y_a\\-Y_b
\end{bmatrix}\cdot{}\frac{1}{Y_a+Y_b}\cdot{}\begin{bmatrix}
-Y_a&-Y_b
\end{bmatrix}\right)
\begin{bmatrix}
\theta_1(t)\\\theta_2(t)
\end{bmatrix}.\]

This process can be generalised to networks with any number transmission lines and nodes. First we form the \(\mathbf{Y}\) matrix for the network, including all floating buses:

\[
\begin{bmatrix}P(t)\\0\end{bmatrix}=
\begin{bmatrix}
\mathbf{Y_{11}}&\mathbf{Y_{12}}\\\mathbf{Y_{21}}&\mathbf{Y_{22}}
\end{bmatrix}
\begin{bmatrix}\theta(t)\\\theta_f(t)\end{bmatrix}\]

We then eliminate the floating buses. This can be done using a matrix operation that is typically referred to as Kron reduction:

\[P(t)=\underbrace{\left(\mathbf{Y_{11}}-\mathbf{Y_{12}}\mathbf{Y_{22}^{\mathrm{-1}}}\mathbf{Y_{21}}\right)}_{\mathbf{Y_{\textrm{red}}}}\theta(t).\]

We are almost done. However in our block diagram the output from our network model was the power transmitted to the components, rather than the power injected into the network at the point of connection. Therefore

\[P_N(t)=-\mathbf{Y_{\textrm{red}}}\theta(t).\]

And this is it! Our structured (but likely not sparse) dynamical model of the transmission network! This equation approximately describes the amount of power that is transmitted to the components (relative the operating point) as a function of the voltage angle deviations (an approximate description of the change in the AC voltage at the component terminals relative to the operating point).

The component dynamics

The components correspond to the consumers and producers of electrical power, for example:

• A fossil fuel power station;
• A wind farm or turbine;
• A factory;
• A house or village.

It may seem strange that some of these components correspond to individual things (e.g. a wind turbine) and others collections of things (e.g. a village). This type of aggregation is quite common in electrical engineering — lets just run with it for now.

As we can see, whether they have been aggregated or not, we need a flexible class of models to describe the many different types of component. To this end, we introduce the following component block diagram:

The component model consists of two dynamical elements, \(G_i\) and \(c_i\). These correspond the parts that depend on the inherent properties of the component, and to the parts that we can design (the control system).

Unlike for the transmission network, we will only specify the class of dynamical models that can be used to describe \(G_i\) and \(c_i\). This is to give our model maximum flexibility, while still imposing some structure. We assume that each \(G_i\) can be described by a state-space model:

\[\begin{aligned}
\dot{x}_{G,i}(t)&=\mathbf{A}_{G,i}x_{G,i}(t)+\mathbf{B}_{G1,i}P_{P,i}(t)+\mathbf{B}_{G2,i}P_{c,i}(t),\;x_{G,i}(0)\in\mathbb{R}^{n_i},\\
\theta_i(t)&=\mathbf{C}_{G1,i}x_{G,i}(t),\\
z_i(t)&=\mathbf{C}_{G2,i}x_{G,i}(t).
\end{aligned}\]

In the above \(x_{G,i}(t)\in\mathbb{R}^{n_i}\) is the state, and \(\mathbf{A}_{G,i},\mathbf{B}_{G1,i},\mathbf{B}_{G2,i},\mathbf{C}_{G1,i},\mathbf{C}_{G2,i}\) are matrices with entries determined by the system parameters and operating point. This is an extremely broad modelling class that can be used to approximate the dynamics of almost any component — a reasonably elementary introduction to this topic, including examples, can be found here.

A subtle but important point here is that this does not imply that we actually have a state-space model of the component. It only implies that there exists a suitable model within the class of state-space models. This gives us the freedom to use other ways to describe the dynamics of the components. For example, we could use a model based on the response of the component to simple experimental inputs, such as a frequency response. This is of particular importance here, because good state-space models may not be available for certain types of component (for example, I have never seen a state-space model of the power-voltage characteristics of a village). The modular design tools that we will develop can actually be applied to an even broader class than the state-space models — the full details are in the paper.

Similarly, we assume that the control system — the part of the component that we can design — also has dynamics described by a state-space model:

\[\begin{aligned}
\dot{x}_{c,i}(t)&=\mathbf{A}_{c,i}x_{c,i}(t)+\mathbf{B}_{c,i}z_i(t),\;x_{c,i}(0)\in\mathbb{R}^{m_i},\\
P_{c,i}(t)&=\mathbf{C}_{c,i}x_c(t)+\mathbf{D}_{c,i}z_i(t).
\end{aligned}\]

Analogously, \(x_{c,i}(t)\in\mathbb{R}^{m_i}\) is the controller state. However this time we can choose the entries of these matrices \(\mathbf{A}_{c,i},\mathbf{B}_{c,i},\mathbf{C}_{c,i},\mathbf{D}_{c,i}\). Remember, our long term goal is to work out how to do this, and this will be the focus of Part V.

These two models can be combined to obtain a single state-space model that describes the full dynamics of the component. This is done by eliminating the internal signals \(z_i(t)\) and \(P_{c,i}(t)\) and yields

\[\begin{aligned}
\underbrace{\begin{bmatrix}\dot{x}_{G,i}\\\dot{x}_{c,i}\end{bmatrix}}_{\dot{x}_i(t)}&=\underbrace{\begin{bmatrix}\mathbf{A}_{G,i}+\mathbf{B}_{G2,i}\mathbf{D}_{c,i}\mathbf{C}_{G2,i}&\mathbf{B}_{G2,i}\mathbf{C}_{c,i}\\
\mathbf{B}_{c,i}\mathbf{C}_{G2,i}&\mathbf{A}_{c,i}\end{bmatrix}}_{\mathbf{A}_i}\underbrace{\begin{bmatrix}{x}_{G,i}\\{x}_{c,i}\end{bmatrix}}_{{x}_i(t)}+\underbrace{\begin{bmatrix}\mathbf{B}_{G1,i}\\0\end{bmatrix}}_{\mathbf{B}_i}P_{P,i}(t),\;x(0)\in\mathbb{R}^{n_i+m_i}\\
\theta_i(t)&=\underbrace{\begin{bmatrix}\mathbf{C}_{G1,i}&0\end{bmatrix}}_{\mathbf{C}_i}x_i(t).
\end{aligned}\]

The structured model

We are now ready to assemble the state-space model for the full power system. First we stack up the state space models for the individual components into one giant system:

\[\begin{aligned}
\underbrace{\begin{bmatrix}\dot{x}_1(t)\\\vdots{}\\\dot{x}_n(t)\end{bmatrix}}_{\dot{x}(t)}
&=\begin{bmatrix}\mathbf{A}_1\\&\ddots{}\\&&\mathbf{A}_n\end{bmatrix}
\underbrace{\begin{bmatrix}{x}_1(t)\\\vdots{}\\{x}_n(t)\end{bmatrix}}_{{x}(t)}+\underbrace{\begin{bmatrix}\mathbf{B}_1\\&\ddots{}\\&&\mathbf{B}_n\end{bmatrix}}_{\mathbf{B}}P_P(t),\;x(0)=\begin{bmatrix}x_1(0)\\\vdots{}\\x_n(0)\end{bmatrix}\\
\theta(t)&=\underbrace{\begin{bmatrix}\mathbf{C}_1\\&\ddots{}\\&&\mathbf{C}_n\end{bmatrix}}_{\mathbf{C}}x(t).
\end{aligned}\]

Then we eliminate \(P_P(t)\) using the equation for conservation of power

\[P_P(t)=P_d(t)+P_N(t)\]

and the transmission network model. This gives
\[
\dot{x}(t)=\underbrace{\left(\begin{bmatrix}\mathbf{A}_1\\&\ddots{}\\&&\mathbf{A}_n\end{bmatrix}-\mathbf{B}\mathbf{Y}_{\textrm{red}}\mathbf{C}\right)}_{\mathbf{A}}x(t)+\mathbf{B}P_d(t)
,\;x(0)=\begin{bmatrix}x_1(0)\\\vdots{}\\x_n(0)\end{bmatrix}.
\]

And that is it! A highly structured, and modularly constructed, dynamical model of a power system. Thank you for reading this far. Next time we will derive a decentralised stability criterion that can be applied to a class of dynamical models with a closely related structure. Here again is a link to the paper for those who can’t wait, and also a link to Part I for those who wish to go back!