Chapter 5 Linear Systems PDF

Chapter 5 Linear Systems One of the key issues in systems behaviour is that of coupling: Nonlinear Linear Linear Uncoupled Coupled Uncoupled Coupled Simple Chaotic Sinusoidal...

Chapter 5 Linear Systems One of the key issues in systems behaviour is that of coupling: Nonlinear Linear Linear Uncoupled Coupled Uncoupled Coupled Simple Chaotic Sinusoidal Non-Periodic A single pendulum exhibits simple, uninteresting periodic motion, however coupling two pendula is a classic example of chaotic behaviour. Even in linear systems, a single mass attached to springs is periodic and sinusoidal, whereas two such masses, although still oscillating, can exhibit motion which is non-periodic. Coupling is central to the study of ecological and social phenomena, since it appears in so many contexts: A town consists of many people, influencing each others’ habits; A city consists of many businesses, who compete in buying and selling; A forest consists of many species, some parasitic, some symbiotic; An ecosystem consists of many ponds, rivers, lakes, streams, with water, nutrients, and pollution flowing between them. On what basis might we be able to take a complex, interacting system and somehow break it into constituent parts? And under what circumstances might decoupling be possible for one type of system and not for another? These are not easy questions to answer, so we need to start with those systems which are relatively simple and mathematically tractable—linear systems. © The Author(s), under exclusive license to Springer Nature Switzerland AG 2021 67 P. Fieguth, An Introduction to Complex Systems, https://doi.org/10.1007/978-3-030-63168-0_5 68 5 Linear Systems Fig. 5.1 Linear Systems: A linear system can refer to several different concepts, depending upon the context. In this chapter we will focus primarily on linear dynamic systems. 5.1 Linearity We begin1 with an arbitrary nth order model ! " z (n) (t) = f z (n−1) (t),... , z (1) (t), z (0) (t), u(t), t (5.1) where z (i) is the i th order time-derivative of z, and u(t) is some external deterministic input. Recall the discussion of state augmentation, allowing higher-order dynamic models to be re-cast as first-order ones. On the basis of that idea, we can narrow our focus to the first-order case: ! " ż(t) = f z(t), u(t), t. (5.2) In this chapter we are considering linear dynamic models, and any linear relationship can be expressed as a matrix-vector product: ż(t) = A(t)z(t) + B(t)u(t) + w(t) (5.3) for some matrices A, B and stochastic (noise) term w. If we assume the model A to be time–stationary (constant over time) and ignore external inputs, then we arrive at the simplest, most fundamental linear-dynamic system: ż(t) = Az(t). (5.4) To be clear, the concept of linear system can mean at least three different, but highly related ideas, illustrated in Figure 5.1. The distinction between these three will be clear from context, as we will not be mixing them. For the time being, we are focusing on linear dynamic systems. In all cases, the key principle of superposition applies: 1The algebra background needed for this chapter is summarized in Appendix A, with the associated notation listed in Appendix D. 5.1 Linearity 69 For a system obeying superposition, the system response to a sum of multiple inputs is just the sum of responses to the individual inputs. Thus, for a linear system of equations Ax = b, If b1 leads to solution x 1 , and b2 leads to solution x 2 Then αb1 + βb2 leads to solution αx 1 + βx 2 For a linear input-output2 system y(t) = h(t) ∗ x(t), If input x1 (t) leads to output y1 (t), and x2 (t) leads to output y2 (t) Then input αx1 (t) + βx2 (t) leads to output αy1 (t) + β y2 (t) For a linear dynamic system2 ż = Az(t), If initial condition z(0) = i 1 leads to time behaviour z 1 (t), and z(0) = i 2 leads to time behaviour z 2 (t), Then initial condition αi 1 + βi 2 leads to behaviour αz 1 (t) + βz 2 (t) The principle of superposition2 is one of the key points of this text and a deep understanding is essential, precisely because it seems so obvious and intuitive. In plain terms, for a linear system the average response is the same as the response to the average: # ! "$ # ! "$ Average System Response Input(t) = System Response Average Input(t) (5.5) Yet virtually all social / ecological / natural systems do not satisfy superposition!, such that (5.5) does not hold. That is, in an ecological (nonlinear) system, two simultaneous disturbances (logging and pollution, for example) which, on their own might have had a modest effect, may together wound or break the system, such that the net response is, in fact, quite different from the sum of the individual ones. For the same reason it is inappropriate to report average pollution levels, such as in a river downstream from a mine, since a constant pollution level τ , versus a level of pollution fluctuating about τ and occasionally surpassing a toxicity threshold, can lead to very different consequences. 2 A dynamic system can have inputs and outputs, and a linear input-output system can be time- dynamic. The distinction is primarily one of representation; that is, whether the system is character- ized by its impulse response h(t), or by its time dynamics A. This text will focus almost entirely on the latter time-dynamic representation, since the concept of impulse response does not generalize well to nonlinear systems. Readers not familiar with convolution (∗) and impulse response h(t) may wish to consult a text on signals and systems, such as. 70 5 Linear Systems 5.2 Modes So, given a first-order linear dynamic system ż = Az (5.6) how can we understand its behaviour? The most beautiful tool in linear algebra is the eigendecomposition. Sadly, most students recall (or were primarily taught) only the mechanics, that the eigendecom- position involves taking a determinant of some matrix and finding the roots of a polynomial, which leaves us with essentially no insight whatsoever. The eigendecomposition takes a system and breaks it into decoupled modes, where the modes are the characteristic or unchanging behaviours of the system. Specifically, suppose we are given the dynamics matrix A in (5.6); if Av = λv (5.7) then we refer to λ as an eigenvalue of A, and v as its corresponding eigenvector or mode. Why does this matter? Observe that if the state at time t = 0 is initialized in the direction of mode v, z(0) = v ż(0) = Az(0) = Av = λv, (5.8) then the change ż also points in the direction of v, so that z continues to point in the direction of mode v. Therefore we can write z(t) = α(t) v and ż(t) = α̇(t) v, (5.9) where α(t) is just a scalar, measuring how far z(t) points in the direction of v. Then, combining (5.8) and (5.9), α̇(t) v = ż(t) = Az(t) = Aα(t)v = α(t)Av = α(t)λv (5.10) from which we find the dynamics of α(t) as α̇(t) = λα(t), (5.11) a simple first-order differential equation, for which the solution follows as α(t) = exp(λt) z(t) = exp(λt) v. (5.12) 5.2 Modes 71 Thus if the state z points in the direction of an eigenvector, then it continues to point in that same direction over time; only the amplitude changes, with the rate of change measured by the time constant τ = 1/λ. In general, however, an n × n matrix A will not have a single mode, rather n of them. If we assume, for the moment, that matrix A is diagonalizable, then we can collect the eigenvalues/eigenvectors     | | λ1   Av i = λi v i , 1 ≤ i ≤ n V = v 1 · · · v n  Λ = ...  (5.13) | | λn where V is invertible and Λ is diagonal. From (5.13) it follows that A V = VΛ (5.14) which, when left-multiplied by V −1 or right-multiplied by V , becomes + A = V Λ V −1 Synthesis / Coupling (5.15) + Λ = V −1 A V. Analysis / Decoupling (5.16) The diagonalization of A in (5.16) refers to the ability of V to transform A to Λ. The power of the eigendecomposition is that, whereas solving ż = Az directly is quite challenging, the eigendecomposition transforms the problem into something simple. Let y(t) represent the amplitudes of the modes of the system, such that y(t) = V −1 z(t) z(t) = V y(t). (5.17) Then ẏ(t) = V −1 ż(t) = V −1 Az(t) = V −1 A V y(t) = Λy(t) (5.18) from which it follows that V −1 ż(t) = Az(t) ẏ(t) = Λy(t) Eig. ∴ ẏi (t) = λi yi (t) (5.19) Easy V z(t) = V y(t) ∴ yi (t) = eλi t yi (0) Eig. 72 5 Linear Systems We can immediately notice a few things: Because Λ is a diagonal matrix, the dynamics for z are decoupled into individual simple scalar dynamics yi (t) for each mode. The dynamics of z(t) need to be initialized, typically by specifying z(0), therefore n , y(0) = V −1 z(0) z(0) = V y(0) = v i yi (0). (5.20) i=1 That is, the weights y1 (0),... , yn (0) of the modes are precisely those weights needed to reconstruct the initial condition z(0) as a weighted sum of mode shapes v1,... , vn. We see that the general solution for z(t) is nothing more than a weighted combi- nation of modes: z(t) = V y(t) = v 1 y1 (t) + · · · + v n yn (t) (5.21) n , = yi (0) · eλi t · v i (5.22) i=1 -./0 -./ 0 -./0 -./0 Combining Weighted Dynamic Modes From (5.22), since z(t) is just a weighted sum of exponentials, broadly speaking there are only two possible behaviours: ! " If Real!λi " < 0 for all i |z(t)| 0 ∴ Stable Dynamics (5.23) If Real λi > 0 for any i |z(t)| ∞ ∴ Unstable Dynamics Therefore all diagonalizable linear systems either decay to zero or grow exponen- tially. Since all natural systems are finite, indefinite exponential growth is not possible, so for any real physical system a linear model is necessarily only an approximation to the actual, bounded, nonlinear system, such as the example in Figure 4.3. 5.3 System Coupling We now know that the eigendecomposition of diagonalizable n × n dynamics A transforms a coupled system of n variables to n decoupled first-order systems: 5.3 System Coupling 73 Example 5.1: The Exponential Function “Anyone who believes exponential growth can go on forever in a finite world is either a madman or an economist.” – Kenneth Boulding The characteristic behaviour (or modes or eigenfunctions) of all linear systems consists of the family of complex exponentials. As a result, with the exception of degenerate cases, all linear systems produce either exponentially shrinking or exponentially growing output over time. In general it is, of course, the exponentially growing systems which concern us. Since people produce more people, and dollars (with interest) produce more dollars, in both cases ż(t) ∝ z(t) (5.24) the solution to which is exponential, which is why such behaviour appears so commonly in societal and economic time series, such as global debt, deforestation, or energy use. So what makes an exponential behaviour so dangerous? In a finite, closed system indefinite exponential growth is not possible, so at some point a limit or carrying capacity is reached, as illustrated in Figure 4.3 and discussed further in Example 6.1. An exponential function must eventually encounter a limit, and history suggests that such limits are encountered in difficult and painful ways. Many people appear to be lulled or complacent of phrases like “modest growth” of 1% or a “low inflation band” of 1% – 3%. However a compounding growth of 1% is still exponential and will ultimately encounter the same system limits, albeit more slowly than under rapid growth. Indeed, our economic system is premised on capital being paid, with interest, in the future; however for capital and interest able to be repaid necessarily and dangerously implies that the money supply must expand exponentially, and indefinitely. The key challenge with exponentials is that they make it exceptionally difficult to learn from history. An exponential function has a finite doubling time: ln 2 eα(t+TDouble ) = 2eαt Doubling time TDouble = , (5.25) α meaning that the time from which a resource is only half used, to being completely exhausted, is just TDouble. Similarly, if the rate at which a resource is used grows expo- nentially as eαt , then since 1 T 1 T +TDouble 1 αT 1 αT eαt dt = e and eαt dt = e (5.26) −∞ α T α it follows that the same amount of a resource was used, for all of recorded (and pre- recorded) history up to time T , as is used in the finite time interval from T to T + TDouble. Further Reading: You may wish to compare the challenges raised here with those of Example 9.6. Problem 5.14 lists a number of follow-up readings and videos related to exponential functions. 74 5 Linear Systems A λ1 λ2 λ3 λ4 Eig. v1 v2 v3 v4 (5.27) -./ 0 -./ 0 Uncoupled Coupled If A is not diagonalizable, this essentially implies that there is a coupling in the system which cannot be undone. One example would be a second-order system, such as the position of a car subject to a random acceleration ẍ(t) = w(t) (5.28) for a random (white noise) process w. Although this dynamic can be converted to first order via state augmentation (with further examples in Appendix C), 2 3 2 3 2 3 ẋ(t) 00 1 z(t) = ż(t) = z(t) + w(t) (5.29) x(t) 10 0 the resulting dynamics matrix cannot be diagonalized since the underlying dynamic really is second order, including teλt terms, which cannot be represented as a weighted sum of first-order exponentials. A dynamic A which is not diagonalizable can, at best, be reduced to a standardized coupling, known as the Jordan form, something like   λ1 Eig.  λ2 1  A   (5.30)  λ2  λ3 graphically corresponding to a maximally reduced system of the form λ1 λ2 λ2 λ3 v1 v2 v2 v3 In general a coupled problem has a variety of possible ways in which it could be decoupled, as illustrated in Figure 5.2. However the intuition remains unchanged from before: An eigendecomposition decomposes a coupled system into uncoupled modes, however some of the modes may be higher order q > 1 and have a characteristic behaviour based on some linear combination of 5.3 System Coupling 75 Fig. 5.2 System Decoupling: Given a fourth-order state, top, there are five possible forms that the decoupling can take, right. eλt teλt... t q−1 eλt (5.31) That is, the simplification / decomposition of (5.19) still applies, however the trans- formed states y i are no longer all necessarily scalar, rather some of them will be vectors, based on the order of the coupling. 5.4 Dynamics Since the preceding discussion looked at the decoupling of coupled systems, we can now focus on decoupled systems, where the simplest possible case is the first-order, scalar, time-invariant one: 76 5 Linear Systems Example 5.2: Coupled Systems Suppose we have two ponds, possibly polluted, with water volumes S1 and S2 , with a water flow of F from one to the other: State z 1 (t) State z 2 (t) F F F S1 S2 (Clean) Clearly there is a coupling between the ponds, since any pollution from Pond 1 flows into Pond 2. Let z 1 (t) and z 2 (t) be the concentration, over time, of a pollutant in the ponds. The system dynamics are 2 3 F F# $ −F/ 0 ż 1 (t) = − z 1 (t) ż 2 (t) = − z 2 (t) − z 1 (t) A = F S1 −F S1 S2 /S2 /S2 Taking the eigendecomposition of A, we find that S1 = S2 A single eigenvalue λ = −F/S , a second-order mode S1 '= S2 Two eigenvalues, λ1 = −F/S1 , λ2 = −F/S2 Thus the second-order Jordan form appears only when the two ponds are identical in volume, which is very unlikely, thus we refer to the second-order case as degenerate. Nevertheless the coupled / Jordan concept has value. Suppose that F = 1, S1 = 10, S2 = 11 z 1 (0) = 5, z 2 (0) = 0 such that the ponds are similar in volume. The eigendecomposition of A gives two first- order modes y1 , y2 , and from (5.17) the transformation is z 1 (t) = 20y1 (t) z 2 (t) = 200y2 (t) − 199y1 (t) Thus the state of the second pond is determined by the small difference between two large numbers y1 , y2 , which does not offer much insight. If we numerically simulate the dynamics, we obtain 0.3 5 Pollution Concentration Pond 1 4 Mode Amplitude 0.2 3 Pond 2 Mode 2 2 0.1 1 Mode 1 0.0 0 0 2 4 6 8 10 0 2 4 6 8 10 Time Time Note that the two modes (left) are, indeed, both simple first-order exponentials, but without intuitive meaning (in this case). On the other hand, the actual concentration states (right) look very nearly like z 1 ( e−λ1 t z 2 ( te−λ2 t where from the Jordan form (5.31) we understand the tet term of z 2 to be due to the fact that we essentially have a coupled second-order system. 5.4 Dynamics 77 Fig. 5.3 Graphical interpretation of a simple linear dynamic system: Wherever ż > 0 we know that z is increasing, and similarly that ż < 0 implies that z is decreasing. Therefore the system to the left, with negative slope α, will tend to return to z̄ if disturbed, and is stable, whereas the system to the right, with positive slope α, will tend to diverge and is unstable. ż = αz + β (5.32) Recall from Section 4.1 the concept of a fixed point, which is trivial to evaluate for the system of (5.32): β 0 = αz̄ + β z̄ = − ż = α(z − z̄) (5.33) α The dynamic “matrix” A is just the scalar α, having eigenvalue λ = α, for which we very well know the behaviour of the system as z(t) = z̄ + eαt · (z(0) − z̄), (5.34) thus, consistent with (5.23), asymptotically the dynamics evolve as α 0 |z| ∞ (unstable) Indeed, Figure 5.3 shows the same behaviour graphically. The graphical interpreta- tion is important, since for many nonlinear systems an exact analytical understanding will be exceptionally challenging, whereas a qualitative graphical interpretation will be relatively straightforward. With an understanding of one-dimensional dynamics, we wish to generalize this understanding to multiple dimensions. Given a dynamic system ẋ = Ax + b, (5.36) as with (5.32) there is a single fixed point, easily evaluated as 0 = A x̄ + b x̄ = −A−1 b. (5.37) 78 5 Linear Systems Suppose the system experiences a disturbance at time t = 0, such that x(0) '= x̄. We assume the fixed point to be stable, so the system dynamics will return the disturbed state to the fixed point. We can consider two difference influences x(0) = x̄ + δ x δ (t) lim x δ (t) = x̄, t ∞ (5.38) x(0) = x̄ + & x & (t) lim x & (t) = x̄. t ∞ If the system is subject to both influences, the system clearly still converges to the fixed point: x(0) = x̄ + δ + & x δ& (t) lim x δ& (t) = x̄ (5.39) t ∞ meaning that x δ& (t) '= x δ (t) + x & (t), (5.40) appearing to contradict superposition. However it is the response to disturbance that superposes, not the constant fixed point. That is, superposition must be understood relative to the fixed point. As a result we will, in general, prefer to translate a linear dynamic problem to shift its fixed point to the origin, so we will study the simplified state z = x + A−1 b such that (5.36) becomes ż = Az. (5.41) That is, z has the same dynamic behaviour as x, just a shifted version placing its fixed point at the origin. Assuming the dynamics to be diagonalizable, which we now know to mean that the individual modes are all first order, then given the decomposition A = VΛ V −1 , from (5.19) we know that the dynamics of z can be written as   eλ1 t  | | ..  −1 z(t) = v 1 · · · v n  .  V z(0) (5.42) | | λn t -./ 0 e = vector c   eλ1 t c  | | 1 ..  = v 1 · · · v n  .  (5.43) | | λn t e cn n , z(t) = v i eλi t ci (5.44) i=1 5.4 Dynamics 79 That is, over time, z(t) is nothing more than a weighted sum of growing or shrinking eigenvectors, where the weights ci are determined by the system initial condition z(0). For the original system of (5.36), having a fixed point at x̄, (5.42) generalizes to n , x(t) = x̄ + v i eλi t ci. (5.45) i=1 That is, precisely the same eigendecomposition and dynamic behaviour of (5.36), but now converging to fixed point x̄ rather than to the origin. Next, suppose we are given the dynamic system 2 3 0 -1 ż = z (5.46) -1 0 for which the eigendecomposition is 2 3 2 3 2 3 0 -1 1 1 λ1 = −1 v 1 = λ2 = +1 v 2 = (5.47) -1 0 1 -1 From (5.42) we therefore know the general form of the solution to this system to be 2 3 2 3 −t 1 t 1 z(t) = c1 e + c2 e. (5.48) 1 -1 Thus the system tends to grow in one direction and shrink in another, known as a saddle, illustrated in the phase portrait of Figure 5.4: A phase plot or phase portrait is a plot of state trajectories of a dynamic system. Once we understand the saddle example well, the concept is relatively easy to gen- eralize to other dynamic matrices, as illustrated in Figure 5.5: for a two dimensional problem we will always have two eigenvectors, where the associated eigenvalues may be positive (unstable), negative (stable), or zero (a degenerate3 case). The eigende- composition is so effective at characterizing a system, that we refer to the pattern or distribution of eigenvalues as the spectrum, however we will see (Section 5.6) that there are circumstances where the spectrum does not quite tell the whole story. It is possible for the eigenvalues to be complex, but because the dynamics matrix is real the eigenvalues and eigenvectors must be present in conjugate pairs 3 Degeneracy will be defined more carefully in Chapter 6; it essentially means that a circumstance (e.g., zero eigenvalue) is highly sensitive to perturbation. 80 5 Linear Systems Fig. 5.4 Phase Plot of a Saddle: If a linear system has one stable mode (v 1 ) and one unstable mode (v 2 ) then the resulting mixture of convergent and divergent behaviour is known as a saddle. Fig. 5.5 Two illustrations of linear dynamic systems: stable (left) and unstable (right). In these examples the eigenvalues are not equal. For example, in the stable system λ1 < λ2 , therefore the v 1 component of the state converges more rapidly to the fixed point than the v 2 component, causing curved trajectories. λ1 = a + ib λ2 = a − ib v 1 = α + iβ v 2 = α − iβ (5.49) Consequently from (5.42) the associated system dynamics must have a form some- thing like ! " ! " eλ1 t = exp (a + ib)t = eat cos bt + i sin bt , (5.50) except that the actual dynamics must clearly be real (no complex i), since the given dynamic model ż = Az is real; the complex term in (5.50) disappears because of the 5.4 Dynamics 81 Fig. 5.6 Systems with Complex Eigenvalues: Complex eigenvalues lead to oscillatory behaviour, where the real part of the eigenvalue determines whether the system is stable (right), a marginally stable centre (left), or unstable (not shown). complex conjugacy of the eigenvalues and eigenvectors. So from (5.50), we observe a spiral / circular behaviour, with the spiral growth / shrinkage controlled by a, and the angular velocity controlled by b. Given a 2D problem ż = Az, if λ = a + ib is a complex eigenvalue of A, with associated complex eigenvector v = α + iβ, then the (real-valued) solution is ! " z(t) = eat cos(bt)α − sin(bt)β. (5.51) Two examples of this case are shown in Figure 5.6. 5.5 Control of Dynamic Systems Although our focus is on the actual dynamics governing a system, most systems do also admit an input, some sort of deliberate influence. Obvious examples include the gas pedal and brake in your car, or the volume control on a radio. 82 5 Linear Systems Fig. 5.7 Feedback Control: Given some dynamic system, we can attempt to control or regulate it by observing the unknown system state z, indirectly, via measurements m, and producing a control input u. In most cases the system will be subject not only to our deliberate signal u, but also to other disturbances and influences, whether human or otherwise. Manipulating a system input in order to regulate the behaviour of a system is known as control theory, a vast topic of which we can only begin to discuss here. There may be two sorts of human-driven “input” to the system: Deliberate feedback control in which the output of the system is observed, and that some function of that output is fed back into the input in order to produce the desired regulation. These inputs are explicitly designed to control the system, such as the re-introduction of native species or the strengthening of wetland protection. Accidental inputs, usually as a side-effect of other human activity, such as the leaching of fertilizers into waterways or the emission of CO2 into the atmosphere. These inputs clearly influence the system, however normally to perturb or imbal- ance the system rather than to regulate it. We will not be designing control strategies in this book, however we will definitely want to understand the human influence on systems, the manner in which systems respond to inputs, whether deliberate or accidental, as illustrated in Figure 5.7. Thus given a linear dynamic system, ż(t) = A(t)z(t) (5.52) as studied throughout this chapter, really there are control inputs u(t) and other accidental / disturbance inputs w(t) present, ż(t) = A(t)z(t) + Bc (t)u(t) + Bo (t)w(t), (5.53) where Bc , Bo describe the manner in which inputs u and w, respectively, affect the system state. Some control problems are not difficult to solve, such as balancing a tennis racket upside-down on your hand, as in Figure 5.8, or regulating the temperature of a 5.5 Control of Dynamic Systems 83 Fig. 5.8 Control of an Unstable System: An inverted pendulum, or tennis racket, has an unstable fixed point when held upright, however the unstable system can be stabilized given an appropriate feedback control input u, based on measurements m of the system state z. house4 in the winter. In contrast, most socio-environmental systems are challenging to control for one of four reasons: 1. Bounded Inputs: For many large-scale systems the size of human input which can be asserted is tiny, such that the input signal is simply too small to control the system. In contrast to accidental inputs, such as deforestation, which can be very large, deliberate inputs need to be undertaken on top of or in opposition to other human activity, which therefore requires time, political will, money, and/or energy, and may be limited by any or all of these reasons. 2. Long Time Constants: In classical control, systems having long time constants (i.e., those with slow dynamics or small eigenvalues) are actually rather easy to control: the system is slow, making it quite simple to design a sufficiently fast controller. However most socio-environmental systems rely on political will or human urgency, which appears only once the system has actually begun to manifest noticeable changes (e.g., global warming, ozone hole), at which point the long time constants ensure that further change will continue to take place for decades. 3. Observability: The control of a system depends on the state of the system being observable or inferable in some way, a topic to be discussed further in Chapter 11. In some cases we may have an exceptionally hard time acquiring measurements, such as in the deep ocean, or we may just not know what to measure, such as the unanticipated atmospheric chemistry surrounding the arctic ozone holes. 4Although, for practical reasons, home temperature control turns out to have nonlinearities present, as will be discussed in Example 6.4. 84 5 Linear Systems 4. Time Lags: The longer the time lags in a system, the more difficult it is to control. System lags may be present in terms of control inputs, or in terms of the state observation. For example, the Covid-19 pandemic (Chapter 13) is challenging in large part because of state–observation lags, the relatively long time between a person acquiring the illness (change of state) and the appearance of symptoms (resulting observation). Similarly, although most environmental systems do not exhibit observational lags, there can be significant lags on the control side—the length of time for people to reach a consensus, to formulate policy, to deal with court challenges etc. all represent delays in affecting a response. These four issues, when combined, can lead to terribly challenging predicaments. In the context of global warming, for example, the time constant for atmospheric CO2 is on the order of decades, so the current rate of temperature rise will continue for many years; at the same time, limiting the temperature rise would require very large inputs, beyond the scope of what humans are willing or politically able to do; it is difficult to measure global CO2 , and it isn’t necessarily obvious what to even measure (methane?, temperature?, permafrost?, arctic ice?); finally, as is painfully obvious to anyone concerned with environmental matters, the repeated failures and delayed attempts to negotiate effective responses to global CO2 has led to extensive delays. We will briefly return to the control problem in the context of nonlinear dynamic systems in Section 7.4, and the Covid-19 challenge in Chapter 13. 5.6 Non-Normal Systems From the general solution to a linear system in (5.42), we know that for a two- dimensional system the dynamic behaviour behaves as z(t) = c1 eλ1 t v 1 + c2 eλ2 t v 2. (5.54) For a stable system, where both eigenvalues have negative real parts, it is therefore natural to assume that a given state converges smoothly to the fixed point at the origin. It is, therefore, perhaps surprising that even if stable, (5.54) can diverge before con- verging to the fixed point. Such a situation can occur if the system is non-normal, meaning that the eigenvectors are not at right angles to one another. There is nothing invalid about non-orthogonal eigenvectors; indeed, we have seen such examples, as in Figure 5.5. The orthogonal case is so commonly studied (Problem 5.9) that we sub- consciously assume that all systems must be so; however, perversely, the orthogonal 5.6 Non-Normal Systems 85 Initial Transient Long Term 250 300 Distance from Fixed Point Distance from Fixed Point 250 200 200 150 150 100 100 50 50 0 0 0 0.1 0.2 0.3 0.4 0.5 0 1 2 3 4 5 Time Time Fig. 5.9 Non- normal Systems: Examining the initial behaviours of two systems, left, which is stable and which is unstable? If they were normal systems, both red and blue would be considered unstable. However after a longer period of time, right, we observe that the non-normal red system of (5.56) is, indeed, stable. It is not possible to say whether the blue system is stable or not, and would require an even longer period of observation, related to the baseline question of Chapter 4. case is actually degenerate, and can therefore really only be expected in contrived circumstances! The issue boils down to taking differences of large numbers. If z = (1001) − (1000) (5.55) then z is small, of course, even though there are large numbers involved. However, if we now add differing dynamics, z(t) = 1001 · e−t − 1000 · e−2t (5.56) then z(t) starts small and grows to a large value before converging back to zero, as illustrated by the red curve in Figure 5.9. The mechanics of the behaviour are examined in detail in Figure 5.10. So the distribution of eigenvalues, the so-called spectrum of a system, clearly reveals whether a system is stable, but does not tell us whether there may be transients present. Testing whether a system is normal or not requires us to also know the eigenvectors; if the eigenvectors are stacked in V , as in (5.13), then system normality is easily checked by taking dot-products of the eigenvectors, as Q diagonal System normal Let Q = V T V, then (5.57) Q not diagonal System non-normal The lesson, discussed in Figure 5.9, is that we need to be quite cautious in assessing a system on the basis of observed behaviour. Similar to the baseline question of 86 5 Linear Systems z2 v1 z2 v1 z2 v 1 , λ1 z2 v1 c1 z(δ) z(0) z(0) z(0) z(0) z1 z1 z1 z1 λ2 < λ1 c2 v2 v2 v 2 , λ2 v2 Fig. 5.10 Non- normal Systems: Suppose we are given non-orthogonal eigenvectors v 1 , v 2 and a system state z(0), for which we can find the components c1 , c2 by projecting the state onto the eigenvectors. Even if both eigenvalues are stable, if they are sufficiently different (c) the mode components will shrink at very different rates, causing the state trajectory initially to grow (d). Chapter 4, asking when to start measuring, we now have a related question: how long to measure. The issue of non-normality plays a significant role in the assessment of resilience (see Example 5.3). Resilience is very much an empirical concept, in that we wish to observe a given system, such as a forest, and then determine the impact of a disturbance, such as building a road through the forest or a dam nearby. Because the system needs to be assessed on the basis of observation, over a finite period of time, it is quite important to appreciate the possibility of non-normality. There are many examples involving the resilience of non-normal systems, including aquatic food chains, predator-prey models, and nutrient flows in tropical systems , where the last of these is explored in some detail in Problem 5.13. Case Study 5: System Decoupling Problem: Describe, in general, the analytical motion of the following spring–mass system: k0 m1 k1 m2 k2 kn + 1 ··· Case Study 5: System Decoupling 87 Example 5.3: Resilience I—Linear Systems Resilience measures how robust or tolerant a given social or ecological system is to some sort of disturbance, of which there are many: Ecological Disturbances: acid rain, logging, pollution, invasive species Social Disturbances: terrorism, electrical blackouts, recessions Asymptotically, after a long time, it is the slowest mode (least negative eigenvalue) which limits how slowly a linear system returns to the fixed point after a disturbance. Therefore the classic measure of system resilience is 4 ! "+ Resilience(A) = min −Real λi (A) (5.58) i The more rapid the return to the stable fixed point, the larger the resilience. It is only the real part of the eigenvalue which is needed here, since it controls the exponential decay; the complex part controls only the frequency of any oscillation. However, as we know from Section 5.6, for non-normal systems there can be significant transient behaviour which is not measured in the asymptotic t ∞ assessment of (5.58), therefore other measures have been proposed to characterize the transient. The reactivity measures the greatest rate of divergence at time t = 0 of the transient relative to the size of perturbation: 5 6 67 8 9 : # ;< 1 d 6z 6 77 A + AT $ Reactivit y(A) = max 6 6 6z 6 d t 7 7 = max Real λi (5.59) zo o i 2 t=0 A reactivity greater than zero implies a non-normal further growth in the initial state beyond the initial perturbation. Whereas the reactivity is an instantaneous assessment at t = 0, to assess the system over time an alternative would be the maximum amplification, 56 68 6z(t)6 κmax = max κ(t) where κ(t) = max 6 6. 6z 6 (5.60) t zo o Given a definition of resilience, there are at least three socio-environmental contexts of interest: Ecological Resilience: Given a list of alternatives, we wish to select a human action which minimizes ecological disruption. That is, to which disturbance is the ecological system most robust? Social Resilience: What strategies can a government put in place to make a society resilient to social disruptions and other challenges? That is, we wish to transform the societal dynamics to increase resilience. Social policy: An effective policy needs to maximize its impact for a given amount perturbation (money or energy). That is, we seek the opposite of resilience: what timing or direction of perturbations maximizes the response? Thus, although (5.58) – (5.60) may seem abstract, the study of resilience is of considerable importance, to be followed up in Example 7.3. Further Reading: There are many good references on the topic of resilience: [3, 5, 6, 9, 12]. 88 5 Linear Systems Not easy!! It is easy enough, however, to write down the system dynamics. It is most convenient for the state z i to represent the deviation of mass i from its rest position, thus in the preceding figure z i = 0. Applying Newton’s F = ma to each mass: ! " ! " m i z̈ i = −ki z i − z i−1 + ki+1 z i+1 − z i (5.61) At the end points the equations look a bit different, since the springs are attached to fixed walls, rather than to a movable mass: ! " ! " m 1 z̈ 1 = −k1 z 1 − 0 + k2 z 2 − z 1 ! " ! " (5.62) m n z̈ n = −kn z n − z n−1 + kn+1 0 − z n The problem, of course, in trying to solve these equations is that n masses can oscillate in all kinds of complicated, interacting ways. The key idea is that the system is exceptionally easy to solve when we have only one mass. Let us first collect equations (5.61) and (5.62) in matrix–vector form:       m1 z̈ 1 (−k1 − k2 ) k2 z1 ..  ..  ..  ..  . .=.... . mn z̈ n kn (−kn − kn+1 ) zn (5.63) -./ 0 -./ 0 -./ 0 -./ 0 M z̈ = K z thus z̈ = M −1 K z. (5.64) The matrix (M −1 K ) describes the interactions between the masses; we can find the eigendecomposition of (M −1 K ) (M −1 K ) = V Λ V −1 (5.65) to rewrite the system of equations: z̈ = M −1 K z = V Λ V −1 z. (5.66) If we then let y = V −1 z, we have ÿ = V −1 z̈ = V −1 V Λ V −1 z = Λ V −1 z = Λ y. (5.67) However our new, transformed state variables y are non-interacting, since the spring– mass system is diagonalizable, meaning that Λ is a diagonal matrix, giving us an Case Study 5: System Decoupling 89 ensemble of simple, periodic single mass systems ÿi = λi yi (5.68) Interestingly, our decoupled system has normalized unit mass with a spring constant controlled by the eigenvalue: −λ m=1 Let us now consider a simplified version of the problem, with equal masses and spring constants: k m k m k m k The dynamic system from (5.61) and (5.62) can neatly be written as   -2k k 0 m z̈ =  k -2k k  z (5.69) 0 k -2k allowing k to be factored out, such that the interactions can be expressed as a simple integer matrix for which the eigendecomposition is    1 √1 1    1 √1 1 −1 -2 1 0 2 2 2 −0.6 2 2 2  1 -2 1  =  1 1   √2 0 − √2   −2  1  √2 0 − √2  1  (5.70) 0 1 -2 1 2 − √1 1 2 −3.4 1 2 − √1 1 2 -./ 0 -./2 0-./ 0 2 MatrixK Eigenvect. in Columns Eigenvalues on Diagonal Recall that the eigendecomposition tells us the modes of the system, those vibrations which will persist over time. For example, the vibration of m 1 on its own is not a mode, since vibrating m 1 will quickly cause m 2 and m 3 to vibrate. However all three weights sliding back and forth — that’s a mode. The eigenvectors in (5.70) actually identify the three modes for us: 90 5 Linear Systems Mode Eigenvector Motion Pattern Spring Constant # $ Slow 1 √1 1 I. 2 2 2 k = 0.6 # $ Medium II. √1 0 − √12 2 k=2 # $ Fast 1 III. − √12 1 2 2 k = 3.4 However our goal was to identify a closed-form solution to the coupled system. How can we use the above decomposition to do this? Suppose we are given an initial condition; let’s deliberately select something that doesn’t look like a mode, such as a disturbance of m 1 on its own: z 1 (0) = 1 z 2 (0) = 0 z 3 (0) = 0 ż 1 (0) = 0 ż 2 (0) = 0 ż 3 (0) = 0 Then, as in (5.19), we have V −1 z(0), ż(0) y(0) = V −1 z(0) Eig. ẏ(0) = V −1 ż(0) Easy (5.71) !√ " V yi (t) = yi (0) cos!√ λi t" + z(t) = V y(t) Eig. ẏi (0) sin λi t The initial conditions determine the constant scalar weights yi , ẏi ; in our particular case         1 1/2 0 0 √ y(0) = V −1 0 = 1/ 2 ẏ(0) = V −1 0 = 0. (5.72) 0 1/2 0 0 Our general solution for z(t) can be written analytically, as in (5.22), consisting of three sinusoids at different frequencies controlled by λi , each multiplying a fixed mode pattern v i : Case Study 5: System Decoupling 91    √  1/2 1/ 2 !√ " √ !√ " z(t) = y1 · cos 0.6t · / 2 + y2 · cos 2t · 0   1   1/2 1/√2 -./0 -./ 0 -./ 0 -./0 -./ 0 -./ 0 weight dynamic shape weight dynamic shape  1/2 !√ " √ + y3 · cos 3.4t · −1/ 2 (5.73) 1/2 -./0 -./ 0 -./ 0 weight dynamic shape The eigendecomposition only needs to be taken once; changing the initial conditions only changes the relative contributions (weights yi ) of the different modes, but does not actually in any way change the mode shapes themselves. As was mentioned on page 67 at the beginning of this chapter, although the individual cosine components are clearly √ periodic, if the relative frequencies of the components are an irrational ratio, such as 2/0.6, then z(t) itself will never actually repeat. The principles explored in this case study, including this particular spring-mass exam- ple, are explored in significantly more detail in Appendix C. Further Reading The references may be found at the end of each chapter. Also note that the textbook further reading page maintains updated references and links. Wikipedia Links: Linear system, Linear Time Invariant Systems Theory, Eigendecomposition, Jordan normal form, Control Theory The reader interested in principles of system decoupling applied to a larger variety of examples should refer to Appendix C. There are a great many excellent texts on linear dynamic systems, particularly those of Gajic and Scheinerman. The field of control theory is vast, and there are many undergraduate texts to consider as next steps, two of which are the books by Friedland and Nise. The topic of non-normal systems is relatively specialized. The book by Trefethen is really the definitive treatment of the subject, but at an advanced level. If the reader is interested, Problems 5.12 and 5.13 develop concrete numerical examples to try. 92 5 Linear Systems Sample Problems (Partial Solutions in Appendix E) Problem 5.1: Short Questions (a) For a given n × n matrix A, how do we define the eigenvalues and eigenvec- tors of A? (b) What do we mean by the mode of a linear dynamic system? (c) What does it mean for a given n × n dynamics matrix A to not be diagonal- izable? Problem 5.2: Analytical—Phase Plot For each of the following systems the associated eigendecomposition is given: 2 3 2 32 3 2 3 2 3 2 3 ẋ 02 x 3 -1 1 λ1 = −2 = + v1 = v2 = ẏ 20 y -1 1 1 λ2 = 2 2 3 2 32 3 2 3 2 3 2 3 ẋ 21 x 2 1 1 λ1 = 2 = + v1 = v2 = ẏ 03 y 3 0 1 λ2 = 3 2 3 2 32 3 2 3 2 3 2 3 ẋ -3 4 x 3 1 -1 = + v= ± i λ = 1 ± 4i ẏ -8 5 y 8 2 0 For each system, (a) Draw the phase plot. (b) What sort of dynamic behaviour does this system possess? Problem 5.3: Analytical—Linear Dynamics The step from (5.33) to (5.34) was just stated, without further argument. Derive, mathematically, that ż = α(z − z̄) z(t) = z̄ + eαt · (z(0) − z̄). (5.74) Problem 5.4: Analytical—Linear Systems A centre is very easy to construct in polar coordinates: ṙ = 0 θ̇ = 1 (5.75) That is, indefinite orbiting behaviour (θ̇ = 1) with no change in radius (ṙ = 0). Sample Problems (Partial Solutions in Appendix E) 93 (a) Derive the dynamics for ẋ, ẏ in Cartesian coordinates. (b) Find the system matrix A corresponding to state [x y]T (c) Either numerically or by hand, find the eigendecomposition of A to validate that a centre dynamic is present. Problem 5.5: Analytical—System Analysis We are given a system and its corresponding eigendecomposition: 2 3 2 32 3 2 3 2 3 ẋ -3 1 x 1 1 λ1 = −4 = v1 = v2 = (5.76) ẏ 1 -3 y -1 1 λ2 = −2 which we now wish to analyze in three different ways: (a) Graphically: Sketch the phase plot. (b) Analytically: Starting from! state (xo ," yo ), use the strategy of (5.19) to ana- lytically find the solution x(t), y(t). (c) Numerically: Apply a simple Forward–Euler discretization (see (8.22)) to (5.76) and simulate the system, starting from several initial points, such as ! " # $ xo , yo = 5 · cos(kπ/4), 5 · sin(kπ/4) 0 ≤ k < 8. Problem 5.6: Analytical—Diagonalization Suppose that A, B are both n × n dynamics matrices, ż = Az(t) ẏ = B y(t) where B is diagonalizable but A is not. Briefly describe the qualitative (concep- tual, high level) differences between the modes of z and those of y. Problem 5.7: Analytical—Discrete Time Let A be the dynamic matrix for a discrete-time linear system z(t + 1) = Az(t). Then we know that z(t) = At z(0). Similar to the continuous-time case discussed in the text, how does having the eigendecomposition of A simplify the compu- tation of At ? 94 5 Linear Systems Problem 5.8: Conceptual—Jordan Form The eigendecomposition tries to decouple the dynamics of a system. Suppose that the eigendecomposition is not able to diagonalize a 6 × 6 matrix A, such that we are left with the partially-decoupled system λ1 λ1 λ1 λ2 λ3 λ3 (a) Give one possible Jordan-form matrix corresponding to this system. (b) Sketch typical impulse responses for each of the three subsystems. Problem 5.9: Analytical—Eigendecompositions In this text we are generally not concerning ourselves too much with proving things mathematically. However there are a few quite elegant results that are relatively easy to derive: (a) For real symmetric matrix A, prove that all of its eigenvalues must be real. As a result, real-symmetric dynamics can never give rise to oscillating / spiral dynamics. (b) For real symmetric matrix A, prove that all of its eigenvectors must be orthog- onal to one another. As a result, real-symmetric dynamics can never give rise to the non-normal dynamics of Section 5.6. Problem 5.10: Numerical/Computational—Eigendecompositions Let’s investigate eigendecompositions numerically. In some numerical mathe- matics program (Matlab, Octave) answer the following: (a) Suppose we create random matrices A + AT A = randn(20); B= 2 Both A and B are random, however B is symmetric whereas A is not. Com- pute the eigenvalues of A and B; what do you observe? (b) Are most matrices diagonalizable? Take a random matrix A = randn(20); find its eigendecomposition, apply the eigenvectors to compute D, the diago- nalization of A, and let dmax be the largest off-diagonal element in D. Perform this test 10000 times and plot the distribution of dmax. Sample Problems (Partial Solutions in Appendix E) 95 (c) It turns out that non-diagonalizability is a degenerate property, meaning that if A is not diagonalizable, A + & A almost certainly is, for nearly any infinites- imal perturbation & A. Look carefully at the Jordan form in (5.30) and also at the eigenvalues you generated in part (a). Develop a strong argument why nearly all matrices must be diagonalizable. Problem 5.11: Numerical/Computational—Non-normal Systems Non-normal systems are not particularly mysterious. Any system that has momentum or inertia, for example, will continue to diverge for a period of time, given an initial kick. Specifically, suppose we are given a point mass, having position x(t) and velocity v(t) in the state 2 3 x(t) z(t) =. v(t) The mass is subject to a weak return force, pulling its position to zero, and is also subject to friction, which pulls its velocity to zero. The dynamics could look something like ẋ ≡ v 2 3 2 32 3 ẋ(t) 0 1 x(t) ż(t) = = v̇(t) -0.01 -0.02 v(t) Return Friction For the given dynamics, (a) Compute the eigendecomposition and demonstrate that the system is (i) Stable, (ii) Oscillatory, and (iii) Non-Normal. (b) Determine that the system can exhibit a transient divergence by computing Reactivit y(A) of (5.59). 96 5 Linear Systems Problem 5.12: Numerical/Computational—Non-normal Systems It is relatively easy to construct stable non-normal systems if we directly specify the eigendecomposition. Let 2 3 2 3 cos(θ) cos(−θ) v1 = v2 = 0 < θ < π/2 sin(θ) sin(−θ) then for θ = π/4 the two eigenvectors are normal (orthogonal), and for other values of θ the eigenvectors are non-normal. Aside from θ, the other aspect of our system which affects the dynamics are the system eigenvalues. The main question concerns the relative values of the eigenvalues, so without any loss of generality we will set λ1 = −1, and then test values λ2 < 0. Since both λ1 , λ2 are negative, the system is clearly stable. 2 3 1 Suppose we always initialize our system as z(0) = 0 As in (5.60), let us measure the degree of amplification κmax in the trajectory z(t) as the furthest extent from the origin: = κmax = max z 1 (t)2 + z 2 (t)2 t where z 1 (t) and z 2 (t) are evaluated at discrete points in time t using (5.42). Numerically evaluate κmax (θ, λ2 ) for a range of values for θ and λ2 : (a) Produce a contour plot for κmax , plotting the contour which separates κmax = 1 (no transient growth) from κmax > 1 (non-normal growth in the transient). (b) Discuss your observations. What are the circumstances which lead to signif- icant transient behaviour? Problem 5.13: Numerical/Computational—Non-normal Systems The paper by Neubert and Caswell has a nice example of a non-normal ecological system, for which a simplified version of the dynamics A appears in Table 5.1. We begin with a look at the values in A: (a) Why are all of the diagonal elements negative? Why should the presence of leaves imply a negative growth rate for leaves, for example? (b) Note that most of the columns of A sum to exactly zero. What is the ecological explanation of this observation? (c) What is the ecological significance of the choice of one column which does not sum to zero? Sample Problems (Partial Solutions in Appendix E) 97 Table 5.1 The dynamics matrix A for Problem 5.13, describing the dynamic ecology of nutrients in a tropical rainforest, modified from. Influence of … Leaves Stems Litter Soil Roots Fruits Detr. Herb. Carn. Leaves -1.56 0.669 0 0 0 0 0 0 0 …on the rate of growth of Stems 0 -0.712 0 0 2.56 0 0 0 0 Litter 1.46 0.036 -6.41 0 0 1.14 0 55.8 17.3 Soil 0 0 0 -0.022 0 0 316 0 0 Roots 0 0 0 0.020 -2.56 0 0 0 0 Fruits 0 0.007 0 0 0 -2.03 0 0 0 Detritivores 0 0 6.41 0 0 0 -316 0 0 Herbivores 0.10 0 0 0 0 0.89 0 -62.6 0 Carnivores 0 0 0 0 0 0 0 6.8 -17.3 (d) Following on the two preceding questions (b,c), prove that if all of the columns of A sum to zero the system cannot be stable. Next, we want to understand the resilience of this system, per the discussion in Example 5.3, in terms of eigendecompositions: (e) Find the eigenvalues of A and determine that the system is stable. From these eigenvalues, compute Resilience(A) as in (5.58). (f) Find the eigenvectors of A and determine that the system is non-normal. (g) Compute Reactivit y(A) in (5.59) from the eigenvalues of the symmetric sum (A + A T )/2. If you examine all of the eigenvalues you will observe that there are two positive ones, meaning that the non-normal transient may have a combined behaviour of two time constants. Problem 5.14: Reading Given that linear systems are characterized by exponentials, it is worth better understanding the properties of the exponential function and the connections to human thinking and behaviour. Read or watch5 one of Chapter 5: Dangerous Exponentials in Crash Course Chapter 3: Exponential Growth by Chris Martenson Arithmetic, Population, and Energy by Albert Bartlett 5 The book and videos may be found online; links are available from the textbook reading questions page. 98 5 Linear Systems and, building on Example 5.1, prepare a half-page discussion summarizing the key attributes of exponential functions and why they pose such a difficulty for people to understand. Problem 5.15: Policy As discussed in Example 5.1, we hear regularly about references to modest eco- nomic growth or monetary inflation of 1% or 2%, however any fixed percentage growth rate, regardless how modest, is still exponential over time. Following up on Problem 5.14, discuss the public policy implications of a political system which appears to be premised, or at least relying upon, indefinite growth. References 1. B. Friedland, Control System Design: An Introduction to State-Space Methods (Dover, 2005) 2. Z. Gajic, Linear Dynamic Systems and Signals (Prentice Hall, 2003) 3. L. Gunderson, C. Holling, Panarchy: Understanding Transformations in Human and Natural Systems (Island Press, 2012) 4. C. Martenson, The Crash Course: The Unsustainable Future of our Economy, Energy, and Environment (Wiley, 2011) 5. D. Meadows, J. Randers, and D. Meadows. Limits to Growth: The 30-Year Update. Chelsea Green, 2004 6. M. Neubert, H. Caswell, Alternatives to resilience for measuring the responses of ecological systems to perturbations. Ecology 78(3), (1997) 7. N. Nise, Control Systems Engineering (Wiley, 2015) 8. A. Oppenheim, A. Willsky, H. Nawab, Signals & Systems (Prentice Hall, 1997) 9. M. Scheffer, Critical transitions in nature and society (Princeton University Press, 2009) 10. E. Scheinerman, Invitation to Dynamical Systems (Prentice Hall, 1995) 11. L. Trefethen, M. Embree, Spectra and Pseudospectra: The Behaviour of Nonnormal Matrices and Operators (Princeton, 2005) 12. B. Walker, D. Salt, Resilience Thinking: Sustaining Ecosystems and People in a Changing World (Island Press, 2006)