User manual

ManualsBrandsNational Instruments ManualsOtherCalculator 370760B-01

MATRIXx

Xmath

Xµ Manual

MATRIXx Xmath Basics

The MATRIXx products and related items have been purchased from Wind

River Systems, Inc. (formerly Integrated Systems, Inc.). These reformatted user

materials may contain references to those entities. Any trademark or copyright

notices to those entities are no longer valid and any references to those entities

as the licensor to the MATRIXx products and related items should now be

considered as referring to National Instruments Corporation.

National Instruments did not acquire RealSim hardware (AC-1000, AC-104,

PCI Pro) and does not plan to further develop or support RealSim software.

NI is directing users who wish to continue to use RealSim software and hardware

to third parties. The list of NI Alliance Members (third parties) that can provide

RealSim support and the parts list for RealSim hardware are available in our

online KnowledgeBase. You can access the KnowledgeBase at

www.ni.com/support.

NI plans to make it easy for customers to target NI software and hardware,

including LabVIEW real-time and PXI, with MATRIXx in the future.

For information regarding NI real-time products, please visit

www.ni.com/realtime or contact us at matrixx@ni.com.

April 2004 Edition

Part Number 370760B-01

Summary of content (384 pages)

PAGE 1
MATRIXx TM Xmath Xµ Manual TM MATRIXx Xmath Basics The MATRIXx products and related items have been purchased from Wind River Systems, Inc. (formerly Integrated Systems, Inc.). These reformatted user materials may contain references to those entities. Any trademark or copyright notices to those entities are no longer valid and any references to those entities as the licensor to the MATRIXx products and related items should now be considered as referring to National Instruments Corporation.
PAGE 2
Support Worldwide Technical Support and Product Information ni.
PAGE 3
Important Information Warranty The media on which you receive National Instruments software are warranted not to fail to execute programming instructions, due to defects in materials and workmanship, for a period of 90 days from date of shipment, as evidenced by receipts or other documentation. National Instruments will, at its option, repair or replace software media that do not execute programming instructions if National Instruments receives notice of such defects during the warranty period.
PAGE 4
Contents 1 Introduction 1 1.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Manual Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 How to avoid really reading this Manual . . . . . . . . . . . . . . . . . . . 3 2 Overview of the Underlying Theory 2.1 2.2 5 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
PAGE 5
CONTENTS iv 2.2.5 2.3 2.4 2.5 Obtaining Robust Control Models for Physical Systems . . . . . . 28 H∞ and H2 Design Methodologies . . . . . . . . . . . . . . . . . . . . . . 29 2.3.1 H∞ Design Overview . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.3.2 Assumptions for the H∞ Design Problem . . . . . . . . . . . . . . 32 2.3.3 A Brief Review of the Algebraic Riccati Equation . . . . . . . . . . 33 2.3.4 Solving the H∞ Design Problem for a Special Case . . . . . . . . . 36 2.3.
PAGE 6
CONTENTS 2.6 v Model Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 2.6.1 Truncation and Residualization . . . . . . . . . . . . . . . . . . . . 65 2.6.2 Balanced Truncation . . . . . . . . . . . . . . . . . . . . . . . . . . 65 2.6.3 Hankel Norm Approximation . . . . . . . . . . . . . . . . . . . . . 68 3 Functional Description of Xµ 71 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.2 Data Objects . . . . . . . .
PAGE 7
CONTENTS vi 3.5 System Interconnection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 3.6 H2 and H∞ Analysis and Synthesis . . . . . . . . . . . . . . . . . . . . . . 95 3.6.1 Controller Synthesis . . . . . . . . . . . . . . . . . . . . . . . . . . 95 3.6.2 System Norm Calculations 3.7 3.8 . . . . . . . . . . . . . . . . . . . . . . 105 Structured Singular Value (µ) Analysis and Synthesis . . . . . . . . . . . 107 3.7.1 Calculation of µ . . . . . . . . . . . . . . . . . . . .
PAGE 8
CONTENTS 4.2 vii 4.1.5 µ Analysis of the H∞ Controller . . . . . . . . . . . . . . . . . . . 138 4.1.6 Fitting D-scales for the D-K Iteration . . . . . . . . . . . . . . . . 140 4.1.7 Design Iteration #2 . . . . . . . . . . . . . . . . . . . . . . . . . . 143 4.1.8 Simulation Comparison with a Loopshaping Controller . . . . . . . 146 A Simple Flexible Structure Example . . . . . . . . . . . . . . . . . . . . 153 4.2.1 The Control Design Problem . . . . . . . . . . . . . . . . . . . . .
PAGE 9
CONTENTS viii A.3 System Response Functions . . . . . . . . . . . . . . . . . . . . . . 398 A.4 System Interconnection . . . . . . . . . . . . . . . . . . . . . . . . 399 A.5 Model Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 399 A.6 H2 and H∞ Analysis and Synthesis . . . . . . . . . . . . . . . . . 399 A.7 Structured Singular Value (µ) Analysis and Synthesis . . . . . . .
PAGE 10
Chapter 1 Introduction Xµ is a suite of Xmath functions for the modeling, analysis and synthesis of linear robust control systems. Robust control theory has developed rapidly during the last decade to the point where a useful set of computational tools can be used to solve a wide range of control problems. This theory has already been applied to a wide range of practical problems. This manual describes the Xµ functions and presents a demonstration of their application.
PAGE 11
CHAPTER 1. INTRODUCTION 2 Notation pdm Dynamic System Meaning Xmath parameter dependent matrix data object Xmath dynamic system data object Code examples and function names are set in typewriter font to distinguish them from narrative text. 1.2 Manual Outline Chapter 2 outlines the applicable robust control theory. Perturbation models and linear fractional transformations form the basis of the modeling framework.
PAGE 12
1.3. HOW TO AVOID REALLY READING THIS MANUAL 1.3 3 How to avoid really reading this Manual The layout of the manual proceeds from introduction to background to syntax detail to application descriptions. This may be tediously theoretical for some. If you are one of those that considers reading the manual as the option of last resort1 then go directly to the applications (Chapter 4). If you have no prior Xmath experience then skimming through Chapter 3 is essential.
PAGE 13
Chapter 2 Overview of the Underlying Theory 2.1 Introduction The material covered here is taken from a variety of sources. The basic approach is described by Doyle [1, 2], and further elaborated upon by Packard [3]. Summaries have also appeared in work by Smith [4] and others. Motivating background can be found in the early paper by Doyle and Stein [5]. An overview of the robust control approach, particularly for process control systems, is given by Morari and Zafiriou [6].
PAGE 14
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 6 several studies involving process control applications, particularly high purity distillation columns. These are detailed by Skogestad and Morari in [15, 16, 17, 18] Section 2.2 introduces robust control perturbation models and linear fractional transformations. Weighted H∞ design is covered in Section 2.3. The analysis of closed loop systems with the structured singular value (µ) is overviewed in Section 2.4. Section 2.
PAGE 15
2.1. INTRODUCTION 7 - z y P11 P12 P21 P22 v u Figure 2.1: The generic robust control model structure Pn Trace(M ) trace of M ( i=1 Mii ) Block diagrams will be used to represent interconnections of systems. Consider the example feedback interconnection shown in Fig. 2.1. Notice that P has been partitioned into four parts. This diagram represents the equations, z y = = P11 v + P12 u P21 v + P22 u v = ∆z.
PAGE 16
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 8 is interpreted to mean that the signal y is the sum of the response of system P21 to input signal v and system P22 to input signal u. In general, we will not be specific about the representation of the system P . If we do need to be more specific about P , then P (s) is the Laplace representation and p(t) is the impulse response. Note that Figure 2.1 is drawn from right to left.
PAGE 17
2.1. INTRODUCTION 9 Euclidean norm. Given,   x1   x =  ...  , xn the Euclidean (or 2-norm) of x, denoted by kxk, is defined by, kxk = n X !1/2 |xi | . i=1 Many other norms are also options; more detail on the easily calculated norms can be found in the on-line help for the norm function. The term spatial-norm is often applied when we are looking at norms over the components of a vector. Now consider a vector valued signal,   x1 (t)   x(t) =  ...  .
PAGE 18
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 10 For persistent signals, where the above norm is unbounded, we can define a power norm, kxi (t)k = 1 lim T →∞ 2T Z !1/2 T |xi (t)| dt 2 −T . (2.1) The above norms have been defined in terms of a single component, xi (t), of a vector valued signal, x(t). The choice of spatial norm determines how we combine these components to calculate kx(t)k. We can mix and match the spatial and time parts of the norm of a signal.
PAGE 19
2.1. INTRODUCTION 11 signals. Strictly speaking, signals in H2 or H2⊥ are not defined on the ω axis. However we usually consider them to be by taking a limit as we approach the axis. A slightly more specialized set is RL2 , the set of real rational functions in L2 . These are strictly proper functions with no poles on the imaginary axis. Similarly we can consider RH2 as strictly proper stable functions and RH⊥ 2 as strictly proper functions with no poles in Re(s)< 0.
PAGE 20
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 12 where σmax denotes the maximum singular value. Not all matrix norms are induced from vector norms. The Froebenius norm (square root of the sum of the squares of all matrix elements) is one such example. Now consider the case where P (s) is a dynamic system and we define an induced norm from L2 to L2 as follows. In this case, y(s) is the output of P (s)u(s) and kP (s)k = max u(s)∈L2 ky(s)k2 .
PAGE 21
2.2. MODELING UNCERTAIN SYSTEMS 13 The set of all systems with bounded ∞-norm is denoted by L∞ . We can again split this into stable and unstable parts. H∞ denotes the stable part; those systems with |P (s)| finite for all Re(s)> 0. This is where the name “H∞ control theory” originates, and we often call this norm the H∞ -norm. Again we can restrict ourselves to real rational functions, so RL∞ is the set of proper transfer functions with no poles on the ω axis.
PAGE 22
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 14 robust control model is therefore a set description and we hope that some members of this set capture some of the uncertain or unmodeled aspects of our physical system. For example, consider the “uncertain” model illustrated in Figure 2.2. This picture is equivalent to the input-output relationship, y = [ (I + ∆Wm )Pnom ] u. y (2.2) Wm ? +m s Pnom u Figure 2.
PAGE 23
2.2. MODELING UNCERTAIN SYSTEMS 15 as specifying a maximum percentage error between Pnom and every other element of P. The system Pnom (s) is the element of P that comes from ∆ = 0 and is called the nominal system. In otherwords, for ∆ = 0, the input-output relationship is y(s) = Pnom (s) u(s). As ∆ deviates from zero (but remains bounded in size), the nominal system is multiplied by (I + ∆Wm (s)).
PAGE 24
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 16 1 Imaginary 0.5 0 -0.5 -1 -0.5 0 0.5 1 1.5 Real Figure 2.
PAGE 25
2.2. MODELING UNCERTAIN SYSTEMS y u +?j 17 Wa P0 u u +j r 6, Figure 2.4: Unity gain negative feedback for the example system, P0 + ∆Wa 2.2.2 Linear Fractional Transformations A model is considered to be an interconnection of lumped components and perturbation blocks. In this discussion we will denote the input to the model by u, which can be a vector valued signal representing input signals such as control inputs, disturbances, and noise.
PAGE 26
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 18 - 1 . z .. m P11 P12 y P21 P22 v u Figure 2.5: Generic LFT model structure including perturbations,∆ A generic model structure, referred to as a linear fractional transformation (LFT), overcomes the difficulties outlined above. The LFT model is equivalent to the relationship, y = P21 ∆(I − P11 ∆)−1 P12 + P22 u, (2.6) where the ∆ is the norm bounded perturbation. Figure 2.
PAGE 27
2.2. MODELING UNCERTAIN SYSTEMS 19 in an LFT format. The open-loop system is described by, y = Fu (Polp , ∆)u, where Polp = 0 Wa . I P0 The unity gain, negative feedback configuration, illustrated in Figure 2.4 (and given in Equation 2.5) can be described by, y = Fu (Gclp , ∆)r, where Gclp = −Wa (I + P0 )−1 Wa (I + P0 )−1 (I + P0 )−1 P0 (I + P0 )−1 Figure 2.5 also shows the perturbation, ∆ as block structured. In otherwords, ∆ = diag(∆1 , . . . , ∆m ). (2.
PAGE 28
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 20 The issue of the invertibility of (I − P11 ∆) is fundamental to the study of the stability of a system under perturbations. We will return to this question in much more detail in Section 2.4. It forms the basis of the µ analysis approach. Note that Equation 2.7 indicates that we have m blocks, ∆i , in our model. For notational purposes we will assume that each of these blocks is square.
PAGE 29
2.2. MODELING UNCERTAIN SYSTEMS 21 w ? Wn j +j? +? y Wu u Pnom u Figure 2.6: Example model: multiplicative output perturbation with weighted output noise References to a robust control model will imply a description of the form given in Equation 2.8. As a example, consider one of the most common perturbation model descriptions, illustrated in Figure 2.6. This model represents a perturbed system with bounded noise at the output.
PAGE 30
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 22 Robust control models are therefore set descriptions. In the analysis of such models it is also assumed that the unknown inputs belong to some bounded set. Several choices of set for the unknown signals can be made, leading to different mathematical problems for the analysis. Unfortunately not all of them are tractable. The following section discusses the assumptions typically applied to the robust control models. 2.2.
PAGE 31
2.2. MODELING UNCERTAIN SYSTEMS 23 Packard [19] discuss the implications of this assumption on robust control theory and we briefly touch upon this in Section 2.4.6. The most common assumption is that ∆ is an unknown, norm-bounded, linear time-invariant system. Systems often do not fall neatly into one of the usual choices of ∆ discussed above. Consider a nonlinear system linearized about an operating point.
PAGE 32
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 24 systems. We will now look at other possible perturbation structures. For more detail on these structures (in the complex case) refer to Packard and Doyle [20]. Consider a blocks which are of the form scalar × identity, where the scalar is unknown. In the following we will include q of these blocks in ∆. The definition of ∆ is therefore modified to be, o n ∆ = diag(δ1 I1 , . . . , δq Iq , ∆1 , . . . , ∆m ) dim(Ij ) = lj × lj , dim(∆i ) = ki × ki .(2.
PAGE 33
2.2. MODELING UNCERTAIN SYSTEMS = Cz −1 (I − z −1 A)−1 B + D = Fu (Pss , z −1 I), 25 where Pss is the real valued matrix, Pss = A B , C D and the scalar × identity, z −1 I, has dimension equal to the state dimension of P (z). This is now in the form of an LFT model with a single scalar × identity element in the upper loop. One possible use of this is suggested by the following. Define, ∆ = δInx δ ∈ C , where nx is the state dimension.
PAGE 34
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 26 where the transport delay, Td , and the combustion lag, Tc , are approximately, Td = 252 v and Tc = 202 . v For the purposes of our example we want to design an air/fuel ratio controller that works for all engine speeds in the range 2,000 to 6,000 rpm. We will use a first order Padé approximation for the delay and express the Td and Tc relationships in an LFT form with a normalized speed deviation, δv .
PAGE 35
2.2. MODELING UNCERTAIN SYSTEMS 27 Putting all the pieces together gives an engine model in the following fractional form. P (s) = Fu (Pmod , ∆), where  Pmod 7.14(s − 19.8) s2 −9.9 s −0.9(s − 15.8)(s − 19.8) s3 −15.87 s    = 0   −27.75 s  141.5(1 + 1.006s) s(s + 1)    ,   9.9 −17.82(s − 15.8)(1 + 1.006s) s(s + 1) and ∆ ∈ B∆, with the structure defined as, ∆= δv I2 δv ∈ R .
PAGE 36
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 28 and ∆ ∈ B∆, with the structure defined as, ∆= diag(δv I2 , ∆1 ) δv ∈ R, ∆1 ∈ C . Note that this is an LFT with a repeated real-valued parameter, δv (|δv | ≤ 1), and a complex perturbation, ∆1 (|∆1 | ≤ 1). Note that as R ⊂ C, assuming that ∆ ∈ C n×n , always covers the case where some of the ∆i (or δj ) are more appropriately modeled as real-valued.
PAGE 37
2.3. H∞ AND H2 DESIGN METHODOLOGIES 29 An area of work, known as identification in H∞ , looks at experimental identification techniques which minimize the worst case H∞ error between the physical system and the model. The following works address this issue: [27, 28, 29, 30, 31, 32, 33, 34, 35]. Applying the more standard, probabilistically based, identification techniques to uncertain systems is also receiving attention.
PAGE 38
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 30 e Ps w ( ) y -K s u ( ) Figure 2.7: LFT configuration for controller synthesis, G(s) = Fl [P (s), K(s)] Note that the interconnection structure, P (s), given here, differs from that discussed in the previous section. Here we set up P (s) so that the input, w, is the unknown signals entering our system. Typical examples would be sensor noise, plant disturbances or tracking commands.
PAGE 39
2.3. H∞ AND H2 DESIGN METHODOLOGIES 31 extend these approaches to the case where P (s) is replaced by Fu (P (s), ∆), ∆ ∈ B∆. 2.3.1 H∞ Design Overview Again, recall from Section 2.1.2, the H∞ is norm of G(s) is, kG(s)k∞ = sup σmax [G(ω)]. ω The H∞ norm is the induced L2 to L2 norm. Therefore minimizing the H∞ norm of G(s) will have the effect of minimizing the worst-case energy of e over all bounded energy inputs at w. Consider γ(K) to be the closed loop H∞ norm achieved for a particular controller K.
PAGE 40
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 32 2.3.2 Assumptions for the H∞ Design Problem There are several assumptions required in order to achieve a well-posed design problem. The DGKF paper gives a state-space solution to the H∞ design problem and we will use a similar notation here. Consider the open loop state-space representation of P (s), partitioned according to the signals shown in Figure 2.7,  A P (s) =  C1 C2 B1 D11 D21  B2 D12 . D22 (2.
PAGE 41
2.3. H∞ AND H2 DESIGN METHODOLOGIES 33 that the effect of all disturbances, w, at every frequency, can be measured by the controller. If either of these conditions are not met then the problem could be ill-posed. It is possible to violate these conditions by using pure integrators as design weights. While this could still give a meaningful design problem, solution via the state-space H∞ approach requires that an approximation be used for the integrator weight.
PAGE 42
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 34 We have partitioned the matrix into two n × n blocks, X1 and X2 . If X1 is invertible, then X = X2 X1−1 , is the unique, stabilizing solution to the ARE. The ability to form X doesn’t depend on the particular choice of X1 and X2 . Given a Hamiltonian, H, we say that H ∈ dom(Ric) if H has no ω axis eigenvalues and the associated X1 matrix is invertible. Therefore, if H ∈ dom(Ric), we can obtain a unique stabilizing solution, X.
PAGE 43
2.3. H∞ AND H2 DESIGN METHODOLOGIES 35 Choose γ > 0 and form the following Hamiltonian matrix, H= A γ −2 BB T . −C T C −AT The following lemma gives a means of checking whether or not kP (s)k∞ < γ. A proof of this lemma is given in DGKF although it is based on the work of Anderson [60], Willems [61] and Boyd et al. [62].
PAGE 44
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 36 any with a zero real part. In practice we must use a tolerance to determine what is considered as a zero real part. Finding a basis for the stable subspace of H involves either an eigenvalue or Schur decomposition. Numerical errors will be introduced at this stage. In most cases using an eigenvalue decomposition is faster and less accurate than using a Schur decomposition. Similarly, forming X = X2 X1−1 will also introduce numerical errors.
PAGE 45
2.3. H∞ AND H2 DESIGN METHODOLOGIES 37 and the components, C1 x and D12 u are orthogonal. D12 is also assumed to be normalized. This essentially means that there is no cross-weighting between the state and input penalties. Assumption (iv ) is the dual of this; the input and unknown input (disturbance and noise) affect the measurement, y, orthogonally, with the weight on the unknown input being unity.
PAGE 46
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 38 L∞ = Z∞ Â∞ = = −Y∞ C2T (I − γ −2 Y∞ X∞ )−1 A + γ −2 B1 B1T X∞ + B2 F∞ + Z∞ L∞ C2 . Actually, the above formulation can be used to parametrize all stabilizing controllers which satisfy, kG(s)k∞ < γ. This can be expressed as an LFT. All such controllers are given by, K∞ = Fl (M∞ , Q), where,  M∞ Â∞ =  F∞ −C2 −Z∞ L∞ 0 I  Z∞ B2 , I 0 and Q satisfies: Q ∈ RH∞ , kQk∞ < γ. Note that if Q = 0 we get back the controller given in Theorem 3.
PAGE 47
2.3. H∞ AND H2 DESIGN METHODOLOGIES 39 a) Choose γ ≥ γopt b) Form H∞ and J∞ c) Check that H∞ ∈ dom(Ric) and J∞ ∈ dom(Ric). d) Calculate X∞ = Ric(H∞ ) and Y∞ = Ric(J∞ ) e) Check that X∞ ≥ 0 and Y∞ ≥ 0 f) Check that ρ(X∞ Y∞ ) < γ 2 g) Reduce γ and go to step b). The value of γ can be reduced until one of the checks at steps c), e) or f) fails. In this case, γ < γopt and we use the X∞ and Y∞ of the lowest successful γ calculation to generate the controller.
PAGE 48
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 40 eigenvalue of X∞ (and Y∞ ) is displayed. The ultimate test of the software is to form the closed loop system and check both its stability and norm. We strongly suggest that the user always perform this step. The numerical issues discussed above are very unlikely to arise in low order systems. Experience has shown that systems with very lightly damped modes are more susceptible to numerical problems than those with more heavily damped modes.
PAGE 49
2.3. H∞ AND H2 DESIGN METHODOLOGIES 41 framework. We again assume the simplifying assumptions used in Section 2.3.4 The H2 design solution is obtained (at least conceptually) from the H∞ design procedure by setting γ = ∞ and using the resulting central controller. It is interesting to compare the H∞ solution, given above, and the H2 solution given below. Define two Hamiltonians, H2 and J2 , by, A −B2 B2T , −C1T C1 −AT −C2T C2 AT .
PAGE 50
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 42 G(s) e Ps w ( ) -K s ( ) Figure 2.8: Closed loop system, G(s), for performance analysis initial choice of γ for the H∞ design procedure. We will see later (Section 2.5) that it can also be used to initialize the D-K iteration procedure when an open-loop H∞ design is poorly conditioned. 2.4 2.4.1 µ Analysis Measures of Performance Section 2.3 presented design performance objectives in terms of the norm (H2 or H∞ ) of a closed loop system.
PAGE 51
2.4. µ ANALYSIS 43 descriptions are considered, where B again denotes the unit ball. ) Z T 1 | w(t) |2 dt ≤ 1 w lim T →∞ 2T −T ( ) Z ∞ 2 2 w kwk2 = | w(t) | dt ≤ 1 ( Power : BP = Energy : BL2 = ( Magnitude : BL∞ = −∞ (2.12) (2.13) ) w kwk∞ = ess sup |w(t)| ≤ 1 (2.14) t These norms are defined for scalar signals for clarity. The choice of w and e as the above sets defines the performance criterion. The performance can be considered as a test on the corresponding induced norm of the system.
PAGE 52
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 44 z e 1 . .. m G11 G12 G21 G22 v w Figure 2.9: Perturbed closed loop system for stability analysis the induced norms, the reader is referred to Doyle [2]. The major advantage of choosing BP or BL2 is that the test for the performance can be considered in terms of the same norm as stability. This has significant advantages when we are considering performance and stability in the presence of perturbations, ∆. 2.4.
PAGE 53
2.4. µ ANALYSIS 45 Consider the case where the model has only one full ∆ block (m = 1 and q = 0 in Equation 2.9). This is often referred to as unstructured, and the well known result (refer to Zames [67] and Doyle and Stein [5]) is given in the following lemma. Lemma 6 (Robust Stability, Unstructured) Fu (G(s), ∆) is stable for all ∆, k∆k∞ ≤ 1, if and only if kG11 (s)k∞ < 1. A generalization of the above is required in order to handle Fu (G(s), ∆) models with more than one full ∆ block.
PAGE 54
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 46 Fu (G(s), ∆) stable for all ∆ ∈ B∆ if and only if kµ(G11 (s))k∞ < 1. where kµ(G11 (s))k∞ = sup µ[G11 (ω)]. ω The use of this notation masks the fact that µ is also a function of the perturbation structure, ∆. The above definition of µ applies to the more general block structure given in Section 2.2.4. We can even consider the some of the blocks to be real valued, rather than complex valued.
PAGE 55
2.4. µ ANALYSIS if and only if 47 kµ(G(s))k∞ < 1, b where µ is taken with respect to an augmented structure ∆, o n b = diag(∆, ∆) ˆ ∆ ∈ ∆, ∆ ˆ = C dim(w)×dim(e) . ∆ ˆ can be thought of as a “performance block” The additional perturbation block, ∆ appended to the ∆ blocks used to model system uncertainty. This result is the major benefit of the choice of input and output signal norms; the norm test for performance is the same as that for stability.
PAGE 56
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 48 For the other extreme consider a single full block (∆ = {∆ | ∆ ∈ C n×n }); the definition of µ is now the same as that for the maximum singular value, ∆ = {∆ | ∆ ∈ C n×n } ⇒ µ(M ) = σmax (M ). Observe that every possible block structure, ∆, contains {λI | λ ∈ C} as a perturbation; and every possible block structure, ∆, is contained in C n×n . These particular block structures are the boundary cases.
PAGE 57
2.4. µ ANALYSIS 49 Actually, the lower bound is always equal to µ but the implied optimization has local maxima which are not global. For the upper bound Safonov and Doyle [72], have shown that finding the infimum is a convex problem and hence more easily solved. However the bound is equal to µ only in certain special cases. Here we use the infimum rather that the miminum because D may have a element which goes to zero as the maximum singular value decreases.
PAGE 58
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 50 and two block structures, ∆1 (compatible with M11 ) and ∆2 (compatible with M22 ). There are two perturbed subsystems that we can study here: Fu (M, ∆1 ), where ∆1 is closed in a feedback loop around M11 ; and Fl (M, ∆2 ), where ∆2 is closed in a feedback loop around M22 . We have already seen that in the case of a dynamic system, the robust stability of Fu (M, ∆1 ) is analyzed by checking that µ1 (M11 ) < 1.
PAGE 59
2.4. µ ANALYSIS 2.4.6 51 State-space Robustness Analysis Tests We will look at some more advanced state-space approaches to the analysis of robust performance. Most users of the software will concentrate on the more common frequency domain analysis methods covered in Section 2.4.3. The analysis tests given here can be implemented with the Xµ functions and the more advanced user can use these to study the robustness of systems with respect to time-varying and nonlinear perturbations.
PAGE 60
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 52 Therefore µ1 (A) < 1 is equivalent to our system being stable. Furthermore, the maximum modulus theorem for a stable system tells us that, kP (z)k∞ = sup σmax (P (z)) |z|≥1 = sup σmax (Fu (Pss , z −1 I) |z −1 |≤1 = sup ∆1 ∈B∆1 µ2 (Fu (Pss , ∆1 ), if ∆2 is defined as a single full block of dimensions matching the input-output dimensions of P (z). The main loop theorem (Theorem 9) immediately suggests the following result.
PAGE 61
2.4. µ ANALYSIS 53 - z,1 I x(k +1) e Gss z - x(k ) w v Figure 2.
PAGE 62
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 54 Note that the nominal system is given by, Gnom (z) = Fu A C1 B1 D11 ,z −1 I , and the perturbed system is, G(z) = Fu (Fl (G, ∆), z −1 I). We assume that ∆ is an element of a unity norm bounded block structure, ∆ ∈ B∆. For the µ analysis we will define a block structure corresponding to Gss , ∆s = diag(δ1 Inx , ∆2 , ∆) δ1 ∈ C, ∆2 ∈ C dim(w)×dim(e) , ∆ ∈ ∆ .
PAGE 63
2.4. µ ANALYSIS 55 iii) There exists a constant β ∈ [0, 1] such that for each fixed ∆ ∈ B∆, G(z) is stable (robust and for zero initial state response, e satisfies kek2 ≤ β kwk2 performance). The frequency domain µ test is implemented by searching for the maximum value of µ over a user specified frequency grid. Theorem 11 shows that this is equivalent to a single, larger, µ test. There are subtle distinctions between the two tests. As we would expect, calculation of the larger µ test is more difficult.
PAGE 64
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 56 i) State-space upper bound: inf σmax [Ds Gss Ds−1 ] < 1; Ds ∈Ds ii) Frequency domain, constant D, upper bound: inf max σmax [Dp Fu (Gss , eω Inx )Dp−1 ] < 1; Dp ∈Dp ω∈[0,2π] iii) Frequency domain upper bound: max inf σmax [Dp Fu (Gss , eω Inx )Dp−1 ] < 1. ω∈[0,2π] Dp ∈Dp In both the state-space and frequency domain, constant D, upper bound, a single D scale is selected to guarantee robust performance over all frequencies.
PAGE 65
2.4. µ ANALYSIS 57 The gap between the state-space (or constant D) upper bound and the frequency domain upper bound is more significant. In the state-space upper bound, a single D scale is selected. This gives robust performance for all ∆ satisfying, kvk ≤ kzk for all e ∈ L2 . This can be satisfied for linear time-varying perturbations or non-linear cone bounded perturbations. The formal result is given in the following theorem (given in Packard and Doyle [20]).
PAGE 66
58 2.4.7 CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY Analysis with both Real and Complex Perturbations The above results only apply to the case where ∆ is considered as a constant complex valued matrix at each frequency. In many engineering applications restricting certain of the ∆ blocks to be real valued may result in a less conservative model. Analysis with such restrictions is referred to as the “mixed” µ problem.
PAGE 67
2.5. µ SYNTHESIS AND D-K ITERATION z e y - P (s) - K (s) 59 v w u Figure 2.11: The generic interconnection structure for synthesis kFu (Fl (P (s), K(s)), ∆)k∞ ≤ 1. This is equivalent to K(s) satisfying µ[Fl (P (s), K(s))] < 1. In other words, the closed loop system satisfies the robust performance specification. Unfortunately this problem has not yet been solved, except is a few special cases.
PAGE 68
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 60 Recall that this is an upper bound for the µ problem of interest, implying that, µ[Fl (P (s), K(s)] ≤ 1, as required. However the upper bound may be conservative, meaning that in order to guarantee that µ[Fl (P (s), K(s)] ≤ 1, we have had to back off on the performance and/or the stability margins. With the appropriate choice of D scalings the upper bound will be much closer to µ.
PAGE 69
2.5. µ SYNTHESIS AND D-K ITERATION 61 dynamic system. This requires fitting an approximation to the upper bound D-scale in the iteration. We will now look at this issue more closely. The D-K iteration procedure is illustrated schematically in Figure 2.12. It can be summarized as follows: i) Initialize procedure with K0 (s): H∞ (or other) controller for P (s). ii) Calculate resulting closed loop: Fl (P (s), K(s)). iii) Calculate D scales for µ upper bound: inf D(ω)∈D σmax [D(ω)Fl (P (s), K(s))D(ω)−1 ].
PAGE 70
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 62 z e a) P (s) y v w u -K 0 - z b) G(s) e v w - c) d) e e D(!) z e D^ (s) z e y G(|!) v D(!), w 1 v D^ , (s) w P (s) 1 - w w u K Figure 2.12: D-K iteration procedure: a) Design H∞ (or other) controller: K0 (s) [step i)]. b) Closed loop perturbed system for µ analysis [step ii)]. c) Frequency by frequency upper bound D(ω) scale approximation to µ analysis [step iii)].
PAGE 71
2.5. µ SYNTHESIS AND D-K ITERATION 63 frequency we would have,    D=   d1 I1 ..   ,  . dm Im Ie where the identity Ie is of dimensions dim(e) × dim(e). The calculation of a new H∞ controller requires a state-space realization of D(ω). For each di in D(ω) we must fit a transfer function approximation, which we will denote by dˆi (s). This is denoted by D̂(s) in the above discussion.
PAGE 72
64 CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY Several aspects of this procedure are worth noting. For the µ analysis and D scale calculation, a frequency grid must be chosen. The range and resolution of this grid is a matter of engineering judgement. The µ analysis can require a fine grid in the vicinity of the lightly damped modes. The order of the initial controller, K0 (s), is the same as the interconnection structure, G(s).
PAGE 73
2.6. MODEL REDUCTION 2.6.1 65 Truncation and Residualization The simplest form of model reduction is state truncation. Consider a system, P (s), with a partitioned state matrix,  A11 P (s) =  A21 C1 A12 A22 C2  B1 B2  . D Truncating the states associated with A22 results in, Ptrun (s) = A11 B1 . C1 D In any practical application we would order the states so that those truncated do not significantly affect the system response.
PAGE 74
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 66 The controllability grammian, Y is defined as, Z ∞ Y = T eAt BB T eA t dt, 0 and the observability grammian, X, is defined as Z ∞ X= T eA t C T CeAt dt. 0 The grammians, X and Y , satisfy the Lyapunov equations, AY + Y AT + BB T = 0 AT X + XA + C T C = 0, and this is typically how they are calculated. We can also see from the definitions that X ≥ 0 and Y ≥ 0.
PAGE 75
2.6. MODEL REDUCTION 67 We will now look at a particular choice of transformation. For a minimal realization, we can always find a transformation that gives, Ŷ = T Y T T = Σ, and X̂ = T −T XT −1 = Σ, where Σ = diag(σ1 , . . . , σn ) and σi ≥ 0, i = 1, . . . , n. This realization, where the grammians are equal, is called a balanced realization. Each mode of the balanced system can be thought of as equally controllable and observable. Balanced realization was first introduced by Moore [77].
PAGE 76
CHAPTER 2. OVERVIEW OF THE UNDERLYING THEORY 68 and Glover [79] independently obtained the following bound on the error induced by balanced truncation. Theorem 13 Given a stable, rational, P (s), and Pbal (s), the balanced truncation of order k < n. Then, kP (s) − Pbal (s)k∞ ≤ 2 n X σi i=k+1 and kP (s) − Pbal (s)kH ≤ 2 n X σi .
PAGE 77
2.6. MODEL REDUCTION 69 Consider the problem of finding the stable, order k realization which minimizes the Hankel norm of the error. Define, Phankel (s) as the minimizing system. Then we have, σk+1 ≤ kP (s) − Phankel (s)kH = inf Pk (s) stable kP (s) − Pk (s)kH . This system also satisfies ∞-norm bounds on the error, as illustrated in the following theorem. Theorem 15 Given a stable, rational, P (s), and the optimal kth order Hankel norm approximation, Phankel (s).
PAGE 78
Chapter 3 Functional Description of Xµ 3.1 Introduction This chapter describes the Xµ functions in the context of their intended usage. Chapter 2 provides the reader with an idea of the theoretical basis behind the various analysis and design calculations. Here we outline the software functions available for doing those calculations. Robust control design uses a subset of data objects provided within Xmath. We discuss the details of the most heavily used objects: Dynamic Systems and pdms.
PAGE 79
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 72 3.2.1 Dynamic Systems Xmath has a dynamic system data object which specifies a dynamic system in terms of A, B, C and D matrices. The dynamic equations of the system are, ẋ(t) = Ax(t) + Bu(t), y(t) = Cx(t) + Du(t), in the continuous case, and x(kT + T ) = y(kT ) = Ax(kT ) + Bu(kT ), Cx(kT ) + Du(kT ), in the discrete time case. The discrete time sample period, T , is stored as part of the data object.
PAGE 80
3.2. DATA OBJECTS 73 As above, these polynomials can be specified by their roots or their coefficients. Note that we can specify the variable, and for continuous systems we use “s”. To create a discrete system “z” is used. # Generate the system from the numerator and denominator # coefficients. numerator = makepoly([-.1,19.97,-7.02725],"s"); denominator = makepoly([1,0.3,0.2725],"s"); # We can also do this by specifying the roots of each # polynomial numerator = -0.1*polynomial([199.3475;0.
PAGE 81
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 74 Because the dynamic system is a built-in data object, the label information will be displayed by simply typing the object name (sys in the above) or appending a question mark to the statement. The core function commentof is used to read the comment attached to a specified variable. 3.2.2 pdms A pdm can be thought of as a three-dimensional matrix.
PAGE 82
3.2. DATA OBJECTS 75 Appending and Merging Data Time functions, for creating simulation inputs for example, can be created by combining pdms. Xmath has two such core functions: concatseg and insertSeg The concatseg appends the data of one pdm to another. The domain is recalculated such that it is always regular. The user can specify a new beginning value for the domain and the spacing is determined by a series of prioritized rules. insertSeg inserts the data from one pdm into another.
PAGE 83
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 76 Data can also be extracted by independent variable index number, by providing a scalar argument to the pdm. In the following example the fifth through tenth and the twentieth independent variables are specified for the smaller pdm, pdm2.
PAGE 84
3.2. DATA OBJECTS 77 # index of the pdm the value 100. idxlst = indexlist([4,1,2]) pdm1(idxlst) = 100 Operations on the Independent Variables/Domain The domain of a pdm is readily changed via the pdm command. The following example illustrates a common application; changing a frequency domain in Hertz to radians/second.
PAGE 85
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 78 # Select columns 1, 3, 4 & 5 and rows 2 & 7 from the # pdm: bigpdm. subpdm = bigpdm([1,3:5],[2,7]) This referencing format can also be used to assign data to parts of a larger pdm. This is shown in the following example.
PAGE 86
3.2. DATA OBJECTS 79 and B. Augmentation for pdms simply performs the augmentation at each domain variable. The domains must be the same. Diagonal augmentation can be performed with the Xµ function daug. This is the equivalent of the matrix augmentation: [A, 0; 0, B], except that up to 20 arguments can be augmented in one function call. Algebraic Operations In binary operations (e.g. +, -, *) between a Dynamic System and a matrix, the matrix is assumed to represent the D matrix of a system.
PAGE 87
80 CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ poles, or a zero D term. Generating random systems is useful for simulating systems with unknown subsystems. Specialized Xµ functions are provided for useful manipulations of the state. For example transforming the A matrix to a real valued, 2 × 2 block diagonal form; here referred to as modal form. These state-space functions are tabulated below.
PAGE 88
3.3. MATRIX INFORMATION, DISPLAY AND PLOTTING 81 # N is the decimation ratio. smallpdm = bigpdm([1:N:length(bigpdm)]) 3.2.5 Continuous to Discrete Transformations Xmath has a single function, discretize, for calculating continuous to discrete transformations. Several options are offered including forward and backward difference, Z-transform equivalent, bilinear and pole-zero matching. Xmath also has a makeContinuous function which performs the inverse of the discretize function.
PAGE 89
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 82 3.3.2 Formatted Display Functions It is often useful to consider complex numbers in terms of frequency and damping. This is particularly applicable when studying the poles or zeros of a system. The command rifd provides a real-imaginary-frequency-damping formatted display for complex valued input. The following example illustrates the most common usage.
PAGE 90
3.3. MATRIX INFORMATION, DISPLAY AND PLOTTING 83 g1 = ctrlplot(sys1g,{bode}); g1 = ctrlplot(sys2g,g1,{bode}); g1 = plot({keep=g1,title = "Bode plots",... legend = ["sys1","sys2"]})? Bode plots 10 1 0.1 Magnitude sys1 0.01 sys2 0.001 0.0001 1e-05 0.01 0.1 1 10 1 10 Frequency 0 -50 Phase (degrees) -100 -150 -200 -250 -300 0.01 0.
PAGE 91
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 84 # Nyquist plots = = = = ctrlplot(sys1g,{nyquist}); ctrlplot(sys2g,g2,{nyquist}); ctrlplot(-1,g2,{nyquist,marker=1,line=0}); plot(g2,{projection="orthographic",... legend=["sys1","sys3","critical point"],... title="Nyquist plots"})? Nyquist plots 0.5 0 -0.5 Imaginary g2 g2 g2 g2 sys1 sys3 critical point -1 -1.
PAGE 92
3.4. SYSTEM RESPONSE FUNCTIONS 3.4 System Response Functions 3.4.1 Creating Time Domain Signals 85 The Signal Analysis Module contains several functions which are useful for building time domain signals: gcos and gsin. Xµ provides gstep for the creation of stair-step signals. The example below illustrates the generation of a sine wave and a step function. # A unit step over 10 seconds time = 0:10:.
PAGE 93
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 86 spacing in the input signal pdm. This means that the input pdm must be regularly spaced. Stair-case input functions with relatively few data points will often give erroneous results. The input signal should be interpolated (with interp) before being multiplied with the Dynamic System. The function defTimeRange will calculate a default integration step size and form the required time vector. This can be useful in creating an input pdm or as an input to interp.
PAGE 94
3.4. SYSTEM RESPONSE FUNCTIONS 87 [a,b,c,d] = abcd(sys) sys = system(a,b,c,d) y1 = system(sys,{X0=[-1;0]})*u # Now plot the result g1 g1 g1 g1 = = = = ctrlplot(u,{line style=2}); ctrlplot(y0,g1,{line style=1}); ctrlplot(y1,g1,{line style=4}); plot(g1,{!grid})? 2 1 0 -1 -2 -3 0 2 4 6 8 10 The native Xmath time domain simulation, using the * operation, has been supplemented with an Xµ function: trsp.
PAGE 95
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 88 and intstep are the interpolation order and integration step size. The system must be continuous. The integration order and sample time are prespecified for discrete time systems making the * operator is suitable for such simulations. Xµ provides a sampled data simulation function (sdtrsp) based on the interconnection structure shown in Figure 3.1. v Ms ( ) y w u - Nz ( ) Figure 3.
PAGE 96
3.4. SYSTEM RESPONSE FUNCTIONS 89 maximum. This can be handy when first examining a high order system. An example is given below. # Create a single-input, two-output system. sys = [1/makepoly([1,0.1,1],"s");... makepoly([1,1],"s")/makepoly([1,10],"s")] # The automatic point selection is used. sysg = freq(sys,{Fmin=0.
PAGE 97
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 90 100 10 1 Magnitude 0.1 0.01 0.001 0.0001 1e-05 1e-06 0.01 0.1 1 10 100 10 100 Frequency 100 50 Phase (degrees) 0 -50 -100 -150 -200 0.01 0.1 1 Frequency Note that the default frequency units in Xmath are Hertz. This applies generically to all functions where the user specifies frequency information, although many functions allow radians/second to be specified via a keyword.
PAGE 98
3.5. SYSTEM INTERCONNECTION 91 bigsys sys1 dim2 ````` ```` dim1 sys2 Figure 3.2: Generic Redheffer interconnection structure 3.5 System Interconnection Interconnections of systems are used extensively in the design, analysis and simulation calculations. The most general form of interconnection, and the one used in forming closed loop systems, is the Redheffer (or star) product. The Xµ function starp performs this operation. This Redheffer product operation is illustrated in Figure 3.2.
PAGE 99
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 92 6act dist ? wght t +?l y - +?l - +?l noise(1) noise(2) 6err 7.5 6 p t c t +l 6, ref - 1.6 Figure 3.3: Example interconnection of subsystems y act err clp ref dist noise(1) noise(2) Figure 3.
PAGE 100
3.5. SYSTEM INTERCONNECTION 93 Using sysic to form this interconnection can be considered as four distinct operations. • Specify the individual subsystems. • Name and dimension the input signals. • Specify, algebraically in terms of subsystem outputs or input signals, the output of the interconnected system. • Specify, algebraically in terms of subsystem outputs or input signals, the input to each of the subsystems. • Name the closed loop system.
PAGE 101
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 94 Specify the outputs in terms of the subsystem names and the input names. Note that individual outputs of MIMO systems and simple arithmetic combinations can be specified. Parenthesis specify which output of a MIMO system is to be used. outputs = ["p(1) + wght"; "7.5*c"; "ref - 1.6*p(1) - 1.6*noise(2) - 1.6*wght"] Now specify the inputs to each of the named subsystem. The order is important. Note how the input to a multiple-input system (”c”) is specified.
PAGE 102
3.6. H2 AND H∞ ANALYSIS AND SYNTHESIS -1.0000e+00 -5.0000e-01 -5.0000e-01 -1.0000e+01 0.0000e+00 6.3048e+00 -6.3048e+00 0.0000e+00 95 (rad/sec) ratio 1.0000e+00 6.3246e+00 6.3246e+00 1.0000e+01 1.0000 0.0791 0.0791 1.0000 The order of the names in the systemname variable, must match the order of the rows in the connections variable and the order of the last arguments in the sysic function call. H2 and H∞ Analysis and Synthesis 3.6 This section discusses the synthesis functions available in Xµ.
PAGE 103
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 96 G(s) e P (s) y -Ks w u ( ) Figure 3.5: Interconnection structure for controller synthesis structure contains more than just the open loop plant. It typically also contains frequency dependent weighting functions and specifies the structure of the interconnection between the open loop plant and the controller. The dimensions of y and u are specified by nmeas and ncon.
PAGE 104
3.6. H2 AND H∞ ANALYSIS AND SYNTHESIS 97 equation solution procedure often becomes poorly conditioned. Displaying intermediate calculation results allows the user to fine tune several tolerances if necessary. The intermediate Hamiltonian and Riccati solution details are displayed as the bisection proceeds. The bisection stopping tolerance, Riccati solution tolerances and the Riccati solution method are specified via keywords. Details on the meaning of these tolerances are given in Section 2.3.5.
PAGE 105
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 98 Open loop plant 100 Magnitude 10 1 0.1 0.01 0.001 0.01 0.1 1 10 100 1 10 100 Frequency 0 Phase (degrees) -100 -200 -300 -400 0.001 0.01 0.1 Frequency The desired closed loop configuration is illustrated in Figure 3.6. t plant u k l + 6, Figure 3.
PAGE 106
3.6. H2 AND H∞ ANALYSIS AND SYNTHESIS 99 p Wperf Wact e(1) e(2) y ref t ?l + , plant t u Figure 3.7: Weighted design interconnection structure: p In order to set up the design problem, we consider ref as an unknown input and the tracking error (input to k), and the actuator signal, u, as outputs to be minimized. These outputs are weighted with the weights Wperf and Wact respectively. The weighted interconnection structure for design, p, is illustrated in Figure 3.7.
PAGE 107
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 100 Weighting functions 100 10 Magnitude 1 0.1 Wperf 0.01 Wact 0.001 0.001 0.01 0.1 1 10 100 Frequency # Form the weighted interconnection structure sysnames sysinp = sysout = syscnx = = ["plant";"Wperf";"Wact"] ["ref";"u"] ["Wperf"; "Wact"; "ref-plant"] ["u";... # input to plant "ref-plant";... # input to Wperf "u"] # input to Wact p = sysic(sysnames,sysinp,sysout,syscnx,plant,...
PAGE 108
3.6. H2 AND H∞ ANALYSIS AND SYNTHESIS 12.500 6.250 3.125 1.562 0.781 1.172 5.2e-01 5.2e-01 5.1e-01 5.0e-01 3.9e-01 4.8e-01 1.7e-03 1.7e-03 1.7e-03 1.7e-03 -9.8e+01 1.8e-03 Gamma value achieved: rifd(Kinf) Poles: 1.0e-02 1.0e-02 1.0e-02 1.0e-02 1.0e-02 1.0e-02 101 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.1719 real imaginary frequency (rad/sec) damping ratio -1.0000e-02 -1.0347e+00 -2.6846e+00 -2.6846e+00 -1.2692e+01 0.0000e+00 0.0000e+00 2.
PAGE 109
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 102 (rad/sec) ratio -2.3161e-01 2.3161e-01 3.4754e+00 -3.4754e+00 0.0000e+00 2.7087e-01 2.7087e-01 3.8203e+00 3.8203e+00 5.2060e+00 0.5186 0.5186 0.4152 0.4152 1.0000 real imaginary frequency (rad/sec) damping ratio -5.0000e-02 -5.0000e-02 -5.0000e+00 -2.0000e+01 -3.1225e-01 3.1225e-01 0.0000e+00 0.0000e+00 3.1623e-01 3.1623e-01 5.0000e+00 2.0000e+01 0.1581 0.1581 1.0000 1.0000 -1.4046e-01 -1.4046e-01 -1.5863e+00 -1.5863e+00 -5.
PAGE 110
3.6. H2 AND H∞ ANALYSIS AND SYNTHESIS 103 Controllers 10 1 Magnitude 0.1 0.01 Kinf K2 0.001 0.0001 0.001 0.01 0.1 1 10 Frequency # Examine sensitivity functions sensinf = inv(1 + plant*Kinf) sensinfg = freq(sensinf,omega) sens2 = inv(1 + plant*K2) sens2g = freq(sens2,omega) g2 = ctrlplot(sensinfg,{logmagplot}); g2 = ctrlplot(sens2g,g2,{logmagplot,line style=4}); g2 = plot(g2,{title="Sensitivity functions",...
PAGE 111
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 104 Sensitivity functions 10 Magnitude 1 Kinf K2 0.1 0.01 0.001 0.01 0.1 1 10 100 Frequency Note that with the interconnection structure shown in Figure 3.7, the closed loop transfer function from ref to e(1) is simply Wperf*sensinf. If the H∞ of the resulting closed loop system was less than one, this would guarantee that sensinf was less than inv(Wperf) at every frequency.
PAGE 112
3.6. H2 AND H∞ ANALYSIS AND SYNTHESIS 105 g3 = ctrlplot(step,g3,{line style=2}); g3 = plot(g3,{title="Step responses",... legend=["Kinf";"K2";"input"]})? Step responses 1.2 1 0.8 0.6 Kinf K2 input 0.4 0.2 0 0 3.6.2 2 4 6 8 10 System Norm Calculations Functions are provided for calculating the H2 and H∞ norms of Dynamic Systems. In the H2 case, this involves the solution of a Lyapunov equation.
PAGE 113
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 106 Bounds on the H∞ norm are returned as out. An estimate of the frequency where the norm is achieved is returned as omega. Further control of the iteration is available via keywords. The following example calculates the H2 and H∞ norms of each of the closed loop systems arising from the previous example. Notice that G2 has the minimum H2 norm and Ginf has the minimum H∞ norm.
PAGE 114
3.7. STRUCTURED SINGULAR VALUE (µ) ANALYSIS AND SYNTHESIS 3.7 107 Structured Singular Value (µ) Analysis and Synthesis This section covers the functions used in the D-K iteration procedure. The primary functions are the calculation of the controller (already discussed), the calculation of µ and the fitting of rational D scales. Several subroutines used for D scale fitting are useful in their own regard and are also discussed here.
PAGE 115
108 CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ The outputs of the µ function are: the upper and lower bounds for µ; the D matrix for the upper bound; the Q matrix for the lower bounds; and a sensitivity estimate for the part of the D matrix corresponding to each block in ∆. The sensitivity estimate is essentially the gradient of the upper bound value with respect to the value of the D scale. It is useful in weighting the D scale fitting procedure. The function syntax is shown below.
PAGE 116
3.7. STRUCTURED SINGULAR VALUE (µ) ANALYSIS AND SYNTHESIS c = 1/sqrt(gamma); d = -sqrt(beta/gamma) f = (1+jay)*sqrt(1/(gamma*beta)) psi1 = -pi/2; psi2 = pi U = V = scl M = [a,0; b,b; c,jay*c; d,f] [0,a; b,-b; c,-jay*c; f*exp(jay*psi1), d*exp(jay*psi2)] = diagonal(random(4,1)+0.1*ones(4,1)) scl*U*V*’*inv(scl) Consider four 1×1 blocks. In this example mu is approximately 0.87. blk1 = [1,1; 1,1; 1,1; 1,1] [mubnds1,D1,Dinv1,Delta1] = mu(M,blk1) max(svd(M))? # a very crude upper bound ans (a scalar) = 3.
PAGE 117
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 110 mubnds2? mubnds2 (a column vector) = 3.17155 3.17155 det(eye(4,4) - M*Delta2)? ans (a scalar) = -2.62055e-16 + 5.82345e-17 j 3.7.2 The D-K Iteration Recall from Section 2.5 that the D-K iteration is used as an approximation to µ synthesis. This section discusses how Xµ implements this procedure. The D-K iteration procedure is as follows. The weighted design interconnection structure is referred to as P. The successive controllers are K i, i = 1,. . .
PAGE 118
3.7. STRUCTURED SINGULAR VALUE (µ) ANALYSIS AND SYNTHESIS 111 (a) Controller and closed loop are satisfactory so stop the iteration. (b) The iteration has converged and the controller and closed loop are not satisfactory. In this case the weighted design problem must be reformulated. (c) The iteration has not yet converged. Continue with step 6. 6. Fit rational approximations to D i and Dinv i. The function to do this, musynfit, is described in more detail in Section 3.7.3.
PAGE 119
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 112 PD Dsys i P Dinvsys i - K j, (j=i+1) Figure 3.9: H∞ controller design. Step 9 in the enumerated procedure There is actually another possibility at step 5; numerical problems cause the iteration to diverge. As γ approaches its optimal value, the numerical properties of the calculation deteriorate. This may lead to mu(G i) increasing as i is increased. This problem is observed more often in systems with very lightly damped modes.
PAGE 120
3.7. STRUCTURED SINGULAR VALUE (µ) ANALYSIS AND SYNTHESIS 113 Both the D and D−1 systems (Dsys and Dinvsys) are returned. The D scale (Dmagdata) comes from a µ calculation on a closed loop system. However, Dsys and Dinvsys are required to multiply the open loop system. They must therefore contain the identity matrices for the inputs and outputs which correspond to the measurements and controls.
PAGE 121
114 CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ creating weights from data and simple system identification. Xµ provides two user callable functions for fitting SISO transfer functions to data. The first is mkphase calculates the phase corresponding to a minimum-phase stable system from magnitude data. This uses the complex cepstrum method described by Oppenheim and Schafer [80, p. 501] to generate the desired frequency response. The syntax of the function is given below.
PAGE 122
3.7. STRUCTURED SINGULAR VALUE (µ) ANALYSIS AND SYNTHESIS plant = 1/makepoly([1,0,-0.01],"s") Wm = makepoly([1,20],"s")/makepoly([1,200],"s") omega = logspace(0.01,100,25) plantg = freq(plant,omega) Wmg = freq(Wm,omega) g1 = ctrlplot(plantg,{logmagplot}); g1 = ctrlplot(Wmg,g1,{logmagplot,line style=4}); g1 = plot(g1,{!grid,legend=["Open loop plant";... "multiplicative perturbation weight"]})? 100 10 1 Magnitude 0.1 0.01 0.001 0.0001 Open loop plant multiplicative perturbation weight 1e-05 1e-06 0.
PAGE 123
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 116 Wperf = makepoly([0.01,1],"s")/makepoly([1,0.01],"s") Wact = 0.1* makepoly([1,1],"s")/makepoly([0.05,1],"s") Wnoise = 0.01 Wref = makepoly([0.005,1],"s")/makepoly([0.
PAGE 124
3.7. STRUCTURED SINGULAR VALUE (µ) ANALYSIS AND SYNTHESIS nms = ["plant";"Wm";"Wperf";"Wact";"Wnoise";"Wref"] inp = ["delt";"ref";"noise";"u"] outp = ["100*Wm"; "Wperf"; "Wact"; "Wref-0.01*delt-plant-Wnoise"] cnx = ["u"; "plant"; "Wref-0.01*delt-plant"; "u"; "noise"; "ref"] p = sysic(nms,inp,outp,cnx,plant,Wm,Wperf,Wact,Wnoise,Wref) The interconnection structure, p, is illustrated in Figure 3.11. p 100*Wm Wperf Wact @@,, ,@ ? +l , +l t 6 t ? +l , 0.
PAGE 125
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 118 blk = [1,1; 2,2] [rpbnds1,D1,Dinv1,Delta1,sens1] = mu(Gg,blk) g3 g3 g3 g3 = = = = ctrlplot(np,{log}); # plot on a log-linear scale ctrlplot(rs,g3,{log,line style=3}); ctrlplot(rpbnds1,g3,{log,line style=[4,5]}); plot(g3,{!grid,title="mu analysis",y lab="Magnitude",... legend=["nominal perf.";"robust stab.";... "robust perf.(upper)";"robust perf. (lower)"]})? mu analysis 4 nominal perf. robust stab. 3 robust perf.(upper) Magnitude robust perf.
PAGE 126
3.7. STRUCTURED SINGULAR VALUE (µ) ANALYSIS AND SYNTHESIS D scale fit, block: 1 100 Magnitude data Previous fit, order: 2 New fit, order: 4 Magnitude 10 1 0.1 0.01 0.01 0.1 1 10 100 Frequency (Hz) # Apply the D scales to another H infinity design Kmu = hinfsyn(Ds*p*Dinvs,nmeas,ncntrls,[0;10]) # Close the loop around the weighted interconnection # structure. Gmu = starp(p,Kmu) # weighted closed loop (2nd it.) omega = logspace(0.
PAGE 127
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 120 mu analysis of robust performance 4 3.5 Kmu Kinf 3 Magnitude 2.5 2 1.5 1 0.5 0 0.01 3.7.4 0.1 1 10 100 Constructing Rational Perturbations For simulation purposes it is useful to be able to construct a rational approximation to the ∆ returned by the µ calculation.
PAGE 128
3.8. MODEL REDUCTION 121 pert = randpert(blk,{sys,sfreq,complex,pnorm}) The user can specify whether the perturbation is a dynamic system or matrix, and whether it is real or complex valued, in addition to specifying the norm. 3.7.5 Block Structured Norm Calculations It is possible to get an idea of the input/output combinations which are limiting the robust stability or robust performance of a system by examining the elements of the product, DM D−1 at the critical frequency.
PAGE 129
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 122 A greater range of model reduction functions is available in the Model Reduction Module. Some of the functions described here are cross-licensed with that module. Section 2.6 describes the theory behind these functions. 3.8.1 Truncation and Residualization Truncation is provided by the truncate function (cross-licensed with the Model Reduction Module). Residualization is performed with the Xµ function: sresidualize.
PAGE 130
3.8. MODEL REDUCTION g1 g1 g1 g1 = = = = 123 ctrlplot(sysout2g,g1,{logmagplot,line style=4}); ctrlplot(residerror,g1,{logmagplot,line style=5}); ctrlplot(truncerror,g1,{logmagplot,line style=6}); plot(g1,{!grid,legend=["original system";... "residualized system";"truncated system";... "residualization error";"truncation error"]})? 0.1 Magnitude 0.01 0.001 original system residualized system truncated system residualization error 0.0001 1e-05 0.01 truncation error 0.1 1 10 100 Frequency 3.8.
PAGE 131
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 124 # Displaying the Hankel singular values shows which # states are close to unobservable and uncontrollable. hsv? hsv (a column vector) = 0.0741834 0.0726887 0.0264105 0.000146401 2.7699e-07 # Compare to the errors from the previous example.
PAGE 132
3.8. MODEL REDUCTION 3.8.3 125 Hankel Singular Value Approximation The function ophank (also cross-licensed from the Model Reduction Module) is used to perform optimal Hankel norm approximation. Recall from Section 2.6.3 that there is an astable system achieving the lower bound. The unstable part of this system is returned as the second argument. The following applies the Hankel norm approximation technique to the previous example.
PAGE 133
CHAPTER 3. FUNCTIONAL DESCRIPTION OF Xµ 126 0.1 0.01 Magnitude 0.001 0.0001 original system Hankel norm approximation error: ophank 1e-05 error: balmoore error: sresidualize 1e-06 0.01 0.
PAGE 134
Chapter 4 Demonstration Examples 4.1 The Himat Example The following demo can be run by executing the following Xmath command: execute file = "$XMATH/demos/xMu/himatdemo" 4.1.1 Problem Description The Himat is a small scale remotely piloted aircraft built to investigate high maneuverability fighter aircraft design. The vehicle was flight tested in the late 1970s. The example studied here considers control of only the longitudinal dynamics.
PAGE 135
CHAPTER 4. DEMONSTRATION EXAMPLES 128 δv Perturbations along the velocity vector. α Angle of attack. I.e. angle between the velocity vector and the aircraft’s longitudinal axis. q Rate-of-change of aircraft attitude angle. θ Aircraft attitude angle. Control can be exerted via the elevon and canard, denoted by δe and δc respectively. The angle of attack (α) and attitude angle (θ) are available as direct measurements. A weighted output disturbance rejection problem will be considered.
PAGE 136
4.1. THE HIMAT EXAMPLE w1 w2 e1 e2 129 clp pertin dist Figure 4.2: Interconnection structure for the himat design example The inputs are the elevon position and the canard position. The outputs that are to be kept small are angle-of-attack (α) and pitch angle (θ). The commands required to enter the state-space description are simply matrix assignments for each of A, B, C and D. Note that the states, inputs and outputs are named. a = [-0.0226, -36.6000, -18.9000, -32.1000;...
PAGE 137
CHAPTER 4. DEMONSTRATION EXAMPLES 130 himat (a state space system) = A -0.0226 0 0.0123 0 -36.6 -1.9 -11.7 0 -18.9 0.983 -2.63 1 -32.1 0 0 0 B 0 -0.414 -77.8 0 C 0 0 57.3 0 D 0 0 0 0 0 0 22.4 0 0 0 0 57.
PAGE 138
4.1. THE HIMAT EXAMPLE 131 angle-of-attack pitch angle System is continuous 4.1.3 Creating a Weighted Interconnection Structure for Design The multiplicative input perturbation weight, Wdel , is constructed as a transfer function. The weight is, Wdel = 50(s + 100) . (s + 10000) The output error weight, Wp , is created as, Wp = 0.5(s + 3) . (s + 0.03) This will be used as the performance weight. The appropriate commands are: wdel = makepoly([50,5000])/makepoly([1,10000]) wp = makepoly([0.5,1.
PAGE 139
CHAPTER 4. DEMONSTRATION EXAMPLES 132 Weights for HIMAT 100 Perturbation weight Performance weight Magnitude 10 1 0.1 0.001 0.01 0.1 1 10 100 1000 10000 Frequency The perturbation weight, Wdel , should actually be two-input, two-output. This is also true of the performance weight, Wp . In this example, we are weighting each performance channel identically.
PAGE 140
4.1. THE HIMAT EXAMPLE 133 sysn = ["himat";"wdel";"wp"] in = ["pert(2)";"dist(2)";"control(2)"] out = ["wdel";"wp";"himat + dist"] inter = ["control + pert"; "control"; "himat + dist"] himat ic = sysic(sysn,in,out,inter,himat,wdel,wp) comment himat ic "Himat design interconnection structure" 4.1.4 H∞ Design The next step is to design an H∞ control law for Himat. The function hinfsyn designs an H∞ control law based on the interconnection structure provided.
PAGE 141
CHAPTER 4. DEMONSTRATION EXAMPLES 134 gamma 6.000 3.400 2.100 1.450 1.775 1.613 1.694 1.653 Hx eig 2.3e-02 2.3e-02 2.3e-02 2.3e-02 2.3e-02 2.3e-02 2.3e-02 2.3e-02 X eig 5.6e-05 5.7e-05 5.9e-05 6.4e-05 6.1e-05 6.2e-05 6.1e-05 6.2e-05 Gamma value achieved: Hy eig 2.3e-02 2.3e-02 2.3e-02 2.3e-02 2.3e-02 2.3e-02 2.3e-02 2.3e-02 Y eig 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 nrho xy 0.0626 0.2020 0.5798 1.4678 0.8652 1.1028 0.9725 1.0343 p/f p p p f p f p f 1.
PAGE 142
4.1. THE HIMAT EXAMPLE 135 Zeros: real imaginary frequency (rad/sec) damping ratio -2.2516e-02 -1.7226e+00 -3.0272e+00 -3.1034e+01 -1.0000e+04 -1.0000e+04 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 2.2516e-02 1.7226e+00 3.0272e+00 3.1034e+01 1.0000e+04 1.0000e+04 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 Next, a magnitude plot of the frequency response of k1 is plotted to check that it looks reasonable. om2 = logspace(0.
PAGE 143
CHAPTER 4. DEMONSTRATION EXAMPLES 136 Controller: k1 1 Magnitude 0.1 0.01 0.001 0.1 1 10 100 1000 10000 100 1000 10000 Frequency 200 Phase (degrees) 100 0 -100 -200 -300 0.1 1 10 Frequency Onto the closed loop, first checking that it is stable by looking at the pole positions.
PAGE 144
4.1. THE HIMAT EXAMPLE -2.2517e-02 -2.2600e-02 -3.0000e-02 -3.0000e-02 -2.9369e+00 -2.9974e+00 -4.8310e+00 -6.5876e+00 -5.8350e+01 -5.8350e+01 -8.8792e+01 -8.8792e+01 -9.9778e+01 -7.4258e+03 -1.0000e+04 -1.0000e+04 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 -5.6049e+01 5.6049e+01 -4.2881e+01 4.2881e+01 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 137 2.2517e-02 2.2600e-02 3.0000e-02 3.0000e-02 2.9369e+00 2.9974e+00 4.8310e+00 6.5876e+00 8.0909e+01 8.
PAGE 145
CHAPTER 4. DEMONSTRATION EXAMPLES 138 Singular value plot of the closed loop 2 1.5 1 0.5 0 0.1 4.1.5 1 10 100 1000 10000 µ Analysis of the H∞ Controller The H∞ control law can be analyzed using µ-analysis. The closed-loop system, g1, has 4 inputs and 4 outputs. The first two inputs and outputs correspond to the uncertainty block, and the second two correspond to the disturbance rejection block, or performance block.
PAGE 146
4.1. THE HIMAT EXAMPLE 139 complex valued blocks. The upper and lower bounds of the µ function are returned in bnds. Also returned are the scaling matrices, D and Dinv, corresponding to the upper bound. The smallest destabilizing perturbation at each frequency is returned as Delta. The variable sens is the sensitivity of the upper bound to the D and Dinv scaling matrices.
PAGE 147
CHAPTER 4. DEMONSTRATION EXAMPLES 140 Max. singular value and mu comparison 1.8 1.6 1.4 1.2 1 max. singular value 0.8 mu upper bound mu lower bound 0.6 0.4 0.1 1 10 100 1000 10000 Note that µ(g1g) is not less than one at all frequencies — we have not met the design objectives. A D-K iteration will be used to lower µ and improve the robust performance with respect to these objectives. 4.1.
PAGE 148
4.1. THE HIMAT EXAMPLE 141 The second thing to note is that the interconnection structure has the additional control inputs and measurement outputs. The D and Dinv systems must be augmented with identities corresponding to these additional inputs and outputs. The musynfit function performs both of these operations. The syntax of musynfit is: [Dsys,Dinvsys] = musynfit(D,blk,nmeas,nctrls,sens,oldDsys) The variables D and sens come directly from the mu function.
PAGE 149
CHAPTER 4. DEMONSTRATION EXAMPLES 142 D scale fit, block: 1 10 Magnitude data Previous fit, order: 0 New fit, order: 3 Magnitude 1 0.1 0.01 0.1 1 10 100 1000 10000 Frequency (Hz) Upper bound comparison, block: 1 10 Data based bound Previous fit bound, order: 0 Magnitude New fit bound, order: 3 1 0.1 0.1 1 10 100 1000 10000 1000 10000 Frequency (Hz) D scale fit weight for block: 1 10 Magnitude 1 0.1 0.01 0.
PAGE 150
4.1. THE HIMAT EXAMPLE 143 comment D1sys "system approx. to D1" comment D1invsys "system approx. to D1inv" 4.1.7 Design Iteration #2 The new D scales can be pre and post multiplied onto the orginal interconnection structure. himat ic2 = D1sys * himat ic * D1invsys [,,nx] = size(himat ic2) display "himat ic2 now has " + string(nx) + " states" himat ic2 now has 20 states comment himat ic2 "interconnection for iteration 2" Note the increase in states due to the inclusion of the D scales.
PAGE 151
CHAPTER 4. DEMONSTRATION EXAMPLES 144 g2 = starp(himat ic,k2) comment g2 "closed loop: iteration 2" [,,nx] = size(k2) display "k2 now also has " + string(nx) + " states" k2 now also has 20 states This design probably resulted in a control law which achieved an infinity norm of approximately 1.1 for the new interconnection structure. k2 has been designed to reduce µ(g2). The singular values may actually get worse. We will see that this is the case here.
PAGE 152
4.1. THE HIMAT EXAMPLE 145 Robustness analysis for the g2 system 1.8 1.6 g2: mu upper bound g2: mu lower bound g1: mu upper bound 1.4 g1: mu lower bound 1.2 1 0.8 0.6 0.4 0.1 comment comment comment comment comment 1 10 100 1000 10000 bnds2 "mu bounds: iteration 2" D2 "D scale: iteration 2" D2inv "D inverse scale: iteration 2" Delta2 "worst case perturbation: iteration 2" sens2 "D scale sensitivity: iteration 2" Here we have only done a single D-K iteration.
PAGE 153
CHAPTER 4. DEMONSTRATION EXAMPLES 146 g3 = starp(himat ic,k3) g3g = freq(g3,omega) [bnds3,D3,D3inv,Delta3,sens3] = mu(g3g,blk) At this point we could do another iteration (to get k4) or perhaps run some simulations to check out k3 more thoroughly. 4.1.8 Simulation Comparison with a Loopshaping Controller A loopshaping design is performed and compared to the H∞ and µ designs. As the loopshaping design procedure is relatively standard, only the final controller is given here (klp). a =[-5.8928e-02,-6.
PAGE 154
4.1. THE HIMAT EXAMPLE 147 d = zeros(2,2) klp = system(a,b,c,d) comment klp "loop shape controller" We will compare the designs, with no error or uncertainty weights, for the nominal case and with a perturbation block of ∆ = [0.1,0;0,-0.1] for the input multiplicative perturbation. The time response will be from 0 to 2 seconds with a sample time of 0.01 seconds. We’ll look at a unit step input into the first channel. The unweighted closed loop system is now formed for each controller.
PAGE 155
CHAPTER 4. DEMONSTRATION EXAMPLES 148 comment gsim mu nom "nominal closed loop sys: mu ctrl" comment gsim hinf nom "nominal closed loop sys: hinf ctrl" comment gsim lp nom "nominal closed loop sys: klp ctrl" A step disturbance is introduced into the first channel. time = 0:2:0.
PAGE 156
4.1. THE HIMAT EXAMPLE 149 Kmu step dist. response (nominal) 1.5 1 0.5 0 -0.5 -1 0 0.5 1 1.5 2 1.5 2 1.5 2 time Kinf step dist. response (nominal) 1.5 1 0.5 0 -0.5 -1 0 0.5 1 time Klp step dist. response: (nominal) 1.5 1 0.5 0 -0.5 -1 0 0.
PAGE 157
CHAPTER 4. DEMONSTRATION EXAMPLES 150 The loopshaping design gives a decoupled response. Both the H∞ and µ designs trade decoupling for speed of response and, as we shall see, robustness with respect to perturbations. The simulation is repeated with a perturbation of size 0.1. Note that this is only 10% of the size perturbation that we were analyzing and designing for in the above. delta = [.1,0;0,-0.
PAGE 158
4.1. THE HIMAT EXAMPLE 151 Kmu step dist. response (perturbed) 1.5 1 0.5 0 -0.5 -1 0 0.5 1 1.5 2 1.5 2 1.5 2 time Kinf step dist. response (perturbed) 1.5 1 0.5 0 -0.5 -1 0 0.5 1 time Klp step dist. response: (perturbed) 1.5 1 0.5 0 -0.5 -1 0 0.
PAGE 159
152 CHAPTER 4. DEMONSTRATION EXAMPLES The loopshaping controller had good nominal performance and very poor robust performance. This was illustrated with a relatively small perturbation. The difference between the µ and H∞ controllers was small in both the nominal and perturbed cases. This may not always be the case for several reasons. Only a single D-K iteration was performed here. Further iterations would further improve the performance of the µ controller.
PAGE 160
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE 4.2 153 A Simple Flexible Structure Example The demonstration script, jplphBdemo.ms, runs through a D-K iteration design for a simple flexible structure problem. The following demo can be run by executing the following Xmath command: execute file = "$XMATH/demos/xMu/jplphBdemo" where $XMATH is the path to your Xmath source location.
PAGE 161
CHAPTER 4. DEMONSTRATION EXAMPLES 154 Mounted on xed optical bench Mounted on exible structure voice coil XXXX XXXXX laser path piezo target mirror - laser system Figure 4.3: Schematic diagram of the JPL Phase B optical testbed design problem problem is the fact that the voice-coil mirror assembly has significant mass and its movement excites a mode in the structure. An identically driven counterbalance effectively makes the piezo-electric actuator reactionless.
PAGE 162
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE 155 static gain. We will include some dynamic uncertainty in the actual design. piezo = 1 4.2.2 Creating the Weighted Design Interconnection Structure The weighted, open-loop design interconnection structure is illustrated in Figure 4.5. For clarity, the two perturbations, ∆1 and ∆2 , have been shown inside the structure. We have included perturbation models for both the actuators, involving the weights Wavoice, Wmvoice and Wmpiezo.
PAGE 163
CHAPTER 4. DEMONSTRATION EXAMPLES 156 The additive weight clearly provides for significant high frequency uncertainty. A multiplicative weight models the low frequency uncertainty. The value selected is somewhat arbitrary and can be considered as a tunable design weight. Wmvoice = 0.1 A multiplicative perturbation is used to model uncertainty in the piezo driver. Experimental observations suggest a 5% deviation from nominal across the frequency range of interest. Wmpiezo = 0.
PAGE 164
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE 157 F = 100 Wactp = makepoly([1/(2*pi*F),1],"s")/... makepoly([1/(200*pi*F),1],"s") Wactp = Wactp*4 The lower frequency response of the voice-coil system means that the Wactv weight should begin rolling up at around 10 Hz. F = 10 Wactv = makepoly([1/(2*pi*F),1],"s")/... makepoly([1/(200*pi*F),1],"s") Wactv = Wactv*0.02 Both of these weights can be adjusted to trade between the relative responses from the voice-coil and the piezo.
PAGE 165
158 CHAPTER 4. DEMONSTRATION EXAMPLES Wdist,Wnoise,Wperf,Wactp,Wactv,piezo) size(P)? ans (a row vector) = 7 6 8 Now select a frequency grid for calculating the frequency responses. Some additional points are included near the oscillatory modes. omega = logspace(1,1000,15)’ omega = sort([omega; [5.275:0.05:5.625]’; [4.5:0.25:6.5]’])’ Examine the frequency response of the open-loop system.
PAGE 166
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE 159 voice coil model 100000 10000 1000 Magnitude 100 10 1 0.1 0.01 0.
PAGE 167
CHAPTER 4. DEMONSTRATION EXAMPLES 160 And examine the design weights. weightsg = freq(weights,omega) gph2 = ctrlplot(weightsg,{logmagplot}); gph2 = plot(gph2,legend=["Wperf";"Wavoice";"Wmvoice";... "Wmpiezo";"Wdist";"Wnoise";"Wactp";"Wactv"],...
PAGE 168
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE 161 Design weights 100 Wperf Wavoice Wmvoice Wmpiezo Wdist Wnoise 10 Wactp Wactv Magnitude 1 0.1 0.01 0.
PAGE 169
CHAPTER 4. DEMONSTRATION EXAMPLES 162 4.2.3 Design of an H∞ Controller An H∞ design is now performed. Recall that we have one interferometer measurement and two controller outputs. nmeas = 1 ncon = 2 glimits = [0;20] [Khinf,gamma] = hinfsyn(P,nmeas,ncon,glimits,{tol=0.25}) Test bounds: 0.0000 < gamma <= 20.0000 gamma Hx eig X eig Hy eig Y eig nrho xy p/f 20.000 2.8e-02 3.1e-07 2.8e-01 -6.5e-17 0.0007 p 10.000 2.8e-02 3.1e-07 2.8e-01 -4.5e-17 0.0026 p 5.000 2.8e-02 3.1e-07 2.8e-01 -1.7e-16 0.0109 p 2.
PAGE 170
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE -3.5543e+02 -3.5543e+02 -6.4598e+02 -6.4598e+02 -3.8358e+02 -3.8358e+02 -1.6934e+03 -2.8107e+03 -6.2832e+03 -6.2832e+04 Zeros: -3.5543e+02 3.5543e+02 -3.3537e+01 3.3537e+01 -5.3777e+02 5.3777e+02 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 5.0265e+02 5.0265e+02 6.4685e+02 6.4685e+02 6.6055e+02 6.6055e+02 1.6934e+03 2.8107e+03 6.2832e+03 6.2832e+04 163 0.7071 0.7071 0.9987 0.9987 0.5807 0.5807 1.0000 1.0000 1.0000 1.0000 We can also look at the controller poles.
PAGE 171
CHAPTER 4. DEMONSTRATION EXAMPLES 164 Controller: Khinf 1 Magnitude 0.1 0.01 0.001 0.
PAGE 172
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE 4.2.4 165 Robustness Analysis The block structure has two perturbation ∆ blocks and a “performance” block. The two voice-coil perturbations are put into a single 1×2 block as they enter the system at the same point. Note that this is not identical to two separate blocks — for example a perturbation in which both blocks have magnitude one is now longer included.
PAGE 173
CHAPTER 4. DEMONSTRATION EXAMPLES 166 gph4 gph4 gph4 gph4 = = = = ctrlplot(npbnds,{log,line style=4}); ctrlplot(rsbnds,gph4,{log,line style=[3,5]}); ctrlplot(rpbnds,gph4,{log,line style=[1,2]}); plot(gph4,{!grid,legend=["Nom perf";"Rob stab (up bnd)";... "Rob stab (lw bnd)"; "Rob perf (up bnd)";...
PAGE 174
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE 167 perturbed closed loop: vc actuator 0.015 Khinf Kmu 0.01 0.005 0 -0.005 -0.01 -0.015 0 0.5 1 1.
PAGE 175
CHAPTER 4. DEMONSTRATION EXAMPLES 168 4.2.5 D-K Iteration We will now perform one D-K iteration to generate the controller Kmu. Significant robustness and performance improvement is achieved with only one iteration. Transfer functions are fit to the D-scales from the previous robust performance µ test. Here we preselect an order of 2 for each D-scale. This has been found to give a satisfactory result.
PAGE 176
4.2.
PAGE 177
CHAPTER 4. DEMONSTRATION EXAMPLES 170 Controller: Kmu 1 Magnitude 0.1 0.01 0.001 0.
PAGE 178
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE 171 We now examine the robustness properties of the new closed loop system. We already know that the robust performance test will be less than the γ value from the H∞ synthesis above (in this case 0.2930). The results are again displayed graphically.
PAGE 179
CHAPTER 4. DEMONSTRATION EXAMPLES 172 mu analysis of Gmu 0.3 0.25 0.2 0.15 Nom perf Rob stab (up bnd) Rob stab (lw bnd) Rob perf (up bnd) 0.1 Rob perf (lw bnd) 0.
PAGE 180
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE 4.2.6 173 A Simulation Study Now the two controllers (Khinf and Kmu) are studied by simulation. An unweighted interconnection is set up with sysic and starp is used to close the loop for each controller. ssnames = ["vcmodel"; "Wavoice"; "Wmvoice"; "Wmpiezo"; "piezo"] inps = ["d1i"; "d2i"; "dist"; "noise"; "vact"; "pact"] ops = ["Wavoice"; "Wmvoice"; "Wmpiezo";... "d1i + vcmodel + d2i + piezo + noise"; "vact"; "pact";...
PAGE 181
CHAPTER 4. DEMONSTRATION EXAMPLES 174 Simulation: disturbance 0.006 0.004 0.002 0 -0.002 -0.004 -0.006 0 0.5 1 1.
PAGE 182
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE gph8 = ctrlplot(u(2,1)); gph8 = plot(gph8,{title="Simulation: noise",...
PAGE 183
CHAPTER 4. DEMONSTRATION EXAMPLES 176 Simulation: noise 0.003 0.002 micrometers 0.001 0 -0.001 -0.002 -0.003 0 0.5 1 1.
PAGE 184
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE 177 A nominal response is calculated by setting ∆ = 0. To get the open-loop simulation model we close the unweighted interconnection structure with a controller equal to zero. deltazero = zeros(2,3) Kzero = zeros(2,1) nomolp = starp(deltazero,Pnom) nomolp = starp(nomolp,Kzero) yolp = nomolp*u gph9 = ctrlplot(yolp(1,1)); gph9 = plot(gph9,{title="open loop beam length",...
PAGE 185
CHAPTER 4. DEMONSTRATION EXAMPLES 178 open loop beam length 10 micrometers 5 0 -5 -10 0 0.5 1 1.
PAGE 186
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE 179 Now we consider the closed-loop nominal response with Khinf and Kmu. The closed-loop system happens to have a large number of high frequency poles which do not contribute significantly to the response. They have the effect of forcing a very fine time discretization in the simulation, resulting in a long calculation time. We remove all poles of frequency greater than 100 rad/sec. by residualization.
PAGE 187
CHAPTER 4. DEMONSTRATION EXAMPLES 180 closed loop beam length 0.4 0.3 Khinf Kmu 0.2 micrometers 0.1 0 -0.1 -0.2 -0.3 0 0.5 1 1.
PAGE 188
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE gph11 = ctrlplot(yclphinf(2,1)); gph11 = ctrlplot(yclpmu(2,1),gph11); gph11 = plot(gph11,{legend=["Khinf";"Kmu"],...
PAGE 189
CHAPTER 4. DEMONSTRATION EXAMPLES 182 closed loop: vc actuator 0.015 Khinf Kmu 0.01 0.005 0 -0.005 -0.01 -0.015 0 0.5 1 1.
PAGE 190
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE gph12 = ctrlplot(yclphinf(3,1)); gph12 = ctrlplot(yclpmu(3,1),gph12); gph12 = plot(gph12,{legend=["Khinf";"Kmu"],...
PAGE 191
CHAPTER 4. DEMONSTRATION EXAMPLES 184 closed loop: piezo actuator 0.03 Khinf Kmu 0.02 0.01 0 -0.01 -0.02 0 0.5 1 1.
PAGE 192
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE 185 Note that Kmu achieves better performance at the expense of greater actuator effort. We will now repeat this simulation for a perturbed system. A bad ∆ is chosen and scaled to have norm 0.5. This is obtained from destabilizing ∆ calculated for the µ lower bound. An all-pass system is fitted to ∆ at one frequency to create a real-rational perturbation. The frequency selected is that where µ is at its maximum.
PAGE 193
186 CHAPTER 4. DEMONSTRATION EXAMPLES gph13 = ctrlplot(ybclpmu(1,1),gph13); gph13 = plot(gph13,{legend=["Khinf";"Kmu"],... title="perturbed closed loop beam length",...
PAGE 194
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE 187 Controller: Khinf 1 Magnitude 0.1 0.01 0.001 0.
PAGE 195
188 CHAPTER 4. DEMONSTRATION EXAMPLES gph14 = ctrlplot(ybclphinf(2,1)); gph14 = ctrlplot(ybclpmu(2,1),gph14); gph14 = plot(gph14,{legend=["Khinf";"Kmu"],...
PAGE 196
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE 189 perturbed closed loop: vc actuator 0.015 Khinf Kmu 0.01 0.005 0 -0.005 -0.01 -0.015 0 0.5 1 1.
PAGE 197
190 CHAPTER 4. DEMONSTRATION EXAMPLES gph15 = ctrlplot(ybclphinf(3,1)); gph15 = ctrlplot(ybclpmu(3,1),gph15); gph15 = plot(gph15,{legend=["Khinf";"Kmu"],...
PAGE 198
4.2. A SIMPLE FLEXIBLE STRUCTURE EXAMPLE 191 perturbed closed loop: piezo actuator 0.03 Khinf Kmu 0.02 0.01 0 -0.01 -0.02 -0.03 0 0.5 1 1.
PAGE 199
Bibliography [1] J. C. Doyle, “Lecture notes on advances in multivariable control.” ONR/Honeywell Workshop, Minneapolis, MN., 1984. [2] J. Doyle, “Structured uncertainty in control system design,” in Proc. IEEE Control Decision Conf., pp. 260–265, 1985. [3] A. K. Packard, What’s new with µ: Structured Uncertainty in Multivariable Control. PhD thesis, University of California, Berkeley, 1988. [4] R. S. Smith, Model Validation for Uncertain Systems. PhD thesis, California Institute of Technology, 1990. [5] J.
PAGE 200
194 BIBLIOGRAPHY [11] G. J. Balas and J. C. Doyle, “Robust control of flexible modes in the controller crossover region,” in Proc. Amer. Control Conf., 1989. [12] G. J. Balas, A. K. Packard, and J. Harduvel, “Application of µ-synthesis techniques to momentum management and attitude control of the Space Station,” in AIAA Guidance, Navigation and Cont. Conf., 1991. [13] R. S. Smith, C.-C. Chu, and J. L.
PAGE 201
BIBLIOGRAPHY 195 [24] D. L. Laughlin, K. G. Jordan, and M. Morari, “Internal model control and process uncertainty: mapping uncertainty regions for SISO controller design,” Int. J. of Control, vol. 44, no. 6, pp. 1675–1698, 1986. [25] R. S. Smith and M. Dahleh, eds., The Modeling of Uncertainty in Control Systems: Proceedings of the 1992 Santa Barbara Workshop. 391 pgs., Springer-Verlag, 1994. [26] M. Gevers, “Connecting identification and robust control: A new challenge,” in Proc. IFAC Symp.
PAGE 202
196 BIBLIOGRAPHY [37] R. Kosut, M. Lau, and S. Boyd, “Parameter set identification of systems with uncertain nonparametric dynamics and disturbances,” in Proc. IEEE Control Decision Conf., vol. 6, pp. 3162–3167, 1990. [38] G. Goodwin, B. Ninness, and M. Salgado, “Quantification of uncertainty in estimation,” in Proc. Amer. Control Conf., pp. 2400–2405, 1990. [39] B. M. Ninness and G. C. Goodwin, “Robust frequency response estimation accounting for noise and undermodeling,” in Proc. Amer. Control Conf.
PAGE 203
BIBLIOGRAPHY 197 [50] J. M. Krause, P. P. Khargonekar, and G. Stein, “Robust parameter adjustment with nonparametric weighted-ball-in-H ∞ uncertainty,” IEEE Trans. Auto. Control, vol. AC-35, pp. 225–229, 1990. [51] R. S. Smith and J. C. Doyle, “Closed loop relay estimation of uncertainty bounds for robust control models,” in Proc. of the 13th IFAC World Congress, vol. 9, pp. 57–60, July 1993. [52] R. J. Schrama and P. M. V.
PAGE 204
198 BIBLIOGRAPHY [63] A. J. Laub, “A Schur method for solving algebraic Riccati equations,” IEEE Trans. Auto. Control, vol. AC-24, pp. 913–921, 1979. [64] T. Pappas, A. J. Laub, and N. R. Sandell, “On the numerical solution of the discrete-time algebraic Riccati equation,” IEEE Trans. Auto. Control, vol. AC-25, pp. 631–641, 1980. [65] W. F. Arnold and A. J. Laub, “Generalized eigenproblem algorithms and software for algebraic Riccati equations,” Proc. IEEE, vol. 72, pp. 1746–1754, 1984. [66] A. J.
PAGE 205
BIBLIOGRAPHY 199 [76] M. Dahleh, A. Tesi, and A. Vicino, “Extremal properties for the parametric robust performance problem,” Tech. Rep. UCSB-ME-91-4, Univ. California, Santa Barbara, Mech. Eng., 1991. also submitted to 30th IEEE CDC. [77] B. C. Moore, “Principal components analysis in linear systems: controllability, observability and model reduction,” IEEE Trans. Auto. Control, vol. AC-26, pp. 17–31, 1981. [78] D. F. Enns, Model Reduction for Control System Design. PhD thesis, Stanford University, 1984.
PAGE 206
Chapter 6 Function Reference 6.1 Xµ Functions The following pages contain descriptions of the Xµ functions. These are also available on-line via the help utility. Each description also gives an illustrative example of the function’s use. The functions are included in alphabetical order. For convenience they are cross-referenced by typical use in the following list. System building and interconnection ...................................................................... 231 ...........................
PAGE 207
CHAPTER 6. FUNCTION REFERENCE 202 rifd ...................................................................... 335 ................................................................. 221 ctrlplot Time response calculations and pdm functions gstep .................................................................... 247 interp ................................................................... 281 mergeseg ................................................................. 285 randpdm ...........
PAGE 208
6.1. Xµ FUNCTIONS 203 Model reduction and state-space functions ................................................................. 205 .............................................................. 293 ................................................................... 315 .............................................................. 319 simtransform ............................................................ 345 sresidualize ............................................................
PAGE 209
CHAPTER 6. FUNCTION REFERENCE 204 Miscellaneous functions conpdm ................................................................... 213 consys ................................................................... 215 ...................................................................... 217 ................................................................ 237 ...................................................................
PAGE 210
205 balmoore balmoore Syntax [SysR,HSV,T] = balmoore(Sys,{nsr,bound}) Parameter List Inputs: Sys Linear, stable, minimal state-space system nsr (optional) If bound is used then a reduction will be performed which meets an error bound specified by the value in nsr, otherwise nsr is the order of the reduced system. If nsr is not specified, the user will be prompted for its value after the Hankel singular values are displayed. Keywords: bound Boolean; meet an upper bound on the error (see nsr).
PAGE 211
Chapter 6. Function Reference 206 Reference B.C. Moore, “Principal Component Analysis in Linear Systems: Controllability, Observability and Model Reduction,” IEEE Trans. Auto. Ctrl., Vol. 26, No. 1, pp. 17–32, Feb. 1981. Example # Create a five state system for reduction. a = daug(-0.891334,[-1.20857,0.799042;-0.799042,-1.20857],... -4.74685,-21.3013) b = [0.0262569;-0.189601;-0.113729;0.211465;-0.538239] c = [0.120725,-0.336942,0.397198,-0.700524,-1.02235] d = 0 sys1 = system(a,b,c,d) fHz = logspace(0.
PAGE 212
balmoore See Also: minimal, ophank.
PAGE 213
209 blknorm blknorm Syntax normM = blknorm(M,blk,p,Frobenius) Parameter List Inputs: M Matrix (or pdm). blk Block structure. See mu section of the manual for a description of the syntax. Scalar valued. Specifies the Holder “p” norm to be used, where 1 ≤ p ≤ inf. Optional. The default is p = 2. p Keywords: Frobenius The Frobenius norm is used. Outputs: normM Matrix (or pdm) norms of each block of M. Description M is partitioned according the input/output partitions determined by blk.
PAGE 214
Chapter 6. Function Reference 210 Examples A = random(3,3)-0.5*ones(3,3)? A (a square matrix) = 0.0618661 0.390622 -0.112622 0.0896177 0.00422128 0.42229 0.185398 -0.150638 0.448818 blknorm(A,[1,1; 1,1; 1,1]) ans (a square matrix) = 0.0618661 0.390622 0.112622 0.0896177 0.00422128 0.42229 0.185398 0.150638 0.448818 blknorm(A,[3,3]) ans (a scalar) = 0.682429 # compare to the following: max(svd(A)) ans (a scalar) = 0.
PAGE 215
211 blknorm 2.23607 7.81025 13.4536 3 7 11 4 8 12 # and compare to norm(B(1,1:2)) ans (a scalar) = 2.
PAGE 216
213 conpdm conpdm Syntax outpdm = conpdm(mat,domain,{skipChks}) Parameter List Inputs: mat constant matrix domain domain over which outsys will be defined. Keywords: skipChks Boolean specifying that syntax checking is to be skipped. Outputs: outpdm pdm Description Creates a pdm data object from a constant matrix. Outpdm represents a constant gain; its value is repeated at every instance of the domain.
PAGE 217
215 consys consys Syntax outsys = consys(mat,{skipChks}) Parameter List Inputs: mat constant matrix Keywords: skipChks Boolean specifying that syntax checking is to be skipped. Outputs: outsys Dynamic System Description Creates a Dynamic System object from a constant matrix. Outsys represents a constant gain; the A, B and C matrices are empty. It is equivalent to, outsys = system([],[],[],mat) and can be useful in specifying constant inputs to freq.
PAGE 218
217 csum csum Syntax [outpdm] = csum(inpdm, {channels}) Parameter List Inputs: inpdm real or complex valued pdm or constant matrix Keywords: channels Sum over channels. outpdm has the same dimensions as inpdm. Outputs: outpdm output pdm Description Perform a cumulative sum over the rows of a matrix or a pdm. If channels is specified then the sum is performed over the domain of the pdm. Examples A = [ones(6,1),random(6,1)]? A (a rectangular matrix) = 1 1 1 1 1 1 0.608453 0.854421 0.0642647 0.
PAGE 219
Chapter 6. Function Reference 218 ans (a rectangular matrix) = 1 2 3 4 5 6 0.608453 1.46287 1.52714 2.35505 3.28128 3.848 pdmA = pdm(A,[1,2,3])? pdmA (a pdm) = domain | Col 1 Col 2 -------+-----------------------1 | Row 1 1 0.608453 | Row 2 1 0.854421 -------+-----------------------2 | Row 1 1 0.0642647 | Row 2 1 0.827908 -------+-----------------------3 | Row 1 1 0.926234 | Row 2 1 0.
PAGE 220
219 csum ans (a pdm) = domain | Col 1 Col 2 -------+----------------------1 | Row 1 1 0.608453 | Row 2 1 0.854421 -------+----------------------2 | Row 1 2 0.672717 | Row 2 2 1.68233 -------+----------------------3 | Row 1 3 1.59895 | Row 2 3 2.
PAGE 221
ctrlplot ctrlplot Syntax graph = ctrlplot(pdm,old graph,{keywords}) Parameter List Inputs: pdm Pdm (or matrix) containing the data to be plotted. old graph (optional) Graphical object to which data is added. Conceptually the same as plot(pdm,{keep=old graph}). Keywords: The following keywords specify the basic plot format. Only one can be selected. timeresp (default) real(pdm) vs. domain.
PAGE 222
Chapter 6. Function Reference 222 linear linear domain. Default = 1 for timeresp keyword. Default = 0 for bode keyword. log logarithmic domain. Default = 1 for bode keyword. Default = 0 for timeresp keyword. Default units can be supplied for the magnitude and phase plots (bode, nichols, logmagplot and phaseplot keywords) with the following keywords. degrees Angles are specified in degrees (default = 1) db log magnitudes are specified in decibels (default = 0).
PAGE 223
ctrlplot 223 Description This function performs some common control system related plotting. The user can use ctrlplot to set up a basic plot and perform some preprocessing of the data. This generates a graphical object containing the data and the user can perform subsequent calls to plot to add things like text, labels, gridding etc. The second argument (optional) is a graphical object to which the pdm will be added. Plots over a domain (bode, timeresp, logmagplot, phaseplot) must be called with a pdm.
PAGE 224
Chapter 6.
PAGE 225
ctrlplot Bode plots 10 1 0.1 Magnitude sys1 0.01 sys2 0.001 0.0001 1e-05 0.01 0.1 1 10 1 10 Frequency 0 -50 Phase (degrees) -100 -150 -200 -250 -300 0.01 0.
PAGE 226
Chapter 6. Function Reference # Nyquist plots g2 = ctrlplot(sys1g,{nyquist}); g2 = ctrlplot(sys2g,g2,{nyquist}); g2 = ctrlplot(-1,g2,{nyquist,marker=1,line=0}); g2 = plot(g2,{projection="orthographic",...
PAGE 227
ctrlplot Nyquist plots 0.5 0 Imaginary -0.5 sys1 sys3 -1 critical point -1.
PAGE 228
Chapter 6. Function Reference 228 # Create a second order lightly damped system to illustrate # time response plotting. The calculation is repeated with # a non-zero initial condition. sys = 5/makepoly([1,1,5],"s") u = gstep([0:0.
PAGE 229
ctrlplot 2 1 0 input -1 x0 = zero non-zero x0 -2 -3 0 2 4 6 8 10
PAGE 230
See Also: plot. Chapter 6.
PAGE 231
daug daug Syntax out = daug (sys1,sys2,...) Parameter List Inputs: sys1 Input systems. These can be dynamical systems and constants, or pdms and constants. .. . Outputs: out “ output system. Description Diagonal augmentation of dynamical system/pdm/constant, matrices.    out =   sys1 0 ... 0 0 sys2 . . . 0 .. .. .. .. . . . . 0 0 . . . sysN      Limitations Only 21 systems can be augmented with a single function invocation.
PAGE 232
Chapter 6. Function Reference 232 ans (a square matrix) = 1 0 0 0 1 0 0 0 0 2 2 0 0 0 0 Inf sys1 = randsys(1,1,2,{stable}) sys1 = system(sys1,{statenames="sys1state"})? sys1 (a state space system) = A -0.886949 B 0.853282 0.012459 C 0.186754 D 0.492058 0.
PAGE 233
daug -1.67106 B 0.579502 C 0.262815 0.436099 D 0.911055 0.808267 X0 1 State Names ----------sys2state System is continuous # Note the effect of the constant in the following daug(sys1,10,sys2) ans (a state space system) = A -0.886949 0 0 -1.67106 B 0.853282 0 0.012459 0 C 0.186754 0 0 0 0 0 0.262815 0.436099 0 0 0 0.
PAGE 234
Chapter 6. Function Reference 234 D 0.492058 0 0 0 0.748961 0 0 0 0 10 0 0 0 0 0.911055 0.808267 X0 0 1 State Names ----------sys1state sys2state System is continuous pdm1 = randpdm(3,2,2)? pdm1 (a pdm) = domain | Col 1 Col 2 -------+--------------------------0 | Row 1 0.810265 0.259043 | Row 2 0.413909 0.359993 -------+--------------------------1 | Row 1 0.691279 0.765686 | Row 2 0.357265 0.76934 -------+--------------------------2 | Row 1 0.547763 0.0962289 | Row 2 0.956117 0.
PAGE 235
daug 0 | 10 -------+----1 | 10 -------+----2 | 10 -------+----daug(pdm1,pdm2) ans (a pdm) = domain | Col 1 Col 2 Col 3 -------+---------------------------------0 | Row 1 0.810265 0.259043 0 | Row 2 0.413909 0.359993 0 | Row 3 0 0 10 -------+---------------------------------1 | Row 1 0.691279 0.765686 0 | Row 2 0.357265 0.76934 0 | Row 3 0 0 10 -------+---------------------------------2 | Row 1 0.547763 0.0962289 0 | Row 2 0.956117 0.
PAGE 236
237 delsubstr delsubstr Syntax [outstr] = delsubstr(str,charstr) Parameter List Inputs: Outputs: str String or vector of strings. charstr String outstr String or vector of strings. Description All occurences of the substring, charstr, within str are deleted. If, by deleting charstr, another occurence of charstr in created, it will not be deleted. Examine the second example closely to see the effect of this. Unless str is a scalar string, deleting a whole string element will cause an error.
PAGE 237
Chapter 6.
PAGE 238
239 fitsys fitsys Syntax [sys] = fitsys(data,npoles,nzeros,weight, {skipchks,Hertz}) Parameter List Inputs: data Complex valued data (pdm). npoles nzeros Order of requested fit. (optional, default = 0). Number of zeros in transfer function. (optional, default = npoles) Weighting function. (scalar, pdm, or Dynamic System) (optional, default = 1). weight Keywords: Outputs: Hertz skipchks Boolean.
PAGE 239
Chapter 6. Function Reference 240 The primary use of this routine is the fitting of D scale weights for mu synthesis iterations. Chebyshev polynomials are used as basis functions for both the numerator and denominator polynomials. WARNING: This routine uses iterative polynomial calculations which are not well conditioned for high order (> 6) fits. Reference For further information see: “Curve Fitter for Pole-Zero Analysis,” J.L. Adcock, Hewlett-Packard Journal, p. 33, January 1987.
PAGE 240
241 fitsys real imaginary frequency (rad/sec) -1.0000e+00 5.0000e-01 5.0000e-01 0.0000e+00 3.1225e+00 -3.1225e+00 1.0000e+00 3.1623e+00 3.1623e+00 damping ratio 1.0000 -0.1581 -0.1581 omega = logspace(0.001,100,200) plantg = freq(plant,omega) # # # # # Use complex cepstrum to fit minimum phase equivalent to the magnitude of the data. One of the principle uses of the fitsys function is fitting approximations to noisy data. To illustrate the concepts, no noise is added here.
PAGE 241
Chapter 6. Function Reference 242 Data and minimum phase fit 100 0 Phase (degrees) -100 -200 original data minimum phase fit -300 -400 0.001 0.01 0.
PAGE 242
243 fitsys # Create fitting weight. # data. 1/s works well for logspaced wght = 1/makepoly([1,0],"s") # Fit new system and compare pole location with # the original. Note that it is minimum phase. nsys = fitsys(cdata,3,3,wght) rifd(nsys) Poles: real imaginary frequency (rad/sec) damping ratio -5.1043e-02 -5.1043e-02 -5.0107e+00 -3.1249e-01 3.1249e-01 0.0000e+00 3.1663e-01 3.1663e-01 5.0107e+00 0.1612 0.1612 1.0000 real imaginary frequency (rad/sec) damping ratio -1.0189e+00 -5.0906e-01 -5.
PAGE 243
Chapter 6. Function Reference 244 Data and minimum phase fit 100 original system Magnitude 10 minimum phase system 1 0.1 0.01 0.001 0.01 0.1 1 10 100 1 10 100 Frequency 100 Phase (degrees) 0 -100 -200 -300 -400 0.001 0.01 0.
PAGE 244
fitsys Limitations Limited to SISO systems.
PAGE 245
247 gstep gstep Syntax gPdm = gstep (ytime,timespec,valspec, {skipChks}) Parameter List Inputs: ytime output time vector (seconds). timespec valspec times for specified step data (optional) value for specified step data (optional) Keywords: skipChks Boolean specifying that syntax checking is to be skipped. Outputs: gPdm pdm containing the step values as a function of time. Description This function creates a PDM over the domain: ytime.
PAGE 246
Chapter 6. Function Reference 248 Data and minimum phase fit 100 original system Magnitude 10 minimum phase system 1 0.1 0.01 0.001 0.01 0.1 1 10 100 1 10 100 Frequency 100 Phase (degrees) 0 -100 -200 -300 -400 0.001 0.01 0.
PAGE 247
gstep See Also randpdm, gcos, gsin, gpulse, gsawtooth, gsquarewave 249
PAGE 248
251 hinfnorm hinfnorm Syntax [out,omega] = hinfnorm(sys,tol,{imag eps,max it}) Parameter List Inputs: sys sys Keywords: imag eps max it Outputs: out omega Dynamic System, frequency response (pdm), or constant gain (matrix). Specifies the relative tolerance of the answer when the input is a Dynamic SystemDefault = 0.001. Epsilon value for determining imaginary eigenvalues of the Hamiltonian. Default = sqrt(eps). Maximum number of iterations. Default = 100. H∞ norm of the input system.
PAGE 249
Chapter 6. Function Reference 252 Stable Dynamic System norms are calculated by an iterative Hamiltonian method. In this case out is a 2×1 vector with upper and lower bounds for the norm. Example # Set up a simple closed loop problem. # This example is given in more detail in the # hinfsyn online help. plant = makepoly([0.1,-0.1,1],"s")*makepoly([1,1],"s")... /(makepoly([1,0.1,.1],"s")*(makepoly([0.2,1],"s"))) # Create weights Wperf = 100/makepoly([100,1],"s") Wact = makepoly([0.5,0.05],"s")/makepoly([0.
PAGE 250
253 hinfnorm 25.000 12.500 6.250 3.125 1.562 0.781 1.172 5.2e-01 5.2e-01 5.2e-01 5.1e-01 5.0e-01 3.9e-01 4.8e-01 Gamma value achieved: 1.7e-03 1.7e-03 1.7e-03 1.7e-03 1.7e-03 -9.8e+01 1.8e-03 1.0e-02 1.0e-02 1.0e-02 1.0e-02 1.0e-02 1.0e-02 1.0e-02 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 1.1719 # form weighted closed loop system clpinf = starp(wghtic,Kinf) # Compare hinfnorm to gamma value. The H infinity # norm should be less than or equal to gamma.
PAGE 251
255 h2norm h2norm Syntax out = h2norm(sys) Parameter List Inputs: sys Continuous time Dynamic System Outputs: out H2 norm of the input system Description The H2 norm of a stable, strictly proper system is calculated. This is given by out = trace(CXC 0 ), where X is the controllability grammian, solving the Lyapunov equation, AX + XA0 + BB 0 = 0. Example # Set up a simple closed loop problem. # This example is also studied in the # hinfsyn on-line help. plant = makepoly([0.1,-0.
PAGE 252
Chapter 6. Function Reference 256 Wperf = 100/makepoly([100,1],"s") Wact = makepoly([0.5,0.05],"s")/makepoly([0.05,1],"s") # Form the weighted interconnection structure sysnames sysinp = sysout = syscnx = = ["plant";"Wperf";"Wact"] ["ref";"control"] ["Wperf"; "Wact"; "ref-plant"] ["control"; ... # input to plant "ref-plant"; ... # input to Wperf "control"] # input to Wact wghtic = sysic(sysnames,sysinp,sysout,syscnx,plant,...
PAGE 253
257 hinfsyn hinfsyn Syntax [k,gfin,stat] = hinfsyn(p,nmeas,ncon,gamma,{keywords}) Parameter List Inputs: p nmeas ncon gamma Keywords: control vector dimension. H∞ norm bound of controller. For a bisection search specify gamma = [gamma min;gamma max]. schur solution real Schur decomposition for Riccati solution (default) eig solution tol eigendecomposition for Riccati solution. tol=value. Specifies the relative tolerance for stopping a bisection fit.
PAGE 254
Chapter 6. Function Reference 258 Description The H∞ (sub)optimal controller for the interconnection, p, is calculated. The resulting closed loop system is illustrated below. z w p y u - k The variables ncon and nmeas are used to specify the dimensions of u and y in the above diagram (ncon = dim(u) and nmeas = dim(y)). The objective is to design a stabilizing controller, k, which minimizes the H2 norm of the closed loop system between w and z.
PAGE 255
259 hinfsyn 1. (a, b2 , c2 ) is stabilizable and detectable 2. d12 and d21 have full rank 3. The matrix [a − jωI, b2 ; c1 , d12 ] has full column rank for all ω 4. The matrix [a − jωI, b1 ; c2 , d21 ] has full row rank for all ω Reference This function uses the state-space formulae given in: “State-space formulae for all stabilizing controllers that satisfy an H∞ norm bound and relations to risk sensitivity,” Keith Glover and John Doyle, Systems & Control Letters 11, pp. 167–172., Oct, 1988.
PAGE 256
Chapter 6. Function Reference 260 wghtic = sysic(sysnames,sysinp,sysout,syscnx,plant,... Wperf,Wact) # Design Hinf controller nctrls = 1 nmeas = 1 gmax = 25 gmin = 0 Kinf = hinfsyn(wghtic,nmeas,nctrls,[gmax;gmin]) Test bounds: 0.0000 < gamma <= gamma Hx eig X eig Hy eig 25.000 5.2e-01 1.7e-03 1.0e-02 12.500 5.2e-01 1.7e-03 1.0e-02 6.250 5.2e-01 1.7e-03 1.0e-02 3.125 5.1e-01 1.7e-03 1.0e-02 1.562 5.0e-01 1.7e-03 1.0e-02 0.781 3.9e-01 -9.8e+01 1.0e-02 1.172 4.8e-01 1.8e-03 1.0e-02 Gamma value achieved: 25.
PAGE 257
261 hinfsyn real imaginary frequency (rad/sec) damping ratio -5.0000e-02 -5.0000e-02 -5.0000e+00 -2.0000e+01 -3.1225e-01 3.1225e-01 0.0000e+00 0.0000e+00 3.1623e-01 3.1623e-01 5.0000e+00 2.0000e+01 0.1581 0.1581 1.0000 1.0000 omega = logspace(0.
PAGE 258
Chapter 6. Function Reference 262 Kinf 10 Magnitude 1 0.1 0.01 0.001 0.001 0.01 0.1 1 10 100 1 10 100 Frequency 40 20 Phase (degrees) 0 -20 -40 -60 -80 -100 0.001 0.01 0.
PAGE 259
263 hinfsyn # Use sysic to create unweighted interconnection ic = sysic("plant",["ref";"ctrl"],["plant";"ref-plant"],... "ctrl",plant) clpinf = starp(ic,Kinf) rifd(clpinf) Poles: real imaginary frequency (rad/sec) damping ratio -5.0000e-02 -5.0000e-02 -9.5105e-01 -9.7350e-01 -9.7350e-01 -5.0000e+00 -5.0901e+00 -1.1812e+01 3.1225e-01 -3.1225e-01 0.0000e+00 -1.2285e+00 1.2285e+00 0.0000e+00 0.0000e+00 0.0000e+00 3.1623e-01 3.1623e-01 9.5105e-01 1.5675e+00 1.5675e+00 5.0000e+00 5.0901e+00 1.
PAGE 260
264 Chapter 6.
PAGE 261
265 hinfsyn Kinf controller: sensitivity function 10 Magnitude 1 0.1 0.01 0.001 0.01 0.
PAGE 262
Chapter 6. Function Reference 266 # Examine step response step = gstep([0:0.
PAGE 263
267 hinfsyn Kinf controller: step response 1.2 1 0.8 0.6 0.4 0.
PAGE 264
268 See also hinfsyn, hinfnorm, h2norm Chapter 6.
PAGE 265
269 h2syn h2syn Syntax k = h2syn(p,nmeas,ncon,{keywords}) Parameter List Inputs: Keywords: Outputs: p nmeas Generalized interconnection structure (Dynamic System) measurement vector dimension. ncon control vector dimension. schur solution real Schur decomposition for Riccati solution (default) eig solution epr eigendecomposition for Riccati solution. epr = value. Tolerance for determining when the Hamiltonian eigenvalues lie on the jω axis. Default = 0.5*sqrt(eps) k H2 optimal controller.
PAGE 266
Chapter 6. Function Reference 270 z w p y u - k The variables ncon and nmeas are used to specify the dimensions of u and y in the above diagram (ncon = dim(u) and nmeas = dim(y)). The objective is to design a stabilizing controller, k, which minimizes the H2 norm of the closed loop system between w and z. p is a state-space system, which can be partitioned with respect to [z; y] and [w; u] in the following way.  a p =  c1 c2 b1 d11 d21  b2 d12  . d22 The following assumptions must hold: 1.
PAGE 267
271 h2syn Reference This function uses the state-space formulae given in: “State-space formulae for all stabilizing controllers that satisfy an H∞ norm bound and relations to risk sensitivity,” Keith Glover and John Doyle, Systems & Control Letters 11, pp. 167–172., Oct, 1988. Example # Set up a simple closed loop problem. # This example is also studied in the # hinfsyn on-line help. plant = makepoly([0.1,-0.1,1],"s")*makepoly([1,1],"s")... /(makepoly([1,0.1,.1],"s")*(makepoly([0.
PAGE 268
Chapter 6. Function Reference 272 K2 = h2syn(wghtic,nmeas,nctrls) rifd(K2) Poles: real imaginary frequency (rad/sec) damping ratio -1.4046e-01 -1.4046e-01 -1.5863e+00 -1.5863e+00 -5.2060e+00 -2.3161e-01 2.3161e-01 3.4754e+00 -3.4754e+00 0.0000e+00 2.7087e-01 2.7087e-01 3.8203e+00 3.8203e+00 5.2060e+00 0.5186 0.5186 0.4152 0.4152 1.0000 real imaginary frequency (rad/sec) damping ratio -5.0000e-02 -5.0000e-02 -5.0000e+00 -2.0000e+01 -3.1225e-01 3.1225e-01 0.0000e+00 0.0000e+00 3.1623e-01 3.
PAGE 269
273 h2syn K2 1 Magnitude 0.1 0.01 0.001 0.0001 0.001 0.01 0.1 1 10 100 1 10 100 Frequency 50 Phase (degrees) 0 -50 -100 -150 0.001 0.01 0.
PAGE 270
Chapter 6. Function Reference 274 # Use sysic to create unweighted interconnection ic = sysic("plant",["ref";"ctrl"],["plant";"ref-plant"],... "ctrl",plant) clp2 = starp(ic,K2) rifd(clp2) Poles: real imaginary frequency (rad/sec) damping ratio -5.0000e-02 -5.0000e-02 -1.8120e-01 -1.8120e-01 -1.6753e+00 -1.6753e+00 -4.9957e+00 -5.0000e+00 3.1225e-01 -3.1225e-01 -4.2333e-01 4.2333e-01 3.4263e+00 -3.4263e+00 0.0000e+00 0.0000e+00 3.1623e-01 3.1623e-01 4.6048e-01 4.6048e-01 3.8139e+00 3.8139e+00 4.
PAGE 271
h2syn sensg = freq(sens,omega) gph2 = ctrlplot(sensg,{logmagplot}); gph2 = plot(gph2,{title="K2 controller: sensitivity function"})? 275
PAGE 272
Chapter 6. Function Reference 276 K2 controller: sensitivity function Magnitude 10 1 0.1 0.001 0.01 0.
PAGE 273
h2syn # Examine step response step = gstep([0:0.
PAGE 274
Chapter 6. Function Reference 278 K2 controller: step response 1 0.8 0.6 0.4 0.
PAGE 275
h2syn See also hinfsyn, h2norm, hinfnorm 279
PAGE 276
281 interp interp Syntax outpdm = interp(inpdm,stepsize,final {keywords}) outpdm = interp(inpdm,domspec, {keywords}) Parameter List Inputs: inpdm stepsize Input pdm. Increment in outpdm domain. final Last value in outpdm domain. Optional: default = max(domain(inpdm)). Regular vector or pdm used to specify the domain of outpdm. domspec Keywords: order Outputs: outpdm Interpolation order.
PAGE 277
Chapter 6. Function Reference 282 This function differs from the interpolate function in that it can handle zero order hold type interpolation and deal with irregularly spaced input pdms. Irregularly spaced output pdms can be generated with the domspec syntax. These features are often useful when dealing with data generated from experiments. Example time = [0:1:5] u = gstep(time,time,time) u1 = interp(u,0.
PAGE 278
283 interp 5 4 1st order interp. 0 order interp.
PAGE 279
284 See Also interpolate Chapter 6.
PAGE 280
285 mergeseg mergeseg Syntax outpdm = mergeseg(pdm1,pdm2, {keywords}) Parameter List Inputs: Keywords: pdm1 input pdm pdm2 input pdm domsort Sort the result of merging the pdms. If !domsort then pdm2 is simply concatenated onto pdm1. Boolean. Default = 1. Sort in increasing order. Boolean. Default = 1. Sort in decreasing order. Boolean. Default = 0. increasing decreasing duplicates tol Outputs: outpdm Leave duplicated domain values in outpdm.
PAGE 281
Chapter 6. Function Reference 286 Example time1 = [0:0.025:1] pdm1 = gsin(time1,{frequency=2}) time2 = [0.8:0.02:1.5] pdm2 = randpdm([],1,1,{dom=time2,zeromean}) outpdm = mergeseg(pdm1,pdm2) gph1 = ctrlplot(outpdm); gph1 = plot(pdm1,gph1,{marker=1,marker style=8,... line=0}); gph1 = plot(pdm2,gph1,{marker=1,marker style=1,...
PAGE 282
287 mergeseg 1 outpdm pdm1 pdm2 0.5 0 -0.5 -1 0 0.5 1 1.
PAGE 283
289 mkpert mkpert Syntax [pertsys] = mkpert(Delta,blk,mubnds,{fselect,pnorm,Hertz}) Parameter List Inputs: Delta blk mubnds Keywords: fselect pnorm Hertz Outputs: pertsys Lower bound perturbation data from mu calculation. (pdm) block structure (refer to mu function documentation). (matrix) Calculated mu bounds (pdm). Scalar valued. Specifies the frequency for interpolation, overriding that obtained from mubnds. Scalar valued.
PAGE 284
290 Chapter 6. Function Reference The bounds, mubnds, are used to determine the ”worst-case” frequency for the interpolation. The norm of pertsys is 1/(norm(Delta(jω))) where ω is the chosen frequency. This is the smallest destabilizing perturbation at that frequency. Both the interpolation frequency and the norm of pertsys can be specified via keywords. The user may be interested only in certain frequency ranges and may have an assumed norm bound on the perturbation (typically unity).
PAGE 285
291 mkphase mkphase Syntax [cdata] = mkphase(magdata, {skipchks,Hertz}) Parameter List Inputs: magdata Magnitude data (pdm) Keywords: skipchks Hertz Boolean. Skip the error checking. (Default = 0) The domain of the pdm is in Hertz. This is the default. !Hertz specifies a domain in rad/sec. Outputs: cdata Complex valued data corresponding to a minimum phase transfer function. Description Fits phase data to magnitude data. The phase is equivalent to that produced by minimum phase system.
PAGE 286
292 Limitations Limited to SISO systems. See Also fitsys, ccepstrum Chapter 6.
PAGE 287
293 modalstate modalstate Syntax outsys = modalstate(sys, {keywords}) Parameter List Inputs: sys Input Dynamic System Keywords: increasing Boolean. Order in terms of increasing magnitude (continuous) or angle (discrete). Default = 1 Boolean. Order in terms of decreasing magnitude (continuous) or angle (discrete). Default = 0 Boolean. Separate the stable from the unstable modes and order separately. Default = 0.
PAGE 288
Chapter 6. Function Reference 294 -8.1602e-01 -8.1602e-01 -4.6142e+00 -4.6142e+00 -2.2478e+01 -3.9525e+01 1.3353e+00 -1.3353e+00 -1.6648e+01 1.6648e+01 0.0000e+00 0.0000e+00 1.5649e+00 1.5649e+00 1.7275e+01 1.7275e+01 2.2478e+01 3.9525e+01 0.5215 0.5215 0.2671 0.2671 1.0000 1.0000 real imaginary frequency (rad/sec) damping ratio -1.0270e+00 -1.0270e+00 -4.5529e+00 -4.5529e+00 -2.4416e+01 -3.9534e+01 -1.5086e+00 1.5086e+00 -1.6621e+01 1.6621e+01 0.0000e+00 0.0000e+00 1.8250e+00 1.8250e+00 1.
PAGE 289
295 mu mu Syntax [mubnds,D,Dinv,Delta,sens] = mu(M,blk) Parameter List Inputs: Outputs: M Matrix or pdm. blk Block structure defined by a matrix of dimension: number of blocks × 2. If the ith block has c outputs and r inputs, then blk(i,:) = [r,c]. The default is equivalent to 1x1 blocks (M must be square). mubnds D,Dinv Upper and lower bounds (in vector form) for mu(M). D-scale matrices giving the calculated upper-bound. mu(M) ≤ msv(D*M*Dinv) Perturbation achieving the lower bound.
PAGE 290
Chapter 6. Function Reference 296 # that mu is not equal to its upper bound for # more than three full blocks. gamma = 3 + sqrt(3); beta = sqrt(3) -1 a = sqrt(2/gamma); b = 1/sqrt(gamma) c = 1/sqrt(gamma); d = -sqrt(beta/gamma) f = (1+jay)*sqrt(1/(gamma*beta)) psi1 = -pi/2; psi2 = pi U = V = scl M = [a,0; b,b; c,jay*c; d,f] [0,a; b,-b; c,-jay*c; f*exp(jay*psi1), d*exp(jay*psi2)] = diagonal(random(4,1)+0.
PAGE 291
297 mu ans (a scalar) = 2.81264 mubnds2? mubnds2 (a column vector) = 2.81264 2.81264 det(eye(4,4) - M*Delta2)? ans (a scalar) = -2.53156e-16 + 4.3828e-17 j For an example of how mu is used for system robustness analysis, refer to the on-line help for musynfit (page 299). Limitations This version of the software cannot handle repeated blocks or real valued blocks.
PAGE 292
299 musynfit musynfit Syntax [Dsys,Dinvsys] = musynfit(Dmag,blk,nmeas,nctrls,.. weight,M,order,{keywords}) Parameter List Inputs: Dmag blk New D matrix from mu calculation (magnitude data only). The domain is assumed to be in Hertz. block structure (refer to mu function documentation). nmeas nctrls scalar: number of measurements scalar: number of controls weight (optional) weighting function for the fit. Typically the sensitivity output of mu is used.
PAGE 293
Chapter 6. Function Reference 300 Keywords: Hertz plotweight fit Boolean. This keyword is mandatory as the function must know whether the domain is in Hertz or radians/second (specified by !Hertz) to fit correctly. Note that the Xmath function freq assumes that the frequency range is specified in Hertz. Boolean, default = 0. This will generate a second plot showing the weighting function. This can be useful is assessing the quality of a given transfer function fit.
PAGE 294
musynfit P = 1/makepoly([1,0,-0.01],"s") W = makepoly([1,20],"s")/makepoly([1,200],"s") # Set up an unweighted interconnection structure # for unity gain negative feedback. We include # the perturbation too. nms = ["P";"W"] inp = ["delt";"ref";"noise";"control"] outp = ["W" ;"ref- delt - P";"control";"ref-delt-P-noise"] cnx = ["control";"P"] ic = sysic(nms,inp,outp,cnx,P,W) # Now include some weights for performance: Wperf = makepoly([0.01,1],"s")/makepoly([1,0.01],"s") Wact = 0.
PAGE 295
Chapter 6. Function Reference 302 1.562 1.406 1.484 1.523 6.8e-01 6.7e-01 6.7e-01 6.7e-01 Gamma value achieved: 8.0e-04 8.1e-04 8.0e-04 8.0e-04 7.7e-03 7.0e-03 7.4e-03 7.5e-03 0.0e+00 -1.1e-18 0.0e+00 -5.0e-19 0.6811 2.0900 1.0394 0.8249 1.5234 G = starp(wghtic,Kinf) omega = logspace(0.
PAGE 296
303 musynfit mu analysis 1.2 1 0.8 nominal perf. robust stab. robust perf.(upper) 0.6 robust perf. (lower) 0.4 0.2 0 0.01 0.
PAGE 297
304 Chapter 6.
PAGE 298
305 musynfit D scale fit, block: 1 10000 Magnitude data Magnitude 1000 Previous fit, order: 0 New fit, order: 3 100 10 1 0.01 0.1 1 10 100 10 100 10 100 Frequency (Hz) Upper bound comparison, block: 1 100 Magnitude 10 1 Data based bound 0.1 Previous fit bound, order: 0 0.01 0.001 0.01 New fit bound, order: 3 0.1 1 Frequency (Hz) D scale fit weight for block: 1 1 Magnitude 0.1 0.01 0.001 0.0001 1e-05 0.01 0.
PAGE 299
Chapter 6. Function Reference 306 # Apply the D scales to another H infinity design Kmu = hinfsyn(Ds*wghtic*Dinvs,nmeas,ncntrls,[gmin;gmax]) Test bounds: 0.0000 < gamma <= 10.0000 X eig Hy eig Y eig nrho xy p/f gamma Hx eig 10.000 6.5e-01 5.6e-07 9.9e-03 -1.0e-15 0.0027 p 5.000 6.4e-01 5.6e-07 9.8e-03 -6.3e-16 0.0113 p 2.500 6.4e-01 5.6e-07 9.2e-03 0.0e+00 0.0517 p 1.250 6.3e-01 5.6e-07 6.0e-03 -7.6e-20 0.4688 p 0.625 5.6e-01 -1.1e+00 7.1e-14 ******* ******* f 0.938 6.1e-01 -2.8e+01 7.
PAGE 300
307 musynfit Kmu & Kinf mu analysis 1.2 1 0.8 0.6 Kmu: robust perf. Kinf: robust perf. 0.4 0.2 0 0.01 0.
PAGE 301
Chapter 6. Function Reference 308 # # # # # Look at the worst case perturbations for each of the Kinf and Kmu controllers. Compare also a random perturbation for each controller. In all cases the perturbation is of size 0.5 and we choose a perturbation which is bad at 1Hz. mupert = mkpert(Delta2,blk,rpbnds2,{fselect=1,pnorm=0.5,Hertz}) mupert = mupert(1,1) # select part to close around top infpert = mkpert(Delta1,blk,rpbnds1,{fselect=1,pnorm=0.
PAGE 302
musynfit ymunom(1,1),ymupert(1,1)]); gph4 = plot(gph4,{legend=["input step";"Kinf nominal";... "Kinf pert.";"Kmu nominal";"Kmu pert.
PAGE 303
Chapter 6. Function Reference 310 1.4 1.2 1 0.8 input step Kinf nominal Kinf pert. Kmu nominal 0.6 Kmu pert. 0.4 0.
PAGE 304
musynfit # Compare with a random perturbation yinfrandp = infrpert*step ymurandp = murpert*step gph5 = ctrlplot([step,yinfnom(1,1),yinfrandp(1,1),... ymunom(1,1),ymurandp(1,1)]); gph5 = plot(gph5,{legend=["input step";"Kinf nominal";... "Kinf rand. pert.";"Kmu nominal";"Kmu rand. pert.
PAGE 305
Chapter 6. Function Reference 312 1.4 1.2 1 0.8 input step Kinf nominal Kinf rand. pert. Kmu nominal Kmu rand. pert. 0.6 0.4 0.
PAGE 306
musynfit See Also: mu, hinfsyn, mkpert, hinfnorm.
PAGE 307
315 ophank ophank Syntax [SysR,SysU,HSV] = ophank(Sys,{nsr,onepass}) Parameter List Inputs: Sys Linear, stable, state-space system (continuous) nsr (optional) Order of the reduced system. If not specified, the user will be prompted for its value after the Hankel singular values are displayed. Keywords: onepass (Boolean) If equal to 1 (true), reduction is calculated in one pass. If false (!onepass or onepass=0), reduction is calculated in (number of states of Sys - nsr) passes. Defaults to 1.
PAGE 308
Chapter 6. Function Reference 316 Uses additional subroutines ophiter, ophred, ophmult and stable. This function is cross-licensed from the Model Reduction Module. Example # Create a five state system for reduction. a = daug(-0.891334,[-1.20857,0.799042;-0.799042,-1.20857],... -4.74685,-21.3013) b = [0.0262569;-0.189601;-0.113729;0.211465;-0.538239] c = [0.120725,-0.336942,0.397198,-0.700524,-1.02235] d = 0 sys1 = system(a,b,c,d) fHz = logspace(0.
PAGE 309
317 ophank 0.1 Magnitude 0.01 0.001 original system reduced: ophank reduced: sresidualize ophank error sresidualize error 0.0001 0.01 0.
PAGE 310
318 See Also minimal, balmoore. Chapter 6.
PAGE 311
319 orderstate orderstate Syntax outsys = orderstate(sys,indx) Parameter List Inputs: Outputs: sys Input Dynamic System indx Lists the desired order of the states in outsys. outsys output dynamic system. Description Reorder the states according to the index argument. The state names and initial condition are also ordered. The input and output names, and period (if discrete) are preserved. Example sys1 = randsys(3,1,1,{stable}) sys1 = system(sys1,{statenames=["s1";"s2";"s3"],... x0=[0.1;0.2;0.
PAGE 312
Chapter 6. Function Reference 320 0.710595 0.688873 C 0.659532 0.181512 0.390497 D 0.15869 X0 0.1 0.2 0.3 State Names ----------s1 s2 s3 System is continuous sys2 = orderstate(sys1,[2,1,3])? sys2 (a state space system) = A -5.41526 3.32295 4.7424 0.65518 -2.27378 -1.2488 9.0643 -8.82005 -13.2781 0.659532 0.390497 B 0.710595 0.0879738 0.688873 C 0.181512 D 0.
PAGE 313
orderstate 0.2 0.1 0.
PAGE 314
323 randpdm randpdm Syntax pdmout = randpdm (ndomain,nrows,ncolumns,{keywords}) Parameter List Inputs: Keywords: ndomain length of the domain nrows ncolumns number of rows in pdmout number of columns in pdmout complex Boolean. A complex valued pdm is generated. Default = 0. Boolean. The values are shifted so that zero is the mean. Default = 0. First value in the domain. This only a bound if !regular. Default = 0. Last value in the domain. Again this is only a bound if !regular.
PAGE 315
Chapter 6. Function Reference 324 Examples pdm0 = randpdm(3,1,2,{zeromean})? pdm0 (a pdm) = domain | Col 1 Col 2 -------+------------------------0 | -0.043917 -0.0789571 -------+------------------------1 | 0.37878 -0.457265 -------+------------------------2 | 0.0466939 -0.34413 -------+------------------------pdm1 = randpdm(4,1,1,{!regular,Dfirst=2,Dlast=19})? pdm1 (a pdm) = domain | ---------+-----------3.1543 | 0.0605523 ---------+-----------4.87235 | 0.900169 ---------+-----------10.0323 | 0.
PAGE 316
325 randpdm 10.0323 | 0.922528 ---------+-----------18.3191 | 0.81113 ---------+-----------pdm3 = randpdm(0,3,2,{complex,zeromean})? pdm3 (a rectangular matrix) = -0.237052 - 0.10845 j 0.097563 + 0.0140395 j 0.0178566 - 0.153032 j -0.250958 - 0.392096 j -0.255142 - 0.419362 j 0.0986898 + 0.
PAGE 317
327 randpert randpert Syntax [pert] = randpert(blk, {sys,sfreq,complex,pnorm}) Parameter List Inputs: blk block structure (refer to mu function documentation). (matrix) Keywords: sys Boolean. Specifies that pert is a dynamic system. Default = !sys, i.e. pert is a matrix. Scalar valued. If a dynamic system is specified, sfreq is the frequency with significant phase. Units are Hertz. Default = 1.0. Boolean. Specifies that pert is complex valued. Default = 1. Scalar valued.
PAGE 318
328 Chapter 6. Function Reference Example The use of randpert is studied in context in the on-line help for musynfit and the manual documentation for musynfit (page 299).
PAGE 319
329 randsys randsys Syntax sys = randsys (nstates,noutputs,ninputs,{keywords}) Parameter List Inputs: Keywords: nstates number of states in sys noutputs ninputs number of outputs in sys number of inputs in sys stable oscillatory Boolean. sys is forced to be stable (default = 1) Boolean. Oscillatory poles are allowed. (default = 1). !oscillatory gives only real eigenvalues (continuous) or positive real eigenvalues (discrete). Boolean. Generate a discrete time system. (default = 0).
PAGE 320
Chapter 6. Function Reference 330 Description A random system, with user specified state, input and output dimension, is generated. Several additional features can be specified by keywords: whether or not the system is stable, where or not it can contain oscillatory modes, whether or not it has a D term, and bounds on the minimum and maximum pole frequencies. If the discrete keyword is specified, a bilinear transformation is performed to get the random discrete system.
PAGE 321
331 randsys real imaginary frequency (rad/sec) damping ratio -1.0408e+00 -1.1390e+00 -2.8050e-01 -2.8050e-01 -1.5337e+00 -5.5566e-01 -5.5566e-01 -3.1454e+00 -5.6798e+00 -6.1039e+00 0.0000e+00 0.0000e+00 -1.4858e+00 1.4858e+00 0.0000e+00 -1.6958e+00 1.6958e+00 0.0000e+00 0.0000e+00 0.0000e+00 1.0408e+00 1.1390e+00 1.5120e+00 1.5120e+00 1.5337e+00 1.7845e+00 1.7845e+00 3.1454e+00 5.6798e+00 6.1039e+00 1.0000 1.0000 0.1855 0.1855 1.0000 0.3114 0.3114 1.0000 1.0000 1.
PAGE 322
Chapter 6. Function Reference 332 10 Magnitude 1 0.1 sys1 (input 1) sys1 (input 2) sys2 0.01 0.001 0.01 0.1 1 10 100 10 100 Frequency 20 0 Phase (degrees) -20 -40 -60 -80 -100 0.01 0.
PAGE 323
333 randsys sys3 = randsys(6,1,1,{discrete,dt=5}) rifd(sys3)? Poles: radius 0.1652 0.4580 0.7872 0.7872 0.9307 0.9307 angle (radians) 0.0000 0.0000 -0.0394 0.0394 -0.2163 0.2163 Zeros: radius 0.1623 0.7981 0.6658 0.8921 0.8921 15.4278 angle (radians) 0.0000 0.0000 0.0000 -0.1981 0.1981 3.
PAGE 324
335 rifd rifd Syntax [stat] = rifd(vec,{discrete,Hertz,degrees}) Parameter List Inputs: vec complex valued vector (or Dynamic System - see below). Keywords: discrete Boolean. Vector is to be interpreted in the z domain, rather than the s domain. Default = 0. Boolean. Display frequency units in Hertz on the s plane. !Hertz gives a display in radians/sec. Default = 0. Boolean. Display angle in degrees for z plane. !degrees gives a display in radians. Default = 0.
PAGE 325
Chapter 6. Function Reference 336 Examples sys1 = randsys(4,3,2,{stable}) rifd(sys1) Poles: real imaginary frequency (rad/sec) damping ratio -1.6847e+00 -2.4383e+00 -8.7457e+00 -1.5041e+01 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 1.6847e+00 2.4383e+00 8.7457e+00 1.5041e+01 1.0000 1.0000 1.0000 1.0000 Zeros: ans (a scalar) = 0 # Compare rifd to A matrix eigenvalues [a,b,c,d] = abcd(sys1) Aeigs = eig(a)? Aeigs (a column vector) = -1.68474 -2.43829 -8.7457 -15.
PAGE 326
337 rifd sys2 = randsys(3,3,2,{stable,discrete}) rifd(sys2) Poles: radius 0.9641 0.2985 0.2985 Zeros: angle (radians) 0.0000 -0.5348 0.
PAGE 327
339 sdtrsp sdtrsp Syntax [v,y,u] = sdtrsp(Sys,dSys,w,tfinal,... {ord,intstep,cdelay}) Parameter List Inputs: Sys dSys w Keywords: tfinal Final time in the simulation (optional). Default = max time specified in w. ord Scalar valued. Specifies order of interpolation for continuous signals in the system. Options are 0 or 1. Default = 0. Scalar valued. Integration step size. If not supplied a default will be calculated based on the system eigenvalues.
PAGE 328
Chapter 6. Function Reference 340 Description Time domain simulation of a sampled data interconnection. The applicable closed loop system is illustrated below. v w Sys y u - dSys This is conceptually the equivalent of: v = starp(Sys, dSys) ∗ w. This function will handle interconnections in which the continuous time signals are not necessarily synchronized with the digital system. Computation delays which are a fraction of the sample period may also be simulated.
PAGE 329
sdtrsp 341 triangle hold equivalent is used. This is the same a linearly connecting the samples at the input to the hold. Some care is needed in the choice of cdelay and intstep. The default for intstep is based on the continuous system eigenvalues and the minimum time spacing in the input vector, w. It will be forced to be an integer divisor of the digital sample time. Intstep must also be a divisor of cdelay.
PAGE 330
Chapter 6. Function Reference 342 inps = "ref" outs = "digP" cnx = ["digC"; ... # input to digP "[ ref;digP]"] # input to digC digclp = sysic(snm,inps,outs,cnx,digP,digC) # Calculate the digital system step response time = [0:T:20*T] step = gstep(time,0,1) digy = digclp*step # # # # # # # # Do a complete sampled data simulation. We will assume that the controller has an input to output calculation delay of 0.1*T.
PAGE 331
343 sdtrsp 4 sampled data calc. discrete calc. 3 2 1 0 -1 -2 0 0.2 0.4 0.6 0.
PAGE 332
344 See Also trsp Chapter 6.
PAGE 333
345 simtransform simtransform Syntax out = simtransform(sys,X) Parameter List Inputs: sys Outputs: X Input system. This may be a state-space system, pdm or constant matrix. Similarity transform. X must be invertible out output system - in the same class as the input. Description Apply a similarity transform to sys. If sys is a matrix or pdm this gives out = inv(X)*sys*X. If sys is a Dynamic System this transform is applied to the state.
PAGE 334
Chapter 6. Function Reference 346 -1.9737e+01 -9.7569e+00 -9.7569e+00 0.0000e+00 -2.9129e+01 2.9129e+01 1.9737e+01 3.0720e+01 3.0720e+01 1.0000 0.3176 0.3176 real imaginary frequency (rad/sec) damping ratio -2.1990e+01 -9.7962e+00 -9.7962e+00 0.0000e+00 -2.9186e+01 2.9186e+01 2.1990e+01 3.0786e+01 3.0786e+01 1.0000 0.3182 0.3182 Zeros: [a,b,c,d] = abcd(sys1) # Do a schur transform on the A matrix # and use for a similarity transform on the system.
PAGE 335
simtransform X0 0 0 0 Input Names ----------Input 1 Output Names -----------Output 1 System is continuous 347
PAGE 336
349 spectrad spectrad Syntax out = spectrad(mat) Parameter List Inputs: mat Square matrix or pdm. Outputs: out spectral radius of the input. Description Calculates the spectral radius (magnitude of the maximum eigenvalue) of the input matrix. Example pdm1 = randpdm(3,2,2,{zeromean})? pdm1 (a pdm) = domain | Col 1 Col 2 -------+-----------------------------0 | Row 1 0.0504605 -0.0914956 | Row 2 0.221744 -0.0231464 -------+-----------------------------1 | Row 1 0.139306 0.496387 | Row 2 -0.342521 0.
PAGE 337
Chapter 6. Function Reference 350 eig(pdm1)? ans (a pdm) = domain | -------+------------------------------0 | Row 1 0.013657 + 0.137601 j | Row 2 0.013657 - 0.137601 j -------+------------------------------1 | Row 1 0.0871876 + 0.409031 j | Row 2 0.0871876 - 0.409031 j -------+------------------------------2 | Row 1 -0.311974 | Row 2 -0.452314 -------+------------------------------spectrad(pdm1)? ans (a pdm) = domain | -------+----------0 | 0.138277 -------+----------1 | 0.41822 -------+----------2 | 0.
PAGE 338
351 sresidualize sresidualize Syntax sysout = sresidualize(sysin,ord) Parameter List Inputs: Outputs: sysin Input Dynamic System ord Order of sysout sysout output Dynamic System Description Residualize the states of sysin to the number specified by ord. The last nx −ord states (nx is the number of states in sysin) are residualized. If ord is greater than the number of states in sysin then sysout = sysin and a warning is displayed.
PAGE 339
Chapter 6. Function Reference 352 sysout1 = sresidualize(sys1,3) sysout2 = truncate(sys1,3) fHz = logspace(0.01,100,100) sys1g = freq(sys1,fHz) sysout1g = freq(sysout1,fHz) sysout2g = freq(sysout2,fHz) residerror = sys1g - sysout1g truncerror = sys1g - sysout2g gph1 = ctrlplot([sys1g,sysout1g,sysout2g,residerror,... truncerror],{logmagplot}); gph1 = plot(gph1,{!grid,legend=["original system";... "residualized system";"truncated system";...
PAGE 340
353 sresidualize 0.1 Magnitude 0.01 0.001 original system residualized system 0.0001 truncated system residualization error truncation error 1e-05 0.01 0.
PAGE 341
354 Chapter 6. Function Reference See also rifd, simtransform, orderstate modalstate, truncate.
PAGE 342
355 starp starp Syntax out = starp (upper,lower,dim1,dim2,skipChks) Parameter List Inputs: upper lower dim1 dim2 Upper object (Dynamic System, constant or pdm) in the interconnection Lower object in the interconnection. The number of outputs of upper to be connected as inputs to lower. (Optional - see description for default) The number of outputs of lower to be connected as inputs to upper. (Optional - see description for default) Keywords: skipChks (Boolean, default = 0).
PAGE 343
Chapter 6. Function Reference 356 specified. This is equivalent to: dim1 = min(upper output dim, lower input dim) and dim2 = min(upper input dim, lower output dim). Examples # Look at a constant matrix, M = random(4,4) lower = [3,4] # The interconnected LFT will have 2 rows and 3 columns result = starp(M,lower)? result (a rectangular matrix) = 0.261939 -0.648616 -0.558202 -1.11273 0.484644 0.
PAGE 344
357 starp # structure for closing control loops. The following is # the structure for a simple unity gain negative feedback # system with plant P. P = 1/makepoly([1,1],"s") M = consys([0,1;1,-1])*daug(1,P) # Test this with a Proportional and PI controller Kp = 10 KpInt = makepoly([10,100],"s")/makepoly([1,0],"s") rifd(starp(M,Kp)) # look at poles and zeros Poles: real imaginary frequency (rad/sec) damping ratio -1.1000e+01 0.0000e+00 1.1000e+01 1.
PAGE 345
Chapter 6. Function Reference 358 -1.0000e+01 0.0000e+00 1.0000e+01 1.0000 # Now consider the system with an output multiplicative # perturbation (of 10%) W = consys(0.1) G = daug(W,1,1)*consys([0,0,1;1,0,1;-1,1,-1])*daug(1,1,P) # The nominal system is constructed by closing the # upper loop with 0. Any thing else is a perturbed # system. Gnom = starp(0,G) Gpert = starp(-1,G) # Now examine the closed loop affects for the PI # controller.
PAGE 346
359 starp Poles: real imaginary frequency (rad/sec) damping ratio -5.0000e+00 -5.0000e+00 -8.0623e+00 8.0623e+00 9.4868e+00 9.4868e+00 0.5270 0.5270 real imaginary frequency (rad/sec) damping ratio -1.0000e+01 0.0000e+00 1.0000e+01 1.
PAGE 347
361 substr substr Syntax littlestring = substr(bigstring,charindex,{skipChks}) Parameter List Inputs: bigstring Input string (a 1×1 string matrix) charindex vector indexing the characters to be returned in the output. Keywords: skipChks Boolean specifying that syntax checking is to be skipped. Outputs: littlestring Output string Description substr takes characters from the string, bigstring and concatentates them to from the output, littlestring.
PAGE 348
363 sysic sysic Syntax [sys] = sysic(sysNames,sysInputs,sysOutputs,sysConnects,... subsys1,subsys2,...) Parameter List Inputs: sysNames sysInputs sysOutputs sysConnects subsys1 .. . Outputs: sys A vector of strings (of the same length as the number of subsystems) naming the subsystems. A vector of strings naming the exogenous inputs to the final system. Each named input must be a different row — however parenthesis can be used to denote vector valued input names.
PAGE 349
Chapter 6. Function Reference 364 Matrices can also be included in the interconnection and are considered to be constant gains when connected with Dynamic Systems and considered to be constant for all domain values when connected with pdms. sysic is able to connect linear systems together. For more complete interconnection and simulation capabilities SystemBuild should be used. Limitations Only 20 subsystems can be interconnected with a single sysic call.
PAGE 350
365 trsp trsp Syntax [y,uint] = trsp(Sys,u,tfinal, ord,intstep) Parameter List Inputs: Sys u tfinal Keywords: ord intstep Outputs: y uint Continuous dynamic system. The initial states are used in the simulation. PDM. Input signal. Final time in the simulation (optional). Default = max time specified in u. Scalar valued. Specifies order of interpolation for continuous signals in the system. Options are 0 or 1. Default = 0. Scalar valued. Integration step size.
PAGE 351
Chapter 6. Function Reference 366 necessary. If intstep is not specified the integration stepsize is determined from the system eigenvalues and the minimum spacing in the input signal, u. A zero or first order discrete equivalent of the systems can be specified. The standard Xmath * operator uses only a zero order discretization. The first order simulation actually uses a triangle hold equivalent. This is equivalent to connecting the samples going into the hold function.
PAGE 352
367 trsp 1.5 1 0.5 0 -0.5 output input -1 -1.5 0 0.5 1 1.
PAGE 353
368 Chapter 6. Function Reference # Compare trsp calculation to standard [ytrsp,uint] = trsp(P,u) gph2 = ctrlplot(ytrsp); gph2 = plot(y,gph2,{line style=4}); gph2 = plot(u,gph2,{line=0,marker=1}); gph2 = plot(uint,gph2,{line style=3,legend=["trsp calc.";... "* calc.";"input";"interpolated input"],title=...
PAGE 354
369 trsp Time response calculation comparisons 1.5 1 0.5 0 -0.5 trsp calc. * calc. input -1 interpolated input -1.5 0 0.5 1 1.
PAGE 355
Chapter 6. Function Reference 370 # Now look at 1st order interpolation [y1trsp,u1int] = trsp(P,u,{ord=1}) gph3 = ctrlplot([y1trsp,u1int]); gph3 = plot(u,gph3,{line=0,marker=1,legend=... ["1st order interpolation";... "interpolated input";"input data"],title=...
PAGE 356
371 trsp Time response calculation comparisons 1.5 1st order interpolation interpolated input input data 1 0.5 0 -0.5 -1 -1.5 0 0.5 1 1.
PAGE 357
373 sresidualize truncate Syntax sysout = truncate(sysin,ord) Parameter List Inputs: Outputs: sysin Input Dynamic System ord Order of the truncated system: sysout sysout Truncated output Dynamic System Description The function truncate is cross-licensed from the model reduction toolbox and has slightly more capabilities than those described here. See the online help for further details. Consider a partitioning of the input Dynamic System, sysin, as follows.
PAGE 358
Chapter 6. Function Reference 374 Example This example is identical to that described for sresidualize and compares the two methods of model reduction. # Create a five state system for reduction. a = daug(-0.891334,[-1.20857,0.799042;-0.799042,-1.20857],... -4.74685,-21.3013) b = [0.0262569;-0.189601;-0.113729;0.211465;-0.538239] c = [0.120725,-0.336942,0.397198,-0.700524,-1.02235] d = 0 sys1 = system(a,b,c,d) # Reduce to a 3 state system by residualization # and truncation.
PAGE 359
375 sresidualize 0.1 Magnitude 0.01 0.001 original system residualized system 0.0001 truncated system residualization error truncation error 1e-05 0.01 0.
PAGE 360
376 Chapter 6. Function Reference See also rifd, simtransform, sresidualize, orderstate modalstate.
PAGE 361
6.2. Xµ SUBROUTINES AND UTILITIES 6.2 377 Xµ Subroutines and Utilities Several subroutines may also be of interest to the user. These subroutines typically perform self contained parts of a calculation. They may be of interest to those developing new robust control algorithms or wishing to study the calculation details of the algorithms given here. Beware of the fact that these subroutines may not contain error checking.
PAGE 362
blkbal 379 blkbal Syntax d = blkbal(M) Description Balances a square matrix assuming only scalar blocks. The Osborne method (growth rate: n2 ) is used for large systems and the Perron method (growth rate: n3 ) for smaller systems. The Perron method will exactly calculate mu for positive matrices.
PAGE 363
hinfcalc hinfcalc Syntax [X,Y,f,h,Ric fail,HX,HY,HXmin,HYmin] = ...
PAGE 364
Chapter 6. Function Reference 382 Parameter List Inputs: p nmeas ncon g epr Generalized interconnection structure (Dynamic System) measurement vector dimension. control vector dimension. H∞ norm of suboptimal controller to be calculated. Referred to in the literature as gamma. Tolerance for determining when the Hamiltonian eigenvalues lie on the jω-axis. Keywords: schur solution eig solution real Schur decomposition for Riccati solution (default) eigendecomposition for Riccati solution.
PAGE 365
hinfcalc Description Form and solve the Riccati equations for the H∞ control problem. X and Y are the resulting Riccati solutions. THIS FUNCTION IS INTENDED ONLY AS A SUBROUTINE CALLED BY THE HINFSYN FUNCTION.
PAGE 366
powermu 385 powermu Syntax [lbnd,delta,errstat] = powermu(M,blk,rp,cp) Description Lower bound power algorithm based on the work of Andy Packard. The vector naming roughly corresponds to that in his thesis.
PAGE 367
387 riccati eig riccati eig Syntax [x1,x2,stat,Heig min] = riccati eig(H,epp) Parameter List Inputs: H epp Hamiltonian matrix. Tolerance for detecting proximity of eigenvalues to the jω axis. Outputs: x1,x2 Basis vectors for stable subspace. See description below. stat Status flag. 0 Stable subspace calculated. 1 Failure to decompose into stable and unstable subspaces. Minimum absolute value of the real part of the eigenvalues of H.
PAGE 368
388 Chapter 6. Function Reference variable, X = x2 x−1 1 , is the stabilizing solution to the Riccati equation. If H has jω axis eigenvalues then no stabilizing solution exists and the function returns a failure status. If any eigenvalue of H is within epp of the jω axis it is considered to lie on the jω axis and no solution is found. This may be due to numerical problems in finding the eigenvalues of poorly conditioned problems even when a stabilizing solution exists.
PAGE 369
389 riccati schur riccati schur Syntax [x1,x2,stat,Heig min] = riccati schur(H,epp) Parameter List Inputs: H epp Hamiltonian matrix. Tolerance for detecting proximity of eigenvalues to the jω axis. Outputs: x1,x2 Basis vectors for stable subspace. See description below. stat Status flag. 0 Stable subspace calculated. 1 Failure to decompose into stable and unstable subspaces. Minimum absolute value of the real part of the eigenvalues of H.
PAGE 370
390 Chapter 6. Function Reference variable, X = x2 x−1 1 , is the stabilizing solution to the Riccati equation. If H has jω axis eigenvalues then no stabilizing solution exists and the function returns a failure status. If any eigenvalue of H is within epp of the jω axis it is considered to lie on the jω axis and no solution is found. This may be due to numerical problems in finding the eigenvalues of poorly conditioned problems even when a stabilizing solution exists.
PAGE 371
Appendices A Translation Between Matlab µ-Tools and Xµ This appendix outlines the functional equivalences between the Matlab µ-Tools and Xmath Xµ. The objective is to provide a smooth transition for users moving from µ-Tools to Xµ. We will assume that the reader is familiar with Matlab µ-Tools and the general operation of Xmath.
PAGE 372
APPENDICES 392 The functionally equivalent commands will be listed, in each sub-section, for convenience. A more detailed discussion is given to illustrate the more subtle differences in the mu and D-K iteration aspects. More importantly, the Himat demo is available in both µ-Tools and Xµ. For a fast start on moving between platforms, study these demos side by side. A.
PAGE 373
A.
PAGE 374
APPENDICES 394 Subblocks: selecting input & outputs In µ-Tools the function sel selects rows and columns from a varying matrix or inputs and outputs from a system matrix. In Xmath these can be obtained by specifying row and column indexes. More flexibility of selecting parts of a pdm can be obtained by using the indexlist function. Augmentation Augmentation is the building of matrices from component pieces. The µ-Tools commands which perform these functions are given in the table below.
PAGE 375
A. TRANSLATION BETWEEN MATLAB µ-TOOLS AND Xµ 395 Note that the transpose and conjugate transpose operators are defined differently for Matlab and Xmath. system/Dynamic System Functions The following functions perform useful manipulations to, or information about, the state of a system or Dynamic System.
PAGE 376
APPENDICES 396 Description peak norm absolute value diagonal matrix round downwards round upwards imaginary part real part complex conjugate norm determinant eigenvalues inverse left division right division pseudo-inverse spectral radius singular values condition number schur decomposition matrix exponential interpolation decimation FFT inverse FFT spectral analysis µ-Tools Function pkvnorm vabs vdiag vfloor vceil vimag vreal vconj vnorm vdet veig vinv or minv vldiv vrdiv vpinv vrho vsvd vrcond vschur ve
PAGE 377
A. TRANSLATION BETWEEN MATLAB µ-TOOLS AND Xµ 397 Miscellaneous Utilities Several utilities are provided in µ-Tools. These are subroutines used by other µ-Tools functions which may be of more general use. Description complex random number fit system to data Eigenvalue based Riccati solution Schur based Riccati solution A.
PAGE 378
APPENDICES 398 A.3 System Response Functions Creating Time Domain Signals The Xmath pdm data object allows the creation of time domain signals via standard and operators. Description cosine waveform sine waveform stair-step waveform general waveform µ-Tools function cos tr sin tr step tr siggen Xmath/Xµ equivalent cos sin gstep standard functions & operators Time Responses Xmath calculates time responses with the * operation.
PAGE 379
A. TRANSLATION BETWEEN MATLAB µ-TOOLS AND Xµ 399 information. A.4 System Interconnection Simple interconnection has already been outlined in the augmentation section above. The more complicated interconnection functions are almost identical in µ-Tools and Xµ. The only real difference is the calling syntax of sysic. Description Redheffer (star) product system interconnection A.
PAGE 380
APPENDICES 400 Description H2 norm calculation H∞ norm calculation H2 controller synthesis H∞ controller synthesis µ-Tools function h2norm hinfnorm h2syn hinfsyn Xmath/Xµ equivalent h2norm hinfnorm h2syn hinfsyn The major syntactical difference is that the Xµ functions do not return the closed loop system. This is easily calculated by a subsequent call to starp. The reason for this is that the D-K iteration changes typically involve a differently weighted closed loop system in subsequent operations.
PAGE 381
A. TRANSLATION BETWEEN MATLAB µ-TOOLS AND Xµ Description structured singular value D scale decoding perturbation decoding block norm calculations rational perturbation random perturbations µ-Tools function mu unwrapd unwrapp blknorm dypert randel 401 Xmath/Xµ equivalent mu not required not required blknorm mkpert randpert In the Xµ mu function, only the default options of the µ-Tools mu calculations are available.
PAGE 382
APPENDICES 402 g1g = freq(g1,omega) [mubnds,Dmagdata] = mu(g1g,blk) From this, frequency domain D-scales are obtained, and a rational approximation is obtained via musynfit. [Dsys,Dinvsys] = musynfit(Dmagdata,blk,nmeas,nctrls,weight,g1g) A difference between the Xµ and µ-Tools implementations of musynfit is that, in the Xµ case, the inverse scaling system, Dinvsys, is generated by musynfit. Notice also that the previous D-scale fit is not passed as an argument.
PAGE 383
A. TRANSLATION BETWEEN MATLAB µ-TOOLS AND Xµ 403 The advantage of this is that in order to restart, or reproduce, an iteration, one need only save the previous controller. The µ-Tools approach requires saving the rational approximation to the previous D-scales. The controller is a more applicable data object to save and the saving of the previous D-scales depends on the quality of the rational approximation.
PAGE 384
Technical Support and Professional Services Visit the following sections of the National Instruments Web site at ni.com for technical support and professional services: • Support—Online technical support resources at ni.