NONLINEAR DYNAMICS OF T H E INVERTED P E N D U L U M UNDER LINEAR-FEEDBACK STABILIZATION By Michael G. Henders Bachelors of Science, Electrical Engineering University of Manitoba A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R S O F A P P L I E D S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S E L E C T R I C A L E N G I N E E R I N G We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A September 1991 © Michael G. Henders, 1991 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Electrical Engineering The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: Abstract Analysis is conducted on a Unear control system used for the stabilization of an inverted pendulum, examining the behaviour of the system throughout its global state space. The dynamic equations of the inverted pendulum are obtained, and a linear state-feedback controller is designed, based on this system model, for its stabilization. The behaviour of the resulting closed-loop system is examined through simulation, in two stages. In the first stage, the feedback controller is modified slightly to decouple the pendulum dynamics from those of its supporting cart, permitting examination of the pendulum dynamics in a two-dimensional state space, using phase plane analysis. In this way, a globally valid description of the pendulum dynamics is obtained, revealing the existence of two state attractors in this modified system. The feedback controller is seen to create a stable point attractor in the phase plane, as intended, but the basin of attraction of this point is of unusual shape and limited extent, due to the nonlinearity of the pendulum dynamics. All initial system states outside this basin are drawn to the second state attractor, which is a stable limit cycle. As the second stage of analysis, the system is examined with its original controller; characterization of the overall dynamics then requires a four-dimensional state space. The continued existence of a point attractor at the state space origin is demonstrated, and the basin of attraction of this point is mapped numerically, using a state space sampling algorithm, which might readily be applied to the analysis of other dynamic systems. This basin is presented as a series of three-dimensional sections of the full four-dimensional region; taken together, these sections outline the entire basin map. The boundaries of the origin-bound basin are seen to be highly irregular, and potentially fractal in nature. ii The existence of an infinite set of such four-dimensional stability regions is demonstrated, and complex, potentially chaotic, behaviour is observed outside the basins of attraction; it is not clear whether or not a four-dimensional counterpart of the limit cycle (from the first stage analysis) exists. in Table of Contents Abstract ii List of Figures vii Acknowledgement ix 1 I N T R O D U C T I O N 1 2 Dynamic Equations And Linear Control 5 2.1 System Equations 6 2.2 Linear Control 7 3 Singularity Analysis: Full-State-Feedback Control 10 3.1 Singularity Location 10 3.2 Linearization 'In-The-Small' 11 3.3 Eigenvalue Analysis 13 4 Singularity Analysis: Reduced-State-Feedback Control 16 4.1 Singularity Location 16 4.2 Linearization 'In-The-Small' 18 4.3 Eigenvalue Analysis 19 5 Global Behaviour: Reduced-State-Feedback Control 22 5.1 Pendulum Dynamics 23 5.2 Cart Dynamics 33 iv 6 State Trajectories: Full-State-Feedback Control 38 6.1 Behaviour Near Singular Points 38 6.2 Unusual Behaviour At Non-Singular Points 43 7 Stability Regions: Full-State-Feedback Control 51 7.1 Algorithm And Software 51 7.2 Discussion Of Results 60 8 Conclusion 80 8.1 Extending The Work 80 8.2 Summary Of Results 81 Bibliography 83 A P P E N D I C E S A Dynamic Equations Derivation 85 B State Feedback Controller Design 89 C Partial Singularity List For Reduced-State Feedback 93 D Linearization Near A General Singularity 94 D.l General Procedure 94 D.2 Linearization Of The x Expression 95 D.3 Linearization Of The 6 Expression 99 D.4 Completing The Linearization 102 E About The Figures 105 v F Simulation Software 107 F . l Using The Program 110 F.2 Ancillary Files 113 F.2.1 Project Control File: SIM-CORE.PRJ 113 F.2.2 Compiler Identification: COMPILER.H 113 F.2.3 Simulation Definitions: SIMUL.H 113 F. 2.4 Keyboard Definitions: KEYSTRKS.H 115 F.3 Driver: SIM-CORE.C 116 F.4 Numerical Integration: NTGRPROC.C 123 F.5 System Definition: DYNAMICS.C 127 F.6 User Interface: DISPPROC.C 133 F. 7 File I/O: FILEPROC.C 149 G Stab i l i t y Region Computation Software 153 . G. l Ancillary Files 156 G. l . l Project Control File: STABVOL.PRJ 156 G.1.2 Volumetric Definitions: VOLUME.H 156 G.2 Driver: STABVOL.C 157 G. 3 Sampling: SMPLPROC.C 165 H Data Manipulation Software 175 H. l 3-D Sectioning: SLIC30F4.C 176 H.2 2-D Sectioning: SLIC20F3.C 184 H.3 3-D Format Conversion: SCATTER3.C 199 vi List of Figures 1.1 Mechanical System: General Arrangement 1 5.1 Origin Region, For 7 = 0.102 24 5.2 Origin Detail, For 7 = 0.102 25 5.3 Origin Region, For 7 = 0.255 27 5.4 Origin Detail, For 7 = 0.255 28 5.5 Origin Attractor And Limit Cycle, For 7 = 0.102 30 5.6 Origin Attractor And Limit Cycle, For 7 = 0.255 31 5.7 Singularity Map With Limit Cycle 32 5.8 Central Unit Of Global Attractor 34 5.9 Detail Of Global Attractor Trajectories 35 5.10 Time History Of Cart Velocity 37 6.1 Origin-Bound State Trajectory: (x, x, 0) Space 40 6.2 Origin-Bound State Trajectory: (x, 0, 9) Space 41 6.3 Odd-n Singularity: Reverse-Time Traces 42 6.4 'Whorl ' In The 9 - 0 State Space 44 6.5 Pseudo-Chaotic Trajectory: Absolute-Sense 0 Value 46 6.6 Pseudo-Chaotic Trajectory: Modulo-360-Degree 0 Value 47 6.7 Pseudo-Chaotic Trajectory: Modulo 0 Value, With x Extension 48 6.8 Difference Of Pseudo-Chaotic Trajectories For 0, Variation 50 vn 7.1 Reference Diagram For STABVOL Algorithm 55 7.2 Origin Stability Region: Section At x, = —150 61 7.3 Origin Stability Region: Section At i:, = -125 62 7.4 Origin Stabibty Region: Section At ii = —100 63 7.5 Origin Stability Region: Section At x; = — 75 64 7.6 Origin Stability Region: Section At i , = —50 65 7.7 Origin Stability Region: Section At i t = —25 66 7.8 Origin Stability Region: Section At x,; = 0 67 7.9 Origin Stability Region: Section At x, = 25 68 7.10 Origin Stability Region: Section At x,; — 50 69 7.11 Origin Stability Region: Section At x,; = 75 70 7.12 Origin Stabibty Region: Section At i , = 100 71 7.13 Origin Stability Region: Section At xt; = 125 72 7.14 Origin Stability Region: Section At Xi = 150 73 7.15 Two-Dimensional Stability Section At i , = 0, 0, = —250 78 7.16 Two-Dimensional Stability Section At x; = 0, 8,\ = 0 79 E . l Sampling Point Representation Using GRAFTOOL 105 viii Acknowledgement I would like to gratefully acknowledge the following people for their contributions to my thesis research: My wife, Cindy, for her patience and support; without her, quite literally, none of this would have happened . . . My advisor, Dr. Avrum Soudack, for his continued encouragement and sup-port throughout the duration of this project. It is greatly appreciated. Dr. Paul Thibault and Mr. David Weiler, for the use of the facilities at Com-bustion Dynamics, and for their patience with my erratic schedule while com-pleting this thesis. The computations required for the work presented here would not have been feasible without the equipment that they provided. This thesis is dedicated to my father, Lloyd Henders, and to the memory of my late mother, Leona Henders (Falk). ix Chapter 1 I N T R O D U C T I O N This thesis reports the results of an investigation into the dynamics of a simple mechanical system, the classic 'inverted pendulum' (Figure 1.1). This system is sometimes used as a tutorial problem in control theory, due to its inherently unstable, highly nonlinear dynamics [Chen] [Dorf]; [Dorf] suggests that the system represents a simple model of a rocket booster on lift-off, although cursory examination of related literature shows this to be a rather extreme simplification of the models actually used by rocket dynamicists. Different aspects of the inverted pendulum have been studied by a number of other researchers, in varying context. [Feng] [Hemami] [Mori] [Watts] T ) Cr X Figure 1.1: Mechanical System: General Arrangement 1 Chapter 1. INTRODUCTION 2 Feng and Yamafuji [Feng] examined a variation on the configuration of Figure 1.1, using a fixed base-pivot for the pendulum, and controlling the pendulum by torquing a reaction arm attached to its upper end; this configuration is based on modelling the balancing of a robotic unicycle. Hemami and Katbab [Hemami] used an inverted pendu-lum model to describe the dynamics of the upper body; in this case, a three-dimensional pendulum was considered, and its angular excursion was limited. Watts [Watts] looked at the basic control problem, as given by Chen or Dorf, but applied optimal regula-tor theory to the design of a controller, and considered several nonlinear effects besides those inherent in the model geometry. Mori, Nishihara, and Furuta [Mori] performed an extensive investigation into the behaviour of the inverted pendulum, which included the construction of an actual mechanical system for testing. Their study permitted the pendulum to hang down from its pivot (at 6 = 180 degrees; see Figure 1.1), and imple-mented a controller capable of 'swinging up' the pendulum into a state where a normal regulator could stabilize it. Surprisingly, there seems to have been little work done on the most basic implemen-tation of control of the inverted pendulum, namely, that accomplished using a simple pole-placement design for a linear controller. Although this system is treated in text-books, discussion of its dynamics seems to be limited to what occurs within a state-space region where linearization is applicable; there is generally no mention of the system be-haviour outside this region, except perhaps to suggest that 'the pendulum falls over'. The work described in this thesis originated with a curiosity about what actually does happen to the inverted pendulum system outside the linear region of its state space, which is relatively small. The system is considered as in [Dorf], but the pendulum itself is assumed to be capable of unlimited rotation (ie, no limits on \9\ ), similar to the implementation of Mori et al, and the excursion of the cart is also unlimited. Control of the pendulum system is accomplished using a linear feedback controller, designed through Chapter 1. INTRODUCTION 3 pole-placement. Two aspects of the controlled system are examined, both of which are concerned particularly with the 'in-the-large' behaviour of the system; ie, its behaviour in those portions of its state space where the controller design linearizations are no longer valid. The earlier work, which has been published in [Hend], examines the dynamics of the system while operating under certain simplifying assumptions; this work makes extensive use of phase plane analysis to characterize the simplified system dynamics. The more recent work comprises a mapping of the stability regions of the controlled system in four-dimensional state space, achieved through computational techniques; these techniques may be found useful in computing stability regions for other dynamic systems. This report comprises a unified presentation of all thesis work conducted on the in-verted pendulum system. The details of the fundamental control problem are presented, deriving the dynamic equations of the physical system and detailing the design of a suit-able stabilization system for it; simplification of the control system is also described, with reasons for such modification. The singular points of the closed-loop dynamic equations are analyzed, in both original and simplified form, as a starting point for investigation of the system dynamics; it is shown that both forms of the equations contain an infinite number of singular points. The results of the work in [Hend] are presented next, giving a globally applicable analysis of the dynamics of the simplified system, and showing that its state is ultimately drawn to one of a pair of attractors; the first of these is the state space origin, and the second is a stable limit cycle. Complex and unexpected state trajectory behaviour is observed. This state space work is then repeated for the unmodified system; although much less comprehensive conclusions are possible in this case, the existence of convoluted, possibly chaotic, state trajectories is demonstrated. This limited study is supplemented by the computation of four-space stability regions for the system, show-ing that these regions are quite extensive, with complex, potentially fractal, boundaries; Chapter 1. INTRODUCTION 4 furthermore, it is shown that the controlled system possesses an infinite number of such regions. Finally, the investigation is summarized, with some general comments on the results, including the suggestion of potential extensions to the work presented. Chapter 2 Dynamic Equations And Linear Control This chapter presents a short-form development of the differential equations describ-ing the dynamics of the inverted pendulum system, and then proceeds with the pole-placement design of a linear state-feedback controller for the stabilization of the system. Full details of these derivations may be found in Appendices A and B. Throughout the reported investigations, the following assumptions are made regarding the system implementation: 1. The system is dissipationless; there is no friction. 2. Unrestricted rotation of the pendulum is permitted. 3. The rotation angle is known in an absolute sense, and is used in this sense by the pendulum control system; ie, 360 degrees rotation is distinct from 0 degrees, or from 720 degrees. 4. Unrestricted lateral motion of the cart is permitted. 5. Unlimited control force is available; no saturation occurs. 6. Where numerical values are required for the system parameters, the following values are used: m = 0.1 kg, M = 1.0 kg, and r = 1.0 m. This is referred to as the 'reference parameter set'. The value of 9.80665 m/s/s is used for gravitational acceleration, g. (See Figure 1.1.) 5 Chapter 2. Dynamic Equations And Linear Control 6 2.1 System Equations The dynamic equations of the inverted pendulum are developed using the technique of Lagrangian analysis, non-dimensionalizing the system masses; the following variables are defined: K = system kinetic energy P = system potential energy L = Lagrangian; L = K — P V = generalized coordinate U = generalized force Lagrange's equation states that ^(dL/dV^J — dL/d*V = U . For the inverted pendu-lum system, the mass-normalized kinetic energy is 2 2 — (x + r sin 9) dt + P ^(rcosfl) (2.1) where p = m/M, x = dx/dt; the normalized potential energy is given by P = figrcosO. Combining this expression with equation 2.1 and taking the indicated time differentials produces the normalized Lagrangian: * X2 ' LL / • \ 2 L = (1 -f n) — + urx9 cos 9 + — (r9) — pgr cos 9 (2.2) The Lagrangian differentiation operations, when complete, result in (1 + n)x + fir (0cos0 - 02sin0) = F x cos 9 + r8 — g sin 9 = 0 (2.3) where F = F/M is the mass-normalized value of the control force. Equations 2.3 are solved for separated expressions for x and 9, producing the differential equations of the Chapter 2. Dynamic Equations And Linear Control 7 inverted pendulum system: pr92 sin 9 + F — pg sin 0 cos 9 x = 1 + p sin2 9 9 = • A g (1 + p) sin 9 — pr92 sin 9 cos 9 — F cos # r (l + p sin2 #) (2.4) In the interests of avoiding possible confusion, it is worth noting that the dynamic equations of the system as developed in [Chen] are incorrect; the solution presented here matches [Dorf], and has been carefully verified on a repeated basis. 2.2 Linear Control The pendulum control system is designed to stabilize the inverted pendulum in the vertical position, while simultaneously bringing the supporting cart to the origin of the x-axis. The result of the control system design is a set of state-variable feedback gains for the controller, which are also used to obtain expressions for the nonlinear differential equations of the closedrloop system, for simulation purposes. To design the controller, the system D. E. s (Equation 2.4) are linearized at the desired r . iT operating point; ie, at x x 9 9 =0, resulting in the following linear system: X 0 1 0 0 X 0 X 0 0 -M 0 X 1 = + 9 0 0 0 1 9 0 9 0 0 [g(l + ri)/r] 0 9 -1 /r (2.5) The general form of this equation is z = [A + BG] z, where Gz = F, and we wish to find values for the components of the G vector, G = [A B C D). For a pole-placement design, desired eigenvalues are specified for [A -f BG]; for the investigations reported Chapter 2. Dynamic Equations And Linear Control 8 here, a single eigenvalue, of multiplicity four, at ( — 1 + JO) is chosen. The required values for G are determined by expanding the characteristic equation for [A + BG] and matching its coefficients with those of the desired characteristic equation; the resulting solution for G is G = [AB C D] = [7 4 7 r [6 + 7 + (1 + p) h \ 4r (1 + 7)] (2.6) where 7 = r/g, giving the state feedback coefficients required for the stabilization of the inverted pendulum system. Substituting the feedback function Gz for F in Equations 2.4 produces the following closed-loop D. E. s for the system: x firO2 sin 9 + 7 (x + 4i) + r [6 + 7 + (1 + p) /j\ 9 + 4r (1 + 7) 9 - jig sin $ cos 9 1 + n sin2 9 (2.7) .. {±±^sin0-//0 2 sin0cos0- {(x + 4.x) /g + [& + 7 + 0 + 4 (1 + 7 ) cost? 9 — • ——~ 1+ fi sin2 9 (2.8) This result marks the point of divergence between the early work, reported in [Hend], on dynamics of a simplified system, and the later work, dealing with the full system as described by Equations 2.7 and 2.8. On examination of these equations, it is evident that the system states are coupled, which complicates investigation of the system, since characterization of the system dynamics requires the use of a four-dimensional state space. It is difficult to visualize such a space in any meaningful way, resulting in the loss of much of the intuition that is the object of the investigation. To deal with this problem, Chapter 2. Dynamic Equations And Linear Control 9 the work in [Hend] (presented in Chapters 4 and 5) considers the closed-loop system with the modification that the feedback gains A and B in the G matrix are arbitrarily set to zero; this is the 'simplified system' referred to. The described system modifications have two main effects: firstly, control over the cart position, x, is lost, as 'its' feedback is eliminated from the control function, and secondly, the 9 dynamics are decoupled from those of x. It should be noted that the x dynamics are still influenced by those of 9; the decoupling is unidirectional. The change is important, however, because, if one is more interested in the action of the pendulum than in that of the cart, its dynamics may now be characterized using only 0 and 9; this reduces the dimension of the state space to two, and permits the use of all the techniques of phase plane analysis. [Stew] The result is a greatly enhanced intuition concerning the behaviour of the controlled pendulum. Chapter 3 Singularity Analysis: Full-State-Feedback Control It is well known that the nature and location of the singular points of a dynamic system define its behaviour; the eigenvalues of a local linearization of the system near the sin-gular points characterize the 'in-the-small' dynamics, and, through global analysis of the singularity types, the 'in-the-large' dynamics may frequently be inferred as well. This analytical power provides motivation for a singularity analysis of the inverted pendulum system, which forms the contents of this and the next chapter. In this chapter, the analy-sis is carried out for the unmodified control system, while Chapter 4 covers the simplified control case. 3.1 Singularity Location The singular points of a dynamic system are those points at which all of the system derivatives are simultaneously equal to zero. [Stew] To find these points for the inverted pendulum system, it suffices to set the x and 6 equations (2.7 and 2.8) equal to zero, with x and 0 also set equal zero, by definition. The result is the following pair of simultaneous nonlinear equations in xp and 0P, the singularity locations: For x: *yxp + r [6 + 7 + (1 + fi) /1] 0P — fig sin 0P cos#p = 0 (3.1) Fore)': [(1 + fi) h]^0p-{(xp/g) +[6+ -r + (l + fi)h] 0P} cos6p = 0 (3.2) 10 Chapter 3. Singularity Analysis: Full-State-Feedback Control 11 The two equations are quite similar; multiplying Equation 3.2 by r, Equation 3.1 by cos dp, and eliminating the resulting common terms gives r(l + p) . sin Op — jig sin $p cos 0P — 0 9 l + p(l- cos2 0P)] sin Bp = 0 (3-3) Given that p must be > 0 for any physically reasonable system, singularity solutions are obtained at 6P = ±mr radians, where n = 0,1,2, • • • ; substituting this solution back into Equation 3.1 gives the result that xp — — [6 + 7 + (1 + fi) / j] gmr metres. The in-verted pendulum system model possesses an infinite number of singular points, regularly spaced along the line in four-dimensional space defined by x = — [6 + 7 + (1 + ft) /f] 9, x = 0,9 = 0. 3.2 Linearization 'In-The-Small' A local linearization of the system dynamics near a singular point is required to carry out eigenvalue analysis, which in turn, forms the basis for identification of the type of the singular point. The linearization is performed by expanding the system equations in a Taylor series about a general singularity, and then eliminating all nonlinear terms from the expansion. The equations prior to expansion are dx _ fir92 sin 9 - ug sin 9 cos 9 + Ax + Bx + C9 + D9 dx x(l + fi sin2 9) de [(1 + n) 11] sin 9 - fi92 sin 9 cos 9- [(Ax + BX + C9 + D9) /r] cos 9 dB ~ 0( l + ^sin20) Chapter 3. Singularity Analysis: Full-State-Feedback Control 12 where the state feedback coefficients are expressed in condensed form as A, B, C, and D, defined in Equation 2.6. The expansion process is lengthy, though not difficult; Appendix D covers the operation in detail. The final result is a pair of equations that are linear in state deviation variables, identified by a d subscript: dxd Axd + Bxd + 6d pr92p cos 9P + C - pg cos (20p)] + 9d (2fir9p sin 9P + D) ~~ (3-6) dxd xd + p (2xp9d sin 9P cos 9P + xd sin2 0P) d9d ~ [(Axd + Bxd) fr] cos 9P - 9d [(D/r) cos 9P + 2fi9p sin 9P cos 0P) ^ ^ dOd 9d + fi (29d9p sin 9P cos 9P + 9d sin2 0P) 9d {\(AXP + Bxp + C9P + D9P) I r] sin 9P - [(C/ r ) - (1 + p) /7] cos 9P - //02 cos (2flp)} 9d + fi (29d9p sin 0P cos 9P + 0d sin2 0P) In these expressions, the (^ -subscripted variables represent small perturbations from the singularity location, which is defined by those variables carrying a p subscript. Equations 3.6 and 3.7 may be cast as a set of linear first-order differential equations, expressed in matrix form as follows: id ^d 9d 9d xd xd 9d 9d 0 1 + fi sin2 0P A B C — pg cos (20p) D 0 0 1 + fi sin2 9P — (A/r) cos 9P -(B/r) cos9p [{Axp + C9p) /r] sin0 p-[(C/r) - ( l + M)/7l COS 9p -(D/r) cos 6P (3.8) Chapter 3. Singularity Analysis: Full-State-Feedback Control 13 and this form may be further reduced to Xd Xd 9d 9d Xd id 9d Od 0 7 0 (-l)"+1/<7 1 4 7 0 4.(-l)n+1/g 0 [<7 + r(6 + 7)] 0 (-l)"+1(6 + 7 ) 0 4 r ( l + 7 ) 1 4(- l )" + 1 ( l + 7) (3.9) using the known locations of the system singularities, as found in Section 3.1. Equa-tion 3.9 describes the response of the system 'in-the-small'; ie, near its singular points. 3.3 Eigenvalue Analysis The exact nature of the system response near each of its singular points depends on the type of the singularity; this categorization is made on the basis of the eigenvalues of the 'in-the-small' linearized system. [Stew] These eigenvalues are obtained by solving det JAI — M T j = 0, where M is the system matrix of the general form of Equation 3.9, z = z T M , and I is a 4 by 4 identity matrix. The resulting determinant is det Al - M r A -1 —7 A — 47 0 0 (-!)"/<? 4(-l)"/<7 0 [g + r(6 + 7)] A 0 - 4 r ( l + 7) -1 ( - l ) » ( 6 + 7 ) A + 4(-l)"(l + 7) (3.10) where n is the singularity number in 0P = ±nn radians, n = 0,1,2, • Chapter 3. Singularity Analysis: Full-State-Feedback Control 14 Expansion of the determinant is accomplished using the method of cofactors [CRC], resulting in the characteristic equation for the linearized system in the vicinity of a general singularity: A4 + 4 [(-1)"(1 + 7) - 7] A3 + [(-1)"(6 + 7) - 7] A2 -f 4(-l)"A + (-1)" = 0 (3.11) Clearly, this equation has two distinct forms: For n even: A4 + 4A3 + 6A2 + 4A + 1 = 0 (3.12) For n odd: A4 - 4(1 + 2 7)A 3 - (6 + 2 7)A 2 - 4A - 1 = 0 For the case of n even, the characteristic equation is identical to that used in designing the feedback control (see Appendix B), so the resulting eigenvalues are those used in the controller design: a single eigenvalue, of multiplicity four, at ( — 1 + jO). These solutions are completely independent of any system parameters. For n odd, characteristic roots depend on the value of 7; in theory, the quartic characteristic equation may be solved symbolically, using a 'resolvant cubic' form [CRC], but in fact, the symbolic solution is too complex to be useful. A numerical root solution [Press], using the 7 value from the reference parameter set (see Chapter 2), yields a pair of real roots, at Ao = —0.4360 and Ai = 5.9716, plus a complex conjugate root pair at A 2, A 3 = —0.3599 ± j0.5045. The four negative real roots of the even-n solution characterize a nodal state attractor in four dimensions; this case includes the state space origin, and confirms that the control system acts according to its design specification, which called for nodal behaviour at the origin. The odd-n eigenvalues are characteristic of a four-dimensional hybrid of a saddle Chapter 3. Singularity Analysis: Full-State-Feedback Control 15 point and a focal point attractor; for these singularities, state trajectories are drawn in directly along one dimension, attracted along spiral paths through two more dimensions, and expelled from the vicinity of the singular point along the fourth dimension of the state space. One particularly interesting point arises from these eigenvalue solutions, namely, that there are two, and only two, sets of eigenvalues which characterize all of the infinite number of singular points found in this system; this is clearly evident in the fact that the eigenvalue solutions are independent of n, the singularity identifier number. This contrasts with the case of reduced-state feedback, discussed next, and has some significant implications, which are discussed in Section 7.2. Unfortunately, there is very little more that can be said about the full-state-feedback case on the basis of its singularity analysis. While it is theoretically possible to infer the complete system dynamics from the singularity analysis, the four-dimensional state space occupied by this system creates profound difficulty in understanding the structure of the separatrices and state flows that must be visualized and integrated to form such a comprehensive view. For this reason, analytical methods for the full-state-feedback system are abandoned at this point, in favor of simulation studies; these investigations are described in Chapters 6 and 7. Chapter 4 Singularity Analysis: Reduced-State-Feedback Control This chapter presents the singularity analysis of the inverted pendulum system, paral-lelling Chapter 3, for the case of modified control feedback, as described in Section 2.2, in which the x and x feedback components are eliminated. 4.1 Singularity Location The singularities of the modified dynamic system are found by solving Equations 3.1 and 3.2 as for the original system, setting the x and x feedback coefficients to zero, since these terms are to be omitted from the simplified feedback function. The resulting equations are For i : r[6 + 7 + (l + p)/7]0p-Ai0sin0pCos0p = O (4.1) For 9: (1 + fx) sm0p - [1 + fi + 7 (6 + 7)] #P COS0p = 0 (4.2) It is worth noting that x does not appear in either of these expressions; one may say either that there are no singularities in the x - x phase space, or alternatively, that any point on the x-axis may act as a singular point, given a singular value for 9. Both expressions specify singular values for 9, and subsequent analysis of the modified system is substantially simplified by concentrating on behaviour in the 9-9 phase space. Singular 16 Chapter 4. Singularity Analysis: Reduced-State-Feedback Control 17 points are found by performing minor reductions on the singularity equations, to obtain For x: sin (28p) = (27///) [6 + 7 + (1 + p) 11] 0P (4.3) For0: tan0p = [l + 7(6 + 7 ) / ( l + M)]^ (4.4) In contrast to the original system case, these modified equations have no direct closed-form solution; it is possible to obtain some general constraints on the singular values, but graphical or numerical solutions are required to obtain the actual singularity locations. By way of constraints, it will be noted that Equation 4.3 specifies the intersection(s) of a sinusoid with a line through the origin; there is only a single such point (at the origin) if the slope of the line exceeds the origin slope of the sinusoid. This is the case if {[1 + fi + 7 (6 + 7)] /fi} > 1, and it is clear that this is always true for any physically reasonable system, for which 7 > 0, fi > 0. Thus, Equation 4.3 specifies one, and only one, singular value. Equation 4.4 specifies a similar intersection problem, involving the tangent function, rather than a sinusoid. Clearly, since the tangent function goes periodically to infinity, these curves exhibit infinitely many intersection points, representing singularities, for any finite line slope; the line slope merely controls the location of the singularities. The spacing of the singular points is not constant, but goes asymptotically to 180 degrees, as \9\ —> 00. There is one potential exception to these statements, within the first ±90 degree excursion from 0 = 0; if the linear slope is less than the tangent curve slope at the origin, no intersections occur in this region, while, for a higher slope, intersections do occur. The origin slope of tan0 is unity, and from Equation 4.4, it is evident that the linear slope must be > 1, given that 7 > 0, p > 0. The case of 'no-intersection-in-this-region' is thus eliminated as a physical impossibility. Chapter 4. Singularity Analysis: Reduced-State-Feedback Control 18 Going beyond these generalities requires the use of numerical methods, which in turn requires consideration of specific parameter sets for the system. This report considers the case of the reference parameter set (see Chapter 2), where m = 0.1 kg, M = 1.0 kg, r = 1.0 m, and g = 9.80665 m/s/s, giving // = 0.100 and 7 = 0.10197; for these parameter values, equations 4.3 and 4.4 reduce to For i : 34.44 9P = sin (29p) For 9: 1.566 6P = tan 9P (4.5) As already discussed, the first of these equations yields a single solution at the origin; the second equation is solved over a range of ±10 revolutions from 9 = 0 using standard numerical root-finding methods [Press], to obtain a list of singularities in this range. The results of this solution are presented as Appendix C. It must be noted that the points listed in Appendix C, other than the phase space origin, are not true singularities, in the sense that the entire system is stationary; in general, 'x is not zero at these points, and thus the supporting cart undergoes acceleration. The absence of x, x feedback, however, means that the motion of the cart has no direct effect on the pendulum dynamics, and thus, the listed points act as singularities in the 9-9 phase space. The pendulum dynamics are conveniently described in the 9-9 phase space by considering these points as 'real' singularities. 4.2 Linearization 'In-The-Small' The local linearization process for the simplified system may be taken as identical to that of the unmodified system, as covered in Section 3.2 and/or Appendix D, with the condi-tion that the reduction of form must stop with Equation 3.8 (Appendix Equation D.38); this restriction is required because the singular points of the reduced system are not lo-cated at multiples of IT radians, like those of the original configuration. Also, the feedback Chapter 4. Singularity Analysis: Reduced-State-Feedback Control 19 coefficients A and B may be set to zero; the linearized result is then given by id id Od Od xd id &d Od 0 0 0 1 + p sin2 0P 0 0 0 C - ug cos (20p) 0 0 D 1+ fi sin2 0P 0 0 (C/r) 0P sin 6p-[(C/r)-(l + fi)h]cosOp -(D/r) cos 0P (4.6) 4.3 Eigenvalue Analysis Singularity eigenvalues are obtained as in Section 3.3, solving det JAI — M T j = 0; the resulting determinant is det Al — M T A + fi sin2 0P) 0 A 0 0 0 0 0 fig cos (20p) — C A - ( C / r ) 0p sin 0P+ [(C/r) - ( l - M / 7 ] cos 0P 0 -D -(l + p sin2 0P) A + (D/r) cos Op (4.7) Substitution of the expanded forms of feedback coefficients C and D, followed by Chapter 4. Singularity Analysis: Reduced-State-Feedback Control 20 expansion of the determinant, results in the characteristic equation for the reduced-feedback system, applicable near a general singular value, 9p: Clearly, this equation gives two eigenvalue solutions at A = 0, indicating neutrally-stable dimensions; these values correspond to the behaviour of the supporting cart. The remaining solutions are obtained by solving the quadratic embedded in Equation 4.8, singular point must be solved individually. Results of such a solution, using the reference parameter set (see Chapter 2), for the singular values computed in Section 4.1, appear in the Appendix C singularity table. Examination of the Appendix C eigenvalue list shows that this system configuration exhibits three classes of singularity; the state space origin acts as a stable focus, while the remaining singular points alternate between saddle points and unstable foci, consider-ing the two system dimensions characterized by the quadratic eigenvalue solution. Also, in contrast to the full-state-feedback case, the eigenvalues are not constant within each singularity class; because the singularity spacing varies, there is a 'phase shift' effect, rel-ative to the trigonometric functions in the singularity equations, resulting in eigenvalues that are different for each individual singularity. The singularity at the state space origin is rather interesting; as noted above, the local system behaviour in its vicinity is that of a focal point, rather than that of a node, although the original system design called for nodal behaviour at the origin. This effect is (4.8) given a singular value 0p. It is readily seen that this quadratic is symmetric about 6P = 0, but a generally applicable solution is not available; the actual eigenvalues of each Chapter 4. Singularity Analysis: Reduced-State-Feedback Control 21 apparently a side effect of the arbitrary modification performed on the designed feedback function, and is perhaps not too surprising. What is interesting, is that at the origin, the characteristic equation, 4.8, reduces to a rather simple form, A2 + 4A(1 + 7) + 6 + 7 = 0; this gives the analytical solution that Ao = —2(1 + 7) ± y/4rf2 + 77 — 2. For the system to be physically meaningful, 7 must be > 0. With this condition, it is easily shown that the state space origin is an unconditional attractor of the system state, relative to 9 and 9. Firstly, for the case where the discriminant of the Ao solution is negative, the eigenvalues are complex conjugate, and it is clear that their real part is negative, indicating a stable origin focus. Secondly, if the discriminant is positive, its value can be no larger than 27, which is the case as 7 —> 00; this forces both of the resulting real eigenvalues to be < 0, corresponding to a stable origin node. There is an important implication to this, namely, that the reduced-feedback control system used in this analysis will always stabilize the inverted pendulum, given an initial state near the origin of the state space, although the details of the local dynamic be-haviour may vary, depending on the value of the parameter 7. The transition from focal to nodal behaviour occurs for 7 = 0.25, and the nature of the origin-local dynamic be-haviour of this system is thus determined by the length of the pendulum strut, assuming g to be fixed. The simple origin structure of the characteristic equation is lost for other singular values; this is unfortunate, since a corresponding loss of insight occurs with regard to the system behaviour near these singularities, due to the intractability of the eigenvalue expression. It is possible that other singular points of the system exhibit type transitions like those seen at the origin, but it is difficult to quantify the conditions that might cause such transitions, and this line of analysis is not pursued any further in this report. Instead, Chapter 5 presents results of simulation work, showing how the singularity analysis done here relates to the global dynamics of the controlled system. Chapter 5 Global Behaviour: Reduced-State-Feedback Control In this chapter, the behaviour of the inverted pendulum system is examined through simulation, and characterized through the use of phase plane diagrams. The simulations use the full nonlinear pendulum model, with the reduced-state-feedback controller, as described in Section 2.2; all simulations are performed using the SIM-CORE program (see Appendix F). This program has been extensively tested through the debugging process inherent in several major software revisions, which effectively provided a cross-check on the overall results, greatly reducing the possibility that 'good' results were actually in error. In addition, several test cases were run on the UBC mainframe, under ACSL [ACSL]; the results were completely comparable to those produced by SIM-CORE, further enhancing confidence in the validity of its solutions. The presentation given here considers first, and most extensively, the behaviour of the pendulum itself, in Section 5.1; the behaviour of the supporting cart is discussed in Section 5.2. It should be noted that none of the Figures in this section, showing phase plane separatrices and/or trajectories, are inferred from analytical results; they are all generated directly from simulations of the controlled system dynamics. Where unstable singularity asymptotes or trajectories are shown, they are generated through reverse-time integration of the system dynamics, which effectively transforms unstable behaviour into stable flow. 22 Chapter 5. Global Behaviour: Reduced-State-Feedback Control 23 5.1 Pendulum Dynamics It is logical to begin investigation of the pendulum dynamics at the origin of the 9-8 phase plane; this is the desired equilibrium state of the control system, and it has al-ready been noted (Section 4.3) that this point is an unconditional attractor for state trajectories. It is not, however, a global attractor. For the reference set of system pa-rameters (given in Chapter 2), the pendulum fails to return to the vertical if its initial static deflection exceeds ±57.541 degrees. The singularity list in Appendix C indicates that this value corresponds to the location of the first non-zero singular points, a pair of saddles; these singularities mark the f?-axis boundary between two distinct modes of system behaviour. The phase plane near the origin is shown in Figure 5.1; the lines shown represent asymptotes of the saddle points, or separatrices, which define regions of behaviour for the system. Trajectories beginning in the area marked 'Stable' are attracted to the origin; ie, the pendulum comes to rest at the vertical. Trajectories beginning outside this area represent those cases for which the control system does not return the pendulum to the vertical; the small arrows on the plot indicate the general direction of state flow. It should be emphasized that for the defined parameter values, the origin is a focal point; thus, the system state oscillates around it, reaching the origin only in infinite time. The enlarged plot in Figure 5.2 demonstrates the way in which trajectories spiral toward the origin; the convergence is quite rapid, indicating that the control system generates a high degree of equivalent damping. Note that the unstable asymptotes identified in this Figure are asymptotes of the saddle points seen in Figure 5.1; they are unstable in the sense that state trajectories in their vicinity are repelled from the saddle points that are their source. Figure 5.3 represents a slight digression; this plot shows the phase plane near the origin o Q c3 •*-> <D .t5 fl o t> >^ frt °^ 4H <L) | a S3 i — i fl <U 600 400 200 0 -200 -400 •600 J I L J I L O r i g i n F o c u s R e g i o n J I I L J I" h I Stable Asymptotes Unstable Asymptotes ' ' I ' ' ' ' I L J L J I I L -100 -75 -50 -25 0 25 50 Pendulum Angle Theta (Degrees) 75 100 - \ O r i g i n F o c u s D e t a i l -- \ - < ^ Focus -\ \ Unstable Asymptotes 0.5 -0.375 -0.25 -0.125 0 0.125 0.25 0.375 0.5 Pendulum Angle Theta (Degrees) Chapter 5. Global Behaviour: Reduced-State-Feedback Control 26 for the case in which parameters are chosen such that the origin is a node, rather than a focus (see Section 4.3). For this case, all other singularities are as listed in Appendix C with regard to type, although their locations are slightly different. The saddle point asymptotes again define operating regions, but one additional set of asymptotes is present: the asymptotes of the node at (0,0). Due to the presence of these asymptotes, there are, in a sense, two types of stability region present in this opposed to a single type for the case of a focus at (0,0). The areas labelled '1-Stable' form one region, in which trajectories approach the origin directly, while, in the remaining areas, labelled '2-Stable', trajectories must reverse their direction, to approach the origin along the stable nodal asymptotes. The physical implication is that, for these parameter values, the control system response may or may not exhibit overshoot, depending on the initial conditions. Regardless, the spiralling of origin-bound trajectories is eliminated, indicating that the system will exhibit (at most) a single-ended overshoot, followed by a monotonic approach to its final state. It may be noted that there are no stable asymptotes identified in the vicinity of the origin in Figure 5.3; this is because the stable origin node asymptotes are effectively coincident with the unstable saddle point asymptotes in this area. Figure 5.4 presents an enlargement of the origin region for the case of a nodal point there, showing the stable asymptotes. The section of the phase plane shown in Figure 5.1 represents the set of states which might be considered if limits were imposed on the angular deviation of the pendulum, but, as such limits have been presumed to be absent for this study, an extended re-gion of the phase plane is now considered. In Figure 5.5, the saddle point asymptotes are integrated to their eventual equilibria. The result is interesting. Firstly, a complex double-spiral region of stability appears, located within the closed figure defined by the stable asymptotes; states in this region have the origin as their ultimate attractor. Sec-ondly, trajectories beginning outside of the origin-bound basin are ultimately drawn to o Q & •*-» 4=1 % o O O n* cn i - H <L> fl PH 750 500 250 0 -250 -500 -750 Stable Asymptotes Unstable Asymptotes O r i g i n N o d e R e g i o n J I I l_ J I I I I I I L J I L J I L J I I L -100 -75 -50 -25 0 25 50 Pendulum Angle Theta (Degrees) 75 100 -0.5 -0.375 -0.25 -0.125 0 0.125 0.25 0.375 0.5 Pendulum Angle Theta (Degrees) Chapter 5. Global Behaviour: Reduced-State-Feedback Control 29 a stable limit cycle, which exhibits an angular excursion of just over one full revolution in either direction from 8 = 0. (Assuming that we begin on the positive-velocity portion of the limit cycle, the pendulum first rotates forward to 8 « 390 degrees, stops, then reverses direction, rotating back to 8 « —390 degrees, stops again, and resumes forward rotation, completing one cycle.) Examination of the time course of angular excursion (not shown) reveals that the limit cycle has a period of 2.40 seconds. A plot corresponding to Figure 5.5, for the case in which the origin is a nodal point, is given in Figure 5.6. It is evident that the phase plane topology for this case is essentially the same as for the case of an origin focus; the change in singularity type at the origin has no effects, beyond those noted previously, relating to overshoot in the vicinity of the origin, and this case is not referred to again. Returning, then, to the case of an origin focus, a map of the system singularities is presented in Figure 5.7, showing the limit cycle as a reference. All trajectories beginning outside the limit cycle in Figure. 5.7 are eventually drawn to the cycle, which acts as a global attractor for system states outside its perimeter. This is an unexpected result, brought about by a complicated inter-weaving of the asymptotes of the repeated saddle point pairs. Global attraction occurs through nesting of a basic repetition unit, consisting of a central attractor, flanked on either side by a saddle point and an unstable focus. For the first such unit, the central attractor is the limit cycle itself; for subsequent units, the outer periphery (if such a thing may be defined) of the previous unit acts as the central attractor. In this case, trajectories are drawn to the central unit, and eventually make their way through its boundary. All such trajectories eventually wind onto the central limit cycle. The first repetition unit of global nesting is detailed in Figure 5.8 and Figure 5.9. The structure at the core of Figure 5.8 is the limit cycle, as seen in Figure 5.5, flanked by saddle points, then unstable foci, as already described. It is easily seen that trajectories 1500 -400 -300 -200 -100 0 100 200 300 400 Pendulum Angle Theta (Degrees) 1200 -400 -300 -200 -100 0 100 200 300 400 Pendulum Angle Theta (Degrees) 1500 o Q +-» <D 4=1 % o 1000 500 00 j=t 0 -500 -1000 •1500 - L i m i l S i n t C y c l e g u l a r i t i a n d es -V 1 r\ v / ex. \/ /z\ v r q v . . . ex. \ / _ Q > s , r\. v \ KJ y\ KJ / \ L Js\ KJ y\ xJ / *\ KJ y\ -i i i , , • , , , , , , • Stable Focus x Saddle Point o Unstable Focus -1200 -800 -400 0 400 800 1200 Pendulum Angle Theta (Degrees) Capter 5. Global Behaviour: Reduced-State-Feedback Control 33 beginning at any point inside the outer closed curve of Figure 5.8 will be drawn to the limit cycle; what is not so obvious, is that trajectories originating outside this curve can pass through it, seeking the limit cycle. Figure 5.9 expands an area near one of the saddle points of Figure 5.8, showing how this is possible. Contrary to appearances, the outer curve is not closed. The unstable asymptote that departs from the bottom left of the left-hand saddle point in Figure 5.8 apparently terminates at the right-hand, detailed saddle point, but this is not the case; it actually passes inside the saddle point, as seen in Figure 5.9, leaving room for incoming trajectories to enter from the outside. The inbound asymptote entering from the top left of Figure 5.9 does not originate from anywhere in the central unit, but rather, comes in from the outside; trajectories approaching it are drawn first to the central unit, then through its boundary, and on to the limit cycle. Subsequent nesting units repeat the basic structure of Figures 5.8 and 5.9, forming a complex basin of attraction centred on the core limit cycle. In summary, the dynamics of the pendulum in the reduced-feedback control system exhibit only two state attractors, as demonstrated by analysis in the 6-0 phase plane, comprising a point attractor and a limit cycle. The phase plane origin attracts all initial states within a region defined by the stable asymptotes of a pair of saddle points, and the limit cycle acts as a global attractor, drawing all other states to itself. 5.2 Cart Dynamics This section comprises a few remarks regarding the behaviour of the cart supporting the inverted pendulum, as observed in the system simulations. The cart exhibits dynamics that are much less interesting than those of the pendulum, described above, although they may be described as being pathological by anyone with a dislike for infinite values. This is probably the most noteworthy feature of the cart dynamics: that the magnitude o Q O u >^ 1—1 a 2000 1500 1000 500 0 -500 -1000 -1500 -2000 -800 -600 -400 -200 0 200 400 600 Pendulum Angle Theta (Degrees) 800 50 o Q +-> 43 •t3 a o o i s . <0 fl 25 0 -25 -50 440 442 444 446 448 Pendulum Angle Theta (Degrees) 450 Chapter 5. Global Behaviour: Reduced-State-Feedback Control 36 of its position coordinate, x, shows a monotonic increase from its initial value, regardless of the initial conditions imposed on the pendulum, so long as the initial cart velocity is zero. The limit of x as time t —» oo is thus undefined, or infinite. The sign of x may go either positive or negative, depending on the overall initial conditions imposed, and the cart velocity, x, may also fluctuate, particularly when the pendulum is operating on its limit cycle, but the velocity does not change sign for such a case. Some typical time-history plots of x are shown in Figure 5.10; the solid-line trace shows the case of an initial conditions vector x, i{ 0, 0, equal to [0.0, 0.0, 57.5, 0.0], while the dashed-line trace is used for the initial state [0.0, 0.0, 57.6, 0.0]. These two ini-tial states bracket the 0 -0 saddle point at 0 = 57.541 degrees, and thereby demonstrate system operation both off and on the 0 limit cycle. It is possible to obtain sign reversals in the cart velocity, if its initial value in non-zero, and happens to have its sign opposite to that of the velocity change imparted by the control system action, but clearly this is a rather artificial case, and cannot be depended on, or even expected to, keep the cart motion bounded. The reason for these difficulties is very simple; the cart moves off to infinity because there is nothing telling it not to. The modifications made to the feedback control system to create the reduced-state-feedback case removed any means of controlling the cart motion, as was noted in Section 2.2; this loss of control represents the price paid for the analytical simplicity obtained for the pendulum dynamics, exploited in Section 5.1. While not detracting from the results obtained by this means, such behaviour on the part of the cart does place the analysis firmly in the realm of the theoretical; a desire for a more potentially 'real' analysis motivates further study of the system under full-state feedback. This work is presented next, comprising Chapters 6 and 7. o • r H <D la OO 150 125 100 75 50 25 0 0 X D o t T i m e H i s t o r i e s J I I I I I I L J I L On Limit Cycle — Off Limit Cycle i i i i i i i i J_ _ l I L 10 15 20 Time (Seconds) 25 30 Chapter 6 State Trajectories: Full-State-Feedback Control The material presented in this chapter comprises a relatively brief look at the behaviour of the inverted pendulum in state space, for the case of full-state-feedback control, as examined through simulation studies. The behaviour of the system in the vicinity of the state space origin is investigated, as is its behaviour near one of the 'odd-n' singularities discussed in Section 3.3; unusual system behaviour in other regions of the state space is also demonstrated. 6.1 Behaviour Near Singular Points Investigation of the system behaviour is initiated at the state space origin, as for the reduced-state-feedback system studies in Chapter 5. This point is a member of the set of 'even-n' singular points for the system, as noted in Section 3.3, meaning that it should exhibit the local behaviour of a four-dimensional nodal point. Obviously, it is difficult to visualize precisely what four-dimensional nodal behaviour should look like, (a common theme throughout this chapter) when viewed in three-dimensional space or two-dimensional projections; the characteristic stable and unstable singularity asymptotes, seen as lines in the Figures of Chapter 5, become four-dimensional attracting and repelling 'sheets' in the state space. Nevertheless, Figures 6.1 and 6.2 show examples of a typical state trajectory through the region near the state space origin. The initial conditions vector for this trajectory is [0.0, 0.0, 51.3, 0.0], which represents a point lying just inside the bounds of stability for 38 Chapter 6. State Trajectories: Full-State-Feedback Control 39 this system; it may be noted that this boundary 6 value is smaller than that found for the reduced-state-feedback case. Apparently, the 0-axis stability bounds are contracted, relative to the reduced-state-feedback case, as a consequence of the addition of x stabi-lization. The Figures show the trajectory itself as a solid line, while projections of the trajectory against the principal state planes are shown as dashed lines; Figure 6.1 plots the trajectory in (x, x, 0) space, ignoring 0, and Figure 6.2 repeats the plot in (x, 0, 0) space, ignoring x. Unfortunately, about all that can be said about this trajectory is that once it gets into the immediate vicinity of the state space origin, it does apparently approach the origin directly, which is consistent with the nodal behaviour expected there, although the gyrations exhibited along the origin-bound path suggest that the singularity asymptote sheets have a rather complex structure in this region. Also, it will be noted that all four system states remain bounded for this case, in contrast to the studies in Chapter 5, in which the supporting cart was seen to move off to infinity in the course of stabilizing the pendulum; under full-state feedback, the cart is returned to x = 0 at the same time as the pendulum is brought to 0 = 0. The situation is somewhat different near the 'odd-n' singular points; Figure 6.3 shows a pair of typical trajectories computed for reverse time-flow from the vicinity of the n = 1 singularity, at [x x 0 9\ = [-520.334, 0.0, 180.0, 0.0]. The left-hand trace in the Figure begins from (x, 0) = (-520.3, 180.0); the right-hand one, from (-520.4, 180.0), and the trajectories are followed through (x, 0, 0) space. It will be noted that near the outer ends of the traces, the x value is increasing very rapidly; in actual fact, the supporting cart experiences an essentially monotonic increase in velocity for these solutions, and quickly departs from the scene. Because these trajectories are the result of a reverse-time solution, they should ac-tually represent incoming state flow, approaching along the stable dimensions of the Figure 6.1: Origin-Bound State Trajectory: (a:, x, 6) Space Figure 6.2: Origin-Bound State Trajectory: (a;, 8, 8) Space Chapter 6. State Trajectories: Full-State-Feedback Control 42 Figure 6.3: Odd-n Singularity: Reverse-Time Traces Chapter 6. State Trajectories: Full-State-Feedback Control 43 singular point; recall that the eigenvalue analysis in Section 3.3 showed that for n odd, the singular points act as a combination of stable focus and saddle point, attracting states along three of their four dimensions. Again, however, the dimensionality of the state space prevents any significant insight into the details of its local structure. Solu-tions run in the forward-time direction from the vicinity of the singular point, like the reverse-time solutions, leave the local space quickly; they subsequently exhibit complex and rather bizarre behaviour, however, which is discussed in Section 6.2. 6.2 Unusual Behaviour At Non-Singular Points This section discusses two different types of unusual effect that have been observed in the state space behaviour of the inverted pendulum system. The first of these is the appear-ance of a 'whorl' in the 9-9 state space; an example of this formation is shown in Fig-ure 6.4. This trajectory resulted as a solution from the initial state [0.0, 0.0, 225.0, 0.0]. The plot is shown in only two dimensions because, like the case of the odd-n sin-gular points, the supporting cart exhibits steadily increasing velocity, and the resulting three-dimensional plot becomes greatly elongated and difficult to interpret. The two-dimensional plot shown, however, exhibits some interesting information; for these initial conditions, the pendulum oscillates over an ever-decreasing angle, at an ever-increasing rate, until shattered by mechanical stress or numerical overflow. Apparently there is some sort of twisting of separatrix hyper-planes in this region which traps the system state, resulting in this type of behaviour. It should be noted that there are no singular points whatsoever in the immediate vicinity of the whorl; it is created by 'in-the-large' interaction of singularity asymptotes in the region. The second type of unusual behaviour noted is a type of highly irregular trajectory, as shown in Figure 6.5; this trajectory is typical of those seen as a forward-time solution 2500 I n i t i a l T h e t a i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 220 230 240 250 260 270 280 290 300 310 320 Pendulum Angle Theta (Degrees) Chapter 6. State Trajectories: Full-State-Feedback Control 45 from the vicinity of odd-n singular points, as well as being the most common type of solution seen for any initial conditions outside the origin-bound basin of attraction. The whorl solution discussed above is rather uncommon; Figure 6.5, beginning from an initial state of [0.0, 0.0, 70.0, 0.0], is much more typical of the unstable solution trajectories. The semi-repetitive form of the Figure 6.5 trajectory immediately suggests that some type of chaotic behaviour is occurring. [Berge] It is not, however, a 'normal' form of chaos, like those discussed in [Berge], because the trajectory, and others like it, does not stay within a bounded region of the state space. It is tempting to plot this trajectory using a modulo-360-degree mapping of 9, as shown in Figure 6.6; in this case, there would seem to be little doubt as to the existence of a chaotic cycle in the system dynamics. The validity of such a plot, however, is questionable, for two reasons. Firstly, although the 9-9 trajectory on the modulo plot appears to remain within a bounded region of the state space, the corresponding trajectory in x - x space is definitely not bounded; Figure 6.7 shows a three-dimensional extension of Figure 6.6 along the x dimension, which continues indefinitely. Extended solution runs also seem to show some expansion of the trajectory space in 9 - 9, although this is not certain, and may be due to accumulated numerical error. Secondly, the fundamental physical basis for all of these plots bes with a control system in which 9 is explicitly not taken in a modulo sense (see the basic research assumptions, in Chapter 2); to suddenly change to a modulo form because it gives a convenient picture, seems somewhat contrived. To expand on this: The feedback controller is explicitly constructed in the form F = Ax + Bx + C9 + D9, as described in Section 2.2; clearly, neglecting the effects of any of the other state variables, F (9) is linearly proportional to 9 itself. This component of the feedback function acts as if it senses the torque exerted by an imaginary torsion spring acting at the pendulum pivot point, and generates a Chapter 6. State Trajectories: Full-State-Feedback Control 46 Figure 6.5: Pseudo-Chaotic Trajectory: Absolute-Sense 6 Value Figure 6.6: Pseudo-Chaotic Trajectory: Modulo-360-Degree 6 Value Chapter 6. State Trajectories: Full-State-Feedback Control 48 Figure 6.7: Pseudo-Chaotic Trajectory: Modulo 9 Value, With x Extension Chapter 6. State Trajectories: Full-State-Feedback Control 49 linear force in proportion to this torque, acting on the supporting cart. Such a torque is not subject to modulo effects; if the spring were a physical entity, it would be under progressively higher stress with each successive revolution of the pendulum, eventually reaching its breaking point. The mathematical fiction has no such limits, and it is difficult to imagine how a modulo function, with its cyclic reset-to-zero behaviour, could be taken as equivalent to a function which increases linearly from 0 to infinity. Conversely, however, if the trajectories are plotted using 9 in its absolute sense, in-terpretation of these results is problematic. The trajectories are complex, pseudo-cyclic, but unbounded, bearing some, but not all, of the usual identifying signs of chaotic be-haviour. These trajectories do exhibit the primary hallmark of chaotic systems, sensitive dependence on initial conditions [Berge], as shown in Figure 6.8. Slight variation in the initial conditions produces trajectories that, while similar in global form, are completely different in evolution; Figure 6.8 shows a plot of the difference between two trajectories having an initial variation of 0.001 degree in the value of 9. It is clear that the solu-tion trajectories are completely divergent. Variation in the error control tolerance of the numerical integrator has a similar effect on trajectories in these state space regions; tra-jectories with identical initial conditions, but slightly different error tolerances on their integration, quickly diverge. The existence of chaotic solutions in the dynamics of this system would certainly not be surprising; chaos is a consequence of the nonlinearity in the system, and of the existence of at least three dimensions in its state space. The inverted pendulum system is highly nonlinear, more so, in fact, than are many of the chaotic systems described in [Berge]. It also possesses the requisite three-plus dimensional state space, but neverthe-less, the interpretation of the results presented in this Section is not at all clear, because of the difficulties described. This aspect of the system dynamics must remain open as a potential subject for further research. Chapter 7 Stabilit y Regions: Full-State-Feedback Control This chapter presents the final set of results obtained for this thesis, comprising a map-ping of the four-space regions within which the linear control system of Section 2.2 is effective in stabilizing the inverted pendulum system. This mapping is achieved through a systematic sampling of initial conditions in the state space, effectively 'growing' a hyper-volume within which the controlled system is known to be stable. The development of the algorithm used for this purpose is discussed first, along with its implementation, followed by presentation and discussion of the stability map itself. 7.1 A l g o r i t h m A n d Software The techniques described in this section arose from the desire to parallel the work in Section 5.1, in which stability regions were obtained under reduced-state feedback, for the case of full-state feedback. It was possible to obtain the Section 5.1 results by theoretical means, aided by simulation, through diagrammatic analysis of the state space topology around the singular points of the system. Plots of the corresponding type of singular-point trajectories for the case of full-state feedback were seen in Chapter 6; it should be evident that the equivalent topological analysis for this case is not feasible. This difficulty led to a progressive series of experiments into the determination of stability bounds for the full-state-feedback system. In the earliest stages, an approximate stability contour in the x - 9 plane (ie, with x, and 0, equal to zero) was obtained by hand, through iterative testing of initial conditions; this contour was found to be roughly 51 Chapter 7. Stability Regions: Full-State-Feedback Control 52 elliptical. Obviously, stabilization must also be feasible for some small non-zero values of ii and/or so it was inferred, without proof, that there existed a nominally ellipsoidal, four-dimensional 'ball' around the state space origin, within which the control system operated as designed. This conclusion naturally led to the consideration of how one might map the actual stability region, to verify that the inference was correct. The 'by hand' technique wias not feasible, due to the large number of samples that would be required. Meanwhile, however, further experimentation with the system simulations produced some results that indicated that the region of stability might be much more complex than expected. These results, in fact, suggested that the stable region, in the three dimensions x, $, and d, might be a sort of helical 'tube', a sort of three-dimensional analogue of the spiral stability region found under reduced-state feedback. The mapping problem then became one of finding a way to locate the boundaries of such a shape in three dimensions; it was decided to ignore x temporarily, with the intention of eventually using it as a parameter, generating a family of three-state contours at different values of the fourth state. It became clear that some form of systematic algorithm would be required to locate the boundaries of such a relatively complex shape; furthermore, it would be necessary for the resulting boundary data be in such a form that it could be presented graphically. It was decided that these goals could be reasonably met by using a slicing algorithm on the stability 'tube', as follows: Beginning from an initial point within the tube, the distance to the stability boundary would be found along each of a series of radial lines in a sampling plane; the resulting set of data points would comprise one data slice. The initial point would then be propagated a short distance along the tube, using negative time-flow integration, and the process repeated, until the entire tube (or a suitable portion of it) had been sampled. The resulting data set would readily lend itself to graphical display as a wireframe plot, or perhaps as a shaded polygon rendering. Chapter 7. Stability Regions: Full-State-Feedback Control 53 A very large amount of time and effort was expended on various adjustments, re-finements, and re-workings of this basic algorithm, but, in the final analysis, it proved to be ineffective, particularly when attempting to map the stability regions for non-zero ii. It did provide confirmation that, at least for = 0, a tubular shape for the three-dimensional stability region, as inferred from earlier simulation runs, was at least approximately correct. Development of the algorithm eventually used (listed in Appendix G) arose from dis-cussions with a colleague, in which he suggested performing a regular grid sampling of the initial state conditions space. This was quickly shown to be infeasible; the hyper-volume of interest was known (or at least suspected) to be rather large, and to contain features having rather high curvature. Achieving a reasonable sampling density, it was anticipated, would require a totally prohibitive amount of computation. (It should per-haps be noted that, when speaking of sampling in this context, each 'sample' actually represents a full simulation run performed on the system equations, beginning from a specified initial state, and continuing until the fate of the pendulum becomes known; ie, does it, or does it not, achieve stability at the vertical? The per-sample computation time is significant.) The suggestion, however, led to the idea that perhaps there might be a way to sample only (or mostly) points within the stability region, stopping on reaching its boundary(s) . The result of this conceptualization, after the addressing of many technical details, is the program STABVOL, listed in Appendix G. It uses the system dynamics and numeri-cal integration routines developed for the simulation program SIM-CORE (Appendix F), and drives them over a systematic sampling pattern, in which the space of initial con-ditions is tested along each state axis until unstable states are encountered. A record of tested stable states is kept, which ultimately becomes the map list for the stability region. This algorithm is essentially a four-dimensional extension of, and variation on, a Chapter 7. Stability Regions: Full-State-Feedback Control 54 'seed-fill' algorithm, well-known in computer graphics circles [Glass], although this fact was not recognized until well after STABVOL was completely designed and operating. Before examining the details of the STABVOL algorithm, a few words regarding the use of the terms 'stable' and 'stability' are in order. These terms are used somewhat loosely in the following discussions, referring variously to several related sets of conditions. In general, 'stable' is used to refer to the inverted pendulum system being in a state such that it is ultimately driven to a position at the x-axis origin, with the pendulum standing vertically upright. In this sense, the terminology effectively refers to a set of initial state conditions for the system; it is in some sense, really the initial state vector which is stable, and the term is also used in this respect, ie, as a 'stable state'. Later discussion considers a case where the pendulum system comes to rest at a state other than the state space origin; this state, and the initial conditions leading to it, are also considered as 'stable', although it is clear that the stability implied is different in some significant ways from that previously discussed. As to the actual determination of stability for the system, or more specifically, for a given initial state vector of the system: This is accomplished through use of a set of heuristic condition tests performed on the system state vector. There is no analyti-cal basis for the tests used, and hence, their accuracy cannot be proven, but extensive experimentation has not demonstrated any exceptions to the implied stability 'rules'. The rules are very simple: Any state vector having 19\ > 400 degrees, or \9 \ > 1500 degrees/second, is considered to be unstable. Any state vector such that the weighted 'state radius' ^(0.020 x)2 + (0.094 xf + (0.020 9)2 + (0.006 #)2 is less than 0.1 is consid-ered to be stable, and all other states are held as indeterminate. These heuristics are encapsulated in the function StablePoint, in Appendix G.3. Within StablePoint, a given initial conditions vector is used as the starting point for a simulation run; the sim-ulation is continued until one of the three heuristic rules is triggered, at which point it is Chapter 7. Stability Regions: Full-State-Feedback Control 55 terminated, with a classification of the initial point passed back as 'Stable' or 'Unstable'. The detailed operation of the STABVOL algorithm is best explained with reference to a diagram; see Figure 7.1, below. The complete operation sequence of the algorithm will be described, as it would occur in processing this example. Q, V7 D3 D2 DO DI V4 V6 VO V3 V2 V5 Figure 7.1: Reference Diagram For STABVOL Algorithm For the purposes of the example, we consider a system with two state variables, P and Q; Figure 7.1 represents the space of initial conditions for this system, and the region enclosed by the dashed line comprises all those initial conditions that are stable, according to the criteria already discussed. Sampling points in the space are indicated by circles, filled if the sample point is stable, and open if not, and four reference directions, DO - D3, are defined in the space. The sample points are grouped together; each group comprises a 'voxel'. 'Voxel' is a contraction of 'volume(tric) element/cell', usually used to denote a three-dimensional entity, but used here in a generalized sense to denote an elemental part of space, in any number of dimensions. Each voxel is identified by a code number, which normally represents its position in the state space, although the voxel Chapter 7. Stability Regions: Full-State-Feedback Control 56 numbers seen in Figure 7.1 have instead been chosen to show the sequence in which the voxels are visited by the STABVOL algorithm. The voxel dimensions along each state axis are set implicitly by the required spatial sampling resolution. The STABVOL algorithm requires a seed point, ie, a state vector that is known to be within the region of stability for the system under consideration. (This point is currently hard-coded into the program; further details may be found in Appendix G.) For the example, we will assume that an initial state vector (0,0) is supplied; examination of Figure 7.1 shows that this state is in fact stable, as required. The algorithm then constructs an initial voxel based at the supplied point by adding offsets to the initial state values in a binary pattern, which results in the generation of state vectors corresponding to the corner points of the voxel; the offsets are set by the size of the sampling grid required. The resulting voxel is placed in a processing list, for subsequent analysis. It will be noted that this type of construction is readily extended to higher dimensions; where Figure 7.1 shows two-dimensional 'area' voxels, true volumes would be created in three dimensions, and hypercubes of appropriate dimensionality would result from further extension. The number of voxel vertices is always 2", where n is the number of dimensions of the state space under consideration. To map the region of system stability, the algorithm indexes through its processing list; for each voxel found in the list, it executes a test that checks each of the voxel vertex states for stability. If none of the vertex states are stable, the voxel is immediately removed from the processing list, but if any of its vertices are stable, the voxel is retained, and a record of its stable vertices is kept. Continuing with the example: There is initially a single voxel, VO, in the list; testing all of its vertices shows that all are stable, so VO is retained in the list, with a code showing that stability information. Immediately after the vertex stability check, the algorithm enters a propagation phase; having analyzed its current operating location, it builds new voxels adjacent to its current Chapter 7. Stability Regions: Full-State-Feedback Control 57 location, along each of the reference directions. This propagation is conditional, however; a new voxel is constructed in a particular direction only if at least one of the vertices on the current voxel face nearest the proposed new voxel is stable. The reason for this is quite simple; if there are no stable vertices on a given face, then the boundary of the stability region must pass through the current voxel, behind the given face, and there is no reason to search further in the indicated direction. In addition, before creating a new voxel, the algorithm checks all existing list entries to see if the proposed voxel already exists. Each successive new voxel is added to the processing list after the last existing voxel. In the example, the algorithm now builds voxels VI, V2, V3, and V4, since all four faces of the initial voxel possess stable vertices, and none of these voxels are yet present in the processing list; these voxels are added in the order indicated, after voxel VO. On completion of the propagation phase, the algorithm returns to the stability testing phase, after incrementing its position in the processing list; the result is that, in the example, voxel VI comes up for consideration. It happens that VI has no stable vertices, so it is removed from the processing list, and the algorithm cycles to voxel V2. This cell also is completely unstable, so it, too, is removed from further consideration, and V3 is loaded. On stability testing, voxel V3 proves to have two stable vertices, so it is retained for processing, and the algorithm re-enters the propagation phase. Propagation is attempted first in direction DO, but the resulting voxel would overlap voxel VO, so it is not created; the bottom face of voxel V3 has one stable vertex, so propagation next occurs along direction DI, creating voxel V5. The V3 face corresponding to direction D2 has no stable vertices, so there is no D2 propagation, and voxel V6 is then created by propagation in direction D3, completing the processing of voxel V3. Voxel V4 is next in line for processing, and its treatment is very similar to that of V3; stable vertices are present, so V4 is retained in the processing list, and its propagation subsequently creates voxel V7. (DO propagation builds V7, DI propagation overlaps VO, Chapter 7. Stability Regions: Full-State-Feedback Control 58 D2 propagation overlaps V6, and D3 propagation is blocked by lack of stable face ver-tices.) On completion of processing for voxel V4, voxel V5 is considered, and immediately removed from the list as being completely unstable; voxel V6 is next, and is retained for its single stable vertex, but generates no new voxels. Finally, voxel V7 is tested, and removed from the list as unstable. Once the complete processing bst has been scanned, as described, the voxels remaining in the list comprise the set of all points lying within the region of system stability. The associated vertex stability codes are also required, since not all vertices of a given voxel are necessarily stable. For the example given, only voxels VO, V3, V4, and V6 remain in the list on completion of processing, and it is clear from examination of Figure 7.1 that these are the only voxels having vertices that fall within the stability region. It is equally clear that not all of the vertices identified by specification of only the voxels, are stable; the vertex stability codes are still needed. This, then, is the basic algorithm used by the program STABVOL; it effectively 'grows' a region of system stability in the space of initial conditions, beginning from a specified seed point, and stopping when the boundaries of the region are reached. Recall that the primary reason for rejection of a simple full-grid-based sampling algo-rithm was the expectation of prohibitive computation time. It quickly became evident, on testing STABVOL, that this would have been the case; STABVOL itself is extremely compute-intensive, as might be expected. Testing and production runs of this program were ultimately performed on a PC-hosted coprocessor board, based on a 40 MHz Intel i860 chip; this chip is capable of 80 MFLOPS, when operating on fully-vectorized code, and achieves approximately three times the speed of a 33 MHz 80386/80387-based PC, when operating in scalar mode. (The board and host system were made available by the principals of Combustion Dynamics Ltd., as noted in the Acknowledgment; their contribution to this work is again gratefully noted.) The code for STABVOL is scalar; Chapter 7. Stability Regions: Full-State-Feedback Control 59 a complete solution run requires approximately 76 hours of execution time on this sys-tem. Examination of the solution set shows that, if a rectangular grid were constructed, encompassing all of the states found to be in the solution set, then approximately 2% of the overall space of initial conditions is actually in the solution set; ie, full-grid sam-pling would require approximately 50 times the computation, to no better result. The 'grow-from-within' algorithm is clearly a superior approach for this problem. It is quite possible that the efficiency of STABVOL might be improved by a redesign based on known seed-fill algorithms, since a great deal of time has been spent on them by a number of researchers, tuning their efficiency. On the other hand, however, little benefit might be derived from such an effort, since the time cost of certain elements of STABVOL is completely different from the time cost of the equivalent elements of a normal seed fill. In particular, the testing of a point for 'inside-ness' is very quick for a graphics fill, but is probably the most computationally expensive operation required by STABVOL. Such tuning efforts would require a carefully considered approach to achieve success. The basic technique used for the STABVOL program is a powerful one, potentially applicable to the study of other dynamic systems. Furthermore, it might be used in ways different from those presented here; as an example, for the inverted pendulum system, one might extend the technique from four-space to eight-space, and map the stability regions of the system dynamics across the space of initial conditions and state feedback coefficients simultaneously. This sort of work represents a potentially fertile area for further research, although the computational cost is likely to be very high. Chapter 7. Stability Regions: Full-State-Feedback Control 60 7.2 Discussion Of Results In this section, the results of a STABVOL solution run on the inverted pendulum dynam-ics are presented. The results comprise a map of those regions within the space of initial conditions of the system for which a simple linear-feedback control system (Section 2.2) stabilizes the pendulum in the vertical position. As previously noted, this map occupies a four-dimensional space, creating certain difficulties in visualization and/or display, which are addressed by the presentation used here. The technique used for presentation is one of sectioning; at each of a set of discrete values of one particular state variable (x is used here), a three-dimensional stablity region is plotted. The set of all such sections, comprising Figures 7.2 to 7.14, taken together and interpolated, represents the complete four-space stability region. The three-dimensional plots themselves are necessarily reduced to two-dimensions for publication; this is ac-complished through the usual projection techniques. [Hearn] [Glass] The projection from three-space to two-space is handled automatically by the GRAFTOOL software ([GRAF], see also Appendix E), used to create these Figures. It is clearly evident that the stability region is quite extensive, and rather complex in structure. Before proceding further, it should also be noted that there is every reason to believe that the stability map presented here shows only one of an infinite set of such regions that exist for this system. This statement has its basis in the eigenvalue analysis from Section 3.3, in which the existence of an infinite set of singular points, identical in local behaviour to the state space origin, was shown. The state space origin lies at the core of the stabibty map shown here (and, in fact, was used as the seed point in its generation); if the local behaviour of other system singularities is identical to that of this point, there must exist a finite region around each of them in which the system is also stable. The implication is clear; each of the even-ra singular points discussed Figure 7.2: Origin Stabil i ty Region: Section A t x, = —150 Figure 7.3: Origin Stability Region: Section At i , = —125 Figure 7.4: Origin Stability Region: Section At ii = —100 Figure 7.5: Origin Stability Region: Section At i , = —75 Chapter 7. Stability Regions: Full-State-Feedback Control 65 Figure 7.6: Origin Stability Region: Section At i ; = -50 Chapter 7. Stability Regions: Full-State-Feedback Control 66 Figure 7.7: Origin Stability Region: Section At i , = —25 Figure 7.8: Origin Stability Region: Section At i , = 0 Figure 7.9: Origin Stability Region: Section At i:,- = 25 Figure 7.10: Origin Stability Region: Section At i ; = 50 Chapter 7. Stability Regions: Full-State-Feedback Control 70 Figure 7.11: Origin Stability Region: Section At x; = 75 Chapter 7. Stability Regions: Full-State-Feedback Control 71 Figure 7.12: Origin Stability Region: Section At i , = 100 Figure 7.13: Origin Stability Region: Section At i : , = 125 Chapter 7. Stability Regions: Full-State-Feedback Control 73 Figure 7.14: Origin Stability Region: Section At ii = 150 Chapter 7. Stability Regions: Full-State-Feedback Control 74 in Section 3.3, of which the origin is one, possesses a surrounding stability region which may be mapped by STABVOL. It is not known, whether the structure of all the resulting stability regions is identical; due to the computing time required for each such a solution, only the origin-centred case has been tested. There is also, perhaps, some justification for treating only this case, since the state space origin is the intended state attractor for the system, as defined during the design of its controller. The Figures actually present two sets of information; the light blue plots represent the actual stability region, in three dimensions, while the red plots are projections of the three-dimensional region, taken against the principal space planes. The combination of these plots makes their interpretation substantially easier. (In the interests of acces-sibility, it should be noted that, when photocopied, the red color, used for the plane projections, becomes a significantly lighter shade of grey than that of the blue three-dimensional plots. This is easy to see on the dense plots, taken near i , = 0, but becomes more difficult to notice as the point density of the plots becomes lower.) The plots shown are not, as one might infer from the appearance of some of them, any sort of shaded rendering; the appearance of solidity is simply an artifact of a higher density of stable points in certain regions of the space. Each point is represented on the plots as a small three-dimensional cross, centred at the sample point itself; this is detailed in Appendix E. The stability region computed by STABVOL is substantially more complex than was anticipated. While its structure is generally that of a helical tube, as originally suspected, at least at i , = 0, such a description is grossly inadequate to describe the actual results. The existence of highly irregular boundaries, forming tendrils and islands in the state space, will be noted; such a structure strongly suggests that the bounds of the stability region may be fractal. [Mand] It is quite possible that the actual degree of complexity is even higher than that seen here; since the stability mapping algorithm used by STABVOL is inherently a sampling process, it is subject to all of the usual problems Chapter 7. Stability Regions: Full-State-Feedback Control 75 relating to Nyquist criterion, bandwidth, and aliasing, in the spatial domain. Should the actual stability region boundaries be fractal, they would, by definition, contain high spatial frequencies, which would be aliased by sampling. There are a number of techniques which might be used to control this effect, or to test for it, including adaptive sampling [Bloom] and spatial filtering [Glass], but each adds significant complication to the overall stability mapping process. It should be noted that, contrary to appearances, the 'islands' are not actually iso-lated from the main body of the stability region. The truth of this statement may be demonstrated in two ways; first, on the basis of the design and operation of the sampling algorithm used in STABVOL, and second, on the basis of an argument through dynam-ics. On the basis of the sampling algorithm, the only way for an island to be mapped would be for the algorithm to step from a stable point on the main body of stability, over a small region of instability, hitting a stable point on the island on its next sample. This is guaranteed by the fact that the algorithm samples in a particular direction only until it encounters an unstable initial point. Such a mis-cue would be impossible to detect, however, because, by definition, there would be no missed sample points between the main body and the island; they would appear to be connected. This, then, is not the source of the 'obviously' isolated regions seen in Figures 7.2 to 7.14. The argument through dynamics is even stronger; consider the following: Given an isolated island of stability, it is clear, by definition, that state trajectories originating within the island have, as their ultimate destination, the state space origin. (This ignores the possibility that their ultimate attractor might be another even-n singular point; these are also stable, but are readily distinguished on the basis of the different end-point state of the trajectory.) Such trajectories are assumed to be continuous; there thus exists at least one series of state vectors defining a continuous path from within the island to within a known region of stability, ie, to the vicinity of the origin. This, however, is a Chapter 7. Stability Regions: Full-State-Feedback Control 76 logical contradiction; if the island is isolated, such a path cannot exist, but if the path does not exist, then points on the island cannot be stable. It is thus demonstrated that any and all points identified as stable by the mapping algorithm are connected parts of one overall stability region. It is logical to question, at this point, exactly what is happening in the plots, if there are really no isolated regions. It seems likely that the appearance of 'islands' is an artifact of the sectioning procedure used in producing these plots. This is not to imply that the technique is faulty; rather, the sectioning of this data set inherently produces an apparent fragmentation, as seen. By way of analogy, one might consider taking sections at different altitudes through a particularly rugged area of topography. Viewed in isolation, the higher-altitude sections would exhibit exactly the type of fragmentation seen here, as data slices were taken closer and closer to the mountain peaks; the inherent continuity of the data set would be lost. This is the reason that, at the beginning of this section, interpolation of the presented data sections was specified; it is necessary to recognize that there is a significant amount of data missing from the picture given. For this reason, it might be preferable to work directly with the full four-dimensional data set, perhaps by performing a projection from four-space into three-space, followed by polygonization and rendering of the resulting three-dimensional data object; such a technique was considered, but eventually rejected as lying beyond the scope of this thesis. The projection operation is a simple extension of that used to go from three-space into two-space, but the polygonization of the point data set is non-trivial [Glass] [Loren] [Bloom], and likely to be computationally expensive, due to the large number of points involved. (The complete four-dimensional data set comprises approximately 215,000 sample points.) Rendering the resulting object, while well-understood [Hearn] [Glass], also promised to be computationally expensive, as the number of polygons, like the number of data points, would be large. Chapter 7. Stability Regions: Full-State-Feedback Control 77 The major advantage of such a technique would be that all of the available data would be presented in one comprehensive picture; it should then be easier to visualize the overall configuration of the stability region, and should also tend to avoid false impressions like that created by the stability 'islands' seen in the section plots. Conversely, however, it would almost certainly be more difficult to see what sort of stability bounds exist within a local region of the state space; even with the three-dimensional plots, this is somewhat difficult. This problem has been addressed by the creation of some data manipulation software, listed in Section H.2, which permits the extraction of two-dimensional sections from any of the three-dimensional data sets. Examples of the resulting data sections are given in Figures 7.15 and 7.16; Figure 7.15 is taken from the x, = 0 data set, in the plane 0, = 250, and Figure 7.16 is taken from the same data set, at 0, = 0. The utility of such sections may perhaps be questioned, but they do provide information that is not readily available from the higher-dimension data set; the ovoid shape of the x, - 0; origin contour, in Figure 7.16, for example, is not easily discerned from the three-dimensional plot in Figure 7.8. Whatever technique is used to present or view the data set, one conclusion becomes quite unavoidable: there is tremendous complexity inherent in the behaviour of this 'simple' control system. o Q cd +-» <D r f l H O fl o o <u co > cd < 2 . H 0 0 fl PH Id • r H •*-» • r H fl 900 700 500 300 100 -100 -300 -500 -700 -900 -----_ -1 1 1 1 I I I I 1 1 1 _ l . . . . I I I I i . i i i i i i i 300 200 100 0 -100 -200 Initial Support Position X (Metres) -300 cd +-> <D 43 H <u 6 a •80, a P H 'cd -400 -300 -200 •100 0 100 200 300 400 •** J L__l I L__L J I I I I L *********** J L * t J L •** J L 300 200 100 0 -100 -200 -300 Initial Support Position X (Metres) Chapter 8 Conclusion 8.1 Extending The Work There is a certain temptation to make this into one of the larger sections of this thesis, not because the work reported here is trivial or incomplete, but because it opens the door to so many different questions that might be investigated. Some of these questions have been noted in the course of the presentation of results; this section expands somewhat on those, and offers some other suggestions for further research. One item noted in the text is the question of chaotic behaviour, or pseudo-chaotic behaviour, in this dynamic system. As noted, the dynamics exhibited possess some, but not all, of the characteristics normally associated with dynamic chaos, and it would be interesting to determine whether or not such chaos is actually present, and if so, why its structure is so unconventional. This question is likely to be intimately related to a more general question, namely, an adequate description of the actual topology of the four-dimensional state space of the system. Should such a description be obtained, it is likely that it would provide the same sort of global insight given by the phase plane work presented in Chapter 5, potentially explaining both the 'chaos' and the 'whorl' behaviours seen in Chapter 6. Another obvious extension to the work presented would be to perform further mapping of stability regions. As noted in the text, it has not been demonstrated whether or not all of the infinite number of stability regions of the inverted pendulum system possess 80 Chapter 8. Conclusion 81 the same structure. One might also work in the direction of refining the information available on the region presented here, through the use of higher spatial resolution or some type of adaptive sampling of the space. Should this particular dynamic system not be of sufficient interest, such research might be directed instead towards similar study of a different set of dynamics. A word of warning might be in order; if changing to the study of a new system, it would be advisable for the researcher to have some reasonable idea of the structure of the stability region to be mapped; the mapping program, already slow, might run for a very long time indeed, should one or more dimensions of the stability region prove to be unbounded! Finally, relating to the study of different systems, one might consider modifications to the inverted pendulum system itself, perhaps in the nature of 'de-idealizing'it through the consideration of more realistic physical effects. Obvious additions to the dynamics include various types of friction and damping, control force saturation, and limited excursions of the mechanical elements. 8.2 S u m m a r y Of R e s u l t s This thesis has reported the results of a study of the dynamics of an inverted pendu-lum, stabilized by a simple linear-feedback control system. The controlled system was examined from two perspectives; firstly, by means of modifications to the control func-tion, insight on the global behaviour of the system was obtained through an essentially geometrical analysis of its state space, and secondly, similar insight was obtained for the unmodified system through numerical techniques. The latter work comprised a mapping of regions of system stability in four-dimensional state space. Unusual results were obtained in both aspects of the study. The modified-system work culminated in a demonstration of the existence of two, and only two, state attractors in Chapter 8. Conclusion 82 the system, comprising a stable point attractor at the state space origin, and a stable limit cycle. It was further demonstrated that the basin of attraction of the origin is limited in extent, and of complex shape, and that all initial states outside this basin are drawn to the limit cycle, frequently following convoluted trajectories to reach it. For the case of the unmodified system, the existence of an origin-bound basin of attraction was again demonstrated; this basin was mapped numerically, showing very complex structure, which suggested that its boundaries might be fractal. It was shown that the basin actually mapped is only one of an infinite number of such regions in the state space of the system. Finally, state trajectories originating outside the stable basins were seen to exhibit extremely convoluted paths, having sensitive dependence on initial conditions, suggesting that such states may evolve through chaotic dynamics. Although the stabilized inverted pendulum is commonly treated as an almost trivial tutorial problem in control theory, the results presented here show the existence of sur-prisingly complex behaviour in its dynamics, demonstrating the unexpected consequences of nonlinearity in its various forms. Bibliography [ACSL] Simulation program: ACSL - Advanced Continuous Simulation Language; ©Mitchell &; Gauthier Associates, Concord, Massachusetts; Mainframe version, 1987. [Berge] Berge, P., Pomeau, Y., and Vidal, C ; Order Within Chaos; Wiley Interscience; 1984. [Bloom] Bloomenthal, J.; Polygonization of implicit surfaces; Computer Aided Geomet-ric Design; Vol. 5, pg. 341-355; 1988. [BOR] C compiler package: Turbo C\ ©Borland International, Inc., 1800 Green Hills Road, P.O. Box 660001, Scotts Valley, CA 95066-0001; Version 2.0, 1988. [Chapra] Chapra, S.C. and Canale, R.P.; Numerical Methods For Engineers - Second Edition; McGraw-Hill Book Company; 1988. [Chen] Chen, O; Linear System Theory and Design; Holt, Rinehart and Winston, Inc.; 1984. [CRC] Beyer, W.H., editor; CRC Standard Mathematical Tables — 25th Edition; CRC Press Inc.; 1979. [Dorf] Dorf, R.O; Modern Control Systems - Third Edition; Addison-Wesley Publishing Company, Inc.; 1980. [Feng] Feng, Q. and Yamafuji, K.; Design and simulation of control systems of an inverted pendulum; Robotica; Vol. 6, No. 3; pg. 235-241; July 1988. [Glass] Glassner, A.S., editor; Graphics Gems; Academic Press, Inc.; 1990. [GRAF] Scientific graphics software package: GRAFTOOL Graphical Analysis System; ©3-D Visions Corporation, 412 S. Pacific Highway, Second Floor, Redondo Beach, CA 90277; Version 3.3, 1990. [Hearn] Hearn, D. and Baker, M.P.; Computer Graphics; Prentice-Hall Inc.; 1986. [Hemami] Hemami, H. and Katbab, A.; Constrained Inverted Pendulum Model For Eval-uating Upright Postural Stability; Trans, of the ASME; Journal of Dynamic Systems, Measurement, and Control; Vol. 104, pg. 343-349; December 1982. 83 Bibliography 84 [Hend] Henders, M.G. and Soudack, A.C.; 'In-the-Large' Behaviour of an Inverted Pen-dulum With Linear Stabilization; International Journal of Non-Linear Mechanics; in press. [Loren] Lorenson, W.E. and Cline, H.E.; Marching Cubes: A High Resolution 3D Surface Construction Algorithm; Computer Graphics; Vol. 21, No. 4; pg. 163-169; July 1987. [Mand] Mandelbrot, B.B.; Fractals: Form, Chance, and Dimension; Freeman Press; 1977. [Mori] Mori, S., Nishihara, H., and Furuta, K.; Control of unstable mechanical system Control of pendulum; International Journal of Control; Vol. 23, No. 5; pg. 673-692; 1976. [NDP] C compiler package: NDP-C860; ©MicroWay, Inc., P.O. Box 79, Kingston, Mas-sachusetts 02364; Version 4.0b, for DOS, 1990-1991. [Press] Press, W.H., Flannery, B.P., Teukolsky, S.A., and Vetterling, W.T.; Numerical Recipes in C - The Art of Scientific Computing; Cambridge University Press; 1988. [Smith] Smith, J.M.; Mathematical Modeling And Digital Simulation For Engineers And Scientists - Second Edition; John Wiley h Sons Ltd.; 1987. [Soud] Soudack, A.C.; Course Notes for UBC-ELEC557: Nonlinear Systems; unpub-lished; 1988. [Stew] Stewart, H.B., and Thompson, J.M.; Nonlinear Dynamics and Chaos; John Wiley k Sons Ltd.; 1986. [Watts] Watts, J.W.; Control of an Inverted Pendulum; American Society For Engineer-ing Education; 92nd ASEE Annual Conference Proceedings (Vol. 2): Engineering Education: Preparation For Life; June 24-28, 1984. Appendix A Dynamic Equations Derivation This Appendix details the development of a set of differential equations describing an idealized implementation of an inverted pendulum system, as sketched in Fig. 1.1. The D. E. s are derived through Lagrangian analysis, non-dimensionalizing the system masses, and ignoring all frictional effects to simplify subsequent analysis. The following variables are defined: K = system kinetic energy V = generalized coordinate P = system potential energy U = generalized force L = Lagrangian; L = K — P Lagrange's equation states that ^ (dL/dV) — dL/dV = U . For the inverted pendu-lum system, the kinetic energy is (A.l) K = M rn — (x + r sin. 8) dt + m dt (r cos 8) and normalization of this expression with respect to the cart mass, M , defining fi = m/M, results in 2 2 d_ dt (x + r sinf?) + 1 dt (r cos 8) (A.2) 85 Appendix A. Dynamic Equations Derivation 86 Performing the indicated differentiation with respect to time gives K= + t U + r0CosoY + £ (-r9sinO)'' (A.3) K = — + £ 2 2 x 2 + 2rx0 cos 0 + (rfl) 2 cos2 0 + (rfl) 2 sin2 0 (A.4) K = — + ^-2 2 z 2 + (rfl)2 + 2rx0cos0 (A.5) A" = (i + AO Y + ^ r x ^ cos e + £ (re)2 (A.6) and the potential energy is given by P = mgr cos 0, where g = gravitational acceleration; this normalizes to P = figr cos 0. Combining this expression with equation A.6 produces the normalized Lagrangian: L = (1 + p) — + firx9 cos 6 + — (re) — figr cos 0 (A.7) Preceding with the Lagrangian differentiation operations, we see that dL/dx is iden-tically zero, by inspection, while for dL/dx, the result is d_ dt dL dx = Y [(1 + p) x + firO cos 6 = (1 + fi) £ + fir (-02 sin 9 + 9 cos 9) (A.8) dt Examination of the dL/89 and dL/89 cases shows that —• = —firx9 sin 9 + figr sin 9 u9 (A.9) Appendix A. Dynamic Equations Derivation 87 and d_ dt 8L 38 = Y [prxcosB + fir2 8^ = fir (-i0sin0 + xcosfl) + fir28 (A.10) dt Identifying the normalized value of applied force F (Fig. 1.1) as F and combining the derived terms for Lagrange's equation results in (l + p)x + pr (r)'cos0 - 02sin0) = F x cos 8 + r0 — g sin 6 = 0 (A.ll) where a common factor of fir has been extracted from the second equation in the pair. To find separated expressions for x and 8, equations A . l l are cast in matrix form as (A.12) (1 + fi) fir cos 8 X fir82 sin8 + F cos 8 r 8 gsin8 whereupon inversion of the matrix multiplying #j leads to a solution as follows X 8 r —fir cos 8 cos8 (1 + fi) fir82 sin0 + F g sin 8 (A.13) r (1 + fi) — fir cos2 8 Appendix A. Dynamic Equations Derivation 88 Finally, the D. E. s for the inverted pendulum system are obtained by performing the indi-cated matrix multiplication, and applying the trigonometric identity cos2 0 = 1 — sin2 0: x = firO2 sin 0 + F — fig sin 0 cos 0 1 + fi sin2 0 g (1 + ft) sin 0 — firO2 sin 0 cos 0 — F cos 0 r (l + fi sin2 0) (A.H) Appendix B State Feedback Controller Design This Appendix details the design of a linear state-feedback controller, intended to sta-bilize the inverted pendulum in the vertical position, while simultaneously bringing the supporting cart to the origin of the x-axis. The end result is a set of state-variable feed-back gains for the controller, which are also used to obtain expressions for the nonlinear differential equations of the closed-loop system, for simulation purposes. Design of the controller begins with linearization of the system D. E. s (Equation A.14) at the desired operating point; ie, at [a: x 8 #j = 0. The approximations sin 9 s=a 9, cos 9 « 1 are used, resulting in the following linear system: x 0 1 0 0 a : 0 x 0 0 — fig O x 1 9 0 0 0 1 0 0 9 0 0 [g(l + fi)/r] 0 9 -1 /r The general form of this equation is z = Az + B F , where, for state feedback control, the control force, F, is formed as a linear combination of the state variable values, ie, F = Gz, with G = [A B C D], giving z = [A + BG] z. The linear controller design process consists of specifying the values comprising the G vector. We specify these values indirectly, by specifying a desired set of eigenvalues for the controlled (linearized) system; it is worth emphasizing that these eigenvalues apply to the linear system specified above, and, for the actual nonlinear system, are valid only to the extent that the linearization 89 F (B. l ) Appendix B. State Feedback Controller Design 90 itself is valid. We arbitrarily choose a single eigenvalue, of multiplicity four, at (—1 + jO), for the controlled system, and determine the required values for G by matching coefficients in the resulting characteristic equations. The characteristic equation for the desired eigenvalue(s) is (A + l) 4 = 0 A4 + 4A3 + 6A2 + 4A + 1 = 0 (B.2) Coefficient matching requires expansion of det {AI — [A -f BG]} = 0, where I is a 4 by 4 identity matrix. The matrix [A + BG] is 0 A 0 1 B 0 0 (C - fig) 0 0 D 1 (~A/r) (-B/r) [(-C/r) + (1 + fi) /7] (-D/r) (B.3) where the parameter 7 (= r/g) is introduced to reduce the dependence of the equations on definite physical parameters of the pendulum system. It should be noted that 7 is not a dimensionless parameter; it has dimensions of sec2, representing, in some sense, the period of the pendulum. Incorporating the AI term produces the determinant for the characteristic equation: det{AI-[A + BG]} = A - l 0 0 -A (A - B) (pg-C) -D 0 0 A -1 A/r B/r {(C/r) - (1 + p) /7] (A + D/r) (BA) Appendix B. State Feedback Controller Design 91 The expansion of the determinant is tedious, and is not included here; it is read-ily accomplished using Laplace's method (method of cofactors). [CRC] The resulting characteristic equation is A 4 + A 3 + C (1 + p) r 7 A 2 + 7 r BX + + _ M_ 7 r A = 0 (B.5) Matching coefficients with the desired characteristic equation (Equation B.2), a set of equations for G is obtained: 0 - 1 0 1/r -1 0 1/r 0 0 [(-P9/r) + (l + p)h] 0 0 [(-pg/r) + (l + p)h] 0 0 0 A B C D e + (i + p)h 4 1 (B.6) This system is readily solved (recalling that r/g = 7), giving G = [A BCD] = [7 4 7 r[6 + 7 + (l + ^)/7] 4r(l + 7)] (B.7) These, then, are the state feedback coefficients required for the stabilization of the inverted pendulum system. Substituting this feedback function for F in Equations A.14 Appendix B. State Feedback Controller Design 92 produces the full nonlinear, closed-loop, system D. E. s: firB2 sin 9 + 7 (x + 4.x) + r [6 + 7 + (1 + ft) fa] 9 + 4r (1 + 7) 9 - fig sin 9 cos 9 X = -• -------1 + fi sin2 9 (B.8) .. <i±*ilsin0 - fi92 sin0cos0 - {{x + Ax) jg + [6 + 7 + ^ ] 0 + 4 (1 + 7) &} COS9 B = 1 + p sin2 B (B.9) Note that the r factor in the denominator of the 8 expression in Equations A. 14 has been divided out; the result of this operation is to put the new B expression in a form containing only the parameters g, 7, and fi, thereby eliminating direct dependence on the physical parameters of the pendulum system. Unfortunately, it is not possible to do this for the x expression; it retains a factor of r in the first numerator term. Appendix C Partial Singularity List For Reduced-State Feedback T Y P E 9 VALUE EIGENVALUES Stable focus Saddle point Unstable focus 0.000 ±57.541 ±262.049 -2.203 ± j l . l l 5 +2.453, -4.819 +0.304 ± J9.110 Saddle point Unstable focus Saddle point ±445.301 ±626.657 ±807.404 +11.791, -12.152 +0.128 ± J14.227 +16.063, -16.262 Unstable focus Saddle point Unstable focus ±987.878 ±1168.206 ±1348.445 +0.081 ± J17.885 +19.384, -19.522 +0.059 ± J20.903 Saddle point Unstable focus Saddle point ±1528.628 ±1708.773 ±1888.890 +22.205, -22.311 +0.047 + J23.535 +24.702, -24.788 Unstable focus Saddle point Unstable focus ±2068.986 ±2249.068 ±2429.137 +0.039 ± J25.899 +26.967, -27.039 +0.033 ± J28.064 Saddle point Unstable focus Saddle point ±2609.196 ±2789.248 ±2969.294 +29.055, -29.117 +0.028 ± J30.073 +31.002, -31.056 Unstable focus Saddle point Unstable focus ±3149.334 ±3329.370 ±3509.402 +0.025 ± J31.956 +32.833, -32.882 +0.023 + J33.734 93 Appendix D Linearization Near A General Singularity D . l General Procedure This Appendix details the derivation of linearized state equations in the vicinity of a general singularity; these equations may be used to compute eigenvalues at the singularity, to characterize the system behaviour in the corresponding region of the state space. The linearization is accomplished using Taylor series expansion, with a state vector given by V = V p + Vjj; ie, a pole location plus a deviation. The following identities and approximations are used in the derivation: sin(o: ± f3) = sin a cos f3 ± cos a sin f3 cos(a ± /3) = cos a cos f3 + sin a sin (3 cos2 a — sin2 a = cos (2a) sin a « a — (a3/3!) cos a w 1 - (a2/2!) sin2 a w a 2 - (a4/3) + (a6/36) cos2 a « 1 - a 2 + ("4/4) sin a cos a « a — (2a3/3) + (a5/12) The linearization is based on the nonlinear state equations, derived in Appendix B, for the inverted pendulum system with state-variable feedback; the ultimate result of that derivation is a pair of equations, giving expressions for x and 8. The initial step of the linearization is to express these equations as x • (dx/dx) and 8 • (d8/d8), respectively; the resulting expressions for dx/dx and d8/d8 are then linearized. 94 0 Appendix D. Linearization Near A General Singularity 95 D.2 Linearization Of The x Expression The expression for x (Equation B.8) is dx_ _ pr§2 sin 9 - fig sin $ cos 0 + Ax + Bx + C9 + D0 dx x(l + fi sin2 #) (D.l) where the variables A, B, C, and D are used to condense the state feedback coefficient expressions. (See Equation B.7.) Substituting V , as defined in Section D.l , into this expression gives dxd fir(0p + 9d)2 sm(9p + 9d) - fig sin(0p + 9d) cos(0p + 9d) dxA (xp + xd) 1 + fi sin2(0p + 0d)] A(xp + xd) + B(xp + xd) + C(9P + 9d) + D(9P + 9d) (ip + xd) [l + fi sin2(9p + 9d)] + (D.2) We consider the numerator (V) and denominator (Q) of this expression separately, be-ginning with the denominator: Q = (xp + xd) [l + fi sin2 (9P + 9dj\ Q = (xp + xd) 1 + fi(s\n 9P cos 9d + cos 9P sin 9d)Q (D.3) Q = (xp + xd) 1 + fi (sin2 9P cos2 9d + 2 sin 9P cos 9P sin 9d cos 9d + cos2 9P sin2 0^ j (D.4) Appendix D. Linearization Near A General Singularity 96 Q « ( i p 4- xd) sin2 0P I 1 - 92d + -j- ) + 2sin0pcos0p I 0d - =^ + ^ ) + 2^ , 6d 4 y ' " " V " 3 12, cos 0 P k , - - + - (D.5) Neglecting all terms containing powers (greater than one) of elements of V d gives Q « {xp + id) 1 + A4 ( s m 2 #p + 20<* sin 0P cos 0P) ] (D.6) and the remaining terms expand as follows: Q « xp + xd +' p (xp sin2 Op + 2xp0d sin 0P cos #p + xd sin2 0P + 2xd0d sin 0P cos 0P) (D-7) Finally, dropping terms containing only elements of V p , or containing more than one element of V d , results in Q ~ xd + fi (2xp0d sin 0P cos 0P + xd sin2 0P) (D.8) It is important to note that a linear term is present in this result for each element of V d that is eliminated by either the power or the cross-product test. This constraint is required to ensure that singularities are correctly characterized by subsequent eigenvalue analysis using the linearized equations. [Soud] Appendix D. Linearization Near A General Singularity 97 The numerator, V, of Equation D.2 comprises three parts: V = fir (ep + ed)2sin(ep + ed)-fig sin (0P + 9d) cos (9P + 9d) + A (xp + xd) + B (ip + id) + C (9P + 9d) + D (9P + 9d) (D.9) Expansion of Part 1 requires trigonometric identities and approximations as follows: Vi = pr (o2p + 19p9d + 92d) (sin 9P cos 9d + cos 9P sin 9d) (D.10) Vx « fir (92p + 29p9d + 92) The resulting higher powers of 9d and 9d are immediately dropped, and the remaining terms expanded, followed by elimination of the pure V p term and the V d cross-product term: Vi « pr (92p sin 9P + 02p9d cos 9P + 29p9d sin 9P + 29p9d9d cos 9P) (D.12) s i n 0 p ( l - ^ J + cos0. ( D . l l ) Vi « /ir (020d cos 0p + 20p0d sin 9P) (D.13) The expansion of Part 2 is as follows: V2 = pg (sin 0P cos 9d + cos #p sin 9d) (cos 0P cos 9d — sin 0P sin 0d) (D.14) Appendix D. Linearization Near A General Singularity 98 T>2 = ng |sin 9P cos 0p (cos2 9d — sin2 9d) -f sin 9d cos 9d (cos2 9p — sin2 0P)] (D.15) ? 2 « w I sin 9P cos 9P 2 , vd \ I o2 °d , Vd 3 36 + 0, - 2 ^ + S ) ( c o s 2*p- s i n 2*p)} ( ° - 1 6 ) Elimination of the pure V p term, and of all the nonlinear powers of the elements of Vd, results in V2 « i>.g9d (cos2 9P - sin2 0P) 7>2 « rig9dcos(29p) (D.17) Part 3 is already linear, requiring only the elimination of the V p terms, to leave V3 * Axd + Bxd + C9d + D9d (D.18) Note that V3 contains a linear term for each component of Vd, justifying the elimination of all nonlinear terms in V. This completes the reduction of V; the final result is obtained by combining equations D.8, D.13, D.18, and D.17 to produce the final 'in-the-small' equation, applicable near a general singularity (xp,xp, 9p,9p): dxd Axd + Bxd + 9d \pr92p cos 9P + C - ug cos (20p)] + 8d (2ur9p sin 9P + Z>) dxd xd + n (2xp9d sin 9P cos 9P + xd sin2 6P) As intended, this equation is linear in the Vd state deviation variables. (D.19) Appendix D. Linearization Near A General Singularity 99 D.3 Linearization Of The 8 Expression Completion of the linearization procedure requires the 8 expression (Equation B.9), which (using condensed feedback coefficients, as in Section D.2), is dB [(1 + p) h] sin 8 - fi82 sin 8 cos 8- [(Ax + Bx + C8 + D§) /r] cos 8 ~d~8 = 8 (l+ fi sin2 8) ( D ' 2 0 ) The state vector V is as defined in Section D.l , and the numerator and denominator of Equation D.20 are defined as TZ and <5, respectively; noting that S differs from Q (Section D.2) only by the presence of the variable 8 instead of x, the S result is immedi-ately obtained, by comparison with Equation D.8, as S » 8d + p (28p8d sin 8P cos 8P + 8d sin2 8P) P-21) The numerator 1Z of Equation D.20 comprises three parts: [(l + p)h]on(ep + ed)- (D.22) P (0P + Od)2 sin (8P + 8d) cos [8P + 8d) -{[A {Xp + xd) + B (xp + xd) + C (8P + 8d) + D ($p + 0d)} / r} cos (8P + 8d) Expansion of Part 1, with the required trigonometric substitutions, gives Hi = [(1 + fi) 11] (sin 8P cos 8d + cos 8P sin 8d) (D.23) Appendix D. Linearization Near A General Singularity 100 + ft] sin0p 1 1 - ^ 1 + c o s 6 p l 6 d - - £ (D.24) and dropping the pure 9P term, along with the nonlinear powers of 0d, the result is Hi ™0d[(l + p)ft] cos6p (D.25) The expansion of Part 2 is rather messy, although the final result is quite compact. Beginning with TZ2 = ft (^ p + 29pdd + @d) ( s l n c o s @d + cos 0P sin 6d) (cos 0P cos 0d — sin 0P sin 0d) (D.26) the expansion of the trigonometric terms parallels the development of Equations D.14 D.16, leading to 7t2 ss p (§1 + 20p0d + 02d) | sin 0P cos 0P 201 03 + 0, ~Jd + (cos20p-sin20p)j (D.27) At this point, dropping all the nonlinear powers of 0d and 0d results in 1Z2 » p (02p + 20p0d) [sin 0P cos 0P + 0d (cos2 0P - sin2 0P)] (D.28) Appendix D. Linearization Near A General Singularity 101 The (cos2 0P — sin2 9p) term is replaced by cos(20p), the expression is expanded, and the terms containing only elements of V p or only Vd cross-product terms are dropped, obtaining the final result for Part 2: Tl2 « ft [bl sin dp cos ep + 9d92p cos (20p) + 29d9p sin 6P cos 9P + 29p9d9d cos (20p)] (D.29) K2 « p [9d92p cos (29p) + 29d9p sin 9P cos0P] (D.30) Expansion of Part 3, with the appropriate trigonometric approximation substitutions, gives 1l3 w {[A(xp + xd) +B(xp + xd) + C(9p + 6d) +D (9p + 0d)}/r} 92A . „ / „ 03 c o s 0 p l l - ^ j -sm9pl9d (D.31) All the linear terms are present, so the 9d and 9d terms are dropped, leading to K3 « (cos 9P - 9d sin 9P) [A (Xp + xd) + B (xp + xd) + C (9P + 9d) + D (9p + 9d)] / r (D.32) Dropping the pure V p terms, and the Vd cross-product terms (including a new 9d term) produces the Part 3 linearization: K3 w \(Axd + Bxd + C9d + D9d) /r] cos 9P - 9d [(Axp + Bxp + C9P + D9P) /r] sin9P (D.33) Appendix D. Linearization Near A General Singularity 102 The final result for this section is obtained by combining equations D.21, D.25, D.30, and D.33, arriving at a linear, 'in-the-small' equation for dd/dB, as follows: ddd - [(Axd + Bxd) fr] cos 9P - 6d \(D/r) cos 0P + 2p6p sin 9P cos 9P) "77T = : 7 = : \ + (D.34) d&d 6d + p \26dBp sin6Pcos Qp + 6d sin2 6P) Bd{[(Axp + Bxp + C6P + Ddp) I r] sinBp - [(C/r) - (1 + p) /7] cos9P - pd2pcos (2BP)} 9d + p (2Bdep sin 9P cos 6P + Sd sin2 Bp) D.4 Completing The Linearization The expressions in Equations D.19 and D.34, taken together, constitute a set of linearized state equations for the controlled inverted pendulum system, valid at points in the system state space near a singularity defined by (xp, xp, 0P, 0p). This may be made explicit by splitting the two derived expressions into four first-order differential equations: dxd _ dx^ dx^ _ V_ dO^ _ dB^ dfa 11 dxd dt 1 dt ~ Q d8d~ dt ' dt ~ S ' dxd dxrt dBd dt dBj. dt [Q v s nf (D.36) dt dt The expressions for V, Q,l^, and S are still relatively complex, but are greatly sim-plified by the fact that, for the inverted pendulum system, xp and 8P are always zero at singular points, by definition. Using this constraint, the following set of 'in-the-small' Appendix D. Linearization Near A General Singularity 103 equations is obtained: dxdIdt didIdt dOd/dt dBd/dt xd (l + H sin2 0pJ Axd + Bxd + 0d[C- pg cos (29p)] + D9d 0 d ( l + /*sin20p) — [(AE<J COS 9P + Bid cos 0P + D9d cos 9P) / r Od {[(C/r) - (1 + p) ft) cos9P - [{Axp + C9P) jr] sin0p} (D.37) Extraction of the state deviation vector from this set puts the equations in a standard matrix form as follows: id id 9d 9d 1 + p sin2 9P Xd id 9d 9d A 0 B 0 C — pg cos (29p) D 1 + p sin2 9P — (A/r) cos 9P -(B/r) cos 9P [(Axp + C9p) /r] sin9p-[(C/r)-(l + p)/j]cos9p -{D jr) cos 9P (D.38) This expression may be further simplified by noting that all singular points of the system are given by 9P = irax radians, with n = 0,1, 2, • • • , as derived in Section 3.1; for these singular values, sin#p = 0, cos#p = (—1)", and cos(2#p) = 1. Application of these Appendix D. Linearization Near A General Singularity 104 simplifications gives 0 A 0 ( - l ) n + 1 (A/r ) 1 B 0 (-l)n+1(B/r) 0 (C-fig) 0 ( _ i ) » + i [ ( C 7 / r ) - ( l + M)/7] 0 Z> 1 (-l)"+1(/J>/r) (D.39) Finally, substitution of the expanded forms for the state feedback coefficients produces Xd Bd 0d Xd id 9d Od xd Xd 0d 0d xd xd Bd Od 0 7 0 (-l)"+1/<7 1 4 7 0 4{-l)n+1/g 0 [g + r(6 + j)} 0 (-l)"+1(6 + 7 ) 0 4r(l + 7) 1 4 ( - l ) n + 1 ( l + 7) (D.40) which expresses the 'in-the-small' response of the closed-loop system near its singular points. Appendix E About The Figures The plots for this report were created using GRAFTOOL (Version 3.3), a powerful scientifically-oriented plotting program for IBM PC-compatible computers, sold by 3-D Visions Corporation [GRAF]; output was performed on an HP LaserJet II printer, for the monochrome plots, and on an HP PaintJet for those that required color. These details are included because the software created for this project is designed specifically around the capabilities and limitations of GRAFTOOL. In particular, the three-dimensional stability region plots in Chapter 7 are explicitly based on its 'Trajec-tory' plotting mode, which permits (with some work) the plotting of each data point as a set of crossed lines, representing the region of space spanned by that sample point. (See Figure E . l , below.) Point representation x Sample points Figure E . l : Sampling Point Representation Using GRAFTOOL 105 Appendix E. About The Figures 106 The data files generated by the program SCATTER3 (Section H.3) are formatted so as to draw this type of point marker when read by GRAFTOOL using the / T option; the program SLIC20F3 (Section H.2) uses a two-dimensional variation of the same idea to create files for two-dimensional stability region slice plots, such as Figure 7.15. Should one wish to pursue similar research using this software, one should be aware that creating the more complex three-dimensional plots requires a long time and a lot of disk space; approximately 7 hours on a 33 MHz 80386/80387 machine, for example, and at least 10 megabytes of free disk space, after the input file is loaded, is needed to generate Figure 7.8. The results, however, are well worth the wait. Appendix F Simulation Software This Appendix comprises listings of all of the software used for dynamics simulation during the investigation of the inverted pendulum system. The code listings are direct copies of the actual source files, but some editing has been performed to fit the listings into the required page margins for this report; errors may have been introduced during this process. Also, certain long string constants have been arbitrarily split across lines, again for reasons of page formatting; these lines must be reconstructed before the code will compile. All simulations for the work reported in this thesis were performed using a specially developed dynamics simulation program, SIM-CORE, written in C. There are certainly a number of different commercial simulation packages that might have been used for this work, based on either personal computers or mainframes, but the PC-based tools were eliminated from consideration early in the course of research, largely on the basis of cost. Another major consideration was that the results be immediately available in graphical form, as a screen plot; working with lists of trajectory data, or waiting for a plotting program to generate hardcopy was deemed unacceptable, effectively eliminating the mainframe-based programs. The remaining choice led to the development of SIM-CORE; some of its component routines were later adapted for use in computing stabibty regions, as described in Chapter 7. 107 Appendix F. Simulation Software 108 SIM-CORE provides the following capabilities: 1. High-quality simulation, using error-controlled Runge-Kutta integrators, with user-specified error tolerances (absolute and relative). Integrator solutions may be run in either the forward or the reverse direction. 2. Immediate X-Y display of any two state variables of the simulated system, permit-ting time-history or phase-plane plots. 3. Dumping of complete state trajectory information to a disk file, for post-processing. 4. Output file format is compatible with the GRAFTOOL [GRAF] software package, facilitating the creation of report-quality hardcopy plots. 5. Interactive, menu-driven operation, conducive to 'quick-look' and 'what-if investi-gation of system dynamics. The original version of SIM-CORE was developed in PASCAL, and used a 'tuneable integrator' [Smith] as its main solution tool; a fourth-order Runge-Kutta algorithm was later added, based on the examples given in [Chapra]. The Smith integrators have the advantage of being quite fast, which is always important, but they have no inherent error control mechanism; although [Smith] does briefly discuss error control, it was not obvious how this should be implemented. When the SIM-CORE code was converted to C, the Smith integrator was eliminated, and the Runge-Kutta implementation was updated, based on the algorithms found in [Press]; it remains a fourth-order implementation. Press et al also discuss other integration algorithms which might be useful for a program like SIM-CORE, but none of these have been implemented as of this writing. The simulation code presented here comprises the C 'translation' of the original PAS-CAL code, which was motivated by a need for certain features not available with the Appendix F. Simulation Software 109 original implementation. The overall structure of the code was also changed substan-tially in the translation process, particularly as regards the user interface, with a view to future enhancements. The simulations produced by this implementation have been verified against the original code and found to be correct, and the code appears to be stable, but testing has not been as extensive as for the original version, and there may still be occasional problems with the new code. The original code is not presented in this report, as the changes necessitated by some of the reported research are so extensive as to render it irrelevant. This code should be viewed as a set of research tools; it is not of production quality, and, while I have attempted to keep the internal documentation and structure consistent, there may be instances where this is not the case. Caveat researcher! The code for the main simulation program, SIM-CORE.C, plus its user interface code, DISPPROC.C, use certain functions specific to Borland International's Turbo C compiler [BOR]; in particular, various graphics functions are used, and the program, as given here, requires the presence of the Borland graphics library, on directory path C : \ T C \ . The complete program SIM-CORE comprises five program modules, plus three special header files: Header files: SIMUL.H, KEYSTRKS.H, COMPILER.H Source files: SIM-CORE.C, NTGRPROC.C, DYNAMICS.C, DISPPROC.C FILEPROC.C The header files are presented in Section F.2, and Sections F.3 to F.7 present the source code listings. As may be inferred from the previous comment regarding graphics functions, this program has been used solely under Turbo C; a project control file for this compiler Appendix F. Simulation Software 110 is included in Section F.2. Some of the modules used by SIM-CORE (NTGRPROC, DYNAMICS, and FILEPROC) are also used by the program STABVOL, presented in Appendix G; since that program has been ported to run under NDP-C860 [NDP], these modules contain conditionally-compiled elements where necessary for portability. F . l Using The Program There has been a significant amount of effort put into making the SIM-CORE program reasonably easy to use, and, hopefully, simple to modify for use on other dynamic system problems. The outer levels of the program are menu-driven; the original intent of the C-language rewrite, mentioned in the previous section, was to structure the whole pro-gram this way, using pop-up dialog boxes and similar techniques. Unfortunately, time constraints prevented completion of this work, which was essentially peripheral to the actual research subject of this thesis. As a result, the inner program levels retain their original interface, using text-mode question-and-answer interaction with the user. A few notes should suffice to explain the overall program operation. The program comes up in graphics mode, showing a blank phase plane display, with a list of command choices presented along the bottom of the screen. The choices are: Start, Resume, Adjust, Clear, Print, and Exit; each one lists an associated control key, which initiates the command. 'Start' and 'Resume' are very similar; both cause the program to begin executing an integration run on the system dynamics. The difference is that 'Start' loads the initial system state from its internal initial conditions vector, which may be specified by the user, while 'Resume' leaves the state vector unchanged. While an integration run is executing, the message Press any key to interrupt. is displayed; the 'Resume' command allows one to use this feature to temporarily pause a run, then continue it without disturbing Appendix F. Simulation Software 111 the system state. The 'Clear' command erases the graphics screen; this is useful for cases that enter a periodic cycle, to determine that the integration is actually still executing, because, by erasing the old data, it becomes possible to see whether or not a fresh trajectory is drawn. 'Print', at the moment, does nothing; the original implementation had a screen print function that was accessed through this command, but it is not implemented in this version, due to technical problems beyond the scope of this discussion. 'Exit' is self-explanatory; hitting the Escape key brings up a prompt for confirmation, then leaves the program if the exit request is validated. 'Adjust' is the access command for the next level of interaction. Typing A at the top program level switches the system into text mode, and prints the following message: P = Phys, S - ICs, I = Integ, D = Disp, <Esc> = Quit, <Ret> = Keep Changes. This is the menu for the second command level, and it provides access to system Physics, Initial Conditions, Integrator parameters, and Display parameters. The 'hot-key' for each of these commands initiates a short question-answer session with the user, to enter new parameter values. The '<Esc>' and '<Ret>' entries are rather important. Normally, one exits this com-mand level, returning to the top level, by hitting the Enter or Return key; in this case, any changes that have been made to the system parameters are copied into the internal system variables, and the original values are lost. If, however, one modifies the param-eters, and then decides that the modified values are incorrect, using the Escape key instead of Enter dumps the modified values, restoring the values that were in effect on entering the 'Adjust' mode. The question-answer dialogs for parameter changes are self-explanatory, but a few comments may be helpful. When a question requests new parameter values, it shows a list of the current values; the number of parameters in the new list must match the Appendix F. Simulation Software 112 number in the original list, or the question is repeated. Parameter values are entered in standard numeric format; exponential notation may be used, and the individual values must be separated by spaces. There is no way, at present, to go back up a question list and change previous entries; if an incorrect value is entered, simply complete the data entry dialog, and then re-enter it from the top. The 'Display' dialog requires the entry of variable names for display; these text entries are case-sensitive, and entries that do not match the list of available choices are rejected. This dialog also provides the option of writing data to a disk file. If this option is used, data is written to the specified file until the Display option is again selected, and the 'write-to-file' option refused, or a new file specified. Note, however, that the file is rewound each time the 'Start' command is executed; it is not possible to write data from multiple runs to a single file. The user should be aware that the files created by this option can quickly become quite large; the program speed is also degraded. Regarding use of SIM-CORE as a general-purpose simulation platform: The intent of the structure used for this program was to encapsulate all of the system dynamics-specific information in the DYNAMICS module (Section F.5). To date, there has been no testing done to determine whether or not this goal was achieved, but, in theory, one need only change the initial definitions in SysModulelnfo, and rewrite RunConstants and Derivs to reflect new system dynamics, to begin work with a different dynamic system simulation. This does, of course, require that a C compiler be available to compile and bnk in the new module, but this should not prove to be a major limitation to anyone working on such a project based on this document. Appendix F. Simulation Software 113 F.2 A n c i l l a r y Files F.2.1 Project Control Fil e : S I M-CORE.PRJ sim-core.c (keystrks.h, simul.h) ntgrproc.c dynamics.c (simul.h) dispproc.c (keystrks.h, simul.h) f i l e p r o c . c (keystrks.h) F.2.2 Compiler Identification: C O M P I L E R . H This file is used to identify the compiler being used, to control conditional compilation of the code. At present, two compilers are supported: Borland International's Turbo C [BOR], and NDP-C860 (Numeric Data Processing C for the i860), from Micro Way, Inc. [NDP] The file contains a single line: For Turbo C: #define Turbo For NDP-C860: #define NDP F.2.3 Simulation Definitions: SIMUL.H /* Header f i l e defining data for project SIM-CORE */ typedef struct { int Num; /* Number of parameters */ char **Name; /* Ptr to 1st of array of ptrs, to 1st of strings */ double *Val;} /* Ptr to 1st elem of array of parameters */ PFrame; /* The physical system under simulation */ /* endtypedef */ typedef struct { int Num; /* Number of states */ char **Name; /* Ptr to 1st of array of ptrs, to 1st of strings */ double *Val; /* Ptr to 1st elem of array of states */ double *Mul;} /* Ptr to 1st elem of array of scaling factors */ IFrame; /* I n i t i a l conditions for the system states */ /* endtypedef */ Appendix F. Simulation Software 114 typedef struct { int Integr; double TStep; double EndT; double *ValO; double *Vall;} MFrame; /* endtypedef */ /* Integration routine to use */ /* Step-size for the integration */ /* End-point for the variable of integration */ /* Ptr to 1st elem of array of parameters */ /* Ptr to 1st elem of array of parameters */ /* Parameters for the integration routines */ typedef struct { int int double *LoLim; double •HiLim; FILE *0File;} DFrame; /* endtypedef */ Num; /* Number of output variables */ •Index; /* Ptr to 1st elem of array of indices, to output values */ /* Ptr to 1st of array of lower scaling l i m i t s */ /* Ptr to 1st of array of upper scaling l i m i t s */ /* Pointer for output data f i l e */ /* Parameters for the output display */ typedef struct { PFrame IFrame MFrame DFrame SysFrame; /* endtypedef */ •Physics; *ICs; •Math; •Display;} typedef struct { int GraphDriver; int GraphMode; int CursorX; int CursorY; char •TextScreen; int MaxX; int MaxY; int NumDimen;} GraphParms; /• endtypedef •/ typedef struct { int XULHC; int YULHC; /• Identifier for graphics driver used •/ /• Identifier for graphics mode used •/ /• Text-mode cursor location: X coordinate •/ /• Text-mode cursor location: Y coordinate •/ /• Pointer to text-mode screen-save buffer •/ /• Maximum X coordinate value: graphics mode •/ /• Maximum Y coordinate value: graphics mode •/ /• Number of dimensions available on display •/ /• Upper left-hand corner of bitmap; X coord. •/ /• Upper left-hand corner of bitmap; Y coord. •/ Appendix F. Simulation Software 115 int *Map;} / * Pointer to saved background bitmap * / BFrame; / * endtypedef * / typedef struct { int int BFrame char Message; / * endtypedef * / LocX; / * Centre of pop-up window; X coordinate * / LocY; / * Centre of pop-up window; Y coordinate * / *BkGnd; / * Pointer to background bitmap frame * / •Msg;} / * Pointer to message text string * / / * Self-contained unit for graphics-mode pop-up text * / enum {RK4 = 0, Smith = 1}; F.2.4 Keyboard Definitions: KEYSTRKS.H / * Header f i l e , with declarations for key-stroke handling * / #include <bios.h> #define GetScanCode(x) (bioskey (x) >> 8) / * These are key-code definitions, as returned by GetScanCode * / enum {Key_A = 30, Key.B = 48, Key.C = 46, Key.D = 32, Key.E = 18}; enum {Key.F = 33, Key.G = 34, Key_H = 35, Key.I = 23, Key.J = 36}; enum {Key_K = 37, KeyJL = 38, Key_M = 50, Key.N = 49, Key.O = 24}; enum {Key.P = 25, Key.Q = 16, Key_R = 19, Key.S = 31, Key.T = 20}; enum {Key.U = 22, Key.V = 47, Key.W = 17, Key.X = 45, Key.Y = 21}; enum {Key.Z = 44, Key.O =11, Key.I = 2, Key_2 = 3, Key_3 = 4}; enum {Key_4 = 5, Key_5 = 6, Key_6 = 7, Key_7 = 8, Key_8 = 9}; enum {Key_9 = 10}; enum {Key.Tilde = 41, Key.ULine = 12, Key.Equals = 13, Key.BkSlsh = 43}; enum {Key.FCNl = 59, Key_FCN2 = 60, Key_FCN3 = 61, Key_FCN4 = 62}; enum {Key_FCN5 = 63, Key_FCN6 = 64, Key_FCN7 = 65, Key_FCN8 = 66}; enum {Key_FCN9 = 67, Key.FCNIO = 68, Key.Spc = 57, Key.Ent = 28}; enum {Key.BkSpc = 14, Key_Tab = 15, Key.Esc = 1, Key_LtBrc = 26}; enum {Key.RtBrc = 27, Key_Coln = 39, Key.Quo = 40, Key.Coma = 51}; enum {Key.Perd = 52, Key_Que = 53, Key_PrtSc = 55, Key.Minus = 74}; enum {Key.Home = 71, Key.Up = 72, Key.PgUp = 73, Key_Lft = 75}; enum {Key_Rt = 77, Key_End = 79, Key.Dn = 80, Key.PgDn = 81}; enum {Key_Ins = 82, Key.Del = 83, Key.Plus = 78, Key_NumLk5 = 76}; / * Commonly-used definitions, in more descriptive form * / enum {Yes = Key.Y, No = Key.N, Enter = Key.Ent, Esc = Key.Esc}; Appendix F. Simulation Software 116 F.3 Driver: S I M - C O R E . C /* */ /* SIM-CORE.C — Main program module f o r dynamic systems s i m u l a t i o n . */ /* */ /* Design & coding: M.G. Henders */ /* Last modi f ied : 25 Ju ly 1991 */ /* */ /* Module type: Ex terna l l y -suppor ted MAIN */ /* Module contents : Main */ /* */ /* Pro ject f i l e : SIM-CORE.PRJ */ /* Support modules: NTGRPROC.C, DYNAMICS.C, DISPPROC.C */ #define Main main /* */ #include <stdio.h> #include <alloc.h> #include <bios.h> #include <ctype.h> #include <math.h> #include <assert.h> #include <io.h> /* S p e c i a l data dec la ra t ions */ #include <simul.h> #include <keystrks.h> vo id Main() { GraphParms GraphFrame; SysFrame SystemFrame; Message TopMenu = {0, 0 , NULL, "S = S t a r t , R = Resume, A = Ad jus t , C = C l e a r , P = P r i n t , <Esc> = E x i t " } ; Message AdjMenuO = {0, 0 , NULL, "P = Phys, S = ICs , I = In teg , D = D isp , <Esc> = Qui t , <Ret> = Keep Changes"}; Message AdjMenul = {0, 0 , NULL, "S = S t a t e s , I = In tegrator , D = D isp lay , <Esc> = Qui t , <Ret> = Keep Changes"}; Message P r i n t i n g = {0, 0 , NULL, "Screen p r i n t i n g i s not a v a i l a b l e . Press any k e y . . . " } ; Message Running = {0, 0 , NULL, "Running system i n t e g r a t i o n . Press any key to i n t e r r u p t . " } ; Message DoneRun = {0, 0 , NULL, " In tegrat ion run complete. Press any key to c o n t i n u e . . . " } ; Message RunError = {0, 0 , NULL, " Integrator e r r o r ; run terminated. Press any k e y . . . " } ; Message Leaving = {0, 0 , NULL, "Prepared to e x i t SIM-CORE; press Y to c o n f i r m . . . " } ; Appendix F. Simulation Software 117 Message •Menu; double *SysState = NULL, *SysRate = NULL, *Scale = NULL, TStep; double T L e f t , DidStep, NextStep, Tiny = 1 .0e -30 ; long PointCount; i n t SysS ize , ConfirmedExit = 0 , KeyCode = 0 , i ; enum {Start = Key .S , Resume = Key.R, Adjust = Key.A, C lear = Key.C, P r i n t = Key.P}; /* Funct ion prototypes */ vo id MayDump ( in t Code, const char *Func); i n t GraphEnter (GraphParms *GData); i n t SysModulelnfo (SysFrame *SData); i n t DrawScreenFrame (const GraphParms *GData); i n t GraphExit (GraphParms *GData); vo id Der ivs (double •OutRate, const double *InState) ; i n t PopUpGraphMsg (Message *Msg); i n t ClearGraphMsg (Message *Msg); i n t EraseGraphics ( vo id ) ; i n t InteractiveChange (SysFrame *SData, const GraphParms *GData); i n t DrawAxes (const SysFrame *SData, const GraphParms *GData); i n t P l o t R e s u l t s (const double *InState, const SysFrame *SData, const GraphParms *GData); v o i d RKQC (double *0ut, i n t N, double hTry, double Eps, const double *State , const double *Rate, const double *Scale , double *hDid, double *hNext); vo id RunConstants (const double *PVals) ; /* I n i t i a l i z e the graphics system. SIM-CORE cannot run without /* g raph ics , so abort the program i f the i n i t i a l i z a t i o n f a i l s . */ i f (GraphEnter (ftGraphFrame) ) { puts ("Graphics system i n i t i a l i z a t i o n f a i l e d ; can ' t run SIM-CORE."); r e t u r n ; } /* endif */ /* Adjust the screen s i z e ; leave room f o r menu i tems, and draw screen frame. •/ GraphFrame.MaxY -= 14; DrawScreenFrame (ftGraphFrame); /•Create the skeleton of the system parameter s t r u c t u r e . •/ i f ( ( (SystemFrame.Physics = malloc (s izeof (PFrame) ) ) == NULL) I I ( (SystemFrame.ICs = malloc (s i zeof (IFrame) ) ) == NULL) I I ( (SystemFrame.Math = malloc (s i zeof (MFrame) ) ) == NULL) I I ( (SystemFrame.Display = malloc (s izeof (DFrame) ) ) == NULL) ) { GraphExit (ftGraphFrame); Appendix F. Simulation Software puts ("Memory allocation failure in i t ia l iz ing SystemFrame; program terminated."); return; } / * endif * / / * Get parameters from the Dynamics module; abort i f ca l l f a i l s . * / SystemFrame.Display->Num = GraphFrame.NumDimen; MayDump (SysModulelnfo (ftSystemFrame), "SysModulelnfo"); SysSize = SystemFrame.ICs->Num; / * Allocate memory for SysState, SysRate, and Scale, . . . i f ( ( (SysState = malloc (SysSize * sizeof (double) ) ) ( (SysRate = malloc (SysSize * sizeof (double) ) ) ( (Scale = malloc (SysSize * sizeof (double) ) ) == { MayDump (-1, "Memory allocation in Main"); } / * endif * / / * ...copy i n i t i a l state information from SystemFrame, . . . * / for (i = 0; i < SysSize; i++) { SysState[i] = SystemFrame.ICs->Val[i] / SystemFrame.ICs->Mul[i]; } / * endfor * / / * . . .and generate internal constants. * / RunConstants (SystemFrame.Physics->Val); / * Init ial ize screen coordinates for pop-up messages. * / TopMenu.LocX = (1 + GraphFrame.MaxX) / 2; TopMenu.LocY = GraphFrame.MaxY + 7; AdjMenuO.LocX = (1 + GraphFrame.MaxX) / 2; AdjMenuO.LocY = GraphFrame.MaxY + 7; AdjMenu1.LocX = (1 + GraphFrame.MaxX) / 2; AdjMenul.LocY = GraphFrame.MaxY + 7; Printing.LocX = (1 + GraphFrame.MaxX) / 2; Printing.LocY = (1 + GraphFrame.MaxY) / 2; Running.LocX = (1 + GraphFrame.MaxX) / 2; Running.LocY = GraphFrame.MaxY + 7; DoneRun.LocX = (1 + GraphFrame.MaxX) / 2; DoneRun.LocY = (1 + GraphFrame.MaxY) / 2; RunError.LocX = (1 + GraphFrame.MaxX) / 2; RunError.LocY = (1 + GraphFrame.MaxY) / 2; */ == NULL) | | == NULL) | | NULL) ) Appendix F. Simulation Software 119 Leaving.LocX = (1 + GraphFrame.MaxX) / 2; Leaving.LocY = (1 + GraphFrame.MaxY) / 2; I* Put up the main menu bar . */ MayDump (PopUpGraphMsg (ftTopMenu), "PopUpGraphMsg"); /* Select the proper menu f o r the 'Adjust Parameters' case. */ i f (SystemFrame.Physics->Num != 0) Menu = ftAdjMenuO; e l s e Menu = ftAdjMenul; /* Enter main program loop ; e x i t only on command. */ whi le (! ConfirmedExit) { /* Look f o r a keyst roke; the user wants to change something. */ f f l u s h ( s t d i n ) ; KeyCode = GetScanCode (0) ; swi tch (KeyCode) { case Ad just : /* Put up the menu f o r i n t e r a c t i v e adjustment. */ MayDump (ClearGraphMsg (ftTopMenu), "ClearGraphMsg"); MayDump (PopUpGraphMsg (Menu), "PopUpGraphMsg"); /* Talk to the user , f o r changes. */ MayDump (InteractiveChange (ftSystemFrame, ftGraphFrame), " Interact iveChange") ; /* Put back the t o p - l e v e l menu. */ MayDump (ClearGraphMsg (Menu), "ClearGraphMsg"); MayDump (PopUpGraphMsg (ftTopMenu), "PopUpGraphMsg"); DrawScreenFrame (ftGraphFrame); break; case C l e a r : /* Erase the graphics screen, and re-draw the coordinate axes. */ MayDump (EraseGraphicsO, "EraseGraphics") ; MayDump (DrawAxes (ftSystemFrame, ftGraphFrame), "DrawAxes"); break; case P r i n t : /* Put up a response box. */ MayDump (PopUpGraphMsg ( f tPr in t ing) , "PopUpGraphMsg"); bioskey (0) ; MayDump (ClearGraphMsg ( f tPr int ing) , "ClearGraphMsg"); break; case S t a r t : Appendix F. Simulation Software 120 /* Copy i n i t i a l s tate condi t ions i n t o system s t a t e , . . . * / f o r ( i = 0; i < SysSize ; i++) {SysState[ i ] = SystemFrame.ICs->Val[i] / SystemFrame.ICs->Mul[i] ; } /* endfor */ /* . . . r e s e t the in tegra to r s t e p - s i z e , . . . */ i f ( (TLeft = SystemFrame.Hath->EndT - SysState[0]) < 0) { TStep = - SystemFrame.Math->TStep; } e l s e { TStep = SystemFrame.Math->TStep; } /* endif */ TLeft = fabs (TLef t ) ; case Resume: /* . . . a n d procede i n t o the system i n t e g r a t i o n . */ MayDump (ClearGraphMsg (ftTopMenu), "ClearGraphMsg"); MayDump (PopUpGraphMsg (ftRunning), "PopUpGraphMsg"); /* Disp lay the frame, axes, and i n i t i a l p o i n t . */ MayDump (DrawScreenFrame (ftGraphFrame), "DrawScreenFrame"); MayDump (DrawAxes (ftSystemFrame, &GraphFrame), "DrawAxes"); MayDump (P lotResul ts (SysState, ftSystemFrame, ftGraphFrame), " P l o t R e s u l t s " ) ; /* Update i n t e r n a l constants . */ RunConstants (SystemFrame.Physics->Val); /* Star t f i l e output , i f requested. */ i f (SystemFrame.Display->0File != NULL) { PointCount = 0; chs ize ( f i l e n o (SystemFrame.Display->0File), 0); f p r i n t f (SystemFrame.Display->0File, 11 \n") ; } /* endif */ /* Loop u n t i l i n t e r r u p t e d , or u n t i l complete. */ f o r ( ; ; PointCount++) { /* Check f o r in te r rupt key. */ i f ( (KeyCode = bioskey ( l ) ) != 0) break; /* Check f o r completion of the i n t e g r a t i o n run . */ i f (TLeft <= 0) { MayDump (PopUpGraphMsg (ftDoneRun), "PopUpGraphMsg"); bioskey (0); MayDump (ClearGraphMsg (ftDoneRun), "ClearGraphMsg"); break; } /* endif */ Appendix F. Simulation Software 121 / • D o the integration step. * / switch (SystemFrame.Math->Integr) { case RK4: / * Compute system derivatives, . . . * / Derivs (SysRate, SysState); / * . . . i n i t i a l i z e tolerances, . . . * / Scale[0] = 0; for (i = 1; i < SysSize; i++) { Scale[i] = fabs (SysState[i]) + fabs (TStep * SysRate[i]) + Tiny; } / * endfor * / / * ...adjust step-size; avoid overshooting endpoint, . . . * / TLeft = fabs (SystemFrame.Math->EndT) - fabs (SysState[0]); i f (TLeft < fabs (TStep) ) { i f (TStep > 0) TStep = TLeft; else TStep = - TLeft; } / * endif * / / * .. .and perform one integration step. * / RKQC (SysState, SysSize, TStep, SystemFrame.Math->ValO[0], SysState, SysRate, Scale, ftDidStep, ftNextStep); •/* Check completion conditions. * / i f ( (DidStep == 0) && (NextStep == 0) ) { i f (SysState[0] != SystemFrame.Math->EndT) { MayDump (PopUpGraphMsg (ftRunError), "PopUpGraphMsg"); bioskey (0); MayDump (ClearGraphMsg (ftRunError), "ClearGraphMsg"); break; } / * endif * / } / * endif * / else { / * Recycle with new step-size. * / TStep = NextStep; } / * endelse * / break; case Smith: MayDump (-1, "Smith integrator not implemented"); break; default: / * Nothing * / ; } / * endswitch * / Appendix F. Simulation Software / * Display the results, . . . * / MayDump (PlotResults (SysState, ftSystemFrame, ftGraphFrame), "PlotResults"); / * . . .writing to output f i l e i f requested. * / i f (SystemFrame.Display->OFile != NULL) { fprintf (SystemFrame.Display->OFile, "•/.f y.f y.f y.f Xfw, s y s state[o], SystemFrame.ICs->Mul[l] * SysState[1], SystemFrame.ICs->Mul[2] * SysState[2], SystemFrame.ICs->Mul[3] * SysState[3], SystemFrame.ICs->Mul[4] * SysState[4]); } / * endif * / } / * endfor * / / * Write the number of data points into the f i l e . * / i f (SystemFrame.Display->OFile != NULL) { rewind (SystemFrame.Display->OFile); fprintf (SystemFrame .Display->OFile, "'/.ld 5\n" , PointCount) ; } / * endif * / MayDump (ClearGraphMsg (ftRunning), "ClearGraphMsg"); MayDump (PopUpGraphMsg (ftTopMenu), "PopUpGraphMsg"); break; case Esc: / * Do you really want to do this? * / MayDump (PopUpGraphMsg (ftLeaving), "PopUpGraphMsg"); fflush (stdin); i f (GetScanCode (0) == Yes) ConfirmedExit = 1; MayDump (ClearGraphMsg (ftLeaving), "ClearGraphMsg"); break; default: / * Nothing * / ; } / * endswitch * / } / * endwhile * / / * Close output f i l e and restore original screen before leaving. * / i f (SystemFrame.Display->OFile != NULL) { fclose (SystemFrame.Display->0File); } / * endif * / GraphExit (ftGraphFrame); } / * end Main * / Appendix F. Simulation Software 123 F .4 Numerical Integration: N T G R P R O C . C /*_ _ */ /* NTGRPROC.C — Numerical i n t e g r a t i o n r o u t i n e s . */ /* */ /* Design ft coding: M.G. Henders */ /* Last modi f ied : 22 A p r i l 1991 */ /* */ /* Nodule type : Support module. */ /* Module contents : RKQC, RK4th */ /* */ /* Support f i l e s : None */ /*_ */ #include <compiler.h> #include <assert.h> ftinclude <math.h> /* Cond i t iona l dec la ra t ions */ #ifdef Turbo /* For Turbo C */ #include <alloc.h> #endif #ifdef NDP /* For NDP C */ #include <stdio.h> #include <stdl ib.h> #endif vo id RKQC (double *0ut, i n t N, double hTry, double Eps, double *State, double *Rate, double *Scale , double *hDid, double *hNext) { /* */ /* Purpose: Stepper rout ine f o r Runge-Kutta i n t e g r a t o r , w i th */ /* v a r i a b l e s t e p - s i z e . See RK4th, below. */ /* */ /* Uses: Out A po in te r to the output (post -step) s ta te vector */ /* N The number of s tate va r iab les i n the system */ /* hTry The s t e p - s i z e to t r y i n i t i a l l y */ /* Eps The s o l u t i o n to lerance requ i red */ /* State A po in te r to the system s tate vector */ /* Rate A po in te r to the system d e r i v a t i v e s vector */ /* Scale A po in te r to the er ror s c a l i n g vector */ /* hDid A po in te r to the s t e p - s i z e a c t u a l l y used */ /* hNext A po in te r to the s t e p - s i z e to t r y next time */ /* */ /* Returns: Nothing */ /* */ s t a t i c double *S1 = NULL, *S2 = NULL, *R1 = NULL, pGrow = - 0 . 2 0 ; Appendix F. Simulation Software s t a t i c double pShrink = - 0 . 2 5 , fCor = 0.0666666666666, Safety = 0 ErrCon = 6 . 0 e - 4 ; double h , H a l f , Temp, ErrMax, TSave, T i ; i n t Done, E r r o r , /* Funct ion prototypes */ vo id RK4th (double *0ut, i n t N, double h , double *S0, double *K1); v o i d Der ivs (double *0utRate, double * Instate ) ; /* A l l o c a t e i n t e r n a l vec to rs , i f not already done. */ i f (SI == NULL) asser t ( (SI = (double*) malloc (N * s i zeof (double) ) ) != NULL) i f (S2 == NULL) asser t ( (S2 = (double*) malloc (N * s i zeof (double) ) ) != NULL) i f (Rl == NULL) asser t ( (Rl = (double*) malloc (N * s i zeof (double) ) ) != NULL) /* I n i t i a l t r i a l s t e p - s i z e va lue . */ h = hTry; E r ro r = 0 ; Done = 0 ; TSave = State [0] ; T i = State [0] ; do { /* S p l i t time step i n h a l f , . . . */ Half - '0.5 * h; /* . . . t a k e two steps at t h i s s t e p - s i z e , . . . */ RK4th (SI , N, H a l f , S t a t e , Rate) ; Der ivs ( R l , S I ) ; RK4th (S2, N, H a l f , S I , R l ) ; /* . . . a d j u s t the t ime , . . . */ T i = T i + h; /* . . . a n d t e s t f o r v a l i d s tep . */ i f (Ti != TSave) { /* Take a step wi th the f u l l s t e p - s i z e , . . . */ RK4th (SI , N, h , S t a t e , Rate) ; /* . . . a n d evaluate accuracy. */ f o r ( i = 1, ErrMax = 0 ; i < N; i++) i. /* Get t runcat ion er ror estimate i n S I , f o r l a t e r use. */ S l [ i ] = S2[ i ] - S l [ i ] ; Temp = fabs ( S l [ i ] / S c a l e [ i ] ) ; i f (Temp > ErrMax) ErrMax = Temp; } /* endfor */ Appendix F. Simulation Software 125 /* Scale to requ i red t o l e r a n c e , . . . */ ErrMax = ErrMax / Eps; /* . . . a n d adjust s t e p - s i z e . */ i f (ErrMax <= 1.0) { /* Step was OK; set f o r next s tep , and set completion f l a g . */ •hDid = h; Done = 1; i f (ErrMax > ErrCon) { *hNext = Safety * h * exp (pGrow * loglO (ErrMax) ) ; } /* endif */ e l s e *hNext = 4 .0 * h; /* Update time output. */ Out[0] = T i ; /* Add i n the t runcat ion e r r o r . */ f o r ( i = 1; i < N; i++) Out[ i ] = S2[ i ] + fCor * S l [ i ] ; } /* endif */ e l se { /* Set new s t e p - s i z e and t r y aga in . */ h = Safety * h * exp (pShrink * loglO (ErrMax) ) ; } /* endelse */ } /* endif */ e l se { /* Step s i z e too s m a l l . */ E r ro r = 1; *hDid = 0; *hNext = 0 ; } /* endelse */ } whi le ( (! Done) && (! Er ror ) ) ; r e t u r n ; } /* end RKQC */ vo id RK4th (double *0ut, i n t N, double h , double *S0, double *K1) { /* */ /* Purpose: Perform 4th -order Runge-Kutta i n t e g r a t i o n s tep . */ /* Reference: W.H. P ress , B .P . F lannery , S .A . Teukolsky, */ /* W.T. V e t t e r l i n g ; 'Numerical Recipes i n C - The Art of S c i e n t i f i c */ /* Computing'; Cambridge U n i v e r s i t y P ress ; 1988; pg 572. */ /* */ /* Uses: Out A po in te r f o r the pos t - s tep s ta te vector */ /* N The number of s tate va r iab les i n the system */ /* h The s t e p - s i z e f o r the i n t e g r a t i o n */ /* SO A po in te r to the system s tate vec to r , p r e - s t e p */ /* KI A po in te r to the system d e r i v a t i v e s v e c t o r , p re - s tep */ /* */ /* Returns: Nothing */ Appendix F. Simulation Software 126 /* */ static double *S1 = NULL, *K2 = NULL, *K3 = NULL, *K4 = NULL; double HalfStep, SixthStep; int i ; / * Function prototypes * / void Derivs (double *0utRate, double *InState); / * Allocate internal vectors, i f not already done. * / i f (SI == NULL) assert ( (SI = (double*) malloc (N * sizeof (double) ) ) i f (K2 == NULL) assert ( (K2 = (double*) malloc (N * sizeof (double) ) ) i f (K3 == NULL) assert ( (K3 = (double*) malloc (N * sizeof (double) ) ) i f (K4 == NULL) assert ( (K4 = (double*) malloc (N * sizeof (double) ) ) NULL) NULL) NULL) NULL) HalfStep = 0.5 * h; SixthStep = 0.16666666666666 * h; / * First step. * / S1[0] = S0[0] + Half Step; for (i = 1; i < N; i++) Sl[i] = S0[i] + HalfStep * KI[i]; Derivs (K2, SI); / * Second step. * / for ( i = 1; i < N; i++) Sl[i] = S0[i] + HalfStep * K2[i]; Derivs (K3, SI); / * Third step. * / S1[0] = S0[0] + h; for (i = 1; i < N; i++) Sl[i] S0[i] + h * K3[i] ; Derivs (K4, SI); / * Final step. * / for (i = 1; i < N; i++) { 0ut[i] = S0[i] + SixthStep * (KI [i] + 2.0 * (K2[i] + K3[i]) + K4[i]); } / * endfor * / return; } / * end RK4th * / Appendix F. Simulation Software 127 F . 5 S y s t e m D e f i n i t i o n : D Y N A M I C S . C /* */ /* DYNAMICS.C — Support module f o r SIM-CORE.C, f o r dynamic systems */ /* s i m u l a t i o n . This module s p e c i f i e s d e t a i l s of the a c t u a l dynamic */ /* system under cons idera t ion . In genera l , to implement a new simu- */ /* l a t i o n , one should rewr i te Der i vs , to provide new system d e r i v a t i v e */ /* f u n c t i o n s , and RunConstants, which i s used simply to pre-compute */ /* constants , to acce lerate the s i m u l a t i o n . Only the i n i t i a l d e c l a r - */ /* a t ions i n SysModuleInfo should be changed; the remaining code i s */ /* intended to handle system changes automat ica l l y . */ /* */ /* Design & coding: M.G. Henders */ /* Last modi f ied : 25 Ju ly 1991 */ /* */ /* Module type: Support module. */ /* Module contents : SysModuleInfo, RunConstants, Der ivs */ /* */ /* Pro jec t f i l e : SIM-CORE.PRJ */ /* Support f i l e s : SIMUL.H */ /* */ #include <compiler.h> #include <stdio.h> #include <string.h> #include <stdl ib.h> #include <f loat f i x .h> #include <math.h> /* S p e c i a l data types */ #include <simul.h> /* Cond i t iona l dec la ra t ions */ #ifdef Turbo /* For Turbo C */ #include <alloc.h> #endif i n t SysModuleInfo (SysFrame *SData) { /* */ /* Purpose: To supply parameters of the dynamic system under s imu- */ /* l a t i o n . This f u n c t i o n should be c a l l e d once, AND ONLY */ /* ONCE, on s t a r t i n g a s imulat ion sess ion . Memory i s a l l o c a t e d */ /* i n t e r n a l l y , and repeated c a l l s may have unexpected e f f e c t s , */ /* i n c l u d i n g the product ion of crashes due to lack of memory. */ /* */ /* Uses: SData, a po in te r to type SysFrame (see f i l e SIMUL.H). */ /* SysModuleInfo w i l l modify the s t ructure po inted to by */ Appendix F. Simulation Software 128 /* SData; i t must be already a l l o c a t e d on entry . */ /* */ /* Returns: The func t ion i t s e l f returns an i n t e g e r , 0 on success fu l */ /* complet ion, < 0 on e r r o r . The f o l l o w i n g e n t r i e s i n the */ /* SData s t ruc ture are updated: */ /* */ /* Physics->Num: Number of adjustable parameters i n the system DEs. */ /* Physics->Name: Po in ters to name s t r i n g s of the adjustable parms. */ /* Phys ics ->Va l : I n i t i a l (defaul t ) values f o r the parameters. */ /* ICs->Num: The number of s ta te v a r i a b l e s i n the system, i n c l u d i n g */ /* the v a r i a b l e of i n t e g r a t i o n (usual ly Time). */ /* ICs->Name: Po in ters to name s t r i n g s f o r the s ta te v a r i a b l e s . */ /* ICs->Val : I n i t i a l (defaul t ) values f o r the system s t a t e s . */ /* ICs->Mul: S c a l i n g f a c t o r s f o r s ta te v a r i a b l e d i s p l a y */ /* Math->Integr: The defau l t numerical i n t e g r a t i o n rout ine to use. */ /* Math->TStep: The i n i t i a l s t e p - s i z e f o r the i n t e g r a t i o n . */ /* Math->EndT: The f i n a l value f o r the v a r i a b l e of i n t e g r a t i o n . */ /* (This may be negat ive , f o r reverse - t ime s o l u t i o n s . ) */ /* Math->ValO: For the RK4 (Runga-Kutta four th -o rder ) i n t e g r a t o r , */ /* the absolute s ta te er ror to le rance . */ /* Math->Val l : For the RK4 i n t e g r a t o r , the r e l a t i v e s ta te e r r o r . */ /* Display->Index: Array ind ices to access s ta te v a r i a b l e s */ /* f o r d i s p l a y . */ /* Display->LoLim: Lower l i m i t value f o r d i s p l a y s c a l i n g . */ /* Display->HiLim: Upper l i m i t value f o r d i s p l a y s c a l i n g . */ /* */ /* NOTE: The ->Name f i e l d s are mult i -e lement po in te r a r r a y s ; space */ /* f o r these arrays i s a l l o c a t e d i n t e r n a l l y to SysModulelnfo, s ince */ /* the s i z e of the arrays i s determined by the ->Num f i e l d e n t r i e s . */ /* S i m i l a r l y , the ->Valx, ->Index, and ->xxLim f i e l d s are m u l t i - */ /* element arrays of numeric va lues , and are a l l o c a t e d i n t e r n a l l y . */ /* The value of the f i e l d Display->Num i s NOT returned by SysModule- */ /* In fo ; on ent ry , i t must be set to the number of d isp layed */ /* v a r i a b l e s f o r the system. */ /* */ /* These dec la ra t ions set up system d e f a u l t s ; change them, and only */ /* them, when implementing a new set of system dynamics. */ i n t PhysNum = 7; s t a t i c char *PhysName = "Mc Mp Lp kX kXDot kTheta kThetaDot"; char *PhysVals = " 1 . 0 0 .1 1.0 1.0 1.0 1.0 1 . 0 " ; i n t ICsNum = 5; s t a t i c char *ICsName = "Time X XDot Theta ThetaDot"; char *ICsVals = " 0 . 0 0 .0 0.0 30.0 0 . 0 " ; char *Scal ing = " 1 . 0 1.0 1.0 57.29577951 57.29577951"; Appendix F. Simulation Software i n t Sys lntegr = RK4; double StepSize = 1 . 0 e - 9 ; double RunEnd = 2000.0; char *AbsErr = " 1 . 0 e - 6 1.0e-6 1.0e-6 1 . 0 e - 6 " ; char *RelErr = " 1 . 0 e - 6 1.0e-6 1.0e-6 1 . 0 e - 6 " ; char *DispName = "Theta ThetaDot"; char *DispRange = " - 4 0 0 . 0 400.0 - 9 0 0 . 0 9 0 0 . 0 " ; /* End of d e c l a r a t i o n of d e f a u l t s , e t c . */ i n t i , j , DSize; char **TName; /* Set up the Physics parameters, . . . */ SData->Physics->Num = PhysNum; i f ( (SData->Physics->Name = (char**) malloc (PhysNum * s i zeo f (char*) ) ) == NULL) re tu rn ( - 1 ) ; i f ( (SData->Physics->Val = (double*) malloc (PhysNum * s i zeo f (double) ) ) == NULL) re tu rn ( - 1 ) ; i f ( (SData->Physics->Name[0] = s t r t o k (PhysName, " ") ) == NULL) r e t u r n ( - 2 ) ; f o r ( i = 1; i < PhysNum; i++) { i f ( (SData->Physics->Name[i] = s t r t o k (NULL, " ") ) == NULL) r e t u r n ( - 2 ) ; } /* endfor */ f o r ( i = 0 ; i < PhysNum; i++, PhysVals++) { SData->Physics->Val[ i ] = s t r t o d (PhysVals, ftPhysVals); } /* endfor */ /* . . . a n d then do the ICs parameters, . . . */ SData->ICs->Num = ICsNum; i f ( (SData->ICs->Name = (char**) malloc (ICsNum * s i zeof (char*) ) ) == NULL) re tu rn ( - 1 ) ; i f ( (SData->ICs->Val = (double*) malloc (ICsNum * s i zeo f (double) ) ) == NULL) re tu rn ( - 1 ) ; i f ( (SData->ICs->Mul = (double*) malloc (ICsNum * s i zeo f (double) ) ) == NULL) re tu rn ( - 1 ) ; i f ( (SData->ICs->Name[0] = s t r t o k (ICsName, " ") ) == NULL) re tu rn ( - 2 ) ; f o r ( i = 1; i < ICsNum; i++) Appendix F. Simulation Software { i f ( (SData->ICs->Name[i] = strtok (NULL, " ") ) == NULL) return (-2); } / * endfor * / for ( i = 0; i < ICsNum; i++, ICsVals++) { SData->ICs->Val[i] = strtod (ICsVals, ftlCsVals); } / * endfor * / for ( i = 0; i < ICsNum; i++, Scaling++) { SData->ICs->Mul[i] = strtod (Scaling, ftScaling); } / * endfor * / / * ...followed by the Math (integrator) parameters, . . . * / SData->Math->Integr = Syslntegr; SData->Math->TStep = StepSize; SData->Math->EndT = RunEnd; i f ( (SData->Math->ValO = (double*) malloc ( (ICsNum - 1) * sizeof (double) ) ) == NULL) return (-1); i f ( (SData->Math->Vall = (double*) malloc ( (ICsNum - 1) * sizeof (double) ) ) == NULL) return (-1); for (i = 0; i < (ICsNum - 1); i++, AbsErr++, RelErr++) { SData->Math->ValO[i] = strtod (AbsErr, ftAbsErr); SData->Math->Vall[i] = strtod (RelErr, ftRelErr); } / * endfor * / / * . . .and, f inal ly , the Display parameters. * / DSize = SData->Display->Num; i f ( (TName = (char**) malloc (DSize * sizeof (char*) ) ) == NULL) return (-1); i f ( (SData->Display->Index = (int*) malloc (DSize * sizeof (int) ) ) == NULL) return (-1); i f ( (SData->Display->LoLim = (double*) malloc (DSize * sizeof (double) ) ) NULL) return (-1); i f ( (SData->Display->HiLim = (double*) malloc (DSize * sizeof (double) ) ) == NULL) return (-1); i f ( (TName[0] • strtok (DispName, " ") ) == NULL) return (-2); for ( i = 1; i < SData->Display->Num; i++) { i f ( (TName[i] = strtok (NULL, " ") ) == NULL) return (-2); } / * endfor * / for (i = 0; i < DSize; i++) { SData->Display->Index[i] = -1; for (j = 0; j < ICsNum; j++) Appendix F. Simulation Software 131 { i f (! strcmp (SData->ICs->Name[j], TName[i]) ) SData->Display->Index[i] = j ; } / * endfor * / i f (SData->Display->Index[i] < 0) return (-3); } / * endfor * / free (TName); for (i = 0; i < DSize; i++, DispRange++) { SData->Display->LoLim[i] = strtod (DispRange, &DispRange); DispRange++; SData->Display->HiLim[i] = strtod (DispRange, ftDispRange); } / * endfor * / SData->Display->OFile = NULL; return (0); } / * end SysModulelnfo * / / * Static variables; pass outside the Derivs function c a l l , for speed. * / double KfX, KfVx, KfA, KfVa, Mc, Mp, Lp, Mpg, MpLp, rLp, Mtg; void RunConstants (double *PVals) { /* _ ^ */ / * Purpose: To calculate secondary parameter values. * / / * * / static double g = 9.80665;/* Gravitational acceleration, (m/sec/sec) * / / * 0 -> Mc, 1 -> Mp, 2 -> Lp, 3 -> kX, 4 -> kXDot, 5-> kTheta, 6 -> kThetaDot * / Mc = PVals[0]; Mp = PVals[1]; Lp = PVals[2]; KfX = Mc * Lp / g; KfVx = 4.0 * Mc * Lp / g; KfA = Mc * Lp * (6.0 + Lp / g) + (Mp + Mc) * g; KfVa = 4.0 * Mc * Lp * (1.0 + Lp / g); KfX *= PVals[3]; KfVx *= PVals[4]; KfA *= PVals[5]; KfVa *= PVals[6]; rLp = 1.0 / Lp; MpLp = Mp * Lp; Appendix F. Simulation Software 132 Mpg = Mp * g ; Mtg = (Mp + Mc) * g ; r e t u r n ; } /* end RunConstants */ vo id Der ivs (double * 0 u t R a t e , double *InState) { /* _ */ /* Purpose: To compute the de r i va t i ves of the dynamic system being */ /* s imulated, given the current system s tate vec to r . */ /* */ /* Uses: Two po in ters to type double, InState and OutRate. */ /* I m p l i c i t l y , InState po in ts to the f i r s t element of the */ /* system s tate vec to r , whi le OutRate po in ts to the f i r s t element */ /* f o r the d e r i v a t i v e s vector . */ double X, Vx, A, Va, U, sqVa, sA, cA, r e c i , sAcA; /* Load l o c a l v a r i a b l e s . */ X = InState [1] ; Vx = InState[2]; A = Instate[3]; Va = InState[4]; /* State feedback c o n t r o l f u n c t i o n . */ U = (KfX * X) + (Kf Vx * Vx) + (KfA * A) + (KfVa * Va) ; /* Ca lcu la te m u l t i p l y - u s e d funct ions once on ly , f o r speed. */ sqVa = Va * Va; sA = s i n (A); cA = cos (A) ; sAcA = sA * cA; r e c i = 1.0 / (Mc + Mp * sA * sA) ; /* Compute the system d e r i v a t i v e s . */ OutRate[1] = Vx; OutRate[2] = r e c i * (U - Mpg * sAcA + MpLp * sA * sqVa); OutRate[3] = Va; OutRate[4] = rLp * r e c i * (Mtg * sA - U * cA - MpLp * sqVa * sAcA); r e t u r n ; } /* end Derivs */ /* /* Returns: Nothing /* */ */ */ Appendix F. Simulation Software 133 F.6 User Interface: D I S P P R O C . C /* _ */ /* DISPPROC.C — Graphics support r o u t i n e s . */ /* */ /* Design ft coding: M.G. Henders */ /* Last modi f ied : 25 Ju ly 1991 */ /* */ /* Module type: Support module. */ /* Module contents : GraphEnter, GraphExit , DrawScreenFrame, DrawAxes,*/ /* EraseGraphics, InteractiveChange, PopUpGraphMsg, */ /* ClearGraphMsg, P l o t R e s u l t s , MayDump */ /* */ /* Pro ject f i l e : SIM-CORE.PRJ */ /* Support f i l e s : SIMUL.H */ /* */ #include <alloc.h> #include <stdio.h> #include <conio.h> #include <graphics.h> #include <bios.h> #include <string.h> #include <stdl ib.h> /* S p e c i a l data dec larat ions */ #include <simul.h> #include <keystrks.h> typedef s t r u c t { char *Label double Va l i n t XULHC i n t YULHC i n t Width i n t Height i n t *BitMap i n t NextRt i n t NextLt i n t NextUp i n t NextDn DataBox; /* endtypedef */ /* Po in ter to i d e n t i f i e r s t r i n g */ /* Data value */ /* Upper l e f t - h a n d corner p o s i t i o n ; X coord */ /* Upper l e f t - h a n d corner p o s i t i o n ; Y coord */ /* Width of box, i n p i x e l s , i n c l u d i n g frame */ /* Height of box, i n p i x e l s , i n c l u d i n g frame */ /* Po in ter to save area f o r the background */ /* Index number of next c u r s o r - r i g h t box */ /* Index number of next c u r s o r - l e f t box */ /* Index number of next cursor -up box */ /* Index number of next cursor-down box */ enum {Up = Key_Up, Dn = Key.Dn, Rt = Key .Rt , Lt = Key_Lft } ; Appendix F. Simulation Software 134 i n t GraphEnter (GraphParms *GData) { /* */ /* Purpose: Sets up the graphics system, i n the highest a v a i l a b l e */ /* r e s o l u t i o n , and passes back informat ion about the */ /* r e s u l t s i n GData. The tex t screen i s saved i n an i n t e r n a l l y - */ /* a l l o c a t e d b u f f e r , and a po in te r to i t i s a l so returned i n GData. */ /* Each c a l l of Graph-Enter overwrites a n y / a l l data i n t h i s b u f f e r . */ /* */ /* Uses: GData, a po in te r to a data s t ructure of type GraphParms. */ /* This s t ructure type i s def ined i n the header f i l e SIMUL.H. */ /* */ /* Returns: An integer va lue ; returns 0 on success fu l complet ion, */ /* or an er ror code < 0 i f the setup procedure f a i l s f o r */ /* any reason. E r ro r messages are p r i n t e d l o c a l l y . */ /* */ s t a t i c char *ScreenBuf = NULL; char *GraphDriverPath = "C:\\TC\\"; i n t Dimensions = 2; /* A l l o c a t e the screen-save b u f f e r , i f not already done. */ i f (ScreenBuf == NULL) { i f ( (ScreenBuf = malloc (80 * 25 * 2) ) == NULL) •{ puts ("NULL po in te r encountered f o r text-mode screen save; GraphEnter t e r m i n a t e d . " ) ; r e t u r n ( - 1 ) ; } /* endif */ } /* endif */ /* Save the text-mode screen s t a t e . */ GData->CursorX = wherex(); GData->CursorY = whereyO; GData->TextScreen = ScreenBuf; i f (! get text (1 , 1, 80, 25, ScreenBuf) ) { puts ("Text-mode screen save f a i l e d ; GraphEnter t e r m i n a t e d . " ) ; re tu rn ( - 2 ) ; } /* endif */ /* I n i t i a l i z e the graphics system. */ GData->GraphDriver = DETECT; initgraph(&(GData->GraphDriver) , &(GData->GraphMode), GraphDriverPath); i f (GData->GraphDriver < 0) { /* Deal w i th graphics i n i t i a l i z a t i o n e r r o r . */ puts (grapherrormsg (GData->GraphDriver) ) ; r e t u r n ( - 3 ) ; Appendix F. Simulation Software 135 } /* endif */ /* I n i t i a l i z e screen v iewport , with c l i p p i n g enabled. */ GData->MaxX = getmaxxO; GData->MaxY = getmaxyO; setviewport (0, 0 , GData->MaxX, GData->MaxY, 1 ) ; i f (graphresult ( ) ) re tu rn ( - 4 ) ; /* Set the number of a v a i l a b l e d i s p l a y dimensions. */ GData->NumDimen = Dimensions; r e t u r n (0) ; } /* end GraphEnter */ i n t GraphExit (GraphParms *GData) { /* */ /* Purpose: Closes the graphics system, and res tores the o r i g i n a l */ /* text-mode screen, us ing informat ion passed i n GData. */ /* */ /* Uses: GData, a po in te r to a data s t ructure of type GraphParms. */ /* This s t ructure type i s def ined i n the header f i l e SIMUL.H. */ /* */ /* Returns: An in teger va lue ; returns 0 on success fu l complet ion, */ /* or an er ror code < 0 i f the e x i t procedure f a i l s f o r */ /* any reason. E r ro r messages are p r i n t e d l o c a l l y . */ /* */ /* F l i p back to graphics page 0 . */ closegraphO ; i f (! put text (1 , 1 , 80, 25, GData->TextScreen) ) { puts ( "Restorat ion of text-mode screen f a i l e d ; GraphExit t e r m i n a t e d . " ) ; re tu rn ( - 1 ) ; } /* endif */ f r e e (GData->TextScreen); GData->TextScreen = NULL; gotoxy (GData->CursorX, GData->CursorY); r e t u r n (0) ; } /* end GraphExit */ i n t DrawScreenFrame (const GraphParms *GData) Appendix F. Simulation Software 136 { /* */ / * Purpose: Draws a rectangular frame around the screen, using * / / * minimum values = 0, maximum values as passed in GData. * / / * Adjusts the viewport so that subsequent _clearviewport calls wi l l * / / * not erase the screen frame. * / / * * / / * Uses: GData, a pointer to a data structure of type GraphParms. * / / * This structure type is defined in the header f i l e SIMUL.H. * / / * * / / * Returns: An integer value; returns 0 on successful completion, * / / * or an error code < 0 i f the procedure fa i l s for any * / / * reason. * / / * * / setviewport (0, 0, GData->MaxX, GData->MaxY, 1); i f (graphresult() ) return (-1); rectangle (0, 0, GData->MaxX, GData->MaxY); i f (graphresult() ) return (-2); setviewport (1, 1, (GData->MaxX - 1), (GData->MaxY - 1), 1); i f (graphresult() ) return (-3); return (0); } / * end DrawScreenFrame * / int DrawAxes (const SysFrame *SData, const GraphParms *GData) { /* */ / * Purpose: Draws coordinate axes, i f visible, based on scaling * / / * ranges and screen dimensions passed in the input * / / * structures. * / /* */ / * Uses: GData, a pointer to a data structure of type GraphParms, * / / * and SData, a pointer to data type SysFrame. These * / / * structure types are defined in the header f i l e SIMUL.H. * / /* */ / * Returns: An integer value; returns 0 on successful completion, * / / * or an error code < 0 i f the procedure fa i l s for any * / / * reason. * / /* */ double Temp; / * Calculate the X-pixel coordinate of the Y axis, . . . * / Temp = SData->Display->HiLim[0] - SData->Display->LoLim[0]; i f (Temp != 0) Temp = - SData->Display->LoLim[0] * GData->MaxX / Temp; else return (-1); Appendix F. Simulation Software 137 /* . . . a n d p l o t the a x i s , i f v i s i b l e . */ i f ( (Temp >= 0) && (Temp <= GData->MaxX) ) { l i n e (Temp, 0, Temp, GData->MaxY); } /* endif */ i f (graphresult ( ) ) re tu rn (-2); /* Then, do the Y - p i x e l coordinate of the X a x i s . */ Temp = SData->Display->HiLim[l] - SData->Display->LoLim[l] ; i f (Temp != 0) Temp = SData->Display->HiLim[l] * GData->MaxY / Temp; e l s e r e t u r n (-3); /* . . . a n d p l o t the a x i s , i f v i s i b l e . */ i f ( (Temp >= 0) && (Temp <= GData->MaxY) ) { l i n e (0, Temp, GData->MaxX, Temp); } /* endif */ i f (graphresult ( ) ) re tu rn (-2); r e t u r n (0); } /* end DrawAxes */ i n t EraseGraphics (void) { /* „/ /* Purpose: C lears the current graphics screen v iewport , and r e - */ /* draws the coordinate axes, i f v i s i b l e . */ /* */ /* Uses: GData, a po in te r to a data s t ructure of type GraphParms, */ /* and SData, a po in te r to a data s t ruc ture of type SysFrame. */ /* These s t ruc ture types are def ined i n the header f i l e SIMUL.H. */ /* */ /* Returns: An integer va lue ; returns 0 on success fu l complet ion, */ /* or an er ror code < 0 i f the procedure f a i l s f o r any */ /* reason. */ /* */ c learv iewpor t ( ) ; i f (graphresult ( ) ) re tu rn (-1); r e t u r n (0); } /* end DrawScreenFrame */ i n t InteractiveChange (SysFrame *SData, const GraphParms *GData) { /* _ */ /* Purpose: Permits the user to adjust the i n t e r n a l system p a r a - */ /* meters, through d i r e c t data ent ry . The parameters are grouped */ /* i n t o four major c l a s s e s : P h y s i c s , S t a t e s , In tegrator , and */ Appendix F. Simulation Software 138 / * Display; any or a l l of these may be selected for editing. * / /* */ / * Uses: SData, a pointer to a data structure of type SysFrame. * / / * This structure type is defined in the header f i l e SIMUL.H. * / / * GData, a pointer to type GraphFrame (see SIMUL.H) is also used. * / /* */ / * Returns: An integer value; returns 0 on successful completion, * / / * or an error code < 0 i f the procedure fa i l s for any * / / * reason. The internal fields accessed by SData wi l l be modified * / / * i f parameters are changed; the function may be exited without * / / * changing any data. * / /* */ int i , ConfirmedExit = 0; enum {Phys = Key_P, States = Key_S, Integ = Key_I, Disp = Key_D}; char *AdjMenuO = "P = Physics, S = States, I = Integrator, D = Display, <Esc> = Quit, <Ret> = Keep Changes"; char *AdjMenul = "S = States, I = Integrator, D = Display, <Esc> = Quit, <Ret> = Keep Changes"; static PFrame *Physics = NULL; static IFrame *ICs = NULL; static MFrame *Math = NULL; static DFrame *Display = NULL; / * Function prototypes * / int UserPhysics (PFrame *Physics); int UserStates (IFrame *ICs); int UserMath (MFrame *Math, const IFrame *ICs); int UserDisp (DFrame *Display, const IFrame *ICs); / * Allocate local parameter frames, i f required, copying data from the input. * / i f (Physics == NULL) { i f ( (Physics = malloc (sizeof (PFrame) ) ) == NULL) return (-1); i f ( (Physics->Val = malloc (SData->Physics->Num * sizeof (double) ) ) == NULL) return (-1); } / * endif * / Physics->Num = SData->Physics->Num; Physics->Name = SData->Physics->Name; for (i = 0; i < Physics->Num; i++) { Physics->Val[i] = SData->Physics->Val[i]; } / * endfor * / i f (ICs == NULL) { i f ( (ICs = malloc (sizeof (IFrame) ) ) == NULL) return (-1); i f ( (ICs->Val = malloc (SData->ICs->Num * sizeof (double) ) ) == NULL) return (-1); Appendix F. Simulation Software i f ( (ICs->Mul = malloc (SData->ICs->Num * sizeof (double) ) ) == NULL) return (-1); } / * endif * / ICs->Num = SData->ICs->Num; ICs->Name = SData->ICs->Name; for ( i = 0; i < ICs->Num; i++) { ICs->Val[i] = SData->ICs->Val[i]; } / * endfor * / for ( i = 0; i < ICs->Num; i++) { ICs->Mul[i] = SData->ICs->Mul[i]; } / * endfor * / i f (Math == NULL) { i f ( (Math = malloc (sizeof (MFrame) ) ) == NULL) return (-1); i f ( (Math->ValO = malloc (SData->ICs->Num * sizeof (double) ) ) == NULL) return (-1); i f ( (Math->Vall = malloc (SData->ICs->Num * sizeof (double) ) ) == NULL) return (-1); } / * endif * / Math->Integr = SData->Math->Integr; Math->TStep = SData->Math->TStep; Math->EndT = SData->Math->EndT; for ( i = 0; i < ICs->Num; i++) { Math->ValO[i] = SData->Math->ValO[i]; Math->Vall[i] = SData->Math->Vall[i]; } / * endfor * / i f (Display == NULL) { i f ( (Display = malloc (sizeof (DFrame) ) ) ™ NULL) return (-1); i f ( (Display->Index = malloc (SData->Display->Num * sizeof (int) ) ) == NULL) return (-1); i f ( (Display->LoLim = malloc (SData->Display->Num * sizeof (double) ) ) == NULL) return (-1); i f ( (Display->HiLim = malloc (SData->Display->Num * sizeof (double) ) ) == NULL) return (-1)j } / * endif * / Display->Num = SData->Display->Num; for (i = 0; i < Display->Num; i++) { Display->Index[i] = SData->Display->Index[i]; Display->LoLim[i] = SData->Display->LoLim[i]; Appendix F. Simulation Software Display->HiLim[i] = SData->Display->HiLim[i]; Display->OFile = SData->Display->OFile; } / * endfor * / / * Fl ip back to text mode to do the user interaction. * / restorecrtmodeO; / * Loop unt i l Quit or Exit. * / while (! ConfirmedExit) { / * Put up the menu. * / i f (Physics->Num != 0) puts (AdjMenuO); else puts (AdjMenul); fflush (stdin); switch (GetScanCode (0) ) { case Phys: i f (Physics->Num != 0) { i f (UserPhysics (Physics) < 0) return (-2); } / * endif * / break; case States: i f (UserStates (ICs) < 0) return (-3); break; case Integ: i f (UserMath (Math, ICs) < 0) return (-4); break; case Disp: i f (UserDisp (Display, ICs) < 0) return (-5); break; case Esc: / * Do you really want to lose changes? * / puts ("Exiting edit mode; parameter changes wi l l be lost. Press Y to confirm..."); fflush (stdin); i f (GetScanCode (0) == Yes) ConfirmedExit = 1; break; case Enter: / * Do you really want to keep changes? * / puts ("Exiting edit mode; parameter values wi l l be updated. Press Y to confirm..."); fflush (stdin); i f (GetScanCode (0) == Yes) ConfirmedExit = 1; Appendix F. Simulation Software 141 / * Copy pointers into the return structure. * / SData->Physics = Physics; SData->ICs = ICs; SData->Math = Math; SData->Display = Display; break; default: / * Nothing * / ; } / * endswitch * / } / * endwhile * / / * Return to graphics mode before exiting. * / setgraphmode (GData->GraphMode); return (0); } / * end InteractiveChange * / int UserPhysics (PFrame *Physics) { int i , j ; char *InString, *TString; / * Display the current values, . . . * / for (i = 0; i < Physics->Num; i++) { printf ("'/.s = '/.g ", Physics->Name[i] , Physics->Val [i] ) ; i f (! ( (i + 1) V. 4) ) puts (""); } / * endfor * / puts ("\n"); / * . . .allocate temporary space, . . . * / i f ( (InString = malloc (255 * sizeof (char) ) ) == NULL) return (-1); TString = InString; / * . . .and query the user for new values. * / do { puts ("Please enter new values for a l l of the Physics parameters."); fflush (stdin); gets (InString); puts (""); / * Parse the input string into blank-separated f ields, keeping a count. * / j = 0; i f (strtok (InString, " ") != NULL) j++; for (i = 1; i < Physics->Num; i++) { i f (strtok (NULL, " ") != NULL) j++; } / * endfor * / / * There must be enough fields available. * / i f (j == Physics->Num) { for (i = 0; i < Physics->Num; i++, InString++) Appendix F. Simulation Software 142 { /* Convert the s t r i n g to numeric va lues . */ Phys ics ->Va l [ i ] = s t r t o d ( InSt r ing , ftlnString); } /* endfor */ } /* endif */ } whi le (j != Physics->Num); f r e e (TSt r ing ) ; r e t u r n (0) ; } /* end UserPhysics */ i n t UserStates (IFrame *ICs) { i n t i , j ; char * InStr ing , *TStr ing; /* Disp lay the current va lues , . . . */ f o r ( i = 0 ; i < ICs->Num; i++) { p r i n t f 07.s = '/.g " , ICs->Name[i] , ICs ->Va l [ i ] ) ; i f (! ( ( i + 1) '/. 4) ) puts (""); } /* endfor */ puts ("\n"); /* . . . a l l o c a t e temporary space, . . . */ i f ( ( InSt r ing = malloc (255 * s i zeof (char) ) ) == NULL) re tu rn ( - 1 ) ; TSt r ing = I n S t r i n g ; /* . . . a n d query the user f o r new va lues . */ do { puts ("Please enter new i n i t i a l values f o r a l l the system s t a t e s . " ) ; f f l u s h ( s t d i n ) ; gets ( InSt r ing ) ; puts (""); /* Parse the input s t r i n g i n t o b lank-separated f i e l d s , keeping a count. */ j = 0 ; i f ( s t r tok ( InSt r ing , " ") != NULL) j++; f o r ( i = 1; i < ICs->Num; i++) { i f ( s t r tok (NULL, " ") != NULL) j++; } /* endfor */ /* There must be enough f i e l d s a v a i l a b l e . */ i f (j == ICs->Num) { f o r ( i = 0; i < ICs->Num; i++, InString++) { /* Convert the s t r i n g to numeric va lues . */ ICs ->Val [ i ] = s t r t o d ( InSt r ing , ftlnString); } /* endfor */ } /* endif */ } whi le (j != ICs->Num); Appendix F. Simulation Software f r e e (TSt r ing ) ; r e t u r n (0) ; } /* end UserStates */ i n t UserMath (MFrame *Math, const IFrame *ICs) { i n t i ; /* Disp lay the current va lues , . . . */ /* puts ("Absolute e r ro r to lerance f o r : " ) ; f o r ( i = 0 ; i < (ICs->Num - 1 ) ; i++) { p r i n t f C7.s = */.g " , ICs->Name[i+l] , Math->ValO[i]) ; i f (! ( ( i + 1) */. 4) ) puts (""); } /* endfor */ puts ("\n"); puts ( "Re lat i ve e r ro r to lerance f o r : " ) ; f o r ( i = 0 ; i < (ICs->Num - 1 ) ; i++) { p r i n t f 07.s = '/.g " , ICs->Name[i+l] , Math->Vall [ i ] ) ; i f (! ( ( i + 1) '/. 4) ) puts (""); } /* endfor */ puts ("\n">; */ p r i n t f ( " Integrator e r ro r to lerance (Absolute and Re la t i ve ) : Math->ValO[0]); /* . . . a n d query the user f o r new va lues . */ do { puts ("Please enter new to lerance v a l u e . " ) ; f f l u s h ( s t d i n ) ; } whi le ( (1 != scanf ("'/.If", &(Math->ValO[0]) ) ) II (Math->Val0[0] < 0) ) ; Math->Vail[0] = Math->Val0[0]; f o r ( i = 1; i < (ICs->Num - 1 ) ; i++) { Math->ValO[i] = Math->Val0[0]; Math->Val l [ i ] = Math->Val0[0]; } /* endfor */ do { p r i n t f ("Current endpoint value i s '/,g; p lease enter new va lue . (May be negat ive)\n" , Math->EndT); f f l u s h ( s t d i n ) ; } whi le (1 != scanf ('"/.If", ft(Math->EndT) ) ) ; re tu rn (0 ) ; } /* end UserMath */ Appendix F. Simulation Software 144 i n t UserDisp (DFrame *Display, const IFrame *ICs) { i n t i , j ; char **TName, * InSt r ing , *TStr ing; FILE *fopenWarn (char *FileName, char *Mode, char *0utWarn); i f ( (TName = mal loc (3 * s i zeof (char*) ) ) == NULL) re tu rn ( - 1 ) ; i f ( ( InSt r ing = malloc (255 * s i zeof (char) ) ) == NULL) re tu rn ( - 1 ) ; TS t r ing = I n S t r i n g ; /* Show current e n t r i e s . */ p u t s ( " A v a i l a b l e va r iab les are Time, X, XDot, Theta, and ThetaDotAn"); puts ("Current ly d isp layed v a r i a b l e s : " ) ; f o r ( i = 0 ; i < Display->Num; i++) { p r i n t f ("'/.s: Range = (*/.g, '/.g)\n", ICs->Name[Display->Index[i]] , D isp lay ->LoL im[ i ] , D isp lay ->HiL im[ i ] ) ; do { puts ("Enter new d isp lay v a r i a b l e , as : Name Lower -Sca le -L imi t U p p e r - S c a l e - L i m i t " ) ; f f l u s h ( s t d i n ) ; gets ( InSt r ing ) ; TName[0] = s t r t o k ( InSt r ing , " " ) ; TName[1] = s t r t o k (NULL, " " ) ; TName[2] = s t r t o k (NULL, " " ) ; Display->Index[ i ] - - 1 ; f o r (j = 0 ; j < ICs->Num; j++) { i f (! strcmp (ICs->Name[j], TName[0]) ) Display->Index[ i ] = j ; } /* endfor */ Display->LoLim[i ] = atof (TName[l]); Display ->HiL im[ i ] = atof (TName[2]); } whi le ( (Display->Index[i ] < 0) I I (TName[l] == NULL) I I (TName [2] == NULL) I I (Display->LoLim[i] == Disp lay ->HiL im[ i ] ) ) ; } /* endfor */ i f (D isp lay ->0Fi le != NULL) f c l o s e (D isp lay ->OFi le ) ; puts ("\nWrite t r a j e c t o r y data to output f i l e ? (Y/N)"); f f l u s h ( s t d i n ) ; do Appendix F. Simulation Software 145 { j = ( i = GetScanCode (0) ) == Yes; } whi le ( ! j && ( i != No) ) ; i f ( j ) { puts ("Enter name f o r output f i l e . " ) ; f f l u s h ( s t d i n ) ; gets (TSt r ing) ; Display->0File = fopenWarn (TSt r ing , "wt" , "WarnOverWrite"); } /* endif */ e lse Display->0File = NULL; f r e e (TName); f r e e (TSt r ing ) ; re tu rn (0); } /* end UserDisp */ i n t PopUpGraphMsg (Message *InMsg) { /* */ /* This f u n c t i o n takes an (X, Y) l o c a t i o n i n LocX and LocY, and a */ /* message s t r i n g , pointed to by Msg. I t determines the s i z e of an */ /* area requ i red to d i sp lay Msg, a l l o c a t e s memory f o r the c o r r e - */ /* sponding b i t -map , saves the current contents of the requ i red area ,*/ /* and d i sp lays the message on-screen. LocX and LocY def ine the */ /* center of the message box; the f u n c t i o n f a i l s i f the r e s u l t i n g */ /* tex t box would go o f f - s c r e e n . The f u n c t i o n returns a po in te r to */ /* the saved screen bit -map on success fu l completion; on f a i l u r e , */ /* NULL i s returned. On e x i t , the f i r s t two words of the area */ /* pointed to by the returned po in te r conta in the (X, Y) coordinates */ /* of the upper l e f t - h a n d corner of the saved area . Use ClearGraph- */ /* Msg, below, to restore the bit -map and d e - a l l o c a t e the b u f f e r . */ /* THESE FUNCTIONS OPERATE IN GRAPHICS MODE ONLY! NOTE: Repeated */ /* c a l l s to PopUpGraphMsg, without c a l l i n g ClearGraphMsg, or us ing */ /* _ f r e e , w i l l eventual ly f i l l a l l a v a i l a b l e memory. WARNING: */ /* DON'T DO IT!! */ s t r u c t tex tse t t ings type ETSet; s t ruc t viewporttype EPSet; i n t MsgCenX, MsgCenY, TextXMin, TextXMax, TextYMin, TextYMax; /* Get the current s e t t i n g s , . . . */ ge t tex tse t t ings (ftETSet); getv iewsett ings (ftEPSet); /* . . .change to new s e t t i n g s , . . . */ s e t t e x t j u s t i f y (CENTER.TEXT, CENTER.TEXT); Appendix F. Simulation Software 146 /* . . . a n d f i g u r e out how b i g the l o c a l viewport should be. */ MsgCenX = (8 + textwidth (InMsg->Msg) ) / 2; MsgCenY = (6 + textheight (InMsg->Msg) ) / 2; TextXMin = InMsg->LocX - MsgCenX; TextXHax = InMsg->LocX + MsgCenX; TextYMin = InMsg->LocY - MsgCenY; TextYMax = InMsg->LocY + MsgCenY; /* Set i t , and check s t a t u s ; e x i t i f tex t window i s i n v a l i d . */ setviewport (TextXMin, TextYMin, TextXMax, TextYMax, 1 ) ; i f (graphresult ( ) ) re tu rn ( - 1 ) ; /* A l l o c a t e the window save a rea , i f necessary, . . . */ i f (InMsg->BkGnd == NULL) { i f ( (InMsg->BkGnd = mal loc (s i zeof (BFrame) ) ) == NULL) re tu rn ( - 2 ) ; i f ( (InMsg->BkGnd->Map = malloc (imagesize (0, 0 , (2 * MsgCenX), (2 * MsgCenY) ) ) ) == NULL) re tu rn ( - 2 ) ; } /* endif */ /* . . . p r e s e r v e the window in format ion , . . . */ InMsg->BkGnd->XULHC = TextXMin; InMsg->BkGnd->YULHC = TextYMin; getimage (0, 0 , (2 * MsgCenX), (2 * MsgCenY), InMsg->BkGnd->Map); /* . . . a n d w r i t e the message s t r i n g to the screen. */ c learv iewpor tO ; rectangle (0, 0 , (2 * MsgCenX), (2 * MsgCenY) ) ; outtextxy ( (MsgCenX + 1 ) , (MsgCenY + 1 ) , InMsg->Msg); /* Put a l l the parameters back the way we found them, and e x i t . */ s e t t e x t j u s t i f y (ETSet .hor iz , ETSet . ve r t ) ; setviewport ( E P S e t . l e f t , E P S e t . t o p , E P S e t . r i g h t , E P S e t . b o t t o m , E P S e t . c l i p ) ; r e t u r n (0) ; } /* end PopUpGraphMsg */ i n t ClearGraphMsg (Message *InMsg) { /* */ /* Used wi th PopUpGraphMsg, above, t h i s f u n c t i o n res to res the saved */ /* screen b i t -map , and d e - a l l o c a t e s the b u f f e r which h e l d i t . The */ /* f u n c t i o n returns 0 f o r successfu l completion, or the value of */ /* . g r a p h r e s u l t , on f a i l u r e . */ s t r u c t viewporttype EPSet; i n t BitMapX, BitMapY, BitMapWidth, BitMapHeight; Appendix F. Simulation Software 147 / * Get the current viewport settings, . . . * / getviewsettings (ftEPSet); / * . . . p u l l the raster coordinates for the bit-map out of the input data, . . . * / BitMapX = InMsg->BkGnd->XULHC; BitMapY = InMsg->BkGnd->YULHC; BitMapWidth = *(InMsg->BkGnd->Map); BitMapHeight = *(InMsg->BkGnd->Map + 1); / * . . .set up the viewport as used for the raster-save, . . . * / setviewport (BitMapX, BitMapY, (BitMapX + BitMapWidth), (BitMapY + BitMapHeight), 1); i f (graphresult() ) return (-11); / * . . .and restore the bit-map; direct screen copy. * / putimage (0, 0, InMsg->BkGnd->Map, COPY.PUT); / * Now, put the parameters back the way we found them, and exit. * / setviewport (EPSet.left,EPSet.top,EPSet.right,EPSet.bottom,EPSet.clip); return (graphresult() ); } / * end ClearGraphMsg * / int PlotResults (const double *InState, const SysFrame *SData, const GraphParms *GData) { double TempO, Tempi; int Color; / * Calculate the X-pixel coordinate, . . . * / TempO = SData->Display->HiLim[0] - SData->Display->LoLim[0]; i f (TempO != 0) { TempO = (Instate[SData->Display->Index[0]] * SData->ICs->Mul[SData->Display->Index[0]] - SData->Display->LoLim[0]) * GData->MaxX / TempO; } else return (-1); / * . . .then, do the Y-pixel coordinate, . . . * / Tempi = SData->Display->HiLim[l] - SData->Display->LoLim[l]; i f (Tempi != 0) { Tempi = (SData->Display->HiLim[l] -Instate[SData->Display->Index[1]] * SData->ICs->Mul[SData->Display->Index[1]]) * GData->MaxY / Tempi; } else return (-2); Appendix F. Simulation Software /* . . . a n d p l o t the p o i n t , i f v i s i b l e . */ i f ( (TempO >= 0) && (TempO <= GData->MaxX) && (Tempi >= 0) (Tempi <= GData->MaxY) ) { Color - getco lor () ; p u t p i x e l (TempO, Tempi, C o l o r ) ; } /* endif */ i f (graphresult ( ) ) re tu rn (-3); r e t u r n (0); } /* end P lo tResu l t s */ vo id MayDump ( in t Code, const char *Func) { i f (Code < 0) { closegraph.0 ; p r i n t f ( "Fa ta l e r r o r , code '/,d, i n y ,s .\n", Code, Func) ; a b o r t ( ) ; } /* endif */ r e t u r n ; } /* end MayDump */ Appendix F. Simulation Software 149 F.7 Fi le I / O : F I L E P R O C . C /* */ /* FILEPROC.C — F i l e handling/manipulation r o u t i n e s . */ /* */ /* Design & coding: M.G. Henders */ /* Last modi f ied : 23 A p r i l 1991 */ /* */ /* Module type: Support module. */ /* Module contents : fopenWarn */ /* */ /* Support f i l e s : None */ /* */ #include <compiler.h> #include <stdio.h> #include <string.h> Wdefine NoAccess access /* Cond i t iona l dec la ra t ions */ #ifdef Turbo /* For Turbo C */ #include <io.h> #include <keystrks.h> #endif #ifdef NDP /* For NDP C */ #include <grex.h> #include <stdl ib.h> #include <keycodes.h> #endif /* WARNING: The NDP C-860 compiler apparently does not recognize " t " */ /* or "b" i n the _fopen mode s t r i n g , l i k e Turbo; t h e r e f o r e , s ince i t */ /* de fau l t s to text-mode f i l e t r a n s l a t i o n , i t i s apparently impossible*/ /* to do binary-mode f i l e I/O wi th t h i s compi ler . The manual r e f e r s */ / • t o _pmode as a v a r i a b l e to s e t , to change the f i l e mode, (def ined */ /* i n s t d i o . h ) , but t h i s i d e n t i f i e r i s not found when attempting to */ /* compile. Examination of s t d i o . h shows i t i s apparently i f d e f ' e d */ /* out f o r the MSDOS compi ler , as i s the setmode f u n c t i o n , normally */ /* used to change the f i l e mode. Consequently, t h i s f u n c t i o n w i l l */ /* probably not work as intended i f used f o r b inary I/O under the NDP */ /* compi ler ; recommend text-mode I/O ONLY. A l s o , the NDP vers ion of */ /* _st rupr i s f a u l t y , and, as i t i s not c r i t i c a l to the f u n c t i o n , i t */ /* i s not used i n the NDP compi la t ion . */ Appendix F. Simulation Software 150 FILE *fopenWarn (char *FileName, char *Mode, char *OutWarn) { /* */ /* Essentially, this function duplicates the standard _fopen */ /* function, but i t has the additional capability of generating */ /* warning messages and prompts under certain conditions. In */ /* particular, a warning w i l l always be generated i f Mode indicates */ /•an open-for-read, and the f i l e cannot be opened. A prompt for a */ /* new filename w i l l then appear. On open-for-write, the warning */ /* mode i s controlled by OutWarn; pass the string "WarnOverWrite" to */ /* generate a warning before over-writing an existing f i l e , and pass */ /* "WarnNotFound" to warn before creating a new output f i l e . */ /• NOTE: If OutWarn does not match one of these strings, the */ /• default open-for-write i s executed, providing no warnings at a l l . */ /* The function returns a f i l e pointer for use i n subsequent f i l e */ /* ops, and the filename f i n a l l y used i s passed back i n FileName. */ /* */ FILE *L F i l e = NULL; int YesKey, Key; /* Check for any of the open-for-output modes. */ i f (strpbrk (Mode, "aw+") != NULL) { /* Some kind of output may take place; check warning mode. */ i f (!strcmp (OutWarn, "WarnOverWrite") ) { /* Test for the existence of the f i l e , ... */ i f (NoAccess (FileName, 0) ) {. /* F i l e does not exist; go ahead and create i t . */ L F i l e = fopen (FileName, Mode); } /* endif */ else { /* The f i l e already exists; check before over-writing i t . */ do { #ifdef Turbo /* For Turbo C */ strupr (FileName); #endif printf ("File '/.s already exists; over-write i t ? (Y/N)\n", FileName); #ifdef Turbo /* For Turbo C */ do { YesKey = (Key = GetScanCode (0) ) == Yes; } while (!YesKey ftft (Key != No) ); #endif #ifdef NDP /* For NDP C */ do { Key = pauseQ; YesKey = (Key == yes) I I (Key == Yes); Appendix F. Simulation Software 151 } whi le (!YesKey k& (Key != no) kk (Key != No) ) ; #endif i f (YesKey) { p r i n t f ( " E x i s t i n g f i l e data w i l l be o v e r - w r i t t e n . \ n " ) ; } /* endif */ e lse { p r i n t f ("Enter new f i l e n a m e . \ n " ) ; scanf 07.s", FileName); p r i n t f ("\n") ; } /* endelse */ } whi le ( (!NoAccess (FileName, 0) ) kk (!YesKey) ) ; /* Do the a c t u a l _fopen c a l l . */ L F i l e = fopen (FileName, Mode); } /* endelse */ } /* endif */ e lse i f (Istrcmp (OutWarn, "WarnNotFound") ) { /* Test f o r the existence of the f i l e , . . . */ i f (NoAccess (FileName, 0) ) { /* F i l e does not e x i s t ; check before c rea t ing i t . */ do { #ifdef Turbo /* For Turbo C */ st rupr (FileName); #endif p r i n t f ( " F i l e */.s does not e x i s t ; create i t ? (Y/N)\n", FileName); #ifdef Turbo /* For Turbo C */ do { YesKey = (Key = GetScanCode (0) ) == Yes; } whi le (!YesKey kk (Key != No) ) ; #endif #ifdef NDP /* For NDP C */ { Key = pause() ; YesKey = (Key == yes) I I (Key == Yes) ; } whi le (!YesKey && (Key != no) kk (Key != No) ) ; #endif i f (!YesKey) { p r i n t f ("Enter new f i l e n a m e . \ n " ) ; scanf 07.s", Fi leName); p r i n t f ("\n") ; } /* endelse */ } whi le ( (NoAccess (FileName, 0) ) kk (!YesKey) ) ; /* Do the a c t u a l _fopen c a l l . */ Appendix F. Simulation Software LFile = fopen (FileName, Mode); } / * endif * / else { / * The f i l e already exists; in this mode, go ahead and over-write. * / #ifdef Turbo / * For Turbo C * / strupr (FileName); #endif printf ("Over-writing data in f i l e '/,s.\n", FileName); LFile = fopen (FileName, Mode); } / * endelse * / } / * endif * / else { / * Normal _fopen ca l l ; no warnings, returns NULL on failure. * / LFile = fopen (FileName, Mode); } / * endelse * / } / * endif * / else { / * Open-for-input ca l l ; bitch to the user i f _fopen doesn't work. #ifdef Turbo / * For Turbo C * / strupr (FileName); #endif for ( ; ( (LFile = fopen (FileName, Mode) ) == NULL); ) {. printf ("Unable to open '/,s; please supply a new filename.\n", FileName); scanf ('"/.s", FileName); #ifdef Turbo / * For Turbo C * / strupr (FileName); #endif printf ("\n"); } / * endfor * / } / * endelse * / return (LFile); } / * end fopenWarn * / Appendix G Stability Region Computation Software This Appendix comprises listings of software used for computation of the stability regions of the inverted pendulum system in four-dimensional state space; the program listed here requires the system dynamics, numerical integration, and file I/O modules that comprise Sections F.5, F.4, and F.7 respectively. The comments made in Appendix F, regarding page formatting and coding structure, apply equally to this code. The main program source file is STABVOL.C (Section G.2); in addition to the mod-ules already mentioned, the header file VOLUME.H and the module SMPLPROC.C are required, and are listed in Sections G.l and G.3 respectively. This code was originally developed under Turbo C [BOR], and a project control file for this compiler is included in Section G.l , but it should be noted that STABVOL will not actually run to completion when compiled under Turbo C, or, most likely, under any DOS real-mode compiler; these compilers cannot access enough memory for the program data. For this reason, as well as for speed, the code was ported to NDP-C860 [NDP], and run on a Number Smasher 860 processor board built by Micro Way, Inc.; this provided adequate memory. The code presented here will compile and run under either compiler, but will simply exit early under Turbo C. The user should also be aware that a run of STABVOL on a 40 MHz i860 requires approximately 76 hours to complete; based on speed comparisons made using STABVOL itself, the same run would require around 230 hours on a 33 MHz 80386/80387, given a suitable compiler. Use of STABVOL is essentially self-explanatory; on startup, the program simply 153 Appendix G. Stability Region Computation Software 154 requests a name for the data file to use for output. All other parameters are currently hard-coded into the program, but are readily found in the initial definition blocks of STABVOL.C. Parameters of particular interest will be the Seed and SpaceSize values, as well as PtsPerAxis. The Seed values (SeedW, SeedX, SeedY, and SeedZ) specify the initial conditions vector for the point at which the stabilty mapping is to be started; note that to actually accomplish mapping, the specified state must lie within the region of stability. The SpaceSize and PtsPerAxis values are related to the fact that the state space, during stability mapping, is sampled at fixed intervals, and the sample points are tracked by a set of index numbers, one for each space coordinate axis. In this implementation of STABVOL, the index numbers are assigned only 8 bits each, to conserve memory, and, in consequence, a maximum of 256 sample values can be identified along each sampling axis. The size of each axis sampling interval is determined by both the SpaceSize value for that axis, and by PtsPerAxis; specifically, GridSize = 2 * SpaceSize / PtsPerAxis. Asa general rule, the SpaceSize values should be set to the approximate widths of the stability region to be mapped, along each mapping axis, and PtsPerAxis should then be adjusted to give the desired sampling resolution for the overall stability region. Note, however, that if PtsPerAxis is set to its maximum allowable value (508, as implemented), any extension of the stability region outside the hypercube defined by the SpaceSize values will create immediate problems with overflowing sample index values. (The allowable maximum PtsPerAxis value is around 512, rather than 256, because each sample index value actually identifies two sample points, not just one.) For this reason, (since the exact extent of the stability region is presumably not known in advance) it is preferable to set the SpaceSize values somewhat larger than the expected size of the stability region, and to use a PtsPerAxis value somewhat lower than the maximum permitted. Appendix G. Stability Region Computation Software 155 STABVOL is a graphics-oriented program; as the stability region is mapped, a rep-resentation of each stable point is plotted on the graphics screen, in a two-dimensional projection of the state space. The scaling of this plot is controlled by the definition of additional values in the STABVOL header, AbsciMinDef, AbsciHaxDef, OrdinMinDef, and OrdinMaxDef; these values set the boundaries of the displayed portion of the state space. Appendix G. Stability Region Computation Software 156 G.l Ancillary Files G.l.l Project Control File: STABVOL.PRJ stabvol.c (keystrks.h, simul.h, volume.h) smplproc.c (simul.h, volume.h) ntgrproc.c dynamics.c (simul.h) fileproc.c (keystrks.h) G.1.2 Volumetric Definitions: VOLUME.H #define p i 3.1415926535 #define MaxIndexVal 127 s t r u c t VStruct {union {char Index[4]; long Key;} IK; /* S p a t i a l l o c a t i o n i n d i c e s */ /* Key value f o r searching */ unsigned SCode;}; /* Vertex s t a b i l i t y code */ typedef s t r u c t VStruct Voxe l ; s t r u c t AStruct { in t i n t i n t i n t DPtr ; D i r ; MaxC; CenC; Min ; Max; Len; Cen; /* /* /* /* /* /* /* /* /* Selector index f o r data value */ Data->Screen conversion (+/- 1) */ Screen coordinate maximum */ Screen coordinate centre */ Data value minimum */ Data value maximum */ Data value range */ Data value centre po int */ Magn i f i ca t ion f a c t o r */ double double double double double Mag;}; typedef s t r u c t AStruct A x i s ; #define sqr(x) ((x)*(x)) Appendix G. Stability Region Computation Software 157 G .2 D r i v e r : S T A B V O L . C /* */ / * STABVOL.C — Main program module for stability volume computation. * / /* */ / * Design ft coding: M.G. Henders * / / * Last modified: 8 June 1991 * / /* */ / * Module type: Externally-supported MAIN * / / * Module contents: Main * / /* */ / * Project f i l e : STABVOL.PRJ * / / * Support modules: NTGRPROC.C, DYNAMICS.C, FILEPROC.C * / #define Main main /* */ #include <compiler.h> #include <stdio.h> #include <assert.h> #include <math.h> / * Special data declarations * / #include <simul.h> #include <volume.h> #define View 1 #if ( (View < 0) I I (View > 2) ) #error Invalid value for View; must be 0 <= View <= 2. #endif #define PtsPerAxis 300 #if ( (PtsPerAxis <= 0) II (PtsPerAxis > (4 * MaxIndexVal) ) ) #error Invalid resolution; use 0 <= PtsPerAxis <= (4 * MaxIndexVal). #endif / * In i t ia l state at which to begin stability region mapping * / #define SeedW 0.0 #define SeedX 0.0 #define SeedY 0.0 tfdefine SeedZ 0.0 / * Nominal dimensions of the space of i n i t i a l conditions that is to be * / / * mapped for stabil ity. These dimensions must be chosen such that no * / / * part of the stability region being mapped extends beyond the walls * / / * of the hypercube of size defined by SpaceSizeW to SpaceSizeZ; the * / / * voxel counting indices wi l l overflow i f this condition is not met, * / / * leading to unpredictable results. * / Appendix G. Stability Region Computation Software 158 #define SpaceSizeV 600.0 #define SpaceSizeX 250.0 #define SpaceSizeY 800.0 tfdefine SpaceSizeZ 1800.0 / * Dimensions of the mapped space, for use in scaling the on-screen * / / * graphics display. These are not particularly c r i t i c a l , but should * / / * generally be chosen to complement the Seed and SpaceSize values in * / / * use for the stability mapping. * / #define AbsciMinDef -375.0 #define AbsciMaxDef 375.0 #define OrdinMinDef -900.0 #define OrdinMaxDef 900.0 / * Conditional declarations * / #ifdef Turbo / * For Turbo C * / #ifndef __HUGE__ #error Requires HUGE memory model. #endif #include <alloc.h> #include <graphics.h> #include <conio.h> #include <dos.h> #include <keystrks.h> extern unsigned _stklen = 0x4000; #endif #ifdef NDP / * For NDP C * / #include <grex.h> #include <stdlib.h> #include <keycodes.h> #endif double GridSize[4]; long MaxVoxListSize; FILE *OutputFile; void Main() { int SysSize, i , GraphMode, BCode, MaxX, MaxY, CenX, CenY; SysFrame SystemData; double xO, yO, *Seed; char *0utputFileName; Axis Absci, Ordin; long HemSize, CurrVox = 0, LastVox = 0; Appendix G. Stability Region Computation Software 159 #ifdef Turbo / * For Turbo C * / int CursorX, CursorY, GraphDriver = DETECT; extern int _8087; char ScreenBuf [80*25*2], *GraphDriverPath = " C : \ \ T C \ \ " ; Voxel huge*VoxList; / * Functions defined in module SMPLPROC.C * / unsigned StableVoxel (Voxel huge*VPtr, int NStates, double *Seed); void DrawVoxelOutline (Voxel huge*Vox, double *Seed, Axis *Absci, Axis *0rdin, int Erase); void SaveVolume (FILE *0File, Voxel huge*VoxList, long LastVox, double *Seed); void BuildVolume (Voxel huge*VoxList, long CurrVox, long *LastVox, int SysSize, double *Seed, Axis *Absci, Axis *0rdin); #endif #ifdef NDP / * For NDP C * / int GraphConfig[4], Dummy; Voxel *VoxList; / * Functions defined in module SMPLPROC.C * / unsigned StableVoxel (Voxel *VPtr, int NStates, double *Seed); void DrawVoxelOutline (Voxel *Vox, double *Seed, Axis *Absci, Axis *0rdin, int Erase); void SaveVolume (FILE *0File, Voxel *VoxList, long LastVox, double *Seed); void BuildVolume (Voxel *VoxList, long CurrVox, long *LastVox, int SysSize, double *Seed, Axis *Absci, Axis *0rdin); #endif / * Functions defined in module DYNAMICS.C * / int SysModuleInfo (SysFrame *SData); void RunConstants (double *PVals); / * Functions defined in module FILEPROC.C * / FILE *fopenWarn (char *FileName, char *Mode, char *0utWarn); / * Create the skeleton of the system parameter structure. * / assert((SystemData.Physics = (PFrame*)malloc(sizeof(PFrame))) != NULL); assert((SystemData.ICs = (IFrame*) malloc (sizeof(IFrame))) != NULL); assert((SystemData.Math = (MFrame*) malloc (sizeof(MFrame))) != NULL); assert((SystemData.Display = (DFrame*)malloc(sizeof(DFrame))) != NULL); / * Get parameters from the Dynamics module; abort i f ca l l f a i l s . * / SystemData.Display->Num = 2; assert (SysModuleInfo (ftSystemData) >= 0); SysSize = SystemData.ICs->Num; Appendix G. Stability Region Computation Software 160 / * Init ial ize secondary parameters. * / RunConstants (SystemData.Physics->Val); / * Allocate state vector and filename space. * / assert ((Seed = (double*)malloc(SysSize * sizeof(double))) != NULL); assert ((OutputFileName = (char*)malloc(255 * sizeof(char))) != NULL); / * Try to maximize the number of voxels that can be created. * / #ifdef Turbo / * For Turbo C * / HemSize = coreleftO; HaxVoxListSize = 0.90 * HemSize / sizeof (Voxel); assert ((VoxList = (Voxel huge*) farmalloc (HaxVoxListSize * sizeof(Voxel))) != NULL); #endif #ifdef NDP / * For NDP C * / HaxVoxListSize = 1000000; do { VoxList = (Voxel*) malloc (MaxVoxListSize * sizeof (Voxel) ); i f (VoxList == NULL) HaxVoxListSize *= 0.99; } while (VoxList == NULL); #endif printf ("There is sufficient memory for '/,ld voxels. \n\n", MaxVoxListSize); / * Init ial ize the state vector and voxel sizes. * / Seed[0] = 0.0; Seed[l] = SeedW; Seed[2] = SeedX; Seed[3] = SeedY; Seed[4] = SeedZ; GridSize[0] = 2.0 * SpaceSizeW / PtsPerAxis; GridSize[l] = 2.0 * SpaceSizeX / PtsPerAxis; GridSize[2] = 2.0 * SpaceSizeY / PtsPerAxis; GridSize[3] = 2.0 * SpaceSizeZ / PtsPerAxis; / * Get an output f i l e for the data, . . . * / printf ("Enter name of output f i le . \n"); scanf ("'/.s", OutputFileName); fflush (stdin); assert ( (OutputFile = (FILE*) fopenWarn (OutputFileName, "w", "WarnOverWrite") ) != NULL); printf ("Data wi l l be written to f i l e '/,s.\n\n", OutputFileName); / * . . . i n i t i a l i z e the voxel l i s t , . . . * / VoxList[0].IK.Key = 0; VoxList[0].SCode = StableVoxel (ftVoxList[0], SysSize, Seed); Appendix G. Stability Region Computation Software /* . . . a n d i n i t i a l i z e d i sp lay parameters. */ Absc i .M in = AbsciHinDef; Absci.Max = AbsciMaxDef; Ordin.Min = OrdinMinDef; Ordin.Max = OrdinMaxDef; Absc i .Len = Absci.Max - A b s c i . M i n ; Ordin.Len = Ordin.Max - Ord in .Min ; Absci .Cen = 0 .5 * (Absci .Min + Absci .Max) ; Ordin.Cen = 0 .5 * (Ordin.Min + Ordin.Max); swi tch (View) { case 0 : A b s c i . D i r = 1; Ordin.DPtr = 0 ; case 1: A b s c i . D i r = 1; Ordin.DPtr = 3 ; case 2: A b s c i . D i r = - 1 ; Ordin.DPtr = 3 ; Ord in .D i r = break; Ord in .D i r = break; Ord in .D i r = break; 1; Absc i .DPt r = 2 - 1 ; Absc i .DPt r = 2 - 1 ; Absc i .DPt r = 0 d e f a u l t : ; /* Nothing */ } /* endswitch */ #ifdef Turbo /* For Turbo C */ /* Save the text-mode s c r e e n s t a t e . */ CursorX = w h e r e x O ; CursorY = whereyO ; assert (gettext (1 , 1 , 80, 25, ScreenBuf) ) ; /* I n i t i a l i z e the graphics system. */ i n i t g r a p h (ftGraphDriver, ftGraphMode, GraphDriverPath); asser t (GraphDriver >= 0 ) ; /* I n i t i a l i z e s c r e e n p a r a m e t e r s . * / MaxX = getmaxxO ; MaxY = getmaxyO ; /* . . .d raw screen boundaries, . . . */ se tco lo r (1 ) ; l i n e (0, 0 , MaxX, 0 ) ; l i n e (MaxX, 0 , MaxX, MaxY); l i n e (MaxX, MaxY, 0 , MaxY); l i n e (0, MaxY, 0 , 0 ) ; #endif Appendix G. Stability Region Computation Software #ifdef NDP / * For NDP C * / / * Init ial ize the graphics system. * / GraphMode = video_configuration (GraphConfig); assert (graphics_mode (GraphMode) ); / * Init ial ize screen parameters. * / get_device_limits (&MaxX, ftMaxY, ftDummy); / * ...draw screen boundaries, . . . * / set.color (1); move (0, 0); draw (MaxX, 0); draw (MaxX, MaxY); draw (0, MaxY); draw (0, #endif / * Finish up the display parameters. * / CenX = MaxX / 2; CenY = MaxY / 2; Absci.MaxC = MaxX; Absci.CenC = CenX; Ordin.MaxC = MaxY; Ordin.CenC = CenY; / * . . .and draw the coordinate axes, i f within the viewport. * / Absci.Mag = MaxX / Absci.Len; Ordin.Mag = MaxY / Ordin.Len; xO = CenX - Absci.Dir * Absci.Mag * Absci.Cen; i f ( (xO >= 0) && (xO <= MaxX) ) { #ifdef Turbo / * For Turbo C * / line ( (int)xO, 0, (int)xO, MaxY); #endif #ifdef NDP / * For NDP C * / move ( (int)xO, 0); draw ( (int)xO, MaxY); #endif } / * endif * / yO = CenY - Ordin.Dir * Ordin.Mag * Ordin.Cen; i f ( (yO >= 0) ftft (yO <= MaxY) ) i #ifdef Turbo / * For Turbo C * / line (0, (int)yO, MaxX, (int)yO); #endif #ifdef NDP / * For NDP C * / move (0, (int)yO); draw (MaxX, (int)yO); #endif } / * endif * / Appendix G. Stability Region Computation Software / * Init ial ize screen viewport, with clipping enabled. * / #ifdef Turbo / * For Turbo C * / setviewport (0, 0, MaxX, MaxY, 1); #endif #ifdef NDP / * For NDP C * / set_clip_limits (0, 0, MaxX, MaxY); #endif / * Map the stability region, . . . * / do { / * Check for request to save accumulated voxel l i s t , or to exit. #ifdef Turbo / * For Turbo C * / i f ( (BCode = GetScanCode (1) ) == Esc) break; else i f (BCode) { SaveVolume (OutputFile, VoxList, LastVox, Seed); do { BCode = GetScanCode (0); } while ( (BCode = GetScanCode (1) ) != 0); } / * endelse * / #endif #ifdef NDP / * For NDP C * / i f ( (BCode = inkey$() ) == Esc) break; else i f (BCode) { SaveVolume (OutputFile, VoxList, LastVox, Seed); BCode = 0; } / * endif * / #endif / * Ignoring completely unstable voxels, . . . * / i f (VoxList[CurrVox].SCode) { / * ...show where we are working, . . . * / DrawVoxelOutline (ftVoxList[CurrVox], Seed, ftAbsci, ftOrdin, 0); / * . . .and expand from the current voxel, . . . * / BuildVolume (VoxList, CurrVox, ftLastVox, SysSize, Seed, ftAbsci, ftOrdin); } / * endif * / / * . . .then move to the next l i s t entry. * / CurrVox++; / * Check for l i s t overflow. * / i f (LastVox >= MaxVoxListSize) BCode = 2; } while (!BCode && (CurrVox <= LastVox) ); Appendix G. Stability Region Computation Software / * . ..and write the whole thing to the output f i l e . * / SaveVolume (OutputFile, VoxList, LastVox, Seed); / * Button up the f i l e , and return to text mode before leaving, fclose (OutputFile); #ifdef Turbo / * For Turbo C * / closegraphO ; ass.ert (puttext (1, 1, 80, 25, ScreenBuf) ); gotoxy (CursorX, CursorY); #endif #ifdef NDP / * For NDP C * / text_mode(); #endif / * Print exit message i f we overflowed the voxel l i s t . * / i f (BCode == 2) printf ("List overflow; too many voxels! The output f i l e is incomplete.\n"); return; } / * end Main * / Appendix G. Stability Region Computation Software 165 G.3 Sampling: S M P L P R O C . C /* */ / * SMPLPROC.C — Sampling routines for stability volume computation. * / /* */ / * Design & coding: M.G. Henders * / / * Last modified: 8 June 1991 * / /* */ / * Module type: Support module * / / * Module contents: * / /* */ / * Project f i l e : STABVOL.PRJ * / / * Support modules: None * / /* */ #include <compiler.h> #include <stdio.h> #include <assert.h> #include <math.h> / * Special data declarations * / #include <simul.h> #include <volume.h> / * Conditional declarations * / #ifdef Turbo / * For Turbo C * / #ifndef __HUGE__ #error Requires HUGE memory model. #endif #include <alloc.h> #include <graphics.h> #include <conio.h> #include <dos.h> #include <keystrks.h> #endif #ifdef NDP / * For NDP C * / #include <grex.h> #include <stdlib.h> #include <keycodes.h> #endif extern double GridSize[4]; long MaxVoxListSize; FILE *OutputFile; Appendix G. Stability Region Computation Software 166 int StablePoint (double *State, int Size) { /* */ / * Run an integration, beginning at State, and return 1 i f the * / / * resulting trajectory goes to the state-space origin; otherwise, * / / * return 0. Note that there is really no way to know that any * / / * given trajectory doesn't EVENTUALLY end up at the origin, after a * / / * great deal of winding about. This function uses some heuristic * / / * condition checks that have been found to give good results for * / / * the dynamics of the inverted pendulum. IN NO WAY CAN THESE * / / * CHECKS BE CONSIDERED TO BE A GENERALLY APPLICABLE TEST FOR * / / * DYNAMIC SYSTEM STABILITY! * / / * * / #define Maxlnteglter 10000L static double *S = NULL, *SDot = NULL, *Scale = NULL; double T = 1.0e-9, Eps = 1.0e-3, DidStep; double MaxAng = 6.981, MaxAngRate = 26.180, StableRad =0.1; double Wtl = 0.020, Wt2 = 0.094, Wt3 = 1.146, Wt4 = 0.344, StateRad; int i , StabIndex; long ItCount = 0; / * Function prototypes * / void Derivs (double *0utRate, double *Instate); void RKQC (double *0ut, int Size, double hTry, double Eps, double *State, double *Rate, double *Scale, double *hDid, double *hNext); int StabilityHeuristics (double *S); / * Allocate local storage, i f not already done, . . . * / i f (S == NULL) { assert ( (S = (double*) malloc (Size * sizeof (double) ) ) != NULL); } / * endif * / i f (SDot == NULL) { assert ((SDot = (double*) malloc (Size * sizeof (double))) != NULL); } / * endif * / i f (Scale == NULL) { assert ((Scale = (double*) malloc (Size * sizeof(double))) != NULL); } / * endif * / / * . . .and transfer data into i t , . . . * / for (i = 0; i < Size; i++) S[i] = State[i]; / * ...converting Theta and Theta-dot to radians. * / S[3] *= pi / 180.0; S[4] *= pi / 180.0; / * Run the system integration, unti l stability can be determined. * / Appendix G. Stability Region Computation Software 167 do { / * Init ial ize derivatives, . . . * / Derivs (SDot, S); / * . . .and tolerances, . . . * / for ( i = 1; i < Size; i++) { Scale[i] = fabs (S[i]) + fabs (T * SDot[i]) + 1.0e-30; } / * endfor * / / * . . .then perform the actual integration, . . . * / RKQC (S, Size, T, Eps, S, SDot, Scale, ftDidStep, &T); / * ...keeping track of how many passes we've done. * / ItCount++; / * Check the stability condition, for possible exit. * / / * F irs t , check for decisive 'out-of-bounds' condition, . . . * / i f ( (fabs (S[3]) > MaxAng) I I (fabs (S[4]) > MaxAngRate) I I (ItCount > Maxlnteglter) ) Stablndex = -1; else { / * If not definitely unstable, compute a weighted state radius, . . . * / StateRad = sqrt (sqr (Wtl * S[l]) + sqr (Wt2 * S[2]) + sqr (Wt3 * S[3]) + sqr (Wt4 * S[4]) ); / * . . .and set 'stable' or 'indeterminate', as appropriate. * / i f (StateRad < StableRad) Stablndex = 1; else Stablndex = 0; } / * endelse * / } while (!Stablndex); i f (Stablndex > 0) return (1); else return (0); } / * end StablePoint * / #ifdef Turbo / * For Turbo C * / unsigned StableVoxel (Voxel huge*VPtr, int NStates, double *Seed) #endif #ifdef NDP / * For NDP C * / unsigned StableVoxel (Voxel *VPtr, int NStates, double *Seed) #endif { /* */ / * Evaluate the stability of sixteen vertices of a 4-space voxel. * / / * Return 0 i f none of the vertices are stable, or, i f any vertex is * / / * stable, return a non-zero value. Bits 0 -> 15 of the return * / / * value are used to indicate the status of vertices 0 -> 15 respec- * / / * tively; each bit is set to 1 i f the corresponding vertex point is * / / * stable. The vertices are numbered as follows: * / Appendix G. Stability Region Computation Software 168 /* */ /* Vertex 0 -> (Mintf, MinX, MinY, MinZ) (0, 0, o, 0) */ /* Vertex 1 -> (MaxW, MinX, MinY, MinZ) (1, o, o, 0) */ /* Vertex 2 -> (MinW, MaxX, MinY, MinZ) (0, 1, o, 0) */ /* Vertex 3 -> (MaxW, MaxX, MinY, MinZ) (1, 1, o, 0) */ /* Vertex 4 -> (MinW, MinX, MaxY, MinZ) (0, o, 1, 0) */ /* Vertex 5 -> (MaxW, MinX, MaxY, MinZ) (1, o, 1, 0) */ /* Vertex 6 -> (MinW, MaxX, MaxY, MinZ) (0, 1, 1, 0) */ /* Vertex 7 -> (MaxW, MaxX, MaxY, MinZ) (1, 1, 1, 0) */ /* Vertex 8 -> (MinW, MinX, MinY, MaxZ) (0, o, o, 1) */ /* Vertex 9 -> (MaxW, MinX, MinY, MaxZ) (1, o, o, 1) */ /* Vertex 10 -> (MinW, MaxX, MinY, MaxZ) (o, 1, o, 1) */ /* Vertex 11 -> (MaxW, MaxX, MinY, MaxZ) (1, 1, o, 1) */ /* Vertex 12 -> (MinW, MinX, MaxY, MaxZ) (0, o, 1, 1) */ /* Vertex 13 -> (MaxW, MinX, MaxY, MaxZ) (1, o, 1, 1) */ /* Vertex 14 -> (MinW, MaxX, MaxY, MaxZ) (0, 1, 1, 1) */ /* Vertex 15 -> (MaxW, MaxX, MaxY, MaxZ) (1, 1, 1, 1) */ /* */ /*— — * / static double *TState = NULL, *UState = NULL; int i , j , k, 1; unsigned FlagBit = 0x0001, RetCode; / * Allocate local storage, i f not already done, . . . * / i f (TState •« NULL) { assert ( (TState = (double*) malloc (NStates * sizeof (double) ) ) != NULL); } / * endif * / i f (UState == NULL) { assert ( (UState = (double*) malloc (NStates * sizeof (double) ) ) != NULL); } / * endif * / / * ...compute the size and reference coordinates for the test voxel, . . . * / for ( i = 0; i < 4; i++) { TState[i+1] = Seed[i+1] + VPtr->IK.Index[i] * GridSize[i]; } / * endfor * / / * . . .and then test the voxel vertices, as offsets from this point. * / RetCode = 0; UState[0] =0.0; for (1 = 0; 1 < 2; 1++) { for (k = 0; k < 2; k++) { for (j = 0; j < 2; j++) { for (i = 0; i < 2; i++) Appendix G. Stability Region Computation Software 169 { / * Set up the test state, . . . * / UState[l] = TState[l] + i * GridSize[0] ; UState[2] = TState[2] + j * GridSize[l]; UState[3] = TState[3] + k * GridSize[2]; UState[4] = TState[4] + 1 * GridSize[3]; / * . . .and do stability test, setting the output code, . . . * / i f (StablePoint(UState, NStates) ) RetCode = RetCode I FlagBit; / * .. .then update the code selector. * / FlagBit = FlagBit « 1; } / * endfor * / } / * endfor * / } / * endfor * / } / * endfor * / / * Pass back the stability code. * / return (RetCode); } / * end StableVoxel * / #ifdef Turbo / * For Turbo C * / void DrawVoxelOutline (Voxel huge*Vox, double *Seed, Axis *Absci, Axis *0rdin, int Erase) #endif #ifdef NDP / * For NDP C * / void DrawVoxelOutline (Voxel *Vox, double *Seed, Axis *Absci, Axis *0rdin, int Erase) #endif { / * Draw 4 lines representing the outline of the voxel passed in Vox. * / static int Color = 1; double DimA, DimO, xO, yO, x l , y l ; / * Cycle the color, so we can see something happening. * / #ifdef Turbo / * For Turbo C * / i f (Erase) setcolor (0); else setcolor (Color); #endif #ifdef NDP / * For NDP C * / i f (Erase) set_color (0); else set_color (Color); #endif Color = (Color + 1) '/, 16; i f (! Color) Color++; / * . . .and load local variables. * / DimA = GridSize[Absci->DPtr] ; Appendix G. Stability Region Computation Software DimO = GridSize[Ordin->DPtr] ; /* Compute the coordinates of the reference p o i n t , . . . */ xO = Seed[l+Absci->DPtr] + Vox->IK.Index[Absci->DPtr] * DimA; yO = Seed[l+Ordin->DPtr] + Vox->IK.Index[Ordin->DPtr] * DimO; /* . . . a n d those of the farthest -removed corner p o i n t , . . . */ x l = xO + DimA; y l = yO + DimO; /* . . . c o n v e r t to screen coord inates , . . . */ xO = A b s c i - >CenC + A b s c i -•>Dir * Absci->Mag * (xO - Absc i - ->Cen) ; yO = Ordin- >CenC + Ordin-•>Dir * Ordin->Mag * (yO - Ordin- ->Cen) ; x l = A b s c i - >CenC + Absc i -•>Dir * Absci->Mag * ( x l - Absc i -•>Cen) ; = Ordin- >CenC + Ordin-•>Dir * Ordin->Mag * ( y l - Ordin- ->Cen) ; /* . . . a n d then draw the voxel edges. */ i f ( (xO >= 0) && (xO <= Absci->MaxC) && (x l >= 0) && ( x l <= Absci->MaxC) ) { i f ( (yO >= 0) && (yO <= Ordin->MaxC) && ( y l >= 0) && ( y l <= Ordin->MaxC) ) { #ifdef Turbo /* For Turbo C */ moveto ( ( int )xO, ( in t )yO) ; l i n e t o ( ( int )xO, ( i n t ) y l ) ; l i n e t o ( ( i n t ) x l , ( i n t ) y l ) ; l i n e t o ( ( i n t ) x l , ( in t )yO) ; l i n e t o ( ( int )xO, ( in t )yO) ; #endif #ifdef NDP /* For NDP C */ move ( ( int )xO, ( in t )yO) ; draw ( ( int )xO, ( i n t ) y l ) ; draw ( ( i n t ) x l , ( i n t ) y l ) ; draw ( ( i n t ) x l , ( in t )yO) ; draw ( ( int )xO, ( in t )yO) ; #endif } /* endif */ } /* endif */ r e t u r n ; } /* end DrawVoxelOutline */ #ifdef Turbo /* For Turbo C */ vo id SaveVolume (FILE *0Fi le , Voxel huge*VoxList, long LastVox, Appendix G. Stability Region Computation Software 171 double *Seed) #endif #ifdef NDP /* For NDP C */ v o i d SaveVolume (FILE * 0 F i l e , Voxel *VoxList , long LastVox, double *Seed) #endif { /* Write the accumulated voxel informat ion out to a f i l e . */ long i , NumVox = 0 ; /* Count the a c t u a l number of a c t i v e voxe ls . */ f o r ( i = 0 ; i <= LastVox; i++) { i f (VoxList [ i ] .SCode) NumVox++; } /* endfor */ rewind ( O F i l e ) ; f p r i n t f (OF i le , "Act i ve Voxels : */.ld of y.ld\n\n", NumVox, (LastVox + 1) ) ; f p r i n t f (OF i le , "Seed S ta te : '/.f '/.f '/.f '/.f y.f\n\n", Seed[0] , Seed[l] , Seed[2], Seed[3] , Seed[4]) ; f p r i n t f (OF i le , "Base l ine Voxel Dimensions: */.f '/.f '/,f 7,f \n\n" , G r i d S i z e [ 0 ] , G r i d S i z e [ l ] , G r i d S i z e [ 2 ] , G r i d S i z e [ 3 ] ) ; f p r i n t f (OF i le , " X XD Th ThD SCode\n\n"); f o r ( i = 0 ; i <= LastVox; i++) { f p r i n t f (OF i le , "y.3hd y.3hd '/.3hd */.3hd y,04x\n", VoxList [ i ] . IK. Index[0] , V o x L i s t [ i ] . I K . I n d e x [ l ] , V o x L i s t [ i ] . I K . I n d e x [ 2 ] , V o x L i s t [ i ] . I K . I n d e x [ 3 ] , VoxL is t [ i ] .SCode) ; } /* endfor */ f f l u s h ( O F i l e ) ; r e t u r n ; } /* end SaveVolume */ unsigned unsigned unsigned unsigned FaceWN = FaceXN = FaceYN = FaceZN = 0x5555, 0x3333, OxOFOF, OxOOFF, FaceWP FaceXP FaceYP FaceZP OxAAAA OxCCCC OxFOFO OxFFOO #ifdef Turbo /* For Turbo long Wi l lOver lap (Voxel Voxel #endif C */ huge*VoxList, long LastVox, huge*Proposed) Appendix G. Stability Region Computation Software 172 #ifdef NDP / * For NDP C * / long WillOverlap (Voxel *VoxList, long LastVox, Voxel *Proposed) #endif { /* */ / * If the Proposed voxel is already l isted, return its l i s t index; * / / * i f the voxel is new, return -1. * / / * * / long i ; / * Scan the voxel l i s t ; break out i f we find an overlap match. * / for ( i = 0; i <= LastVox; i++) { i f (Proposed->IK.Key == VoxList[i].IK.Key) return ( i); } / * endfor • / return (-1); } / * end WillOverlap * / #ifdef Turbo / * For Turbo C * / void AddNewVoxel (Voxel huge*VoxList, long CurrVox, long *LastVox, int Dir, int SysSize, double *Seed, Axis *Absci, Axis *0rdin) #endif #ifdef NDP / * For NDP C * / void AddNewVoxel (Voxel *VoxList, long CurrVox, long *LastVox, int Dir, int SysSize, double *Seed, Axis *Absci, Axis *0rdin) #endif { int PropDir; char Offset; / * Add a new voxel, i f there is room for i t . * / (•LastVox)++; i f (•LastVox < MaxVoxListSize) { PropDir = abs (Dir) - 1; Offset = (char) (2 • (Dir / abs (Dir) ) ); / • Set up proposed voxel indices, . . . • / VoxList[•LastVox].IK.Key = VoxList[CurrVox].IK.Key; VoxList[•LastVox].IK.Index[PropDir] += Offset; / • . . .and check the new voxel for overlap with existing voxels. * / i f (WillOverlap (VoxList, ( (•LastVox) - 1), ftVoxList[•LastVox]) < 0) { / • Determine the stability status of the new voxel, . . . • / VoxList[•LastVox].SCode = StableVoxel (ftVoxList[•LastVox], SysSize, Seed); Appendix G. Stability Region Computation Software 173 /* ...and draw i t , i f i t i s not completely unstable. */ i f (VoxList[*LastVox].SCode) { DrawVoxelOutline (ftVoxList[*LastVox], Seed, Absci, Ordin, 0); } /* endif */ } /* endif */ else (*LastVox)—; /* Reset the index, to re-use this l i s t entry. */ } /* endif */ return; } /* end AddNewVoxel */ #ifdef Turbo /* For Turbo C */ void BuildVolume (Voxel huge*VoxList, long CurrVox, long *LastVox, int SysSize, double *Seed, Axis *Absci, Axis *0rdin) #endif #ifdef NDP /* For NDP C */ void BuildVolume (Voxel *VoxList, long CurrVox, long *LastVox, int SysSize, double *Seed, Axis *Absci, Axis *0rdin) #endif { unsigned CurrSCode; /* Load the s t a b i l i t y state of the current voxel, ... */ CurrSCode = VoxList[CurrVox].SCode; /* ...and propagate; start with the negative-going W face, ... */ i f (CurrSCode ft FaceWN) { AddNewVoxel (VoxList, CurrVox, LastVox, -1, SysSize, Seed, Absci, Ordin); } /* endif */ /* ...then, the positive-going W face, ... */ i f (CurrSCode ft FaceWP) { AddNewVoxel (VoxList, CurrVox, LastVox, 1, SysSize, Seed, Absci, Ordin); } /* endif */ /* ...next, the negative-going X face, ... */ i f (CurrSCode ft FaceXN) { AddNewVoxel (VoxList, CurrVox, LastVox, -2, SysSize, Seed, Absci, Ordin); } /* endif */ /* ...the positive-going X face, ... */ i f (CurrSCode ft FaceXP) { AddNewVoxel (VoxList, CurrVox, LastVox, 2, SysSize, Seed, Absci, Ordin); Appendix G. Stability Region Computation Software } / * endif * / / * . . .the negative-going Y face, . . . * / i f (CurrSCode ft FaceYN) { AddNewVoxel (VoxList, CurrVox, LastVox, -3, SysSize, Seed, Absci, Ordin); } / * endif * / / * . . .the positive-going Y face, . . . * / i f (CurrSCode ft FaceYP) { AddNewVoxel (VoxList, CurrVox, LastVox, 3, SysSize, Seed, Absci, Ordin); } / * endif * / / * . . .then the negative-going Z face, . . . * / i f (CurrSCode ft FaceZN) { AddNewVoxel (VoxList, CurrVox, LastVox, -4, SysSize, Seed, Absci, Ordin); } / * endif * / / * . . .and, f inal ly , the positive-going Z face. * / i f (CurrSCode ft FaceZP) { AddNewVoxel (VoxList, CurrVox, LastVox, 4, SysSize, Seed, Absci, Ordin); } / * endif * / return; } / * end BuildVolume * / Appendix H Data Manipulation Software This Appendix comprises listings of software used for manipulating stability region data files, as produced by the region computation program listed in Appendix G. The first program, SLIC30F4.C, (Section H.l) extracts three-dimensional data sections from the four-dimensional data set, while SLIC20F3.C, (Section H.2) the second program, op-erates on three-dimensional data sections, extracting two-dimensional slices from them. This program also plots a display of the data section as it is extracted. SLIC20F3 generates an output file which may be read directly by the GRAFTOOL software package [GRAF] (see Appendix E); the output of SLIC30F4 is not directly usable for GRAFTOOL input, but may be converted to a suitable format using SCAT-TERS. C, which is listed in Section H.3. These programs use FILEPROC.C (Section F.7), but are otherwise self-contained; they are intended for compilation using Turbo C. [BOR] 175 Appendix H. Data Manipulation Software 176 H . l 3-D Sectioning: SLIC30F4.C /* + / /* SLIC30F4.C — Program f o r e x t r a c t i n g 3-dimensional sect ions from a */ /* 4-dimensional data f i l e , as generated by STABVOL.EXE. */ /* */ /* Design ft coding: M.G. Henders */ /* Last modi f ied : 21 June 1991 */ /* */ /* Module type : Ex terna l l y -suppor ted MAIN */ /* Module name: SLIC30F4.C */ /* Module contents : */ /* */ /* Support modules: FILEPROC.C */ #define Main main /* */ #include <stdio.h> #include <string.h> #include <alloc.h> #include <assert.h> #include <keystrks.h> #include <volume.h> vo id Main ( i n t argc , char *argv[]) { char *InputFilename, *OutputFilename, *TStr ing, SectIndex; double *Seed, SectCoord = 0 . 0 , G r i d S i z e [ 4 ] ; FILE * InputF i le , *0utputF i le ; Voxel TVox, MinVox, MaxVox; i n t i , Key = 0 , F l a g , S l i c e A x i s = Key_X, S l i c e A x i s P t r ; long Count; unsigned Mask, P i c k e r , OutCode; unsigned FaceWN = 0x5555, FaceWP = OxAAAA; unsigned FaceXN = 0x3333, FaceXP = OxCCCC; unsigned FaceYN = OxOFOF, FaceYP = OxFOFO; unsigned FaceZN = OxOOFF, FaceZP = OxFFOO; FILE *fopenWarn (char *FileName, char *Mode, char *OutWarn); /* A l l o c a t e s ta te vector and s t r i n g space. */ assert ((Seed = (double*) malloc (5 * s izeof (double) ) ) != NULL); asser t ((InputFilename = (char*) malloc (255 * s i zeo f ( char ) ) ) != NULL); asser t ((OutputFilename = (char*)malloc (255 * s i zeof (char ) ) ) != NULL); assert ( (TStr ing = (char*) malloc (255 * s i zeof (char ) ) ) != NULL); Appendix H. Data Manipulation Software 177 /* Look for conraiand-line arguments; prompt i f not found/no good. */ i f (argc != 3) "This i s a program to extract 3-space data from 4-space f i l e s generated by"); "the program STABVOL. In general, the input and output filenames may be"); "specified on the command l i n e , as SLIC30F4 i n p u t f i l e outputfile. Please"); "enter a name for the input f i l e , now.\n"); { puts puts puts puts gets puts puts gets (OutputFilename); puts (""); } /* endif */ else { strcpy (InputFilename, argv [1]); strcpy (OutputFilename, argv[2]); } /* endelse */ InputFilename); "\nNow enter a name for the output f i l e . If the selected f i l e i s already"); "present, you w i l l be prompted to confirm overwrite.\n"); /* Try to open the input f i l e , bitching to the user i f not found. */ InputFile = fopenWarn (InputFilename, " r t " , "WarnNotFound"); /* Get extremal values from the f i l e , for reference. */ puts ("Scanning input f i l e to determine data ranges. Please wait..."); /* Read f i l e header information, ... */ for ( i = 0; i < 7 for ( i = 0; i < 5 for ( i = 0; i < 3 for ( i = 0; i < 4 for ( i = 0; i < 5 i++) fscanf (InputFile, '"/.s", TString); i++) fscanf (InputFile, '"/.If", ft(Seed[i]) ); i++) fscanf (InputFile, '"/.s", TString); i++) fscanf (InputFile, '"/.If", ft(GridSize[i]) ); i++) fscanf (InputFile, '"/.s", TString); /* ...and then go for the actual data. */ Count = 0; while (fscanf (InputFile, '"/.hd '/.hd V.hd V.hd */.x", ft(TVox.IK.Index[0]), ft(TVox.IK.Index[1]), ft(TVox.IK.Index[2]), ft(TVox.IK.Index[3]), ftTVox.SCode) == 5) { i f (!Count) { MinVox.IK.Key = TVox.IK.Key; MaxVox.IK.Key = TVox.IK.Key; } /* endif */ for ( i = 0; i < 4; i++) Appendix H. Data. Manipulation Software { i f (TVox.IK.Index[i] < MinVox.IK.Index[i]) { MinVox.IK.Index[i] = TVox.IK.Index[i]; } / * endif * / else i f (TVox.IK.Index[i] > MaxVox.IK.Index[i]) { MaxVox.IK.Index[i] = TVox.IK.Index[i]; } / * endelse * / } / * endfor * / Count++; i f ( (Count '/. 1024) == 0) printf ("I"); else i f ( (Count '/. 1024) == 512) printf ("+"); } / * endwhile * / puts (""); / * Outer level program loop, to do multiple output f i l es . * / do { rewind (InputFile); i f (Key) { puts ("\nEnter new output filename."); gets (OutputFilename); puts (""); } / * endif * / / * Try to open the output f i l e , warning i f i t already exists. * / OutputFile = fopenWarn (OutputFilename, "wt", "WarnOverWrite"); / * Get sl icing axis. * / puts ("Change slicing axis? (Y/N)"); do { Flag = (Key = GetScanCode (0) ) == Yes; } while (IFlag && (Key != No) ); i f (Flag) { puts ("Coordinate mapping: X -> W, XDot -> X, Theta -> Y, ThetaDot -> Z"); switch (SliceAxis) { case Key_W: puts ("Now slicing at constant W value. Select new slice normal. (W/X/Y/Z)\n"); break; case Key_X: puts ("Now slicing at constant X value. Select new slice normal. (W/X/Y/Z)\n"); Appendix H. Data Manipulation Software break; case Key_Y: puts ("Now slicing at constant Y value. Select new slice normal. (W/X/Y/Z)\n"); break; case Key_Z: puts ("Now slicing at constant Z value. Select new slice normal. (W/X/Y/Z)\n"); break; default: ; / * Nothing * / } / * endswitch * / do { Flag = (SliceAxis = GetScanCode (0) ) == Key_W; } while (!Flag kk (SliceAxis != Key_X) kk (SliceAxis != Key.Y) (SliceAxis != Key.Z) ); } / * endif * / / * Set up axis selector. * / switch (SliceAxis) { case Key_W: SliceAxisPtr = 0 case Key_X: SliceAxisPtr = 1 case Key_Y: SliceAxisPtr = 2 case Key_Z: SliceAxisPtr = 3 default: ; / * Nothing * / }• / * endswitch * / break; break; break; break; / * Get sl icing coordinate. * / puts ("Change slicing coordinate? (Y/N)"); do { Flag = (Key = GetScanCode (0) ) == Yes; } while (!Flag kk (Key != No) ); i f (Flag) { printf ("Current slicing axis extremals: near <'/,f, y,f>\n", Seed[SliceAxisPtr+l] + MinVox.IK.Index[SliceAxisPtr] * GridSize[SliceAxisPtr], Seed[SliceAxisPtr+l] + MaxVox.IK.Index[SliceAxisPtr] * GridSize[SliceAxisPtr]); Appendix H. Data. Manipulation Software 180 p r i n t f ("Enter coordinate value f o r s l i c i n g p lane ; now = */.f An", SectCoord); do { i = scanf 07.1f", ftSectCoord); i f ( i != 1) puts ( " Inva l id value entered; t r y a g a i n . " ) ; } whi le ( i != 1 ) ; puts (""); } /* endif */ /* Transfer f i l e header information to output f i l e , . . . */ f o r ( i = 0 ; i < 7; i++) fscanf ( I n p u t F i l e , '"/.s", T S t r i n g ) ; f o r ( i = 0 ; i < 5 ; i++) fscanf ( InputF i l e , '"/.If", &(Seed[i ] ) ) ; f o r ( i = 0 ; i < 3 ; i++) fscanf ( I n p u t F i l e , "7,s", T S t r i n g ) ; f o r ( i = 0 ; i < 4; i++) fscanf ( I n p u t F i l e , '"/.If", & (Gr idS i ze [ i ] ) ) ; f o r ( i = 0 ; i < 5; i++) fscanf ( I n p u t F i l e , "'/.s", T S t r i n g ) ; switch (S l i ceAx i s ) { case Key_W: f p r i n t f (OutputF i le , "Sect ion on X at */.f \n\n" , SectCoord); break; case Key_X: f p r i n t f (OutputF i le , "Sect ion on XDot at */f\n\n", SectCoord); break; case Key_Y: f p r i n t f (OutputF i le , "Sect ion on Theta at y.f\n\n", SectCoord); break; case Key_Z: f p r i n t f (OutputF i le , "Sect ion on ThetaDot at '/.f \n\n" , SectCoord); break; d e f a u l t : ; /* Nothing */ } /* endswitch */ f p r i n t f (OutputF i le , "Seed S ta te : */.f */.f */.f */,f */.f \n\n" , Seed[0] , S e e d [ l ] , Seed[2], Seed[3] , Seed[4]) ; f p r i n t f (OutputF i le , "Basel ine Voxel Dimensions: */,f '/.f */.f '/.f\n\n", G r i d S i z e [ 0 ] , G r i d S i z e [ l ] , G r i d S i z e [ 2 ] , G r i d S i z e [ 3 ] ) ; f p r i n t f (OutputF i le , " X XD Th ThD SCode\n\n"); /* . . . a n d then set up f o r data s l i c i n g . */ Appendix H. Data Manipulation Software i f ( (SectCoord - Seed[SliceAxisPtr+l]) >= 0) { Sectlndex = (char) ( (SectCoord - Seed[SliceAxisPtr+l] (GridSize[SliceAxisPtr] / 2) ) / GridSize[SliceAxisPtr]); switch (SliceAxis) { case Key_H: i f (Sectlndex '/. 2) { Mask = FaceWP; Sectlndex—; } / * endif * / else Mask = FaceWN; break; case Key_X: i f (Sectlndex '/. 2) { Mask = FaceXP; Sectlndex—; } / * endif * / else Mask = FaceXN; break; case Key_Y: i f (Sectlndex '/, 2) { Mask = FaceYP; Sectlndex—; } / * endif * / else Mask = FaceYN; break; case Key_Z: i f (Sectlndex */. 2) { Mask = FaceZP; Sectlndex—; } / * endif * / else Mask = FaceZN; break; default: ; / * Nothing * / } / * endswitch * / } / * endif * / else { Sectlndex = (char) ( (SectCoord - Seed[SliceAxisPtr+l] (GridSize [SliceAxisPtr] / 2) ) / GridSize[SliceAxisPtr]); switch (SliceAxis) { case Key_W: i f ( (- Sectlndex) '/. 2) { Mask = FaceWP; Sectlndex—; } / * endif * / Appendix H. Data Manipulation Software 182 else Mask = FaceWN; break; case Key_X: i f ( (- Sectlndex) '/. 2) { Mask = FaceXP; Sectlndex—; } / * endif * / else Mask = FaceXN; break; case Key_Y: i f ( (- Sectlndex) '/. 2) { Mask = FaceYP; Sectlndex--; } / * endif * / else Mask = FaceYN; break; case Key_Z: i f ( (- Sectlndex) 7. 2) { Mask = FaceZP; Sectlndex—; } / * endif * / else Mask = FaceZN; break; default: ; / * Nothing * / } / * endswitch * / } / * endelse * / puts ("Processing data f i l e . Please wait.. ."); Count = 0; while (f scanf (InputFile, "V.hd */.hd */.hd 7.hd 7.x", ft(TVox.IK.Index[0]), ft(TVox.IK.Index[1]), ft(TVox.IK.Index[2]), ft(TVox.IK.Index[3]), ftTVox.SCode) == 5) { / * F irs t , test for a basic index value match, . . . * / i f (TVox.IK.Index[SliceAxisPtr] == Sectlndex) { / * . . .and, i f found, pull out the relevant hyper-face, . . . * / for (i = 0, Picker = 0x0001, OutCode = 0; i < 16; i++, Picker « = 1) { i f (Mask ft Picker) OutCode |= (Picker ft TVox.SCode); else OutCode <<= 1; } / * endfor * / OutCode » = 8; / * .. .then write the resulting 3-voxel to the output f i l e . * / i f (OutCode) Appendix H. Data Manipulation Software 183 { fprintf (OutputFile, '7.3hd TVox.IK.Index[0], TVox.IK.Index[3], } / * endif * / } / * endif * / '/.3hd */.3hd y.3hd y.04x\n" , TVox.IK.Index[1], TVox.IK.Index[2], OutCode); Count++; i f ( (Count */. 1024) == 0) printf ("I"); else i f ( (Count '/. 1024) == 512) printf ("+"); } / * endwh.ile * / puts ("\nCreate another data section? Hit <Esc> to exit, any other key to continue."); Key = GetScanCode (0); fflush (stdin); / * Button up the output f i l e before looping or leaving, . . . * / fclose (OutputFile); } while (Key != Esc); / * . . .and close the input f i l e before f inal exit. * / fclose (InputFile); } / * end Main * / Appendix H. Data, Manipulation Software 184 H.2 2-D Sectioning: SLIC20F3.C /* */ / * SLIC20F3.C — Program for extracting 2-dimensional sections from a * / / * 3-dimensional data f i l e , as extracted from a 4-space * / / * f i l e by SLIC30F4.EXE. 4-space f i les are created by STABVOL.EXE. * / /* */ / * Design ft coding: M.G. Henders * / / * Last modified: 24 June 1991 * / /* */ / * Module type: Externally-supported MAIN * / / * Module name: SLIC20F3.C * / / * Module contents: * / /* */ / * Support modules: FILEPROC.C * / #define Main main /*_ */ #include <stdio.h> #include <graphics.h> #include <string.h> #include <conio.h> #include <bios.h> #include <stdlib.h> #include <alloc.h> #include <assert.h> #include <math.h> #include <keystrks.h> #include <volume.h> #define AbsciMinDef -375.0 #define AbsciMaxDef 375.0 #define OrdinMinDef -900.0 #define OrdinMaxDef 900.0 extern unsigned _stklen = 0x4000; void Main (int argc, char *argv[]) { char ScreenBuf [80*25*2], *GraphDriverPath = " C : \ \ T C \ \ " , Sectlndex; char *InputFilename, *OutputFilename, *TString; long Count; Voxel TVox, MinVox, MaxVox; FILE *InputFile, *0utputFile; int CursorX, CursorY, GraphMode, GraphDriver = DETECT; int MaxX, MaxY, CenX, CenY, Key = 0, Flag, i ; Appendix H. Data Manipulation Software 185 int Slice4Axis, Slice3Axis = Key_X, Slice3AxisPtr; double xO, yO, GridSize[4], *Seed, SectCoord = 0.0; Axis Absci, Ordin; unsigned Mask, Picker, OutCode; unsigned FaceON = 0x0055, FaceOP = OxOOAA; unsigned FacelN = 0x0033, FacelP = OxOOCC; unsigned Face2N = OxOOOF, Face2P = OxOOFO; / * Function prototypes * / FILE *fopenWarn (char *FileName, char *Mode, char *OutWarn); void OutputStablePoints (Voxel *Vox, double *Seed, double *GridSize, Axis *Absci, Axis *0rdin, FILE *0File); / * Allocate state vector and string space. * / assert ((Seed = (double*) malloc (5 * sizeof(double))) != NULL); assert ((InputFilename = (char*) malloc (255 * sizeof(char))) != NULL); assert ((OutputFilename = (char*)malloc(255 * sizeof(char))) != NULL); assert ((TString = (char*) malloc (255 * sizeof(char))) != NULL); / * Save the text-mode screen state. * / CursorX = wherex(); CursorY = w h e r e y O ; assert (gettext (1, 1, 80, 25, ScreenBuf) ); / * Look for command-line arguments, and prompt i f not found/no good. * / i f (argc != 3) "This is a program to extract 2-space data from 3-space f i les generated by"); "the program SLIC30F4. In general, the input and output filenames may be"); "specified on the command l ine, as SLIC20F3 inputfile outputfile. Please"); "enter a name for the input f i l e , now.\n"); { puts puts puts puts gets puts puts gets (OutputFilename); puts (""); } / * endif * / else { strcpy (InputFilename, argv [1]); strcpy (OutputFilename, argv[2]); } / * endelse * / InputFilename); "\nNow enter a name for the output f i l e . If the selected f i l e is already"); "present, you wi l l be prompted to confirm overwrite.\n"); Appendix H. Data Manipulation Software 186 / * Init ial ize display parameters. * / Absci.Min = AbsciMinDef; Absci.Max = AbsciMaxDef; Ordin.Min = OrdinMinDef; Ordin.Max = OrdinMaxDef; / * Outer loop, for different input f i l es . * / do i i f (Key) { puts ("\nEnter new input filename."); gets (InputFilename); puts (""); } / * endif * / / * Try to open the input f i l e , bitching to the user i f not found. * / InputFile = fopenWarn (InputFilename, "rt", "WarnNotFound"); / * Get extremal values from the f i l e , for reference. * / puts ("Scanning input f i l e to get data ranges. Please wait.. ."); / * Read f i l e header information, . . . * / for (i = 0; i < 2; i++) fscanf (InputFile, '"/.s", TString); fscanf (InputFile, "'/,s", TString); i f (Istrcmp (TString, "X") ) Slice4Axis = 0; else i f (Istrcmp (TString, "XDot") ) Slice4Axis = 1; else i f (Istrcmp (TString, "Theta") ) Slice4Axis = 2; else Slice4Axis = 3; for ( i = 0; i < 4; i++) fscanf (InputFile, "'/.s", TString); for ( i = 0; i < 5; i++) fscanf (InputFile, '"/.If", ft (Seed [i]) ); for ( i = 0; i < 3; i++) fscanf (InputFile, '"/.s", TString); for (i = 0; i < 4; i++) fscanf (InputFile, "'/.lf", ft (GridSize [i]) ); for (i = 0; i < 5; i++) fscanf (InputFile, '"/.s", TString); / * . . .and then go for the actual data. * / Count = 0; while (fscanf (InputFile, '7.hd */.hd */.hd */.hd '/.x", ft(TVox.IK.Index[0]), ft(TVox.IK.Index[1]), ft(TVox.IK.Index[2]), ft(TVox.IK.Index[3]), ftTVox.SCode) == 5) { i f (ICount) { MinVox.IK.Key = TVox.IK.Key; MaxVox.IK.Key = TVox.IK.Key; } / * endif * / for (i = 0; i < 4; i++) { i f (TVox.IK.Index[i] < MinVox.IK.Index[i]) { MinVox.IK.Index[i] = TVox.IK.Index[i]; } / * endif * / else i f (TVox.IK.Index[i] > MaxVox.IK.Index[i]) { MaxVox.IK.Index[i] = TVox.IK.Index[i]; } / * endelse * / Appendix H. Data Manipulation Software } /* endfor */ Count++; i f ( (Count V. 1024) == 0) p r i n t f ("I"); e lse i f ( (Count */. 1024) == 512) p r i n t f ("+"); } /* endwhile */ puts (""); /* Outer l e v e l program loop , to do m u l t i p l e output f i l e s . */ do { rewind ( I n p u t F i l e ) ; i f (Key) { puts ("\nEnter new output f i l e n a m e . " ) ; gets (OutputFilename); puts (""); } /* endif */ /* Try to open the output f i l e , warning i f i t a l ready e x i s t s OutputFi le = fopenWarn (OutputFilename, "wt" , "WarnOverWrite /* Loop to set up output p l o t the way you want i t . */ do { rewind ( I n p u t F i l e ) ; rewind (OutputF i le ) ; /* Get s l i c i n g a x i s . */ puts ("Change s l i c i n g ax is? (Y/N)"); do { F lag = (Key = GetScanCode (0) ) == Yes; } whi le ( IF lag kk (Key != No) ) ; i f (Flag) { switch (S l i ce4Ax is ) { case 0 : puts ("Coordinate mapping: XDot -> X, Theta -> Y, ThetaDot -> Z" ) ; break; case 1: puts ("Coordinate mapping: X -> X, Theta -> Y, ThetaDot -> Z" ) ; break; case 2: puts ("Coordinate mapping: X -> X, XDot -> Y, ThetaDot -> Z" ) ; break; Appendix H. Data Manipulation Software 188 case 3: puts ("Coordinate mapping: X -> X, XDot -> Y, Theta -> Z"); break; default: ; / * Nothing * / } / * endswitch * / switch (Slice3Axis) { case Key_X: puts ("Now slicing at constant X value. Select new slice normal. (X/Y/Z)\n"); break; case Key_Y: puts ("Now slicing at constant Y value. Select new slice normal. (X/Y/Z)\nM); break; case Key_Z: puts ("Now slicing at constant Z value. Select new slice normal. (X/Y/Z)\n"); break; default: ; / * Nothing * / } / * endswitch * / do { Flag = (Slice3Axis = GetScanCode (0) ) == Key_X; } while (!Flag kk (Slice3Axis != Key_Y) kk (Slice3Axis != Key_Z) ); } / * endif * / / * Set up axis selector. * / switch (Slice3Axis) { case Key_X: switch (Slice4Axis) case 0: Slice3AxisPtr = 1; break; case 1: Slice3AxisPtr = 0; break; case 2: Slice3AxisPtr = 0; break; case 3: Slice3AxisPtr = 0; break; default: ; / * Nothing * / } / * endswitch * / break; case Key_Y: Appendix H. Data Manipulation Software 189 switch (S l i ce4Ax is ) { case 0 : S l i c e 3 A x i s P t r = 2 case 1: S l i c e 3 A x i s P t r = 2 case 2: S l i c e 3 A x i s P t r = 1 case 3 : S l i c e 3 A x i s P t r = 1 d e f a u l t : ; /* Nothing */ } /* endswitch */ break; case Key_Z: switch (S l i ce4Ax is ) { case 0 : S l i c e 3 A x i s P t r = 3 case 1: S l i c e 3 A x i s P t r = 3 case 2: S l i c e 3 A x i s P t r = 3 case 3 : S l i c e 3 A x i s P t r = 2 d e f a u l t : ; /* Nothing */ } /* endswitch */ break; d e f a u l t : ; /* Nothing */ } /* endswitch */ /* Get s l i c i n g coordinate. */ puts ("Change s l i c i n g coordinate? (Y/N)"); do { F lag = (Key = GetScanCode (0) ) == Yes; } whi le ( IF lag ftft (Key != No) ) ; i f (Flag) { p r i n t f ("Current s l i c i n g ax is extremals: near <*/,f, '/f>\n", Seed[Sl ice3AxisPtr+l ] + MinVox. IK . Index[Sl ice3AxisPtr ] * Gr idS ize [ S l i c e 3 A x i s P t r ] , Seed[Sl ice3AxisPtr+l ] + MaxVox. IK. Index[Sl ice3AxisPtr] * G r i d S i z e [ S l i c e 3 A x i s P t r ] ) ; p r i n t f ("Enter coordinate value f o r s l i c i n g p lane ; now = '/,f An", SectCoord); do { i = scanf ('"/.If", ftSectCoord) ; i f ( i != 1) puts ( " Inva l id value entered; t r y a g a i n . " ) ; } whi le ( i != 1 ) ; puts (""); } /* endif */ /* Get viewing parameters. */ puts ("Modify d i sp lay sca l ing? (Y/N)"); break; break; break; break; break; break; break; break; Appendix H. Data Manipulation Software 190 do { F lag = (Key = GetScanCode (0) ) == Yes; } whi le ( IF lag &ft (Key != No) ) ; i f (Flag) { p r i n t f ("Enter min. and max. values f o r absc issa data ; now = 7,f, y . f .\n" , A b s c i . M i n , Absc i .Max) ; do { i = scanf ("'/.If '/.If", &(Absc i .Min) , ft(Absci.Max) ) ; i = ( i != 2) || (Absci .Min == Absci .Max) ; i f ( i ) { puts ( " Inva l id value entered; t r y a g a i n . " ) ; } /* endif */ } whi le ( i ) ; i f (Absc i .Min > Absci.Max) { xO = A b s c i . M i n ; Absc i .M in = Absci .Max; Absci.Max = xO; } /* endif */ p r i n t f ("Enter min. and max. values f o r ord inate data ; now = V.f, y . f .\n" , Ord in .Min , Ordin.Max); do { i = scanf ("'/.If '/.If", ft(Ordin.Min) , ft(Ordin.Max) ) ; i = ( i != 2) II (Ordin.Min == Ordin.Max); i f ( i ) { puts ( " Inva l id value entered; t r y a g a i n . " ) ; } /* endif */ } whi le ( i ) ; puts (""); i f (Ordin.Min > Ordin.Max) { xO = Ordin .Min ; Ordin.Min = Ordin.Max; Ordin.Max = xO; } /* endif */ } /* endif */ Absc i .Len = Absci.Max - A b s c i . M i n ; Ordin.Len = Ordin.Max - Ord in .Min ; Absci .Cen = 0 .5 * (Absci .Min + Absci .Max) ; Ordin.Cen = 0.5 * (Ordin.Min + Ordin.Max); swi tch (S l i ce3Ax is ) { case Key_X: A b s c i . D i r = 1 case Key_Y: A b s c i . D i r = - 1 case Key_Z: A b s c i . D i r = 1 d e f a u l t : ; /* Nothing */ } /* endswitch */ Ord in .D i r = - 1 ; break; Ord in .D i r = - 1 ; break; Ord in .D i r = 1; break; Appendix if. Data Manipulation Software 191 switch (Slice4Axis) { case 0: switch (Slice3Axis) case Key_X: Absci DPtr = 2; Ordin DPtr = 3; break; case Key.Y: Absci DPtr = 1; Ordin DPtr = 3; break; case Key.Z: Absci DPtr = 2; Ordin DPtr = 1; break; default: ; /* Nothing */ } /* endswitch */ break; case 1: switch (Slice3Axis) case Key.X: Absci .DPtr = 2; Ordin DPtr = 3; break; case Key.Y: Absci .DPtr = 0; Ordin DPtr = 3; break; case Key.Z: Absci DPtr = 2; Ordin DPtr = 0; break; default: ; /* Nothing */ } /* endswitch */ break; case 2: switch (Slice3Axis) case Key.X: Absci .DPtr = 1; Ordin DPtr = 3; break; case Key.Y: Absci .DPtr = 0; Ordin DPtr = 3; break; case Key.Z: Absci .DPtr = 1; Ordin DPtr = 0; break; default: ; /* Nothing */ } /* endswitch */ break; case 3: switch (Slice3Axis) case Key.X case Key.Y case Key.Z default: Absci.DPtr = 1; Absci.DPtr = 0; Absci.DPtr = 1; /* Nothing */ } /* endswitch */ break; default: ; /* Nothing */ } /* endswitch */ Ordin.DPtr Ordin.DPtr Ordin.DPtr = 0 break; break; break; /* Initialize the graphics system. */ if (GraphDriver == DETECT) { initgraph (ftGraphDriver, ftGraphMode, GraphDriverPath); assert (GraphDriver >= 0); HaxX = getmaxx(); MaxY = getmaxyO ; } /* endif */ else Appendix H. Data, Manipulation Software 192 { setgraphmode (GraphMode); } /* endelse */ /* . . .d raw screen boundaries, . . . */ se tco lo r (2) ; l i n e (0, 0 , MaxX, 0 ) ; l i n e (MaxX, 0 , MaxX, MaxY); l i n e (MaxX, MaxY, 0 , MaxY); l i n e (0, MaxY, 0 , 0 ) ; /* F i n i s h up the d i sp lay parameters. */ CenX = MaxX / 2; CenY = MaxY / 2; Absci.MaxC = MaxX; Absci.CenC = CenX; Ordin.MaxC = MaxY; Ordin.CenC = CenY; /* . . . a n d draw the coordinate axes, i f w i t h i n the viewport . */ Absci.Mag = MaxX / Absc i . Len ; Ordin.Mag = MaxY / Ordin .Len; xO = CenX - A b s c i . D i r * Absci.Mag * Absc i .Cen ; i f ( (xO >= 0) && (xO <= MaxX) ) { l i n e ( ( int )xO, 0 , ( int )xO, MaxY); } /* endif */ yO = CenY - Ord in .D i r * Ordin.Mag * Ordin.Cen; i f ( (yO >= 0) && (yO <= MaxY) ) { l i n e (0, ( int )yO, MaxX, ( in t )yO) ; } /* endif */ /* I n i t i a l i z e screen v iewport , wi th c l i p p i n g enabled. */ setviewport (0, 0 , MaxX, MaxY, 1 ) ; /* Transfer f i l e header informat ion to output f i l e , . . . */ f o r ( i = 0 ; i < 4; i++) fscanf ( I n p u t F i l e , '"/.s", T S t r i n g ) ; f scanf ( I n p u t F i l e , '"/.If", ftxO); f o r ( i = 0 ; i < 2; i++) fscanf ( InputF i le , '"/.s", T S t r i n g ) ; f o r ( i = 0 ; i < 5; i++) fscanf ( InputF i le , '"/.If", ft(Seed[i])); f o r ( i = 0 ; i < 3 ; i++) fscanf ( InputF i le , '"/.s", T S t r i n g ) ; f o r ( i = 0 ; i < 4; i++) fscanf ( InputF i le , '"/.If", ft (Gr idSize [ i ] ) ) ; f o r ( i = 0 ; i < 5; i++) fscanf ( InputF i le , '"/.s", T S t r i n g ) ; switch (S l i ce4Ax is ) { case 0 : f p r i n t f (OutputF i le , "\\Sect ion on X at '/.f \ n " , xO) ; switch (S l i ce3Ax is ) { case Key_X: f p r i n t f (OutputF i le , "\\Sect ion on XDot at 7,f \ n " , SectCoord); Appendix H. Data, Manipulation Software 193 break; case Key_Y: fprintf (OutputFile, "\\Section on Theta at '/.f\n", SectCoord); break; case Key_Z: fprintf (OutputFile, "WSection on ThetaDot at '/.f\n", SectCoord); break; default: ; / * Nothing * / } / * endswitch * / break; case 1: fprintf (OutputFile, "WSection on XDot at y.f\n", xO) ; switch (Slice3Axis) { case Key_X: fprintf (OutputFile, "WSection on X at */.f \n", SectCoord); break; case Key_Y: fprintf (OutputFile, "WSection on Theta at '/.f\n", SectCoord); break; case Key_Z: fprintf (OutputFile, "WSection on ThetaDot at 7,f \n", SectCoord); break; default: ; / * Nothing * / } / * endswitch * / break; case 2: fprintf (OutputFile, "WSection on Theta at V.f \n", xO) ; switch (Slice3Axis) { case Key_X: fprintf (OutputFile, "WSection on X at '/,f\n", SectCoord); break; Appendix H. Data, Manipulation Software 194 case Key_Y: fprintf (OutputFile, "WSection on XDot at */.f \n", SectCoord); break; case Key.Z: fprintf (OutputFile, "WSection on ThetaDot at '/.f \n", SectCoord); break; default: ; / * Nothing * / } / * endswitch * / break; case 3: fprintf (OutputFile, "WSection on ThetaDot at 7,f \n", xO); switch (Slice3Axis) { case Key.X: fprintf (OutputFile, "WSection on X at y.f\n", SectCoord); break; case Key.Y: fprintf (OutputFile, "WSection on XDot at '/.f\n", SectCoord); break; case Key.Z: fprintf (OutputFile, "WSection on Theta at '/.f \n", SectCoord); break; default: ; / * Nothing * / } / * endswitch * / break; default: ; / * Nothing * / } / * endswitch * / fprintf (OutputFile, "\\Use GRAFTOOL to display plots.\n\n"); / * . . .and then set up for data sl icing. * / i f ( (SectCoord - Seed[Slice3AxisPtr+l]) >= 0) { Sectlndex = (char) ( (SectCoord - Seed[Slice3AxisPtr+l] + (GridSize[Slice3AxisPtr] / 2) ) / Appendix H. Data Manipulation Software GridSize[Slice3AxisPtr]); switch (Slice3Axis) { case Key_X: i f (Sectlndex '/, 2) { Mask = FaceOP; Sectlndex—; } / * endif * / else Mask = FaceON; break; case Key.Y: i f (Sectlndex */, 2) { Mask = FacelP; Sectlndex—; } / * endif * / else Mask = FacelN; break; case Key.Z: i f (Sectlndex '/, 2) { Mask = Face2P; Sectlndex—; } / * endif * / else Mask = Face2N; break; default: ; / * Nothing * / } / * endswitch * / } / * endif * / else { Sectlndex = (char) ( (SectCoord - Seed[Slice3AxisPtr+l] (GridSize[Slice3AxisPtr] / 2) ) / GridSize[Slice3AxisPtr]); switch (Slice3Axis) { case Key.X: i f ( (- Sectlndex) */. 2) { Mask = FaceOP; Sectlndex--; } / * endif * / else Mask = FaceON; break; case Key.Y: i f ( (- Sectlndex) 2) { Mask = FacelP; Sectlndex—; } / * endif * / else Mask = FacelN; break; Appendix H. Data Manipulation Software case Key_Z: i f ( (- Sectlndex) */. 2) { Mask = Face2P; Sectlndex--; } / * endif * / else Mask = Face2N; break; default: ; / * Nothing * / } / * endswitch * / } / * endelse * / while (fscanf (InputFile, '"/.hd '/.hd */.hd '/.hd 7.x ", ft(TVox.IK.Index[0]), ft(TVox.IK.Index[1]), ft(TVox.IK.Index[2]), ft(TVox.IK.Index[3]), ftTVox.SCode) == 5) { / * F irs t , test for a basic index value match, . . . .*/ i f (TVox.IK.Index[Slice3AxisPtr] == Sectlndex) { / * . . .and, i f found, pul l out the relevant face, . . . * / for ( i = 0, Picker = 0x0001, OutCode =0; i < 8; i++, Picker « = 1) { i f (Mask ft Picker) OutCode |= (Picker ft TVox.SCode); else OutCode <<= 1; } •/* endfor * / OutCode » = 4; TVox.SCode = OutCode; switch (Slice3Axis) { case Key_X: OutputStablePoints (ftTVox, Seed, GridSize, ftAbsci, ftOrdin, OutputFile); break; case Key_Y: OutputStablePoints (ftTVox, Seed, GridSize, ftAbsci, ftOrdin, OutputFile); break; case Key_Z: OutputStablePoints (ftTVox, Seed, GridSize, ftAbsci, ftOrdin, OutputFile); break; default: ; / * Nothing * / } / * endswitch * / } / * endif * / } / * endwhile * / Appendix H. Data, Manipulation Software 197 / * Toggle back to text mode to talk to the user. * / Key = GetScanCode (0); restorecrtmodeO; } while (Key != Esc); puts ("\nCreate another data section? Hit <Esc> to exit, any other key to continue."); Key = GetScanCode (0); fflush (stdin); / * Button up the output f i l e before looping or leaving, . . . * / fclose (OutputFile); } while (Key != Esc); puts ("\nProcess another f i le? Hit <Esc> to exit, any other key to continue."); Key = GetScanCode (0); fflush (stdin); / * . . .and close the input f i l e before new f i l e , or f inal exit. * / fclose (InputFile); } while (Key != Esc); closegraph(); assert (puttext (1, 1, 80, 25, ScreenBuf) ); gotoxy (CursorX, CursorY); } / * end Main * / void OutputStablePoints (Voxel *Vox, double *Seed, double *GridSize, Axis *Absci, Axis *0rdin, FILE *0File) { / * _ * / / * Draw crosses, centred at a l l stable vertices of the supplied * / / * 2-voxel, with axis lengths determined by the voxel size. * / / * */ double MagA, MagO, DimA, DimO, xO, yO, x, y, x l , y l , x2, y2; int i , j ; unsigned Picker = 0x0001; / * . . .and load local variables. * / MagA = Absci->Dir * Absci->Mag; MagO = Ordin->Dir * 0rdin->Mag; DimA = GridSize [Absci->DPtr]; DimO = GridSize[0rdin->DPtr]; / * . . .and find the "minimum coordinates" point on the voxel. * / xO = Seed[l+Absci->DPtr] + Vox->IK.Index[Absci->DPtr] * DimA; yO = Seed[l+0rdin->DPtr] + Vox->IK.Index[0rdin->DPtr] * DimO; / * Now, process the voxel vertices as offsets from this point. * / for (j = 0; j < 2; j++) { for (i = 0; i < 2; i++) { / * Check to see i f we need to draw this point, f i r s t . * / Appendix H. Data Manipulation Software i f (P icker ft Vox->SCode) { /* Need to drav i t ; f i n d endpoints f o r the l o c a t o r x = xO + ( i * DimA); y = yO + (j * DimO); x l = xO + ( i * DimA) - (0.5 * DimA); y l = yO + (j * DimO) - (0.5 * DimO); x2 = xO + ( i * DimA) + (0.5 * DimA); y2 = yO + (j * DimO) + (0.5 * DimO); /* . . . a n d output to the data f i l e , . . . */ f p r i n t f (OF i le , "*/.d y.d\n", 2, 2 ) ; f p r i n t f (OF i le , "V.f %f \ n " , x l , y) ; f p r i n t f (OF i le , "*/.f '/.f \ n " , x2 , y) ; f p r i n t f (OF i le , '7,d y.d\n", 2, 2 ) ; f p r i n t f (OF i le , '7.f 7.f \ n " , x , y l ) ; f p r i n t f (OF i le , "7,f V.f \ n " , x , y2) ; /* . . . t h e n convert to screen coord inates , . . . */ x = Absci->CenC + MagA * (x - Absci->Cen); y = Ordin->CenC + MagO * (y - Ordin->Cen); x l = Absci->CenC + MagA * ( x l - Absci->Cen); y l = Ordin->CenC + MagO * ( y l - Ordin->Cen); x2 = Absci->CenC + MagA * (x2 - Absci->Cen); y2 = Ordin->CenC + MagO * (y2 - Ordin->Cen); /* . . . a n d draw the screen l o c a t o r . */ i f ( ( x l >= 0) ftft ( x l <= Absci->MaxC) ftft (x2 >= 0) ftft (x2 <= Absci->MaxC) ) { i f ( ( y l >= 0) ftft ( y l <= Ordin->MaxC) ftft (y2 >= 0) ftft (y2 <= Ordin->MaxC) ) { l i n e ( ( i n t ) x l , ( i n t ) y , ( i n t ) x 2 , ( i n t ) y ) ; l i n e ( ( i n t ) x , ( i n t ) y l , ( i n t ) x , ( i n t ) y 2 ) ; } /* endif */ } /* endif */ } /* endif */ P i c k e r « = 1; } /* endfor */ } /* endfor */ r e t u r n ; } /* end OutputStablePoints */ Appendix H. Data Manipulation Software 199 H.3 3-D Format Conversion: SCATTER3.C /* */ /* SCATTER3.C — A program for creating 'scatter plot' GRAFTOOL files */ /* from 3-dimensional data files, as extracted from */ /* 4-space files by SLIC30F4.EXE. The 4-space files are created by */ /* STABVOL.EXE. */ /* */ /* Design ft coding: M.G. Henders */ /* Last modified: 24 June 1991 */ /* */ /* Module type: Externally-supported MAIN */ /* Module name: SCATTER3.C */ /* Module contents: */ /* */ /* Support modules: FILEPROC.C */ #define Main main /*- */ #include <stdio.h> #include <string.h> #include <alloc.h> #include <assert.h> #include <keystrks.h> #include <volume.h> double xl , y l , z l , x2, y2, z2; FILE *InputFile, *OutputFile; assert ((InputFilename = (char*) malloc(255 * sizeof(char))) != NULL); assert ((OutputFilename = (char*)malloc(255 * sizeof(char))) != NULL); void Main (int argc, char *argv[]) { char *InputFilename, *OutputFilename, *TString; Voxel TVox; double *Seed, GridSize[4], SectCoord, x, y, z, xO, yO, zO; l l , l , l , , , ; I I t il , 0utputF int i , j , k, Key = 0, SliceAxis, IndX, IndY, IndZ; long Count; unsigned Picker; FILE *fopenWarn (char *FileName, char *Mode, char *0utWarn); /* Allocate state vector and string space. */ assert ((Seed = (double*) malloc (5 * sizeof (double))) != NULL); assert ((InputFilename = (char*) malloc(255 * sizeof(char))) != NUL assert ((OutputFilename = (char*)malloc(255 * sizeof(char))) != NUL assert ((TString = (char*) malloc (255 * sizeof (char))) != NULL); /* Look for command-line arguments, and prompt if not found/no good. */ i f (argc != 3) Appendix H. Data Manipulation Software 200 { puts puts puts puts gets puts puts gets (OutputFilename); puts (""); } /* endif */ e l s e { s t rcpy (InputFilename, argv [1] ) ; } /* endelse */ "This i s a program to convert 3-space data f i l e s i n t o GRAFTOOL f i l e s , f o r " ) ; " s c a t t e r - p l o t t i n g . In genera l , the input and output f i lenames may b e " ) ; " s p e c i f i e d on the command l i n e , as SCATTER3 i n p u t f i l e o u t p u t f i l e . P l e a s e " ) ; "enter a name f o r the input f i l e , now.\n") ; InputFilename); "\nNow enter a name f o r the output f i l e . I f the se lec ted f i l e i s a l r e a d y " ) ; "present , you w i l l be prompted to conf i rm o v e r w r i t e . \ n " ) ; s t rcpy (OutputFilename, argv [2] ) ; /* Outer l e v e l loop , to do m u l t i p l e f i l e s . */ do { i f (Key) { puts ("Enter new input f i l e n a m e . " ) ; gets ( InputFilename); puts ("") puts ("Enter new output f i l ename." ) gets (OutputFilename); puts ("") } /* endif */ /* Try to open the f i l e s , b i t c h i n g to the user i f i t doesn't work. */ InputF i le = fopenWarn (InputFilename, " r t " , "WarnNotFound"); OutputFi le = fopenWarn (OutputFilename, "wt" , "WarnOverWrite"); /* Read f i l e header, . . . */ f o r ( i = 0 ; i < 2; i++) fscanf ( InputF i l e , "7,s", T S t r i n g ) ; f scanf ( I n p u t F i l e , "7.s", TS t r ing ) ; i f (Istrcmp (TSt r ing , "X") ) S l i c e A x i s = 0 ; e l se i f (Istrcmp (TSt r ing , "XDot") ) S l i c e A x i s = 1; e l se i f (Istrcmp (TSt r ing , "Theta") ) S l i c e A x i s = 2; e l se S l i c e A x i s = 3 ; fscanf ( I n p u t F i l e , "7.s '/.If", T S t r i n g , ftSectCoord) ; f o r ( i = 0 ; i < 2 f o r ( i = 0 ; i < 5 f o r ( i = 0; i < 3 f o r ( i = 0; i < 4 f o r ( i = 0; i < 5 i++) fscanf ( I n p u t F i l e , '"/.s", T S t r i n g ) ; i++) fscanf ( I n p u t F i l e , '"/.If", &(Seed[i] ) ) ; i++) fscanf ( I n p u t F i l e , "V.s", T S t r i n g ) ; i++) fscanf ( I n p u t F i l e , '"/.If", & (Gr idSize [ i ] ) ) ; i++) fscanf ( I n p u t F i l e , '"/.s", T S t r i n g ) ; Appendix H. Data Manipulation Software 201 switch (SliceAxis) { case 0: fprintf (OutputFile, "WSection on X at '/.f \n", SectCoord); IndX = 1; IndY = 2; IndZ = 3; break; case 1: fprintf (OutputFile, "WSection on XDot at y.f\n", SectCoord); IndX = 0; IndY = 2; IndZ = 3; break; case 2: fprintf (OutputFile, "WSection on Theta at '/.f\n", SectCoord); IndX = 0; IndY = 1; IndZ = 3; break; case 3: fprintf (OutputFile, "WSection on ThetaDot at V.f \n", SectCoord); IndX = 0; IndY = 1; IndZ = 2; break; default: ; / * Nothing * / } / * endswitch * / fprintf (OutputFile, "\\Use GRAFT00L for display.\n\n"); puts ("Processing data f i l e . Please wait.. ."); Count = 0; while (f scanf (InputFile, '7.hd */.hd '/.hd */.hd 7.x", k(TVox.IK.Index[0]), 4(TVox.IK.Index[1]), k(TVox.IK.Index[2]), &(TVox.IK.Index[3]), ftTVox.SCode) == 5) {. / * Write coords for a 3-space cross at each stable vertex point. * / Picker = 0x0001; / * Voxel coordinates * / x = Seed[l+IndX] + TVox.IK.Index[IndX] * GridSize[IndX] y = Seed[l+IndY] + TVox.IK.Index[IndY] * GridSize[IndY] z = Seed[l+IndZ] + TVox.IK.Index[IndZ] * GridSize[IndZ] for (k = 0; k < 2; k++) { for (j = 0; j < 2; j++) { for (i = 0; i < 2; i++) { / * Check vertex stability. * / Appendix H. Data, Manipulation Software 202 i f (Picker & TVox.SCode) { / * Vertex coordinates * / xO = x + i * GridSize[IndX]; yO = y + j * GridSize[IndY]; zO = z + k * GridSize[IndZ]; / * Locator endpoint coordinates * / xl = xO - 0.5 * GridSize[IndX]; y l = yO - 0.5 * GridSize[IndY]; z l = zO - 0.5 * GridSize[IndZ]; x2 = xO + 0.5 * GridSize[IndX]; y2 = yO + 0.5 * GridSize[IndY]; z2 = zO + 0.5 * GridSize[IndZ]; / * Output locator to f i l e . * / fprintf (OutputFile, '7.d */.d\n", 2, 3); fprintf (OutputFile, '7,f */.f It\n", x l , yO, zO); fprintf (OutputFile, "%f */.f y.f\n", x2, yO, zO) ; fprintf (OutputFile, '7,d '/.d\n", 2, 3); fprintf (OutputFile, '7,f 7.f V.f \n", xO, y l , zO) ; fprintf (OutputFile, '7.f V.f */.f \n", xO, y2, zO) ; fprintf (OutputFile, '7,d y.d\n", 2, 3); fprintf (OutputFile, '7,f 'I.i */.f \n", xO, yO, zl) ; fprintf (OutputFile, "V.f V.f V . f \ n " , xO, yO, z2) ; } / * endif * / Picker <<= 1; } / * endfor * / } / * endfor * / } / * endfor * / Count++; i f ( (Count '/. 1024) == 0) printf ("I"); else i f ( (Count 1024) == 512) printf ("+"); } / * endwhile * / puts ("\nProcess another f i le? Hit <Esc> to exit, any other key to continue."); Key = GetScanCode (0); fflush (stdin); / * Button up the f i les before looping or leaving, . . . * / fclose (InputFile); fclose (OutputFile); } while (Key != Esc); } / * end Main * /
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Nonlinear dynamics of the inverted pendulum under linear-feedback...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Nonlinear dynamics of the inverted pendulum under linear-feedback stabilization Henders, Michael G. 1991
pdf
Page Metadata
Item Metadata
Title | Nonlinear dynamics of the inverted pendulum under linear-feedback stabilization |
Creator |
Henders, Michael G. |
Publisher | University of British Columbia |
Date Issued | 1991 |
Description | Analysis is conducted on a linear control system used for the stabilization of an inverted pendulum, examining the behaviour of the system throughout its global state space. The dynamic equations of the inverted pendulum are obtained, and a linear state-feedback controller is designed, based on this system model, for its stabilization. The behaviour of the resulting closed-loop system is examined through simulation, in two stages. In the first stage, the feedback controller is modified slightly to decouple the pendulum dynamics from those of its supporting cart, permitting examination of the pendulum dynamics in a two-dimensional state space, using phase plane analysis. In this way, a globally valid description of the pendulum dynamics is obtained, revealing the existence of two state attractors in this modified system. The feedback controller is seen to create a stable point attractor in the phase plane, as intended, but the basin of attraction of this point is of unusual shape and limited extent, due to the nonlinearity of the pendulum dynamics. All initial system states outside this basin are drawn to the second state attractor, which is a stable limit cycle. As the second stage of analysis, the system is examined with its original controller; characterization of the overall dynamics then requires a four-dimensional state space. The continued existence of a point attractor at the state space origin is demonstrated, and the basin of attraction of this point is mapped numerically, using a state space sampling algorithm, which might readily be applied to the analysis of other dynamic systems. This basin is presented as a series of three-dimensional sections of the full four-dimensional region; taken together, these sections outline the entire basin map. The boundaries of the origin-bound basin are seen to be highly irregular, and potentially fractal in nature. The existence of an infinite set of such four-dimensional stability regions is demonstrated, and complex, potentially chaotic, behaviour is observed outside the basins of attraction; it is not clear whether or not a four-dimensional counterpart of the limit cycle (from the first stage analysis) exists. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-11-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0098490 |
URI | http://hdl.handle.net/2429/29971 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1991_A7 H46.pdf [ 14.07MB ]
- Metadata
- JSON: 831-1.0098490.json
- JSON-LD: 831-1.0098490-ld.json
- RDF/XML (Pretty): 831-1.0098490-rdf.xml
- RDF/JSON: 831-1.0098490-rdf.json
- Turtle: 831-1.0098490-turtle.txt
- N-Triples: 831-1.0098490-rdf-ntriples.txt
- Original Record: 831-1.0098490-source.json
- Full Text
- 831-1.0098490-fulltext.txt
- Citation
- 831-1.0098490.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0098490/manifest