National Instruments 370760B-01 manual

MA TRIXx TM Xmath TM X µ Manual MA TRIXx Xmath Basics The MA TRIXx products and related items ha ve been purchased from W ind Riv er S ystems, Inc. (formerly Integrated Systems, Inc.). Thes e refo rmatted user materials may contain references to th ose entities. Any trademark or copyright notices to those entities are no longer valid and any refer[...]

Support Worldwide Technical Support and Product Info rmation ni.com National Instruments Corporate Headquarters 11500 North Mopac Expressway Aust in, Texas 78759-3504 USA Tel: 512 683 0100 Worldwide Offices Australia 1800 300 800, Austria 43 0 662 45 79 90 0, Belgium 32 0 2 757 00 20, Brazil 55 11 3262 3599, Canada (Calgary) 403 274 9391, Ca nada ([...]

Important Information Warranty The media on which you receive National In struments software are warranted not t o fail to execute p rogramming instruction s, due to defects in materials and workmanship, for a period of 90 days from date of shipment, as eviden ced by receipts or other documentation. N ational Instruments will , at its option, repai[...]

Con ten ts 1 In tro duction 1 1 . 1 N o t a t i o n ....... ....... .......... ......... .... 1 1.2 Manual Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 How to avoid really r eading this Manual . . . . . . . . . . . . . . . . . . . 3 2 Ov ervi ew o f the Underl ying Theory 5 2 . 1 I n t r o d u c t i o n ..... ....[...]

iv CONTENTS 2.2.5 Obtaining Robust Control Mo dels for Ph ysica l Systems . . . . . . 28 2.3 H ∞ and H 2 D e s i g n M e t h o d o l o g i e s ........... ......... .. 2 9 2.3.1 H ∞ D e s i g n O v e r v i e w ..... ........... ......... 3 1 2.3.2 Assump tions for the H ∞ D e s i g n P r o b l e m .......... .... 3 2 2.3.3 A Brief Review of t[...]

CONTENTS v 2 . 6 M o d e l R e d u c t i o n .. ....... .......... ......... .... 6 4 2.6.1 T runcation and Residualization . . . . . . . . . . . . . . . . . . . . 65 2 . 6 . 2 B a l a n c e d T r u n c a t i o n ... .......... ......... .... 6 5 2 . 6 . 3 H a n k e l N o r m A p p r o x i m a t i o n .......... ......... .. 6 8 3 F unctional Descr[...]

vi CONTENTS 3 . 5 S y s t e m I n t e r c o n n e c t i o n ..... .......... ......... .... 9 1 3.6 H 2 and H ∞ A n a l y s i s a n d S y n t h e s i s ......... ......... .... 9 5 3 . 6 . 1 C o n t r o l l e r S y n t h e s i s ...... ......... ......... .. 9 5 3 . 6 . 2 S y s t e m N o r m C a l c u l a t i o n s ......... ......... .... 1 0 5 [...]

CONTENTS vii 4.1.5 µ An alysis of the H ∞ C o n t r o l l e r ........ ......... .. 1 3 8 4.1.6 Fit ting D -scales for the D - K I t e r a t i o n .......... ...... 1 4 0 4 . 1 . 7 D e s i g n I t e r a t i o n # 2 ... .......... ......... .... 1 4 3 4.1.8 Sim ula tion Comparison with a Lo opshaping Co nt r oller . . . . . . . 146 4.2 A Simple F[...]

viii CONTENTS A . 3 S y s t e m R e s p o n s e F u n c t i o n s .. ......... ....... .... 3 9 8 A . 4 S y s t e m I n t e r c o n n e c t i o n ...... ......... ......... 3 9 9 A . 5 M o d e l R e d u c t i o n ..... .......... ........... .. 3 9 9 A.6 H 2 and H ∞ A n a l y s i s a n d S y n t h e s i s ...... ........... 3 9 9 A.7 Str uctured [...]

Chapter 1 In tro duction X µ is a suite of Xmath functions for the mo deling , ana lysis and syn thesis of linear robust co n trol systems. Robust con tro l theor y has developed rapidly during the last decade to the point where a useful set of computational tools can b e used to solve a wide range of control problems. This theor y has already b e[...]

2 CHAPTER 1. INTR ODUCTION Notatio n Meaning pdm Xmath parameter dependen t matrix data ob ject Dynamic System Xmath dynam ic system data ob ject Co de examples and function names are set in t yp ewriter font to distinguish them from narrativ e text. 1.2 Man ual Outline Chapter 2 outlines the applica ble ro bust control theor y . Perturbation mo de[...]

1.3. HO W TO A VOID RE ALL Y READING THIS MANUAL 3 1.3 Ho w to a v oid really reading this M an ual The layout of the man ual pro ceeds from in tro duction to background t o syn tax detail to application descriptio ns. This may be tediously theoretical for some. If you are one of those that considers reading the manual as the o ption of last resort[...]

Chapter 2 Ov erview of the Underlying Theory 2.1 In tro duction The material cov ered here is taken from a v ariety of sources. The basic approach is described b y Do yle [1, 2], and further elaborated upon b y Pack ard [3]. Summaries ha ve also appeared in work by Smi th [4] and others. Motiv ating background can b e found in the ear ly paper by D[...]

6 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y several studies in volving pro cess control applications, par ticularly high pur it y distillation columns . These are detailed b y Skogestad an d Morari in [15, 16, 17, 18 ] Section 2 .2 introduces robust control p erturba tion mo dels and linear fr actiona l transform ations. W eight ed H ∞ d[...]

2.1. INTRODUCTION 7 Figure 2.1: The generic rob ust cont rol model stru cture T race( M ) trace of M ( P n i =1 M ii ) Block diagrams will be used to represent in terconnections of systems. Consider the example feedback in ter connection shown in Fig . 2.1. Notice that P has b een partitioned in to four parts. This diagram represen ts the equ ation[...]

8 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y is interpreted to mean that the signa l y is the sum o f the resp o nse of system P 21 to input sig nal v and system P 22 to input signal u . In general, w e will not b e specific a b out the represen tation of the system P . If we do need to be more specific ab out P ,t h e n P ( s ) is the La[...]

2.1. INTRODUCTION 9 Euclidean norm. Given, x =    x 1 . . . x n    , the E uclidean (o r 2 -norm) of x , denoted by k x k , is defined by , k x k = n X i =1 | x i | ! 1 / 2 . Man y other norms are also opti ons; more detail on the easily calculated n orms can b e found in the on-line help for the norm function. The term spatial-n[...]

10 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y F or p ersistent signals, where the ab ov e norm is un b ounded, w e ca n define a p ow er nor m, k x i ( t ) k = lim T →∞ 1 2 T Z T − T | x i ( t ) | 2 dt ! 1 / 2 . (2.1) The ab ov e norms hav e been defined in terms of a single comp onent, x i ( t ), of a v ector v alued s ignal, x ( t[...]

2.1. INTRODUCTION 11 signals. Strictly sp eaking, sig nals in H 2 or H ⊥ 2 are not defined on the ω axis. Ho wev er w e usually consi der the m to be b y taking a lim it as we approac h the axis. A slig ht ly mo re sp ecialized set is RL 2 , the set of real rational functions in L 2 . These are strictl y prop er functi ons with no poles on t[...]

12 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y where σ max denotes the maxim um sing ular v alue. No t all matrix no rms are induced from vector norms. The F r o eb enius nor m (square ro ot o f the sum of the squares of all matrix elements) is one suc h example. Now consider the case where P ( s ) is a dynamic system and we define an indu[...]

2.2. MOD ELING UNCER T AIN SYSTEM S 13 The set o f all s ystems with b ounded ∞ -nor m is denoted by L ∞ . W e can again split this into stable and unstable parts. H ∞ denotes the sta ble part; those systems with | P ( s ) | finite for all Re( s ) > 0. This is where the name “ H ∞ control t heory” originates, and we often call t his[...]

14 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y robust control model is therefo re a set descriptio n and we hop e that so me members of this set capture some o f the uncertain or unmo deled as pects o f our physical system. F or exa mple, consider the “uncertain” mo del illustrated in Figure 2.2 . This picture is equiv alen t to the inpu[...]

2.2. MOD ELING UNCER T AIN SYSTEM S 15 as sp ecifying a maximum percen tage error b etw een P nom and eve ry other elem ent of P . The system P nom ( s ) is the elemen t of P that comes fro m ∆ = 0 and is called the nominal system. In otherwords, for ∆ = 0, the input-o utput r elationship is y ( s )= P nom ( s ) u ( s ). As ∆ deviates fro m z[...]

16 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y Real 0 0.5 1 -0.5 1.5 Imaginary -0.5 0 0.5 -1 1 Figure 2.3: Nyquist diagram of the set of systems, P[...]

2.2. MOD ELING UNCER T AIN SYSTEM S 17 Figure 2.4: Unit y gain negativ e feedbac k for the example syst em, P 0 +∆ W a 2.2.2 Linear F ractional T ransformations A mo del is considered to be an interconnection of lumped comp onents and perturbatio n blo cks. In this discussion we will denote the input to the mo del by u , which can be a vector v a[...]

18 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y Figure 2.5 : Generic LFT mo del s tructure including p erturbatio ns,∆ A generic model s tructure, r eferred to a s a linear fractiona l transfor mation (LFT), ov ercomes the difficulties outlined above. The LFT mo del is equiv alen t to the relationship, y =  P 21 ∆( I − P 11 ∆) − [...]

2.2. MOD ELING UNCER T AIN SYSTEM S 19 in an LFT format. The open-lo op system is described by , y = F u ( P olp , ∆) u, where P olp =  0 W a IP 0  . The unity ga in, nega tive feedbac k config uration, illustr ated in Figure 2.4 (and given in Equation 2.5) can be described by , y = F u ( G cl p , ∆) r, where G clp =  − W a ( I + P [...]

20 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y The issue o f the inv er tibilit y o f ( I − P 11 ∆) is fundamental to the study of the stability of a system under per turbations. W e will r eturn to this question in muc h more deta il in Section 2.4. It for ms the basis of the µ analysis app roac h. Note that Eq uation 2.7 indicates tha[...]

2.2. MOD ELING UNCER T AIN SYSTEM S 21 Figure 2.6 : Example mo del: multiplicativ e output p er turbation with weighted output noise References to a robust co ntrol model will imply a description of the fo rm given in Equation 2. 8. As a exa mple, consider one of the most common p erturbation mo del descriptions, illustrated in Figure 2.6. This mo [...]

22 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y Robust control models a re therefor e set descriptions. In the a nalysis of such models it is also assumed that the unknown inputs belo ng to some bo unded set. Several c hoices of set for the unknown signals can be made, leading to different mathematica l problems for the ana lysis. Unfo rtuna[...]

2.2. MOD ELING UNCER T AIN SYSTEM S 23 Pac k ar d [19] discuss the implications of this a ssumption on ro bust con tr ol theor y and we briefly touch upon this in Section 2.4.6 . The most common assumption is that ∆ is an unknown, norm-bounded, linea r time-inv ar iant system. Systems often do not fall neatly in to one of the usual choices o f ?[...]

24 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y systems. W e will now lo ok at other p ossible p erturba tion structures. F or more detail on these struct ures ( in the complex case) refer to P ack ard and Do yle [20]. Conside r a blo cks w hich are of the form scalar × iden tit y , where the scalar is unknown. In the following we will inclu[...]

2.2. MOD ELING UNCER T AIN SYSTEM S 25 = Cz − 1 ( I − z − 1 A ) − 1 B + D = F u ( P ss ,z − 1 I ) , where P ss is the r eal v alued matr ix, P ss =  AB CD  , and the scalar × identit y , z − 1 I , has dimension equal to the state dimension of P ( z ). This is now in the form of an LFT mo del with a sing le scalar × identit y ele[...]

26 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y where the transport dela y , T d , and th e comb ustion lag, T c , are appro ximately , T d = 252 v and T c = 202 v . F or the purp oses of our example we w ant to design an air /fuel ratio controller that w orks for all engi ne sp eeds in the range 2,000 to 6,000 rp m. W e will use a first ord[...]

2.2. MOD ELING UNCER T AIN SYSTEM S 27 Putting all the pieces together gives an engine mo del in the follo wing fractional form. P ( s )= F u ( P mod , ∆) , where P mod =        − 15 . 87 s 7 . 14( s − 19 . 8) s 2 141 . 5(1 + 1 . 006 s ) s ( s +1 ) 0 − 9 . 9 s 9 . 9 − 27 . 75 s − 0 . 9( s − 15 . 8) ( s − 19 . 8)[...]

28 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y and ∆ ∈ B∆ , with the structure defined as, ∆ =  diag( δ v I 2 , ∆ 1 )     δ v ∈R , ∆ 1 ∈C  . Note that this is an LFT with a rep eated real-v alued parameter, δ v ( | δ v |≤ 1), and a complex p erturba tion, ∆ 1 ( | ∆ 1 |≤ 1). Note that a s R⊂C ,[...]

2.3. H ∞ AND H 2 DESI GN METHODOLOGIES 29 An area of work, kno wn as iden tificat ion in H ∞ , lo oks at experimen tal iden tification techniques whic h minimize the worst case H ∞ error bet w een the ph ysical system and the mo del. The fo llowing w o rks address this issue: [27, 28, 29, 30, 31, 32, 33, 34, 35]. Applying the more standard,[...]

30 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y Figure 2.7: LFT confi guration for control ler synth esis, G ( s )= F l [ P ( s ) ,K ( s )] Note that the interconnection structure, P ( s ), given here, diff ers from that discussed in the previous sect ion. Here w e set up P ( s ) so that the input, w , is the unknown signals en tering our s[...]

2.3. H ∞ AND H 2 DESI GN METHODOLOGIES 31 extend these approac hes t o the case where P ( s ) i s replaced b y F u ( P ( s ) , ∆), ∆ ∈ B∆ . 2.3.1 H ∞ Design Ov erview Again, reca ll from Section 2.1.2, the H ∞ is norm of G ( s )i s , k G ( s ) k ∞ =s u p ω σ max [ G ( ω )] . The H ∞ norm is the induced L 2 to L 2 norm. Ther [...]

32 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y 2.3.2 Assumptions for the H ∞ Design Problem There are sev eral assumpt ions required in order to achiev e a w ell-p osed design problem. The DGKF paper giv es a state-space solut ion to the H ∞ design problem and w e will use a similar notation here. Consider th e op en lo op state-space re[...]

2.3. H ∞ AND H 2 DESI GN METHODOLOGIES 33 that the eff ect of all disturban ces, w , at ev ery frequency , can b e measured b y the controller. If either of these conditions are not met then the problem could be ill-p osed. It is p ossible to violate these conditions by using pure in tegrators as design weigh ts. While this could still give a me[...]

34 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y W e hav e partitioned the matr ix in to tw o n × n blocks, X 1 and X 2 .I f X 1 is inv ertible, then X = X 2 X − 1 1 , is the unique, stabilizing so lution to the ARE. The ability to form X doesn’t dep end on the pa rticular choice of X 1 and X 2 . Given a Hamiltonia n, H ,w es a yt h a t H[...]

2.3. H ∞ AND H 2 DESI GN METHODOLOGIES 35 Cho ose γ> 0 and fo rm the following Ha miltonian matr ix, H =  Aγ − 2 BB T − C T C − A T  . The following lemma gives a means o f chec king whether or not k P ( s ) k ∞ <γ .A p r o o f o f this lemma is given in DGKF a lthough it is based on the w o rk of Ander son [60], Willems [6[...]

36 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y an y with a zero real part. In practice w e must use a tolerance t o determine what is considered as a zero real part. Finding a ba sis for the stable subspace of H in volv es either an eigenv alue o r Sch ur decomposition. Numerical errors will be in tro duced at thi s stage. In most cases usin[...]

2.3. H ∞ AND H 2 DESI GN METHODOLOGIES 37 and the component s, C 1 x and D 12 u are orthogonal . D 12 is also assumed to be normalized. This essentially means that there is no cro ss-weigh ting betw een the sta te and input penalties. Assumption ( iv ) is the dual of this; the input a nd unknown input (distu rbance and noise) affect the measurem[...]

38 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y L ∞ = − Y ∞ C T 2 Z ∞ =( I − γ − 2 Y ∞ X ∞ ) − 1 ˆ A ∞ = A + γ − 2 B 1 B T 1 X ∞ + B 2 F ∞ + Z ∞ L ∞ C 2 . Actually , the ab ov e formulat ion can b e used to parametrize all stabilizing con trollers which satisfy , k G ( s ) k ∞ <γ . This can b e expressed[...]

2.3. H ∞ AND H 2 DESI GN METHODOLOGIES 39 a) Choose γ ≥ γ opt b) F orm H ∞ and J ∞ c) Check that H ∞ ∈ dom(Ric) and J ∞ ∈ dom(Ric). d) Calculate X ∞ = Ric( H ∞ )a n d Y ∞ = Ric( J ∞ ) e) Check that X ∞ ≥ 0a n d Y ∞ ≥ 0 f ) Check that ρ ( X ∞ Y ∞ ) <γ 2 g) Reduce γ a n dg ot os t e pb ) . The v alue of γ [...]

40 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y eigen v alue of X ∞ (and Y ∞ ) is displayed. The ultimate test of the softw a re is to form the closed lo op system a nd chec k both its stability and nor m. W e str ongly sug gest that the user alwa ys perform this step . The numerical issues discussed ab ov e ar e very unlikely to ar ise i[...]

2.3. H ∞ AND H 2 DESI GN METHODOLOGIES 41 framew ork. W e again assume the simplif ying assumptions used in Sect ion 2 .3.4 The H 2 design solution is obtained (at lea st conceptually) from the H ∞ design pro cedure by setting γ = ∞ and using the resulting central controller. It is interesting to compa re the H ∞ solution, given ab ov e, a[...]

42 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y Figure 2.8: Closed loop system, G ( s ), for performance analysis initial choice of γ fo r the H ∞ design pr o cedure. W e will see later (Section 2 .5) that it can also b e used to initialize the D - K iter ation pro cedur e when an o p en-lo op H ∞ design is po orly conditioned. 2.4 µ An[...]

2.4. µ ANAL YSIS 43 descriptions are considered, where B again denotes the unit ball. Power : BP = ( w     lim T →∞ 1 2 T Z T − T | w ( t ) | 2 dt ≤ 1 ) (2.12) Energy : B L 2 = ( w     k w k 2 2 = Z ∞ −∞ | w ( t ) | 2 dt ≤ 1 ) (2. 13) Magnitude : B L ∞ = ( w     k w k ∞ = ess sup t | w ( t ) [...]

44 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y Figure 2.9: Perturb ed closed lo op s ystem for stability analysis the induced norms, the reader is referred to Do y le [2]. The ma jor adv ant a ge of c ho o sing BP or B L 2 is that t he test for th e p erformance can b e considered in terms of the same norm as stability . This has sig nifica[...]

2.4. µ ANAL YSIS 45 Consider the case where the model ha s only one full ∆ block ( m =1 a n d q =0 i n Equation 2.9 ). This is often r eferred to as unstructured, a nd the well kno wn r esult (refer to Zames [67] and Doyle and Stein [5 ]) is g iven in the following lemma. Lemma 6 (Robust St abilit y , Unstruct ured) F u ( G ( s ) , ∆) is stabl[...]

46 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y F u ( G ( s ) , ∆) stable fo r al l ∆ ∈ B∆ if and on ly if k µ ( G 11 ( s )) k ∞ < 1 . wher e k µ ( G 11 ( s )) k ∞ =s u p ω µ [ G 11 ( ω )] . The use of this notation masks the fact that µ is al so a function of the p erturbati o n structure, ∆ . The ab ove definition[...]

2.4. µ ANAL YSIS 47 if and on ly if k µ ( G ( s )) k ∞ < 1 , wher e µ is taken wi th r esp e ct to an a ugmente d structur e b ∆ , b ∆ = n diag(∆ , ˆ ∆)    ∆ ∈ ∆ , ˆ ∆= C dim( w ) × dim ( e ) o . The additio nal p erturba tion blo ck, ˆ ∆ can b e thou ght of as a “performance block” app ended to the ∆ bl[...]

48 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y F or the o ther extreme consider a single full blo c k ( ∆ = { ∆ | ∆ ∈C n × n } ); the definition of µ is now the same as tha t for the maximum singular v alue, ∆ = { ∆ | ∆ ∈C n × n }⇒ µ ( M )= σ max ( M ) . Observe that every po ssible blo ck structure, ∆ , contains { [...]

2.4. µ ANAL YSIS 49 Actually , the lower bound is always equal to µ but the implied optimization ha s lo cal maxima which are not glo bal. F or the upper b ound Sa fonov and Doyle [72], hav e shown that finding the infim um is a conv ex problem and hence more easily solved. Howev er the bo und is equal to µ only in certa in sp ecial cases. Her[...]

50 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y and tw o blo ck structures, ∆ 1 (compa tible with M 11 )a n d ∆ 2 (compatible with M 22 ). There are t wo perturbed subsystems th at we can stud y here: F u ( M, ∆ 1 ), where ∆ 1 is closed in a feedbac k lo op aro und M 11 ;a n d F l ( M, ∆ 2 ), where ∆ 2 is clos ed in a feedback loo[...]

2.4. µ ANAL YSIS 51 2.4.6 State-sp ace Robustness Analysis T ests W e will lo ok at some more adv anced state-space approac hes to the analysis of robust performance. Most users of the soft ware wi ll concentrat e on the more common frequency domain analysis metho ds cov ered in Section 2 .4.3. The a nalysis tests given here can b e implement ed w[...]

52 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y Therefore µ 1 ( A ) < 1 is equiv alen t to our system b eing stable. F urthermo re, the maximum modulus theorem for a stable system tells us that, k P ( z ) k ∞ =s u p | z |≥ 1 σ max ( P ( z )) =s u p | z − 1 |≤ 1 σ max ( F u ( P ss ,z − 1 I ) =s u p ∆ 1 ∈ B∆ 1 µ 2 ( F u ([...]

2.4. µ ANAL YSIS 53 Figure 2.10: P erturb ed system for st ate-space robust ness tests[...]

54 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y Note that the nominal system is given by , G nom ( z )= F u  AB 1 C 1 D 11  ,z − 1 I  , and the perturbed system is, G ( z )= F u ( F l ( G, ∆) ,z − 1 I ) . W e a ssume that ∆ is an element of a unity norm b o unded blo ck structure, ∆ ∈ B∆ . Fo r t h e µ analysis w e [...]

2.4. µ ANAL YSIS 55 iii ) There exists a c onstant β ∈ [0 , 1] s uch that for e ach fixe d ∆ ∈ B∆ , G ( z ) is stable and for zero initial st ate r e sp ons e, e satisfies k e k 2 ≤ β k w k 2 (r obust p erformanc e). The frequency domain µ test is implemented by sea rching for the maximum v alue o f µ ov er a user specified freque[...]

56 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y i ) State-spa ce upp er bound: inf D s ∈D s σ max [ D s G ss D − 1 s ] < 1; ii ) F re quency domain, constant D , upper b ound: inf D p ∈D p max ω ∈ [0 , 2 π ] σ max [ D p F u ( G ss , e ω I nx ) D − 1 p ] < 1; iii ) F re quency domain upp er bo und: max ω ∈ [0 , 2 π [...]

2.4. µ ANAL YSIS 57 The gap betw een the state- space (or constant D ) upp er b o und and the frequency domain upp er bo und is mor e significant. In the sta te-space upper b ound, a single D scale is selected . This giv es robust performance for all ∆ satisfy ing, k v k≤k z k for all e ∈L 2 . This can b e satisfied for linea r time-v a ry[...]

58 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y 2.4.7 Analysis with b oth Real and Complex Perturbations The ab ove results o nly apply to the case where ∆ is considered as a co nstant complex v alued matrix a t each frequency . In many engineering applications restr icting certain of the ∆ blo c ks to b e real v alued may result in a les[...]

2.5. µ SYNTHESIS AND D - K ITERA TION 59 Figure 2.11: The generic in terconnecti o n structu re for synthesis k F u ( F l ( P ( s ) ,K ( s )) , ∆) k ∞ ≤ 1. This is equiv alen t to K ( s ) satisfying µ [ F l ( P ( s ) ,K ( s ))] < 1. In other words, the closed lo op system sa tisfies the ro bust performance specification. Unfortunately [...]

60 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y Recall that this is an upp er bo und for the µ pr oblem of in terest, implying that, µ [ F l ( P ( s ) ,K ( s )] ≤ 1 , as required. How ever the upp er b ound may be co nserv ative, meaning that in o rder to guaran tee that µ [ F l ( P ( s ) ,K ( s )] ≤ 1, we ha ve h ad to back off on th[...]

2.5. µ SYNTHESIS AND D - K ITERA TION 61 dynamic system. This requires fitting an appro x imation to the upper b ound D -scale in the itera tion. W e will now look a t this issue mor e closely . The D - K iteration pro cedur e is illustr ated schematically in Figure 2.12 . It can b e summarized as f ollo ws: i ) Initia lize pro cedure with K 0 ( [...]

62 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y Figure 2. 12: D - K iterat ion pro cedure: a) Design H ∞ ( or other) con trolle r: K 0 ( s )[ s t e p i )]. b) Closed lo op p erturbed system for µ analysi s [step ii )]. c) F requency b y frequency upper bo und D ( ω ) s cale a pproximation to µ analysis [step iii )]. d) Scaling of H ∞ d[...]

2.5. µ SYNTHESIS AND D - K ITERA TION 63 frequency we w ould hav e, D =      d 1 I 1 . . . d m I m I e      , where the identit y I e is of dimensions dim( e ) × dim( e ). The calculation of a new H ∞ con tro ller requires a state-space rea lization of D ( ω ). F or eac h d i in D ( ω ) we m ust fit a transfer f[...]

64 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y Sev eral asp ects of t his pro cedure are worth noti ng. F or the µ analysis and D scale calculation, a fr equency grid must b e chosen. The range and reso lution of this g rid is a matter o f engineering judgemen t. The µ analysis ca n require a fine grid in the vicinity of the lightly dampe[...]

2.6. MODEL REDUCTION 65 2.6.1 T runcation and Residualization The simplest fo rm of mo del r eduction is state truncation. Co nsider a system, P ( s ), with a partitio ned state matrix, P ( s )=   A 11 A 12 B 1 A 21 A 22 B 2 C 1 C 2 D   . T runca ting the states a sso cia ted with A 22 re sults in, P tr un ( s )=  A 11 B 1 C 1 D  [...]

66 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y The controllabilit y grammian, Y is defined a s, Y = Z ∞ 0 e At BB T e A T t dt, and the observ abilit y grammian, X , is defined as X = Z ∞ 0 e A T t C T C e At dt. The grammians, X and Y , satisfy the Lyapunov equations, AY + YA T + BB T =0 A T X + XA + C T C =0 , and this is t ypically [...]

2.6. MODEL REDUCTION 67 W e will now lo ok at a particular choice of transfor mation. F or a minimal r ealizatio n, w e can always find a transformation that gives, ˆ Y = TY T T =Σ , and ˆ X = T − T XT − 1= Σ , where Σ = diag( σ 1 ,...,σ n )a n d σ i ≥ 0, i =1 ,...,n . This r ealization, wher e the grammians ar e equal, is called a b[...]

68 CHAPTER 2. OVER VIE W OF THE UNDE RL YING THEO R Y and Glov er [79] independently obtained the following b ound on the error induced by balanced trunca tion. Theorem 13 Given a st able, r at ional, P ( s ) ,a n d P bal ( s ) , the b alanc e d tr uncation of or der k< n . Then, k P ( s ) − P bal ( s ) k ∞ ≤ 2 n X i = k +1 σ i and k P ( [...]

2.6. MODEL REDUCTION 69 Consider the problem of finding the stable, order k realization w hich minimizes the Hank el norm of the error. Defi ne, P hankel ( s ) a s the minimizing system. Then we ha ve, σ k +1 ≤k P ( s ) − P hankel ( s ) k H = inf P k ( s ) stable k P ( s ) − P k ( s ) k H . This system also satisfies ∞ -norm b ounds on [...]

Chapter 3 F unctional Description of X µ 3.1 In tro duction This chapter describes the X µ functions in the co n text of their in tended usage. Chapter 2 provides the reader with an idea of the theo retical basis b ehind the v arious analysis and design ca lculations. Here w e outline the softw are functio ns a v aila ble for doing those ca lcula[...]

72 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ 3.2.1 D ynamic System s Xmath has a dynamic system data ob ject 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 conti nuous case, and x ( kT + T )= Ax ( kT )+ Bu ( kT ) ,[...]

3.2. D A T A OBJECTS 73 As above, these p olynomials can be sp ecified b y their ro ots or their coefficients. Note that w e can sp ecify the v ariable, and for con tin uous systems w e use “s”. T o create a discrete system “z” is used. # Generate the system from the numerator and denominator # coefficients. numerator = makepoly([-.1,19.97[...]

74 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ Because the dynamic system is a built-in data ob ject, the lab el info rmation will be display ed by s imply typing the ob ject name ( sys in the ab ov e) o r app ending a question mark to the statement. The core function commentof is used to read the commen t attached to a sp ecified v ariable. 3.2.2 p[...]

3.2. D A T A OBJECTS 75 App ending and Merging Data Time functions, for cre ating simulation inputs for example, can b e cr eated by com bining pdm s. Xmath has t wo suc h core functi ons: concatseg and insertSeg The concatseg app ends the data of one pdm to a nother. The domain is recalcula ted suc h that it is alwa ys regula r. The user can sp ec[...]

76 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ Data ca n also b e extracted b y indep endent v ariable index num b er , by providing a scalar argumen t to the pdm . In the following example the fifth through tenth and the t wen tieth indep endent v ariables are specified for the smaller pdm , pdm2 . size(pdm1) ans (a row vector) = 1 1 100 pdm2 = pd[...]

3.2. D A T A OBJECTS 77 # index of the pdm the value 100. idxlst = indexlist([4,1,2] ) pdm1(idxlst) = 100 Op erations on the Indep endent V ariables/Domain The domain of a pdm is readily changed via the pdm co mmand. The fo llowing example illustrates a co mmon applica tion; changing a frequency do main in Her tz to radians/s econd. # The following[...]

78 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ # Select columns 1, 3 ,4&5a n d rows 2 & 7 from the # pdm: bigpdm. subpdm = bigpdm([1,3:5],[2 ,7]) This referen cing format can also be used to assign data to parts of a larger pdm . This is shown in the fo llowing example. # Replace the 3,2 block of a pdm with its absolute value pdm1(3,2) = abs([...]

3.2. D A T A OBJECTS 79 and B . Augmentation for pdm s simp ly p erforms the augmen tation at eac h domain v ariable. The domains m ust b e the same. Diagonal augmen tation can be p erformed with the X µ functio n daug . This is the equiv alent of the matrix augmen tation: [A, 0; 0, B] , except t hat up to 20 argumen ts can b e aug ment ed in one [...]

80 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ poles, or a zero D term. Generating rando m systems is useful for simulating systems with unknown subsystems. Specialized X µ functions a re provided for useful manipulations o f the state. F or example transform ing the A matr ix to a r eal v alued, 2 × 2 block diagonal form; here referred t o as mo d[...]

3.3. MA TRIX INFORMA TION, DISPLA Y AND PLOTTING 81 # N is the decimation ratio. smallpdm = bigpdm([1:N:length(bigpdm )]) 3.2.5 Con tinuous to Discrete T ransformations Xmath has a single function, discretize , fo r calculating contin uo us to discrete transforma tions. Several o ptions a re o ffered including forward and ba ckw ard differ ence, [...]

82 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ 3.3.2 F ormatted Displa y F unctions It is often useful to consider complex num b ers in terms of frequency and damping. This is particular ly applicable when studying the p oles o r zeros of a system. The comma nd rifd provides a real-imaginary-fr equency-damping formatted display for complex v alued in[...]

3.3. MA TRIX INFORMA TION, DISPLA Y AND PLOTTING 83 g1 = ctrlplot(sys1g, { bode } ); g1 = ctrlplot(sys2g,g1, { bo de } ); g1 = plot( { keep=g1,title = "Bode plots",... legend = ["sys1","sys2"] } )? Frequency 0.1 1 0.01 10 Magnitude 0.0001 0.001 0.01 0.1 1 1e-05 10 Bode plots sys1 sys2 Frequency 0.1 1 0.01 10 Phase (deg[...]

84 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ # Nyquist plots g2 = ctrlplot(sys1g, { nyqui st } ); g2 = ctrlplot(sys2g,g2, { ny quist } ); g2 = ctrlplot(-1,g2, { nyqui st,marker= 1,line=0 } ); g2 = plot(g2, { projection=" orthograph ic",... legend=["sys1","sys3" ,"critical point"],... title="Nyquist plots[...]

3.4. SYSTEM RESP ONSE FUNCTIONS 85 3.4 System Resp onse F unctions 3.4.1 Creating Time Domain Signals The Signal Analysis Mo dule cont a ins several functions which are useful for building time domain signals: gcos and gsin .X µ provides gstep for the creation of stai r-step signals. The example b elow illustrates the g eneration of a sine wav e a[...]

86 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ spacing in the input signa l pdm . This mea ns that the input pdm must be regularly spaced. Stair-cas e input functions with relatively few data po ints will often give erroneous results. The input signal should b e interpolated (with interp ) before being multiplied with t he D ynamic System . The funct[...]

3.4. SYSTEM RESP ONSE 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 = ctrlplot(u, { line style=2 } ); g1 = ctrlplot(y0,g1, { line style=1 } ); g1 = ctrlplot(y1,g1, { line style=4 } ); g1 = plot(g1, { !grid } )? 2468 0 10 -2 -1 0 1 -3 2 The native Xmath time domain simulation, us[...]

88 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ and intstep are the in terp olation order an d integrat ion step size. The system m ust b e con tinu ous. The int egration order and sample t ime are prespecified for discret e tim e systems making the * op erato r is suita ble for such sim ulatio ns. X µ provides a sampled data sim ula tion function ([...]

3.4. SYSTEM RESP ONSE FUNCTIONS 89 maximum. This can b e ha ndy when first examining a high order system. An example is giv en b elow. # Create a single-input, two-output system. sys = [1/makepoly([1,0.1,1 ],"s");... makepoly([1,1],"s")/m akepoly([1 ,10],"s") ] # The automatic point selection is used. sysg = freq(sys,[...]

90 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ Frequency 0.1 1 10 0.01 100 Magnitude 1e-05 0.0001 0.001 0.01 0.1 1 10 1e-06 100 Frequency 0.1 1 10 0.01 100 Phase (degrees) -150 -100 -50 0 50 -200 100 Note that the defa ult frequency units in Xmath are Hertz. This a pplies generically to all functions where the user specifies frequency informa tion, [...]

3.5. SYSTEM INTERCONNECTION 91 Figure 3.2: Generic Red heffer in terconnect ion structure 3.5 System In terconnection In terconnections of systems are used extensively in the design, analysis and sim ula tion calculati ons. The most general form of in terconnection , and the one used in formin g closed lo op systems, is the Redheffer (or star) pr[...]

92 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ Figure 3.3: Example in terconnecti on of subsystem s Figure 3.4: Example in terconnected syst em[...]

3.5. SYSTEM INTERCONNECTION 93 Using sysic to fo rm this interconnection can be considered a s four distinct o p erations. • Specify the individual subsystems. • Name and dimension the input signals. • Specify , algebra ically in terms of subsystem outputs or input signals, the output of the interconnected system. • Specify , algebra ically[...]

94 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ Specify the o utputs in terms of the subsystem names and the input names. Note that individual o utputs of MIMO systems and s imple ar ithmetic combinations can b e specified. Paren thesis sp ecify which output of a MIMO system is to be used. outputs = ["p(1) + wght"; "7.5*c"; "[...]

3.6. H 2 AND H ∞ ANAL YSIS AND SYNTHESIS 95 (rad/sec) ratio -1.0000e+00 0.0000e+00 1.0000e+00 1.0000 -5.0000e-01 6.3048e+00 6.3246e+00 0.0791 -5.0000e-01 -6.3048e+00 6.3246e+00 0.0791 -1.0000e+01 0.0000e+00 1.0000e+01 1.0000 The order of the names in th e systemname v ariable, m ust matc h the order of the ro ws in the connections v ariabl e and [...]

96 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ Figure 3.5: In terconnection structu re for control ler syn thesis structure contains more than just the op en lo op plant. It t ypica lly also con ta ins frequency dependent w eighting functions and specifies the structur e of t he in terconnection b etw een the ope n lo op pla nt and the controller. T[...]

3.6. H 2 AND H ∞ ANAL YSIS AND SYNTHESIS 97 equation solutio n pro cedure o ften b ecomes po o rly conditio ned. Displaying in termedia te calculation results allows the user to fine tune several tolerances if necessary . The intermediate Hamiltonian a nd Riccati so lution details a re displayed as the bisection pro ceeds. The bisection s toppin[...]

98 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ Frequency 0.01 0.1 1 10 0.001 100 Magnitude 0.1 1 10 0.01 100 Open loop plant Frequency 0.01 0.1 1 10 0.001 100 Phase (degrees) -300 -200 -100 -400 0 The desired clos ed lo op configur ation is illustrated in Figure 3.6. Figure 3.6: Closed loop configuration[...]

3.6. H 2 AND H ∞ ANAL YSIS AND SYNTHESIS 99 Figure 3.7: W eight ed design in terconnect ion structure: p In order to set up the desig n problem, we consider ref as an unknown input a nd the trac king error (input to k ), and the actuator signal, u , as outputs to be minimized. These outputs are weigh ted with the weigh ts Wperf and Wact respectiv[...]

100 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ Frequency 0.01 0.1 1 10 0.001 100 Magnitude 0.01 0.1 1 10 0.001 100 Weighting functions Wperf Wact # Form the weighted interconnection structure sysnames = ["plant";"Wperf";"Wact"] sysinp = ["ref";"u"] sysout = ["Wperf"; "Wact"; "[...]

3.6. H 2 AND H ∞ ANAL YSIS AND SYNTHESIS 101 12.500 5.2e-01 1.7e-03 1.0e-02 0.0e+00 0.0000 p 6.250 5.2e-01 1.7e-03 1.0e-02 0.0e+00 0.0000 p 3.125 5.1e-01 1.7e-03 1.0e-02 0.0e+00 0.0000 p 1.562 5.0e-01 1.7e-03 1.0e-02 0.0e+00 0.0000 p 0.781 3.9e-01 -9.8e+01 1.0e-02 0.0e+00 0.0000 f 1.172 4.8e-01 1.8e-03 1.0e-02 0.0e+00 0.0000 p Gamma value achieve[...]

102 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ (rad/sec) ratio -1.4046e-01 -2.3161e-01 2.7087e-01 0.5186 -1.4046e-01 2.3161e-01 2.7087e-01 0.5186 -1.5863e+00 3.4754e+00 3.8203e+00 0.4152 -1.5863e+00 -3.4754e+00 3.8203e+00 0.4152 -5.2060e+00 0.0000e+00 5.2060e+00 1.0000 Zeros: real imaginary frequency damping (rad/sec) ratio -5.0000e-02 -3.1225e-01 3[...]

3.6. H 2 AND H ∞ ANAL YSIS AND SYNTHESIS 103 Frequency 0.01 0.1 1 10 0.001 100 Magnitude 0.001 0.01 0.1 1 0.0001 10 Controllers Kinf K2 # Examine sensitivity functions sensinf = inv(1 + plant*Kinf) sensinfg = freq(sensinf,omega) sens2 = inv(1 + plant*K2) sens2g = freq(sens2,omega) g2 = ctrlplot(sensinfg, { lo gmagplot } ); g2 = ctrlplot(sens2g,g2[...]

104 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ Frequency 0.01 0.1 1 10 0.001 100 Magnitude 0.1 1 0.01 10 Sensitivity functions Kinf K2 Note tha t with the interconnection structure shown in Figure 3 .7, the clo sed lo op transfer fun ction from ref to e (1) is simply Wperf*sensinf .I f t h e H ∞ of the r esulting closed loop system was less th a n[...]

3.6. H 2 AND H ∞ ANAL YSIS AND SYNTHESIS 105 g3 = ctrlplot(step,g3, { line style=2 } ); g3 = plot(g3, { title="Step responses",... legend=["Kinf";"K2";" input"] } )? 2468 0 10 0.2 0.4 0.6 0.8 1 0 1.2 Step responses Kinf K2 input 3.6.2 System Norm Calculations F unctions a re provided for calculating the H 2[...]

106 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ Bounds on the H ∞ norm are ret urned as out . An estimate of the frequency where the norm is achiev e d is returned a s omega . F ur ther con tro l of the itera tion is av aila ble via k eyw ords. The following example calcula tes the H 2 and H ∞ norms of eac h of the closed lo op systems ari sing f[...]

3.7. STRUCTURED SINGULAR V ALUE ( µ ) ANAL YSIS AND SYNTHE SIS 107 3.7 Structured Singular V alue ( µ ) Analysis and Syn thesis This section cov ers the functions used in the D - K i teration procedure. The primary funct ions are the cal culation of the con troller ( already dis cussed), the cal culation of µ and the fitting of r ational D scal[...]

108 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ The outputs of the µ functio n a re: the upper and low er bounds for µ ;t h e D matrix fo r the upp er b o und; the Q matr ix for the lower b o unds; and a sensitivity estimate for the part of the D matrix corresponding to eac h blo ck in ∆. The sensitivity esti mate is essent ia lly the gradient o [...]

3.7. STRUCTURED SINGULAR V ALUE ( µ ) ANAL YSIS AND SYNTHE SIS 109 c = 1/sqrt(gamma); d = -sqrt(beta/gamma) f = (1+jay)*sqrt(1/(gamma* beta)) psi1 = -pi/2; psi2 = pi U = [a,0; b,b; c,jay*c; d,f] V = [0,a; b,-b; c,-jay*c; f*exp(jay*psi1), d*exp(jay*psi2)] scl = diagonal(random(4,1) +0.1*ones( 4,1)) M = scl*U*V*’*inv(scl) Consider four 1 × 1 bloc[...]

110 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ 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 iter ation is used a s an approximation to µ syn thesis. This section discusses how X µ implemen ts this pro ced[...]

3.7. STRUCTURED SINGULAR V ALUE ( µ ) ANAL YSIS AND SYNTHE SIS 111 (a) Cont roller and closed loop are satisfactory so stop the iterati o n. (b) The iteration has con verged and the controll er and closed lo op are not satisfactory . In this case the weighted design problem m ust be reformulated. (c) The iter ation has not yet con verged. Con tinu[...]

112 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ Figure 3. 9: H ∞ controller design. Step 9 in the en umer ated pro cedure There is a ctually another p ossibility at step 5 ; numerical problems cause the itera tion to diverge. As γ approac hes its optimal v alue, the n umerical properties of the calcul ation deteri orate. This ma y lead to mu(G i) [...]

3.7. STRUCTURED SINGULAR V ALUE ( µ ) ANAL YSIS AND SYNTHE SIS 113 Both the D and D − 1 systems ( Dsys and Dinvsys ) are returned. The D scale ( Dmagdata ) comes from a µ ca lculation on a closed lo op system. Howev er, Dsys and Dinvsys are requir ed to m ultiply the op en lo op system. They must therefore contain the iden tity matrices f o r t[...]

114 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ creating weigh ts fro m data and simple system identification. X µ provides t wo user callable functions for fitting SISO tra nsfer functions to da ta. The first is mkphase calculates the phase corresp onding to a minim um-phase stable system from magnitude data. This uses the complex cepstr um meth[...]

3.7. STRUCTURED SINGULAR V ALUE ( µ ) ANAL YSIS AND SYNTHE SIS 115 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, { logm agplot } ); g1 = ctrlplot(Wmg,g1, { logm agplot,lin e sty[...]

116 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ Wperf = makepoly([0.01,1],"s")/ makepoly([ 1,0.01],"s ") Wact = 0.1* makepoly([1,1],"s")/make poly([0.05 ,1],"s") Wnoise = 0.01 Wref = makepoly([0.005,1],"s")/ makepoly([ 0.05,1],"s ") Wperfg = freq(Wperf,omega) Wactg = freq(Wact,omega) Wnoiseg[...]

3.7. STRUCTURED SINGULAR V ALUE ( µ ) ANAL YSIS AND SYNTHE SIS 117 nms = ["plant";"Wm";"Wperf ";"Wact";" Wnoise";"W ref"] inp = ["delt";"ref";"noise ";"u"] outp = ["100*Wm"; "Wperf"; "Wact"; "Wref-0.01*delt-plant-W[...]

118 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ blk = [1,1; 2,2] [rpbnds1,D1,Dinv1,Delt a1,sens1] = mu(Gg,blk) g3 = ctrlplot(np, { log } ); # plot on a log-linear scale g3 = ctrlplot(rs,g3, { log,l ine style=3 } ); g3 = ctrlplot(rpbnds1,g3, { log,line style=[4,5] } ); g3 = plot(g3, { !grid,title= "mu analysis",y lab="Magnitude",..[...]

3.7. STRUCTURED SINGULAR V ALUE ( µ ) ANAL YSIS AND SYNTHE SIS 119 Frequency (Hz) 0.1 1 10 0.01 100 Magnitude 0.1 1 10 0.01 100 D scale fit, block: 1 Magnitude data Previous fit, order: 2 New fit, order: 4 # Apply the D scales to another H infinity design Kmu = hinfsyn(Ds*p*Dinvs,n meas,ncntr ls,[0;10]) # Close the loop around the weighted interco[...]

120 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ 0.1 1 10 0.01 100 Magnitude 0.5 1 1.5 2 2.5 3 3.5 0 4 mu analysis of robust performance Kmu Kinf 3.7.4 Constructing Rational Perturbations F or sim ula tion pur p oses it is useful to be able to construct a rational a pproximation to the ∆ r eturned by the µ calculat ion. The approac h i s to c ho os[...]

3.8. MODEL REDUCTION 121 pert = randpert(blk, { sys,sfreq, complex,pn orm } ) The user ca n specify whether the p er turbation is a dynamic system o r matrix, and whether it is real or co mplex v alued, in addition to specifying the no rm. 3.7.5 Blo ck Structured Norm Calculations It is possible to get an idea of the input/o utput combinations whic[...]

122 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ A g reater r ange of mo del reductio n functio ns is av a ilable in the Mo del Reduction Mo dule. Some of the functions descr ib ed here are cro ss-licensed with that mo dule. Section 2 .6 describ es the theory b ehind these functions. 3.8.1 T runcation and Residualization T runcatio n is provided b y t[...]

3.8. MODEL REDUCTION 123 g1 = ctrlplot(sysout2g,g1, { logmagplo t,line style=4 } ); g1 = ctrlplot(residerror,g 1, { logmagp lot,line style=5 } ); g1 = ctrlplot(truncerror,g 1, { logmagp lot,line style=6 } ); g1 = plot(g1, { !grid,legend =["origina l system";... "residualized system";"truncated system";... "residua[...]

124 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ # 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. g2 = ctrlplot(sys1g, { logma gplot } ); g2 = ctrlplot(sysout3g,g[...]

3.8. MODEL REDUCTION 125 3.8.3 Hank el Singular V alue Approximation The function ophank (also cr oss-licensed from the Mo del Reductio n Mo dule) is used to perform optimal Hank el norm appro ximation. Recall f rom Section 2.6.3 that there is an astable system achieving the lower bound. The unstable part of this system is retur ned as the second a[...]

126 CHAPTER 3. FUNCTIONAL DESCRIPTION OF X µ Frequency 0.1 1 10 0.01 100 Magnitude 1e-05 0.0001 0.001 0.01 1e-06 0.1 original system Hankel norm approximation error: ophank error: balmoore error: sresidualize[...]

Chapter 4 Demonstration Examples 4.1 The Himat Example The following demo can be run by executing the f ollowing Xmath command: execute file = "$XMATH/demos/xMu/ himatdemo" 4.1.1 Problem Description The Hima t is a small scale remotely piloted a ircraft built to investigate high maneuverabilit y figh ter aircraft design. The v ehicle was[...]

128 CHAPTER 4. DEMONSTRA TION EXAMPLES δv Perturbations along the velo city v ector. α Angle of attack. I .e. angle betw een the veloc it y vector and the aircra ft’s longitudinal axis. q Rate-of -c hange of aircraft attitude angle. θ Aircraft at titude angle. Con trol can be exert ed via the elevon and canard, denot ed b y δ e and δ c resp [...]

4.1. TH E HIMA T EXAMPLE 129 Figure 4.2: In terconnection st ructure for the himat design exam ple The inputs are the elev o n pos ition and the canard po sition. The outputs that ar e to be k ept small are angle-of-attac k ( α ) and pitch angle ( θ ). The commands required to enter the state-space description are simply matrix assignmen ts for e[...]

130 CHAPTER 4. DEMONSTRA TION EXAMPLES himat (a state space system) = A -0.0226 -36.6 -18.9 -32.1 0 -1.9 0.983 0 0.0123 -11.7 -2.63 0 00 1 0 B 00 -0.414 0 -77.8 22.4 00 C 0 57.3 0 0 0 0 0 57.3 D 00 00 X0 0 0 0 0 State Names ----------- forward speed angle-of-attack pitch rate pitch angle Input Names ----------- elevon canard Output Names ----------[...]

4.1. TH E HIMA T EXAMPLE 131 angle-of-attack pitch angle System is continuous 4.1.3 Creating a W eigh ted In terconnection Structure for Design The multiplicativ e input p erturba tion weight, W del , is constructed as a tran sfer functi on. The weigh t is, W del = 50( s + 100) ( s + 10000) . The output error weigh t, W p , is created as, W p = 0 .[...]

132 CHAPTER 4. DEMONSTRA TION EXAMPLES Frequency 0.01 0.1 1 10 100 1000 0.001 10000 Magnitude 1 10 0.1 100 Weights for HIMAT Perturbation weight Performance weight The p erturbat ion w eight, W del , should actually be tw o-input, tw o -output. This is also true of the performance w eigh t, W p . In this example, w e are weigh ting eac h performanc[...]

4.1. TH E HIMA T 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,h imat[...]

134 CHAPTER 4. DEMONSTRA TION EXAMPLES gamma Hx eig X eig Hy eig Y eig nrho xy p/f 6.000 2.3e-02 5.6e-05 2.3e-02 0.0e+00 0.0626 p 3.400 2.3e-02 5.7e-05 2.3e-02 0.0e+00 0.2020 p 2.100 2.3e-02 5.9e-05 2.3e-02 0.0e+00 0.5798 p 1.450 2.3e-02 6.4e-05 2.3e-02 0.0e+00 1.4678 f 1.775 2.3e-02 6.1e-05 2.3e-02 0.0e+00 0.8652 p 1.613 2.3e-02 6.2e-05 2.3e-02 0.[...]

4.1. TH E HIMA T EXAMPLE 135 Zeros: real imaginary frequency damping (rad/sec) ratio -2.2516e-02 0.0000e+00 2.2516e-02 1.0000 -1.7226e+00 0.0000e+00 1.7226e+00 1.0000 -3.0272e+00 0.0000e+00 3.0272e+00 1.0000 -3.1034e+01 0.0000e+00 3.1034e+01 1.0000 -1.0000e+04 0.0000e+00 1.0000e+04 1.0000 -1.0000e+04 0.0000e+00 1.0000e+04 1.0000 Next, a magnitude p[...]

136 CHAPTER 4. DEMONSTRA TION EXAMPLES Frequency 1 10 100 1000 0.1 10000 Magnitude 0.01 0.1 0.001 1 Controller: k1 Frequency 1 10 100 1000 0.1 10000 Phase (degrees) -200 -100 0 100 -300 200 Onto the closed loop, first chec king that it is stable by loo king at the po le po sitions. rifd(g1) Poles: real imaginary frequency damping (rad/sec) ratio[...]

4.1. TH E HIMA T EXAMPLE 137 -2.2517e-02 0.0000e+00 2.2517e-02 1.0000 -2.2600e-02 0.0000e+00 2.2600e-02 1.0000 -3.0000e-02 0.0000e+00 3.0000e-02 1.0000 -3.0000e-02 0.0000e+00 3.0000e-02 1.0000 -2.9369e+00 0.0000e+00 2.9369e+00 1.0000 -2.9974e+00 0.0000e+00 2.9974e+00 1.0000 -4.8310e+00 0.0000e+00 4.8310e+00 1.0000 -6.5876e+00 0.0000e+00 6.5876e+00 [...]

138 CHAPTER 4. DEMONSTRA TION EXAMPLES 1 10 100 1000 0.1 10000 0.5 1 1.5 0 2 Singular value plot of the closed loop 4.1.5 µ A n a l y s i so ft h e H ∞ Con troller The H ∞ con trol law can b e analyzed using µ -analysis. The closed-lo o p system, g1 ,h a s4 inputs and 4 outputs. The first t wo inputs and outputs corresp ond to the uncert a i[...]

4.1. TH E HIMA T EXAMPLE 139 complex v alued blo cks. The upp er and lower bounds of the µ function are return ed in bnds . Also returned ar e the sca ling ma trices, D and Dinv , corresponding to the upper b ound. The smallest destabilizing p ertur bation at each frequency is returned a s Delta . The v ar iable sens is the sensitivity of the upp [...]

140 CHAPTER 4. DEMONSTRA TION EXAMPLES 1 10 100 1000 0.1 10000 0.6 0.8 1 1.2 1.4 1.6 0.4 1.8 Max. singular value and mu comparison max. singular value mu upper bound mu lower bound Note that µ ( g1g ) is not less than one at all frequencies — we ha ve not met the design ob jectiv es. A D - K iteratio n will be used to low er µ and imp rov e the[...]

4.1. TH E HIMA T EXAMPLE 141 The second thing to note is that the in ter connection structure has the additional control inputs and meas urement outputs. The D and Dinv systems must be augmented wit h iden tities corresponding to these additional inputs and outputs. The m usynfit function per forms both of these op eratio ns. The syntax of m usyn?[...]

142 CHAPTER 4. DEMONSTRA TION EXAMPLES Frequency (Hz) 1 10 100 1000 0.1 10000 Magnitude 0.1 1 0.01 10 D scale fit, block: 1 Magnitude data Previous fit, order: 0 New fit, order: 3 Frequency (Hz) 1 10 100 1000 0.1 10000 Magnitude 1 0.1 10 Upper bound comparison, block: 1 Data based bound Previous fit bound, order: 0 New fit bound, order: 3 Frequency[...]

4.1. TH E HIMA T 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 b e pre and p ost m ultiplied o n to the orginal in terconnection structure. himat ic2 = D1sys * himat ic * D1invsys [,,nx] = size(himat ic2) display "himat ic2 now has [...]

144 CHAPTER 4. DEMONSTRA TION EXAMPLES 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 pr obably resulted in a co ntrol la w which achiev ed an infinity norm of approximately 1.1 for the new interc[...]

4.1. TH E HIMA T EXAMPLE 145 1 10 100 1000 0.1 10000 0.6 0.8 1 1.2 1.4 1.6 0.4 1.8 Robustness analysis for the g2 system g2: mu upper bound g2: mu lower bound g1: mu upper bound g1: mu lower bound comment bnds2 "mu bounds: iteration 2" comment D2 "D scale: iteration 2" comment D2inv "D inverse scale: iteration 2" comme[...]

146 CHAPTER 4. DEMONSTRA TION EXAMPLES g3 = starp(himat ic,k3) g3g = freq(g3,omega) [bnds3,D3,D3inv,Delta3 ,sens3] = mu(g3g,blk) A t this point we could do ano ther iteration (to get k4 ) or p erha ps run some simulations to chec k out k3 more thoroughly . 4.1.8 Sim ulation Comparison with a Loopshaping Con troller A lo opsha ping design is p erfor[...]

4.1. TH E HIMA T EXAMPLE 147 d = zeros(2,2) klp = system(a,b,c,d) comment klp "loop shape controller" W e will compare the designs, with no e rror or uncert aint y we igh ts, for th e nominal case and with a p erturba tion blo ck of ∆ = [0.1,0;0,-0 .1] fo r the input multiplicativ e pe rturbatio n. The time resp onse will b e fro m 0 to[...]

148 CHAPTER 4. DEMONSTRA TION EXAMPLES 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 in tr o duced into the first c hannel. time = 0:2:0.01 u = gstep(time) # a unit step is [...]

4.1. TH E HIMA T EXAMPLE 149 time 0.5 1 1.5 0 2 -0.5 0 0.5 1 -1 1.5 Kmu step dist. response (nominal) time 0.5 1 1.5 0 2 -0.5 0 0.5 1 -1 1.5 Kinf step dist. response (nominal) time 0.5 1 1.5 0 2 -0.5 0 0.5 1 -1 1.5 Klp step dist. response: (nominal)[...]

150 CHAPTER 4. DEMONSTRA TION EXAMPLES The lo opsha ping design gives a decoupled resp onse. Both the H ∞ and µ designs trade decoupling for speed of resp onse and, as w e shall see, robustness with resp ect to per turbations. The simulation is rep ea ted with a per turbation o f size 0 .1. Note that this is o nly 10% of the size p ertur bation [...]

4.1. TH E HIMA T EXAMPLE 151 time 0.5 1 1.5 0 2 -0.5 0 0.5 1 -1 1.5 Kmu step dist. response (perturbed) time 0.5 1 1.5 0 2 -0.5 0 0.5 1 -1 1.5 Kinf step dist. response (perturbed) time 0.5 1 1.5 0 2 -0.5 0 0.5 1 -1 1.5 Klp step dist. response: (perturbed)[...]

152 CHAPTER 4. DEMONSTRA TION EXAMPLES The lo opsha ping controller had goo d nominal perfo rmance and very p o or r obust pe rformance. This was illustra ted with a relatively small pertur bation. The difference betw een the µ and H ∞ controllers was small in b oth the nominal a nd p ertur be d cases. This ma y not alw ays b e the case for sev[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 153 4.2 A Simple Flexible Structure Example The demonstration scrip t, jplphBdemo.ms , runs throu gh a D - K iteration desig n for a simple flexible structur e problem. The following demo can be run by executing the f ollowing Xmath command: execute file = "$XMATH/demos/xMu/ jplphBdemo " where $[...]

154 CHAPTER 4. DEMONSTRA TION EXAMPLES Figure 4.3: Sc hematic diagram of t he JPL Phase B opti cal testb ed design proble m problem is the fact that the voice-coil mirror assem bly has significant mass and its mov ement excites a mo de in the structur e. An iden tica lly driven co unterbalance effectiv ely mak es the piezo-electric actuator react[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 155 static gain. W e will include some dynamic uncertaint y in the actual design. piezo = 1 4.2.2 Creating the W eigh ted Design In terconnection Structure The weighted, open-lo op design interconnection structure is illustra ted in Figure 4.5. F or clarit y , the t wo perturbat ions, ∆ 1 and ∆ 2 , hav[...]

156 CHAPTER 4. DEMONSTRA TION EXAMPLES The additive weigh t clearly provides for significant high frequency uncer taint y . A m ultiplicative w eight mo dels the low frequency uncertain ty . The v alue selected is somewhat arb itrary and can be considered as a tunable design w eigh t. Wmvoice = 0.1 A m ultiplica tive per turbation is used to mo de[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 157 F = 100 Wactp = makepoly([1/(2*pi*F),1] ,"s")/... makepoly([1/(200*pi*F) ,1],"s") Wactp = Wactp*4 The low er frequency response of the v o ice-coil system mea ns that the W actv we i g ht should begin rolling up at around 10 Hz. F=1 0 Wactv = makepoly([1/(2*pi*F),1] ,"s")/[...]

158 CHAPTER 4. DEMONSTRA TION EXAMPLES Wdist,Wnoise,Wperf,W actp,Wactv ,piezo) size(P)? ans (a row vector) = 7 6 8 Now select a frequency grid for calculating the frequency resp onses. Some additiona l po ints are included near the oscillato ry mo des. omega = logspace(1,1000,15)’ omega = sort([omega; [5.275:0.05:5.625]’; [4.5:0.25:6.5]’])’[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 159 Frequency 10 100 1 1000 Magnitude 0.01 0.1 1 10 100 1000 10000 0.001 100000 voice coil model Frequency 10 100 1 1000 Phase (degrees) -20 0 20 40 60 80 -40 100[...]

160 CHAPTER 4. DEMONSTRA TION EXAMPLES And examine the design weigh ts. weightsg = freq(weights,omega) gph2 = ctrlplot(weightsg, { logma gplot } ); gph2 = plot(gph2,legend=["Wperf ";"Wavoice ";"Wmvoice ";... "Wmpiezo";"Wdist";"W noise";"Wa ctp";"Wact v"],... title="[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 161 Frequency 10 100 1 1000 Magnitude 0.01 0.1 1 10 0.001 100 Design weights Wperf Wavoice Wmvoice Wmpiezo Wdist Wnoise Wactp Wactv[...]

162 CHAPTER 4. DEMONSTRA TION EXAMPLES 4.2.3 Design of an H ∞ Con troller An H ∞ design is now performed. Recall that we hav e one interferometer measuremen t and tw o co nt r oller outputs. nmeas = 1 ncon = 2 glimits = [0;20] [Khinf,gamma] = hinfsyn(P,nmeas,nc on,glimit s, { tol=0.2 5 } ) Test bounds: 0.0000 < gamma <= 20.0000 gamma Hx e[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 163 -3.5543e+02 -3.5543e+02 5.0265e+02 0.7071 -3.5543e+02 3.5543e+02 5.0265e+02 0.7071 -6.4598e+02 -3.3537e+01 6.4685e+02 0.9987 -6.4598e+02 3.3537e+01 6.4685e+02 0.9987 -3.8358e+02 -5.3777e+02 6.6055e+02 0.5807 -3.8358e+02 5.3777e+02 6.6055e+02 0.5807 -1.6934e+03 0.0000e+00 1.6934e+03 1.0000 -2.8107e+03 0[...]

164 CHAPTER 4. DEMONSTRA TION EXAMPLES Frequency 10 100 1 1000 Magnitude 0.001 0.01 0.1 0.0001 1 Controller: Khinf Frequency 10 100 1 1000 Phase (degrees) -50 0 50 100 150 -100 200[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 165 4.2.4 Robustness Analysis The blo ck str ucture has tw o p ertur bation ∆ blo cks and a “per formance” blo ck. The t wo voice-coil per turbations are put in to a sing le 1 × 2 blo ck as they enter the system at the same poin t. Note that this is not iden tical to t wo sep arate blocks — f o r [...]

166 CHAPTER 4. DEMONSTRA TION EXAMPLES gph4 = ctrlplot(npbnds, { log,lin e style=4 } ); gph4 = ctrlplot(rsbnds,gph4, { lo g,line style=[3,5] } ); gph4 = ctrlplot(rpbnds,gph4, { lo g,line style=[1,2] } ); gph4 = plot(gph4, { !grid,legend= ["Nom perf";"Rob stab (up bnd)";... "Rob stab (lw bnd)"; "Rob perf (up bnd)&q[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 167 0.5 1 1.5 0 2 -0.01 -0.005 0 0.005 0.01 -0.015 0.015 perturbed closed loop: vc actuator Khinf Kmu[...]

168 CHAPTER 4. DEMONSTRA TION EXAMPLES 4.2.5 D - K Iteration W e will now perfo rm one D - K itera tion to generate the con tro ller Kmu . Significant robustness and p e rformance improvemen t is achiev ed with only one iteration. T ransfer f unctions are fit to t he D -scales from the previous robust performance µ test. Here w e preselect an or[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 169 Kmug = freq(Kmu,omega) gph5 = ctrlplot(Kmug, { bode } ); gph5 = plot(gph5, { title="Contro ller: Kmu" } )?[...]

170 CHAPTER 4. DEMONSTRA TION EXAMPLES Frequency 10 100 1 1000 Magnitude 0.001 0.01 0.1 0.0001 1 Controller: Kmu Frequency 10 100 1 1000 Phase (degrees) -50 0 50 100 150 -100 200[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 171 W e now examine the robustness p roperties of the new closed lo op system. W e already kno w that the robust p erformance t est will b e less than t he γ v alue fr om the H ∞ syn thesis ab ove (in this case 0.2930). The results are again display ed graphically . Gmug = freq(Gmu,omega) npbnds = norm([...]

172 CHAPTER 4. DEMONSTRA TION EXAMPLES 10 100 1 1000 0.05 0.1 0.15 0.2 0.25 0 0.3 mu analysis of Gmu Nom perf Rob stab (up bnd) Rob stab (lw bnd) Rob perf (up bnd) Rob perf (lw bnd)[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 173 4.2.6 A Sim ulation Study No w the t w o con trollers ( Khinf and Kmu ) are studied b y simulation. An unw eighted interconnection is set up with sysic and starp is used to clos e the lo op for each con troller. ssnames = ["vcmodel"; "Wavoice"; "Wmvoice"; "Wmpiezo&quo[...]

174 CHAPTER 4. DEMONSTRA TION EXAMPLES 0.5 1 1.5 0 2 -0.004 -0.002 0 0.002 0.004 -0.006 0.006 Simulation: disturbance[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 175 gph8 = ctrlplot(u(2,1)); gph8 = plot(gph8, { title="Simula tion: noise",... y lab="micrometers" } )?[...]

176 CHAPTER 4. DEMONSTRA TION EXAMPLES 0.5 1 1.5 0 2 micrometers -0.002 -0.001 0 0.001 0.002 -0.003 0.003 Simulation: noise[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 177 A nominal r esp onse is calcula ted by setting ∆ = 0. T o get the op en-lo o p simulation mo del we close the unw eighted in ter connection structure with a co ntroller equa l to zero . deltazero = zeros(2,3) Kzero = zeros(2,1) nomolp = starp(deltazero,P nom) nomolp = starp(nomolp,Kzer o) yolp = nomo[...]

178 CHAPTER 4. DEMONSTRA TION EXAMPLES 0.5 1 1.5 0 2 micrometers -5 0 5 -10 10 open loop beam length[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 179 Now w e consider the closed-lo op nominal resp onse with Khinf and Kmu . The closed-lo o p system happens to hav e a large n umber of high frequency p oles whic h do not con tribute significantly to the respo nse. They hav e the effect o f forcing a very fine time discretization in the sim ulation, [...]

180 CHAPTER 4. DEMONSTRA TION EXAMPLES 0.5 1 1.5 0 2 micrometers -0.2 -0.1 0 0.1 0.2 0.3 -0.3 0.4 closed loop beam length Khinf Kmu[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 181 gph11 = ctrlplot(yclphinf(2,1)) ; gph11 = ctrlplot(yclpmu(2,1),gp h11); gph11 = plot(gph11, { legend=["Kh inf";"Kmu" ],... title="closed loop: vc actuator" } )?[...]

182 CHAPTER 4. DEMONSTRA TION EXAMPLES 0.5 1 1.5 0 2 -0.01 -0.005 0 0.005 0.01 -0.015 0.015 closed loop: vc actuator Khinf Kmu[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 183 gph12 = ctrlplot(yclphinf(3,1)) ; gph12 = ctrlplot(yclpmu(3,1),gp h12); gph12 = plot(gph12, { legend=["Kh inf";"Kmu" ],... title="closed loop: piezo actuator" } )?[...]

184 CHAPTER 4. DEMONSTRA TION EXAMPLES 0.5 1 1.5 0 2 -0.01 0 0.01 0.02 -0.02 0.03 closed loop: piezo actuator Khinf Kmu[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 185 Note that Kmu achiev es better performance at the expense of greater actu ator effort. W e will now repeat this sim ula tion for a pertur b ed system. A bad ∆ is chosen and scaled to ha ve norm 0.5. This is obtained fro m destabilizing ∆ c alculated for the µ lower bo und. An all-pa ss system is [...]

186 CHAPTER 4. DEMONSTRA TION EXAMPLES gph13 = ctrlplot(ybclpmu(1,1),g ph13); gph13 = plot(gph13, { legend=["Kh inf";"Kmu" ],... title="perturbed closed loop beam length",... y lab="micrometers" } )?[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 187 Frequency 10 100 1 1000 Magnitude 0.001 0.01 0.1 0.0001 1 Controller: Khinf Frequency 10 100 1 1000 Phase (degrees) -50 0 50 100 150 -100 200[...]

188 CHAPTER 4. DEMONSTRA TION EXAMPLES gph14 = ctrlplot(ybclphinf(2,1) ); gph14 = ctrlplot(ybclpmu(2,1),g ph14); gph14 = plot(gph14, { legend=["Kh inf";"Kmu" ],... title="perturbed closed loop: vc actuator" } )?[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 189 0.5 1 1.5 0 2 -0.01 -0.005 0 0.005 0.01 -0.015 0.015 perturbed closed loop: vc actuator Khinf Kmu[...]

190 CHAPTER 4. DEMONSTRA TION EXAMPLES gph15 = ctrlplot(ybclphinf(3,1) ); gph15 = ctrlplot(ybclpmu(3,1),g ph15); gph15 = plot(gph15, { legend=["Kh inf";"Kmu" ],... title="perturbed closed loop: piezo actuator" } )?[...]

4.2. A SIMPLE FLEX IBLE STR UCTURE EXAMPLE 191 0.5 1 1.5 0 2 -0.02 -0.01 0 0.01 0.02 -0.03 0.03 perturbed closed loop: piezo actuator Khinf Kmu[...]

Bibliograph y [1] J. C. Doyle, “Lecture notes on adv ances in multiv ariable control.” ONR/Honeyw ell W orkshop, Minneapolis, MN., 1984. [2] J. Doyle, “Structured uncertaint y in cont r ol system design,” in Pr o c . IEEE Cont r ol De cision Con f. , pp . 260–265, 1985. [3] A. K. Pac k ard, What’s new wit h µ : Structur e d Unc ertai n[...]

194 BIBLIOGRAPHY [11] G. J. Ba las and J. C. Doyle, “Robust con tr ol of flexible modes in the controller crossov er r egion,” in Pr o c. Amer. Control Conf. , 1989. [12] G. J. Ba las, A. K. Pack ar d, and J. Harduvel, “Application of µ -syn thesis techniques to momentum management and attitude co ntrol of the Space Statio n,” in AIAA Gui[...]

BIBLIOGRAPHY 195 [24] D. L. Laughlin, K. G. Jordan, and M. Morari, “In ternal mo del con trol and process uncertaint y: mapping uncerta int y regions for SISO controller design,” Int. J. of Control , v ol. 44, no. 6, pp . 1675–1698, 1986. [25] R. S. Smith and M. Dahleh, eds., The Mo deling of Uncertainty in Cont r ol Systems: P roceed i n g s[...]

196 BIBLIOGRAPHY [37] R. Kosut, M . Lau, and S. Bo yd, “Param eter set ident ification of systems wit h uncertain non parametric dynam ics and distu rbances,” in Pr o c. I EEE Cont r ol De cision Con f. , v ol. 6, pp. 3162–3167, 1990. [38] G. Go o dwin, B. Ninness, and M. Salgado, “Quantification o f uncertaint y in estimation,” in Pro [...]

BIBLIOGRAPHY 197 [50] J. M. Krause, P . P . Khargonek ar, and G. Stein , “Robust param eter adjust ment with nonparametric w eigh ted-ball-in- H ∞ uncertaint y ,” IEEE T r an s. Auto. Con tr ol , v ol. AC -35, pp. 225–229, 1990. [51] R. S. Smith a nd J. C. Doyle, “Closed lo op r elay estimation of uncertaint y b o unds for robust control [...]

198 BIBLIOGRAPHY [63] A. J. La ub, “A Sch ur method for solving algebr aic Riccati equatio ns,” IEEE T r ans. A uto. Co ntr ol , v ol. A C-24, pp . 913–921, 1979. [64] T. Pappas, A. J. Laub, a nd N. R. Sandell, “On the n umer ical solution o f the discrete-time algebraic Riccati equation,” IEEE T rans. A u to. Cont r ol , v ol. A C-25, pp[...]

BIBLIOGRAPHY 199 [76] M. Dahleh, A . T esi , and A. V icino, “Extremal prop erties for t he parametric robust pe rformance pro blem,” T ech. Rep. UCSB-ME-9 1-4, Univ. Califo rnia, Santa Barbara, Mec h. Eng., 1991. also submitted to 30th IEEE CD C. [77] B. C. Mo ore, “Principal comp onents analysis in linear systems: controllability , observ a[...]

Chapter 6 F unction Reference 6.1 X µ F unctions The following pages contain descriptions of the X µ functions. These are also av ailable on-line via the help utility . Each description also gives an illustrative example of the function’s use. The functio ns are included in alphab etical o rder. F or con venience they a re cross-refere nced by [...]

202 CHAPTER 6. FUNCTION RE FERENCE rifd ....... ........ ......... ........ ....... ........ ......... ........ ...... 3 3 5 ctrlplot ...... ........ ......... ........ ....... ........ ......... ........ .. 2 2 1 Time r esp onse c alculations and pdm functions gstep ..... ......... ........ ......... ...... ......... ........ ......... ..... 2 4 7[...]

6.1. X µ FUNCTIONS 203 Mo del r e duction and state-sp ac e functions balmoore ...... ........ ......... ........ ....... ........ ......... ........ .. 2 0 5 modalstate ....... ........ ......... ........ ....... ........ ......... ...... 2 9 3 ophank ...... ........ ......... ........ ...... ......... ........ ......... .... 3 1 5 orderstate ...[...]

204 CHAPTER 6. FUNCTION RE FERENCE Misc el lane ous functions conpdm ...... ........ ......... ........ ...... ......... ........ ......... .... 2 1 3 consys ...... ........ ......... ........ ...... ......... ........ ......... .... 2 1 5 csum ....... ........ ......... ........ ....... ........ ......... ........ ...... 2 1 7 delsubstr ...... ...[...]

balmoore 205 balmo ore Syn t ax [SysR,HSV,T] = balmoore(Sys, { nsr,bound } ) Pa r a m e t e r Li s t Inputs: Sys Linear, stable, minima l state-spac e system nsr (optional ) If bound is used then a reduction will b e p er- formed whic h meets an error b ound specified by the v alue in nsr ,o t h e r w i s e nsr is the order of the reduced system. [...]

206 Chapter 6. F unctio n Reference Reference B.C. Mo ore, “Prin cip al Comp onen t Analysis in Linear Systems: Control lability, Observabilit y and Mo del R ed uction ,” IEEE T rans. Auto. Ct rl., V ol. 26, N o. 1, pp. 17–32, F eb. 1981. Exam ple # Create a five state system for reduction. a = daug(-0.891334,[-1.208 57,0.79904 2;-0.79904 2,-[...]

balmoore 207 See Also: minimal , ophank .[...]

blknorm 209 blknorm Syn t ax normM = blknorm(M,blk,p,Frobeni us) Pa r a m e t e r Li s t Inputs: M Matrix (or pdm ). blk Blo ck structure . See mu sectio n of the manual for a de- scription of the syn ta x. p Scalar v alued. Specifies the Ho lder “p” norm to b e used, where 1 ≤ p ≤ inf. Optional. The default is p = 2 . Keywords: F ro beniu[...]

210 Chapter 6. F unctio n Reference Examples A = random(3,3)-0.5*ones(3 ,3)? A (a square matrix) = 0.0618661 0.0896177 0.185398 0.390622 0.00422128 -0.150638 -0.112622 0.42229 0.448818 blknorm(A,[1,1; 1,1; 1,1]) ans (a square matrix) = 0.0618661 0.0896177 0.185398 0.390622 0.00422128 0.150638 0.112622 0.42229 0.448818 blknorm(A,[3,3]) ans (a scalar[...]

blknorm 211 2.23607 3 4 7.81025 7 8 13.4536 11 12 # and compare to norm(B(1,1:2)) ans (a scalar) = 2.23607 See Also norm[...]

conpdm 213 conp dm Syn t ax outpdm = conpdm(mat,domain, { skipChk s } ) Pa r a m e t e r Li s t Inputs: mat constant matrix domain domain over which outsys will be defined. Keywords: skipChks Bo olea n specifying that synt a x chec king is to be skipp ed. Outputs: outp dm pdm Description Creates a pdm data ob ject from a constan t matrix. Outp dm [...]

consys 215 consys Syn t ax outsys = consys(mat, { skipChks } ) Pa r a m e t e r Li s t Inputs: mat constant matrix Keywords: skipChks Bo olea n specifying that synt a x chec king is to be skipp ed. Outputs: outsys Dynamic System Description Creates a Dynamic System ob j ect from a constan t matrix. Outsys represen ts a constant gain; the A , B and [...]

csum 217 csum Syn t ax [outpdm] = csum(inpdm, { channels } ) Pa r a m e t e r Li s t Inputs: inpdm real or complex v alued pdm or constant matrix Keywords: channels Sum ov er c ha nnels. outpdm ha s the same dimensions as inpdm. Outputs: outp dm o utput pdm Description Pe rform a cum ulati ve sum o ver the row s of a matri x or a pdm . If c ha nnel[...]

218 Chapter 6. F unctio n Reference ans (a rectangular matrix) = 1 0.608453 2 1.46287 3 1.52714 4 2.35505 5 3.28128 6 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 -------+--------[...]

csum 219 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.24905 -------+-------------- ---------[...]

ctrlplot 221 ctrlplot Syn t ax graph = ctrlplot(pdm,old graph, { keywords } ) Pa r a m e t e r Li s t Inputs: p dm Pdm (or matrix) containing the data to be plotted. old graph (optional) Graphical ob ject to which data is added. Con- ceptually the same as plot(pdm, { keep=old graph } ) . Keyw ords: The foll owing keyw ords specify the basi c plot f[...]

222 Chapter 6. F unctio n Reference linear linear domain. Default = 1 for timeresp k eyword. Defa ult = 0 for bo de k eyword. log logarith mic dom ain. Default = 1 for bo de ke yword. De- fault = 0 for timeresp keyw ord. Default units can b e supplied for the magnitude and phase plots (b o de, nichols, logmagplo t and pha seplot keyw ords) with the[...]

ctrlplot 223 Description This function per forms some co mmon control system related plotting. The user can use ctrlplot to set up a basic plot and per form some pr epro cessing of the data. This generates a graphical ob ject containing the data and the user ca n per form subsequen t calls to plot to add things like text , lab els, gridding etc. Th[...]

224 Chapter 6. F unctio n Reference legend = ["sys1","sys2"] } )?[...]

ctrlplot 225 Frequency 0.1 1 0.01 10 Magnitude 0.0001 0.001 0.01 0.1 1 1e-05 10 Bode plots sys1 sys2 Frequency 0.1 1 0.01 10 Phase (degrees) -250 -200 -150 -100 -50 -300 0[...]

226 Chapter 6. F unctio n Reference # Nyquist plots g2 = ctrlplot(sys1g, { nyqui st } ); g2 = ctrlplot(sys2g,g2, { ny quist } ); g2 = ctrlplot(-1,g2, { nyqui st,marker= 1,line=0 } ); g2 = plot(g2, { projection=" orthograph ic",... legend=["sys1","sys3" ,"critical point"],title="N yquist plots" } )?[...]

ctrlplot 227 Real -1 0 1 -2 2 Imaginary -1.5 -1 -0.5 0 -2 0.5 Nyquist plots sys1 sys3 critical point[...]

228 Chapter 6. F unctio n Reference # 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.05:10],[0;1 ;5],[0;1;- 1]) y0 = sys*u [a,b,c,d] = abcd(sys) sys = system(a,b,c,d) y1 = system(sys, { X0=[-1[...]

ctrlplot 229 2468 0 10 -2 -1 0 1 -3 2 input x0 = zero non-zero x0[...]

230 Chapter 6. F unctio n Reference See Also: plot .[...]

daug 231 daug Syn t ax out = daug (sys1,sys2,...) Pa r a m e t e r Li s t Inputs: sys1 Input systems. These can b e dynamical systems and co n- stants, or pdms and consta nt s. . . .“ Outputs: out output system. Description Diagonal augmen tation of dyn amical system/pdm/constan t, matri ces. out =      sys1 0 ... 0 0 sys2 ... 0 . [...]

232 Chapter 6. F unctio n Reference ans (a square matrix) = 110 0 002 0 002 0 000I n f 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.748961 X0 0 State Names ----------- sys1state System is continuous sys2 = randsy[...]

daug 233 -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.012459 0 0 0 0 0 0.579502 C 0.186754 0 00 0 0.262815 0 0.436099[...]

234 Chapter 6. F unctio n Reference D 0.492058 0.748961 0 0 001 0 0 0 0 0 0.911055 0 0 0 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 -------+-------------- ------[...]

daug 235 0| 1 0 -------+----- 1| 1 0 -------+----- 2| 1 0 -------+----- 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[...]

delsubstr 237 delsubstr Syn t ax [outstr] = delsubstr(str,charstr) Pa r a m e t e r Li s t Inputs: str String or vector of strings. charstr String Outputs: outstr Str ing or vector of strings. Description All o ccur ences of the substring, charstr , w ithin str a re deleted. If, by deleting charstr , another o ccurence of charstr in cr eated, it wi[...]

238 Chapter 6. F unctio n Reference out2 (a column vector of strings) = string one aaa xy # If executed, the following would give an error # out3 = delsubstr(strvec,"a" )?[...]

fitsys 239 fitsys Syn t ax [sys] = fitsys(data,npoles ,nzeros,we ight, { skipchks,Hertz } ) Pa r a m e t e r Li s t Inputs: data Complex v alued data ( pdm ). npo les Order of requested fit. (optional, defa ult = 0). nzeros Num b er of zeros in t ransfer function . (optional, defau lt = npo les) weigh t W eigh ting functio n. (scalar, pd m ,o r D[...]

240 Chapter 6. F unctio n Reference The primary use of this routine is the fitting of D scal e w eight s for mu syn thesis iterations. Cheb yshev polynomials are used as basis funct ions for b oth the n umerator and denominator p o lynomials. W ARNING: This routine uses iterative polyno mial calculations which are not well conditioned for high or [...]

fitsys 241 real imaginary frequency damping (rad/sec) ratio -1.0000e+00 0.0000e+00 1.0000e+00 1.0000 5.0000e-01 3.1225e+00 3.1623e+00 -0.1581 5.0000e-01 -3.1225e+00 3.1623e+00 -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[...]

242 Chapter 6. F unctio n Reference Frequency 0.01 0.1 1 10 0.001 100 Phase (degrees) -300 -200 -100 0 -400 100 Data and minimum phase fit original data minimum phase fit[...]

fitsys 243 # Create fitting weight. 1/s works well for logspaced # data. 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 damping (rad/sec) ratio -5.1043e-02 -3.1249e-01 3.1663e-01 0.1612 -5.1[...]

244 Chapter 6. F unctio n Reference Frequency 0.01 0.1 1 10 0.001 100 Magnitude 0.1 1 10 0.01 100 Data and minimum phase fit original system minimum phase system Frequency 0.01 0.1 1 10 0.001 100 Phase (degrees) -300 -200 -100 0 -400 100[...]

fitsys 245 Limita tions Limited to SISO systems . See Also tfid[...]

gstep 247 gstep Syn t ax gPdm = gstep (ytime,timespec,valspec , { skipChks } ) Pa r a m e t e r Li s t Inputs: ytime output time v ector (seconds). timespec times for specified step data (opti onal) v alspec v alue for specified step dat a (option al) Keywords: skipChks Bo olea n specifying that synt a x chec king is to be skipp ed. Outputs: gPdm[...]

248 Chapter 6. F unctio n Reference Frequency 0.01 0.1 1 10 0.001 100 Magnitude 0.1 1 10 0.01 100 Data and minimum phase fit original system minimum phase system Frequency 0.01 0.1 1 10 0.001 100 Phase (degrees) -300 -200 -100 0 -400 100[...]

gstep 249 See Also randpdm , gcos , gsin , gpulse , gsawtooth , gsquarewave[...]

hinfnorm 251 hinfnorm Syn t ax [out,omega] = hinfnorm(sys,tol, { imag eps,max it } ) Pa r a m e t e r Li s t Inputs: sys Dynamic System , frequency response ( pdm ), or constant gain (matr ix). sys Sp ecifies the relative tolerance of the answ er when the in- put is a Dynamic System Def ault = 0.001. Keyw ords: imag eps E psilon v alue for determi[...]

252 Chapter 6. F unctio n Reference Stable Dynamic System norms a re calculated by an itera tive Hamiltonian metho d. In this case out is a 2 × 1 vector with upp er and low er b o unds for the no rm. Exam ple # 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],&q[...]

hinfnorm 253 25.000 5.2e-01 1.7e-03 1.0e-02 0.0e+00 0.0000 p 12.500 5.2e-01 1.7e-03 1.0e-02 0.0e+00 0.0000 p 6.250 5.2e-01 1.7e-03 1.0e-02 0.0e+00 0.0000 p 3.125 5.1e-01 1.7e-03 1.0e-02 0.0e+00 0.0000 p 1.562 5.0e-01 1.7e-03 1.0e-02 0.0e+00 0.0000 p 0.781 3.9e-01 -9.8e+01 1.0e-02 0.0e+00 0.0000 f 1.172 4.8e-01 1.8e-03 1.0e-02 0.0e+00 0.0000 p Gamma[...]

h2norm 255 h2norm Syn t ax out = h2norm(sys) Pa r a m e t e r Li s t Inputs: sys Co n tin uous time Dynamic System Outputs: out H 2 no rm of the input system Description The H 2 norm of a stable, str ictly prop er system is calcula ted. This is g iven b y out = trace( CXC 0 ) , where X is the cont rollability grammian, solving the Ly apunov equatio[...]

256 Chapter 6. F unctio n Reference Wperf = 100/makepoly([100,1],"s ") Wact = makepoly([0.5,0.05],"s") /makepoly( [0.05,1]," s") # Form the weighted interconnection structure sysnames = ["plant";"Wperf";"Wact"] sysinp = ["ref";"control"] sysout = ["Wperf"; &qu[...]

hinfsyn 257 hinfsyn Syn t ax [k,gfin,stat] = hinfsyn(p,nmeas,ncon,gamma, { keywords } ) Pa r a m e t e r Li s t Inputs: p Generalized interconnection structur e ( Dynamic Sys- tem ) nmeas measuremen t v ector dimension. ncon con trol v ector dimension. gamma H ∞ nor m bound o f cont r oller. F or a bisection search spec- ify gamma = [gamma min;ga[...]

258 Chapter 6. F unctio n Reference Description The H ∞ (sub)optimal controller for the interconnection, p , is ca lculated. The resulting closed lo op system is illustrated b elow. p k - zw y u The v ariab les ncon and nmeas are used to sp ecify the dimensions of u and y in the ab ove diagram ( ncon = dim( u )a n d nmeas = dim( y )). The ob ject[...]

hinfsyn 259 1. ( a, b 2 ,c 2 ) is stabiliza ble and detectable 2. d 12 and d 21 ha ve full rank 3. The matri x [ a − jωI, b 2 ; c 1 ,d 12 ] has fu ll column ran k for all ω 4. The matri x [ a − jωI, b 1 ; c 2 ,d 21 ] has full row rank for all ω Reference This function uses the state-spa ce formulae given in: “Stat e-sp ace formu lae for a[...]

260 Chapter 6. F unctio n Reference wghtic = sysic(sysnames,sy sinp,sysou t,syscnx,p lant,... Wperf,Wact) # Design Hinf controller nctrls = 1 nmeas = 1 gmax = 25 gmin = 0 Kinf = hinfsyn(wghtic,nmeas,nct rls,[gmax; gmin]) Test bounds: 0.0000 < gamma <= 25.0000 gamma Hx eig X eig Hy eig Y eig nrho xy p/f 25.000 5.2e-01 1.7e-03 1.0e-02 0.0e+00 0[...]

hinfsyn 261 real imaginary frequency damping (rad/sec) ratio -5.0000e-02 -3.1225e-01 3.1623e-01 0.1581 -5.0000e-02 3.1225e-01 3.1623e-01 0.1581 -5.0000e+00 0.0000e+00 5.0000e+00 1.0000 -2.0000e+01 0.0000e+00 2.0000e+01 1.0000 omega = logspace(0.001,100,200) Kinfg = freq(Kinf,omega) gph1 = ctrlplot(Kinfg, { bode } ); gph1 = plot(gph1, { title="[...]

262 Chapter 6. F unctio n Reference Frequency 0.01 0.1 1 10 0.001 100 Magnitude 0.01 0.1 1 0.001 10 Kinf Frequency 0.01 0.1 1 10 0.001 100 Phase (degrees) -80 -60 -40 -20 0 20 -100 40[...]

hinfsyn 263 # Use sysic to create unweighted interconnection ic = sysic("plant",["ref"; "ctrl"],[" plant";"re f-plant"] ,... "ctrl",plant) clpinf = starp(ic,Kinf) rifd(clpinf) Poles: real imaginary frequency damping (rad/sec) ratio -5.0000e-02 3.1225e-01 3.1623e-01 0.1581 -5.0000e-02 -3.12[...]

264 Chapter 6. F unctio n Reference gph2 = ctrlplot(sensg, { logmagpl ot } ); gph2 = plot(gph2, { title="Kinf controller: sensitivity function" } )?[...]

hinfsyn 265 Frequency 0.01 0.1 1 10 0.001 100 Magnitude 0.1 1 0.01 10 Kinf controller: sensitivity function[...]

266 Chapter 6. F unctio n Reference # Examine step response step = gstep([0:0.1:10],0,1) y = clpinf*step gph3 = ctrlplot([y,step]); gph3 = plot(gph3, { title="Kinf controller: step response" } )?[...]

hinfsyn 267 2468 0 10 0.2 0.4 0.6 0.8 1 0 1.2 Kinf controller: step response[...]

268 Chapter 6. F unctio n Reference See also hinfsyn , hinfnorm , h2norm[...]

h2syn 269 h2syn Syn t ax k = h2syn(p,nmeas,ncon, { ke ywords } ) Pa r a m e t e r Li s t Inputs: p Generalized interconnection structur e ( Dynamic Sys- tem ) nmeas measuremen t v ector dimension. ncon con trol v ector dimension. Keyw ords: sc hur solution real Sch ur decomp osition for Riccati solution (default) eig solution eigendecomp ositio n f[...]

270 Chapter 6. F unctio n Reference p k - zw y u The v ariab les ncon and nmeas are used to sp ecify the dimensions of u and y in the ab ove diagram ( ncon = dim( u )a n d nmeas = dim( y )). The ob jective is to design a sta bilizing con troller, k , whic h minimizes the H 2 norm of the closed loop system betw een w and z . p is a s tate-space syst[...]

h2syn 271 Reference This function uses the state-spa ce formulae given in: “Stat e-sp ace formu lae for al l s tabilizing contr ol lers that sat isfy an H ∞ norm b ound and r elati ons to risk sensiti vity,” Keith Glov e r and John Doyle, Systems & Co nt r ol Letters 11 , pp. 167–172., Oct, 1988. Exam ple # Set up a simple closed loop p[...]

272 Chapter 6. F unctio n Reference K2 = h2syn(wghtic,nmeas,nc trls) rifd(K2) Poles: real imaginary frequency damping (rad/sec) ratio -1.4046e-01 -2.3161e-01 2.7087e-01 0.5186 -1.4046e-01 2.3161e-01 2.7087e-01 0.5186 -1.5863e+00 3.4754e+00 3.8203e+00 0.4152 -1.5863e+00 -3.4754e+00 3.8203e+00 0.4152 -5.2060e+00 0.0000e+00 5.2060e+00 1.0000 Zeros: re[...]

h2syn 273 Frequency 0.01 0.1 1 10 0.001 100 Magnitude 0.001 0.01 0.1 0.0001 1 K2 Frequency 0.01 0.1 1 10 0.001 100 Phase (degrees) -100 -50 0 -150 50[...]

274 Chapter 6. F unctio n Reference # Use sysic to create unweighted interconnection ic = sysic("plant",["ref"; "ctrl"],[" plant";"re f-plant"] ,... "ctrl",plant) clp2 = starp(ic,K2) rifd(clp2) Poles: real imaginary frequency damping (rad/sec) ratio -5.0000e-02 3.1225e-01 3.1623e-01 0.1581[...]

h2syn 275 sensg = freq(sens,omega) gph2 = ctrlplot(sensg, { logmagpl ot } ); gph2 = plot(gph2, { title="K2 controller: sensitivity function" } )?[...]

276 Chapter 6. F unctio n Reference Frequency 0.01 0.1 1 10 0.001 100 Magnitude 1 0.1 10 K2 controller: sensitivity function[...]

h2syn 277 # Examine step response step = gstep([0:0.1:10],0,1) y = clp2*step gph3 = ctrlplot([y,step]); gph3 = plot(gph3, { title="K2 controller: step response" } )?[...]

278 Chapter 6. F unctio n Reference 2468 0 10 0.2 0.4 0.6 0.8 0 1 K2 controller: step response[...]

h2syn 279 See also hinfsyn , h2norm , hinfnorm[...]

interp 281 in terp Syn t ax outpdm = interp(inpdm,stepsize,final { keywords } ) outpdm = interp(inpdm,domspec, { keywords } ) Pa r a m e t e r Li s t Inputs: inpdm Input pdm. stepsize Incremen t in o utp dm domain. final La st v alue in outp dm doma in. Optional: default = max(domain(inp dm)). domsp ec Regula r vector or pdm used to sp ecify the d[...]

282 Chapter 6. F unctio n Reference This function differs from the interpolate function in tha t it can handle zero order hold type interp olation a nd deal with irregula rly spaced input p dms. Irreg ularly spaced output pdms can be generated wit h the domspec syn tax. These featu r es are often useful when dealing with data g enerated from exp e[...]

interp 283 1234 0 5 1 2 3 4 0 5 1st order interp. 0 order interp. original pdm[...]

284 Chapter 6. F unctio n Reference See Also interpolate[...]

mergeseg 285 mergeseg Syn t ax outpdm = mergeseg(pdm1,pdm2, { keywords } ) Pa r a m e t e r Li s t Inputs: pdm1 input pdm pdm2 input pdm Keyw ords: dom sort Sort the result of merging the pdm s. If !domsort then pdm2 is simply concatenated o nto pdm1 . Bo o lean. De- fault = 1. increasing Sort in increasing order. Bo olea n. Default = 1. decreasing[...]

286 Chapter 6. F unctio n Reference Exam ple time1 = [0:0.025:1] pdm1 = gsin(time1, { frequency= 2 } ) time2 = [0.8:0.02:1.5] pdm2 = randpdm([],1,1, { dom=ti me2,zerom ean } ) outpdm = mergeseg(pdm1,pdm 2) gph1 = ctrlplot(outpdm); gph1 = plot(pdm1,gph1, { marker =1,marker style=8,... line=0 } ); gph1 = plot(pdm2,gph1, { marker =1,marker style=1,...[...]

mergeseg 287 0.5 1 0 1.5 -0.5 0 0.5 -1 1 outpdm pdm1 pdm2[...]

mkpert 289 mkp ert Syn t ax [pertsys] = mkpert(Delta,blk,mubnds, { f select,pn orm,Hertz } ) Pa r a m e t e r Li s t Inputs: Delta Low er b ound p ertur bation da ta fro m mu calculatio n. ( pdm ) blk blo ck structur e (refer to m u function do cumentation). (matri x) mubnds Calculated mu bounds ( pdm ). Keywords: fselect Scalar v alued. Specifies[...]

290 Chapter 6. F unctio n Reference The bo unds, m ubnds, ar e used to determine the ”worst-case” fr equency for the in terpolation. The norm of pertsys is 1/(norm(D elta( jω ))) where ω is the c hosen frequency . This is the smallest destabilizing pe rturbatio n at that frequency . Both the in terpo lation frequ ency and the norm of pertsys [...]

mkphase 291 mkphase Syn t ax [cdata] = mkphase(magdata, { skipchks,Hertz } ) Pa r a m e t e r Li s t Inputs: magdata Magnitude data ( pdm ) Keyw ords: skipch ks Bo olean. Skip the error ch ecking. (D efault = 0) Hertz The domain of the pdm is in Hertz. This is the default. !Hertz sp ecifies a domain in ra d/sec. Outputs: cdata Complex v alued data[...]

292 Chapter 6. F unctio n Reference Limita tions Limited to SISO systems . See Also fitsys , ccepstrum[...]

modalstate 293 mo dalstate Syn t ax outsys = modalstate(sys, { keywords } ) Pa r a m e t e r Li s t Inputs: sys Input D ynamic System Keywords: increasing Bo olean. O rder in terms of incr easing magnitude (contin- uous) o r angle (discr ete). Default = 1 decreasing Bo olean. Order in ter ms o f decrea sing magnitude (contin- uous) o r angle (discr[...]

294 Chapter 6. F unctio n Reference -8.1602e-01 1.3353e+00 1.5649e+00 0.5215 -8.1602e-01 -1.3353e+00 1.5649e+00 0.5215 -4.6142e+00 -1.6648e+01 1.7275e+01 0.2671 -4.6142e+00 1.6648e+01 1.7275e+01 0.2671 -2.2478e+01 0.0000e+00 2.2478e+01 1.0000 -3.9525e+01 0.0000e+00 3.9525e+01 1.0000 Zeros: real imaginary frequency damping (rad/sec) ratio -1.0270e+0[...]

mu 295 mu Syn t ax [mubnds,D,Dinv,Delta, sens] = mu(M,blk) Pa r a m e t e r Li s t Inputs: M Matrix or p dm. blk Blo ck structure defined b y a matrix o f dimension: num b er of blo cks × 2. If the ith blo ck has c outputs and r inputs, then blk(i,:) = [r,c]. The default is equiv alent to 1x1 blocks (M must be square). Outputs: mubnds Upp er a nd[...]

296 Chapter 6. F unctio n Reference # 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 = [a,0; b,b; c,jay*c; d,f] V = [0,a; b,-b; c,-jay*c; f*exp(j[...]

mu 297 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 F or an example of how m u is used for system robustness analysis, refer to the on-line help for m usynfit (p age 299). Limita tions This version of the softw ar e cannot handle r ep eated blo [...]

musynfit 299 m usynfit Syn t ax [Dsys,Dinvsys] = musynfit(Dmag,blk,nmeas,nc trls,.. weight,M,order, { keyw ords } ) Pa r a m e t e r Li s t Inputs: Dmag New D matrix fr om mu calcula tion (magnitude data o nly). The domain is assumed t o be in Hertz. blk blo ck structure (refer to mu function do cumentation). nmeas scalar: n umber of measurements [...]

300 Chapter 6. F unctio n Reference Keyw ords: He rtz Boolean. This ke yword i s mandatory as th e funct ion m ust know whether the do main is in Her tz or ra dians/second (specified by !Hertz) to fit correctly . Note that the Xmath function freq assumes that the frequency range is sp ecified in Hertz. plotw eight Bo olean, default = 0 . This wi[...]

musynfit 301 P = 1/makepoly([1,0,-0.01] ,"s") W = makepoly([1,20],"s")/m akepoly([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 ";"cont[...]

302 Chapter 6. F unctio n Reference 1.562 6.8e-01 8.0e-04 7.7e-03 0.0e+00 0.6811 p 1.406 6.7e-01 8.1e-04 7.0e-03 -1.1e-18 2.0900 f 1.484 6.7e-01 8.0e-04 7.4e-03 0.0e+00 1.0394 f 1.523 6.7e-01 8.0e-04 7.5e-03 -5.0e-19 0.8249 p Gamma value achieved: 1.5234 G = starp(wghtic,Kinf) omega = logspace(0.01,100,25) Gg = freq(G,omega) G11g = Gg(1,1) G22g = G[...]

musynfit 303 0.1 1 10 0.01 100 0.2 0.4 0.6 0.8 1 0 1.2 mu analysis nominal perf. robust stab. robust perf.(upper) robust perf. (lower)[...]

304 Chapter 6. F unctio n Reference # Fit transfer functions to D1 & Dinv1 for a mu # synthesis iteration [Ds,Dinvs] = musynfit(D1,blk,nmeas,ncnt rls,sens1 ,Gg, { Hertz } )[...]

musynfit 305 Frequency (Hz) 0.1 1 10 0.01 100 Magnitude 10 100 1000 1 10000 D scale fit, block: 1 Magnitude data Previous fit, order: 0 New fit, order: 3 Frequency (Hz) 0.1 1 10 0.01 100 Magnitude 0.01 0.1 1 10 0.001 100 Upper bound comparison, block: 1 Data based bound Previous fit bound, order: 0 New fit bound, order: 3 Frequency (Hz) 0.1 1 10 0.[...]

306 Chapter 6. F unctio n Reference # Apply the D scales to another H infinity design Kmu = hinfsyn(Ds*wghtic*Di nvs,nmeas, ncntrls,[g min;gmax] ) Test bounds: 0.0000 < gamma <= 10.0000 gamma Hx eig X eig Hy eig Y eig nrho xy p/f 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.6[...]

musynfit 307 0.1 1 10 0.01 100 0.2 0.4 0.6 0.8 1 0 1.2 Kmu & Kinf mu analysis Kmu: robust perf. Kinf: robust perf.[...]

308 Chapter 6. F unctio n Reference # 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 } ) mup[...]

musynfit 309 ymunom(1,1),ymupert(1 ,1)]); gph4 = plot(gph4, { legend=["inpu t step";"Kinf nominal";... "Kinf pert.";"Kmu nominal";"Kmu pert."] } )?[...]

310 Chapter 6. F unctio n Reference 123 0 4 0.2 0.4 0.6 0.8 1 1.2 0 1.4 input step Kinf nominal Kinf pert. Kmu nominal Kmu pert.[...]

musynfit 311 # Compare with a random perturbation yinfrandp = infrpert*step ymurandp = murpert*step gph5 = ctrlplot([step,yinfnom(1 ,1),yinfra ndp(1,1),. .. ymunom(1,1),ymurandp(1 ,1)]); gph5 = plot(gph5, { legend=["inpu t step";"Kinf nominal";... "Kinf rand. pert.";"Kmu nominal";"Kmu rand. pert."] [...]

312 Chapter 6. F unctio n Reference 123 0 4 0.2 0.4 0.6 0.8 1 1.2 0 1.4 input step Kinf nominal Kinf rand. pert. Kmu nominal Kmu rand. pert.[...]

musynfit 313 See Also: mu , hinfsyn , mkpert , hinfnorm .[...]

ophank 315 ophank Syn t ax [SysR,SysU,HSV] = ophank(Sys, { nsr,o nepass } ) Pa r a m e t e r Li s t Inputs: Sys Linear , stable, state-space system (contin uo us) nsr (option al) Order of the reduced system. I f not sp ecified, the user will b e prompted for its v a lue after the Hankel singular v alues a re display e d. Keywords: onepass (Bo olea[...]

316 Chapter 6. F unctio n Reference Uses additional s ubroutines ophiter , ophred , ophmult and stable . This function is cross-licensed from the Model Reduction Mo dule. Exam ple # Create a five state system for reduction. a = daug(-0.891334,[-1.208 57,0.79904 2;-0.79904 2,-1.2085 7],... -4.74685,-21.3013) b = [0.0262569;-0.189601;- 0.113729;0 .21[...]

ophank 317 Frequency 0.1 1 10 0.01 100 Magnitude 0.001 0.01 0.0001 0.1 original system reduced: ophank reduced: sresidualize ophank error sresidualize error[...]

318 Chapter 6. F unctio n Reference See Also minimal , balmoore .[...]

orderstate 319 orderstate Syn t ax outsys = orderstate(sys,indx) Pa r a m e t e r Li s t Inputs: sys Input D ynamic System indx Lists the desired order of the states in outsys. Outputs: outsys o utput dynamic system. Description Reorder th e s tates according to the index argument . The state names and init ial condition are also ordered. The input[...]

320 Chapter 6. F unctio n Reference 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 0.65518 9.0643 3.32295 -2.27378 -8.82005 4.7424 -1.2488 -13.2781 B 0.710595 0.0879738 0.688873 C 0.181512 0.65953[...]

orderstate 321 0.2 0.1 0.3 State Names ----------- s2 s1 s3 Input Names ----------- Input 1 Output Names ------------ Output 1 System is continuous[...]

randpdm 323 randp dm Syn t ax pdmout = randpdm (ndomain,nrows,ncolumns, { keywords } ) Pa r a m e t e r Li s t Inputs: ndomain length o f the domain nrows num b er of rows in p dmout ncolumns num b er of columns in p dmo ut Keywords: complex Bo olean. A complex v alued pdm is gener ated. Default = 0. zeromean Bo o lean. The v alues are shifted so t[...]

324 Chapter 6. F unctio n Reference 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,[...]

randpdm 325 10.0323 | 0.922528 ---------+------------ 18.3191 | 0.81113 ---------+------------ pdm3 = randpdm(0,3,2, { complex,z eromean } )? pdm3 (a rectangular matrix) = -0.237052 - 0.10845 j -0.250958 - 0.392096 j 0.097563 + 0.0140395 j -0.255142 - 0.419362 j 0.0178566 - 0.153032 j 0.0986898 + 0.375528 j[...]

randpert 327 randp ert Syn t ax [pert] = randpert(blk, { sys,sfreq,complex,pnorm } ) Pa r a m e t e r Li s t Inputs: blk block structur e (refer to mu function do cumentation). (matri x) Keywords: sys Bo olean. Sp ecifies that p ert is a dynamic system. Default = !sys, i.e. per t is a matrix. sfreq Scalar v alued. If a dynamic system is sp ecifie[...]

328 Chapter 6. F unctio n Reference Exam ple The use of randpert is studied in co ntext in the on-line help fo r musynfit and the man ual do cumen tation f or musynfit (page 299).[...]

randsys 329 randsys Syn t ax sys = randsys (nstates,noutputs,ninput s, { keywor ds } ) Pa r a m e t e r Li s t Inputs: nstates n um be r of sta tes in sys noutputs num b er o f outputs in sys ninputs num b er o f inputs in s ys Keywords: stable Bo o lean. sys is forced to b e sta ble (default = 1) oscillator y Bo olean. Os cillatory p o les are all[...]

330 Chapter 6. F unctio n Reference Description A rando m system, with user sp ecified state, input and output dimension, is generated. Several additional features can b e specified b y keyw ords: whether or no t the system is stable, where or not it ca n contain oscillatory mo des, whether or no t it ha s a D term, and b ounds on the minimum and[...]

randsys 331 real imaginary frequency damping (rad/sec) ratio -1.0408e+00 0.0000e+00 1.0408e+00 1.0000 -1.1390e+00 0.0000e+00 1.1390e+00 1.0000 -2.8050e-01 -1.4858e+00 1.5120e+00 0.1855 -2.8050e-01 1.4858e+00 1.5120e+00 0.1855 -1.5337e+00 0.0000e+00 1.5337e+00 1.0000 -5.5566e-01 -1.6958e+00 1.7845e+00 0.3114 -5.5566e-01 1.6958e+00 1.7845e+00 0.3114 [...]

332 Chapter 6. F unctio n Reference Frequency 0.1 1 10 0.01 100 Magnitude 0.01 0.1 1 0.001 10 sys1 (input 1) sys1 (input 2) sys2 Frequency 0.1 1 10 0.01 100 Phase (degrees) -80 -60 -40 -20 0 -100 20[...]

randsys 333 sys3 = randsys(6,1,1, { discrete, dt=5 } ) rifd(sys3)? Poles: radius angle (radians) 0.1652 0.0000 0.4580 0.0000 0.7872 -0.0394 0.7872 0.0394 0.9307 -0.2163 0.9307 0.2163 Zeros: radius angle (radians) 0.1623 0.0000 0.7981 0.0000 0.6658 0.0000 0.8921 -0.1981 0.8921 0.1981 15.4278 3.1416[...]

rifd 335 rifd Syn t ax [stat] = rifd(vec, { discrete,Hertz,d egrees } ) Pa r a m e t e r Li s t Inputs: vec complex v a lued vector (or Dynamic System - see b elow). Keyw ords: discret e Bo olean. V ector is to be in terpreted in the z domain, rather than the s domain. Default = 0. Hertz Bo olean. Display frequency units in Hertz on the s plane. !H[...]

336 Chapter 6. F unctio n Reference Examples sys1 = randsys(4,3,2, { stable } ) rifd(sys1) Poles: real imaginary frequency damping (rad/sec) ratio -1.6847e+00 0.0000e+00 1.6847e+00 1.0000 -2.4383e+00 0.0000e+00 2.4383e+00 1.0000 -8.7457e+00 0.0000e+00 8.7457e+00 1.0000 -1.5041e+01 0.0000e+00 1.5041e+01 1.0000 Zeros: ans (a scalar) = 0 # Compare rif[...]

rifd 337 sys2 = randsys(3,3,2, { stable,di screte } ) rifd(sys2) Poles: radius angle (radians) 0.9641 0.0000 0.2985 -0.5348 0.2985 0.5348 Zeros:[...]

sdtrsp 339 sdtrsp Syn t ax [v,y,u] = sdtrsp(Sys,dSys,w,tfinal,. .. { ord,intstep,cdelay } ) Pa r a m e t e r Li s t Inputs: Sys Con tinuous dynamic system. This is the upper system in the LFT. The initial states are used in the sim ula tion. dSys Digita l dynamic system. Lower system in the LFT. The initial sta tes are used in the simulation. w P D[...]

340 Chapter 6. F unctio n Reference Description Time doma in simulation of a sa mpled data interconnection. The applica ble closed lo o p system is illustra ted be low. Sys dSys - vw y u This is co nceptually the equiv alent of: v = starp ( Sys , dSys ) ∗ w . This functio n will handle interconnections in which the co nt inuous time signals ar e [...]

sdtrsp 341 triangle ho ld equiv alent is used. This is the same a linea rly connecting the samples a t the input to the hold. Some care is needed in the choice of cdelay and intstep. The default for intstep is ba sed on the contin uo us system eigenv alues and the minimum time spa cing in the input vector, w. It will be fo rced to be an integer div[...]

342 Chapter 6. F unctio n Reference inps = "ref" outs = "digP" cnx = ["digC"; ... # input to digP "[ ref;digP]"] # input to digC digclp = sysic(snm,inps,ou ts,cnx,dig P,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 s[...]

sdtrsp 343 0.2 0.4 0.6 0.8 0 1 -1 0 1 2 3 -2 4 sampled data calc. discrete calc.[...]

344 Chapter 6. F unctio n Reference See Also trsp[...]

simtransform 345 sim transform Syn t ax out = simtransform(sys,X) Pa r a m e t e r Li s t Inputs: sys Input system. This ma y b e a state-space system, pdm or constant matrix. X Similarit y transf orm. X must be inv ertible Outputs: out output system - in the sa me class as the input. Description Apply a similari t y tra nsform to sys .I f sys is a[...]

346 Chapter 6. F unctio n Reference -1.9737e+01 0.0000e+00 1.9737e+01 1.0000 -9.7569e+00 -2.9129e+01 3.0720e+01 0.3176 -9.7569e+00 2.9129e+01 3.0720e+01 0.3176 Zeros: real imaginary frequency damping (rad/sec) ratio -2.1990e+01 0.0000e+00 2.1990e+01 1.0000 -9.7962e+00 -2.9186e+01 3.0786e+01 0.3182 -9.7962e+00 2.9186e+01 3.0786e+01 0.3182 [a,b,c,d] [...]

simtransform 347 X0 0 0 0 Input Names ----------- Input 1 Output Names ------------ Output 1 System is continuous[...]

spectrad 349 sp ectrad Syn t ax out = spectrad(mat) Pa r a m e t e r Li s t Inputs: ma t Square matr ix or pdm . Outputs: out sp ectr al radius o f the input. Description Calculates the sp ectral radius (magnitude of the maximum eigen v alue) of the input matrix. Exam ple pdm1 = randpdm(3,2,2, { zeromea n } )? pdm1 (a pdm) = domain | Col 1 Col 2 --[...]

350 Chapter 6. F unctio n Reference 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 [...]

sresidualize 351 sresidualize Syn t ax sysout = sresidualize(sysin,ord) Pa r a m e t e r Li s t Inputs: sysin Input D ynamic System ord Order of sysout Outputs: sysout output Dynamic System Description Residualize the states o f sysin to the num b er sp ecified by ord . The last n x − ord states ( n x is the n umbe r of sta tes in sysin ) are re[...]

352 Chapter 6. F unctio n Reference sysout1 = sresidualize(sys1,3) sysout2 = truncate(sys1,3) fHz = logspace(0.01,100,10 0) 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], { logmagplo t[...]

sresidualize 353 Frequency 0.1 1 10 0.01 100 Magnitude 0.0001 0.001 0.01 1e-05 0.1 original system residualized system truncated system residualization error truncation error[...]

354 Chapter 6. F unctio n Reference See also rifd , simtransform , orderstate modalstate , truncate .[...]

starp 355 starp Syn t ax out = starp (upper,lower,dim1,dim 2,skipChks ) Pa r a m e t e r Li s t Inputs: upper Upper ob ject ( D ynamic System , constan t or pdm )i nt h e interconnection low er Low er ob ject in the in terconnection. dim1 The n umber o f outputs of upp er to b e connected as inputs to lower. (Optiona l - see descr iption for defa u[...]

356 Chapter 6. F unctio n Reference specified. This is equiv alent 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 m[...]

starp 357 # 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])*dau g(1,P) # Test this with a Proportional and PI controller K p=1 0 KpInt = makepoly([10,100],"s")/ makepoly([ 1,0],"s") ri[...]

358 Chapter 6. F unctio n Reference -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. Gno[...]

starp 359 Poles: real imaginary frequency damping (rad/sec) ratio -5.0000e+00 -8.0623e+00 9.4868e+00 0.5270 -5.0000e+00 8.0623e+00 9.4868e+00 0.5270 Zeros: real imaginary frequency damping (rad/sec) ratio -1.0000e+01 0.0000e+00 1.0000e+01 1.0000[...]

substr 361 substr Syn t ax littlestring = substr(bigstring,charinde x, { skipChk s } ) Pa r a m e t e r Li s t Inputs: bigstring Input string (a 1 × 1 s tring ma trix) ch arindex vector indexing the c haracters t o b e returned in the output . Keywords: skipChks Boolean specifying that syntax chec king is to be skipp ed. Outputs: littlestring Outp[...]

sysic 363 sysic Syn t ax [sys] = sysic(sysNames,sysInput s,sysOutp uts,sysCon nects,... subsys1,subsys2,...) Pa r a m e t e r Li s t Inputs: sysNa mes A vector of strings (of the sa me length as the num b er of subsystems) naming the subsystems. sysInputs A vector of strings naming the exo genous inputs to the final system. Each named input m ust [...]

364 Chapter 6. F unctio n Reference Matrices can a lso be included in the interconnection and are consider ed to b e consta nt gains when connected with Dynamic System s and con sidered to be constan t f or all domain v alues when connected with pdm s. sysic is able to connect linear systems to gether. F or more complete interconnection and sim ula[...]

trsp 365 trsp Syn t ax [y,uint] = trsp(Sys,u,tfinal, ord,intstep) Pa r a m e t e r Li s t Inputs: Sys Con tinuous dynamic system. The initial states ar e used in the simulation. u PDM. Input signal. tfinal Final time in the simulation (optiona l). Default = ma x t i m es p e c i fi e di nu . Keyw ords: ord Sc alar v alued. Specifies order of int[...]

366 Chapter 6. F unctio n Reference necessary . If int step is not specified the in tegration stepsize is d etermined from t he system eigen v alues and the minim um spac ing in the input signal, u . A zero or fi rst order discrete equiv alen t of the systems can b e sp ecified. The standard Xmath * op erator uses only a zero order discretizatio[...]

trsp 367 0.5 1 1.5 0 2 -1 -0.5 0 0.5 1 -1.5 1.5 output input[...]

368 Chapter 6. F unctio n 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,mark er=1 } ); gph2 = plot(uint,gph2, { line style=3,legend=["trsp calc.";... "* calc.";"input";"interpola ted input"[...]

trsp 369 0.5 1 1.5 0 2 -1 -0.5 0 0.5 1 -1.5 1.5 Time response calculation comparisons trsp calc. * calc. input interpolated input[...]

370 Chapter 6. F unctio n Reference # Now look at 1st order interpolation [y1trsp,u1int] = trsp(P,u, { ord=1 } ) gph3 = ctrlplot([y1trsp,u1int]) ; gph3 = plot(u,gph3, { line=0,mark er=1,legen d=... ["1st order interpolation";... "interpolated input";"input data"],title=... "Time response calculation comparisons&qu[...]

trsp 371 0.5 1 1.5 0 2 -1 -0.5 0 0.5 1 -1.5 1.5 Time response calculation comparisons 1st order interpolation interpolated input input data[...]

sresidualize 373 truncate Syn t ax sysout = truncate(sysin,ord) Pa r a m e t e r Li s t Inputs: sysin Input D ynamic System ord Or der of the truncated system: s ysout Outputs: sysout T runcated output Dynamic System Description The function truncate is cross-licensed from the mo del reduction to olb ox and has slightly more capabilities than those[...]

374 Chapter 6. F unctio n Reference Exam ple This example is iden tica l to tha t describ ed for sresidualize and compares the tw o m et hod s of m ode l re du ct i on . # Create a five state system for reduction. a = daug(-0.891334,[-1.208 57,0.79904 2;-0.79904 2,-1.2085 7],... -4.74685,-21.3013) b = [0.0262569;-0.189601;- 0.113729;0 .211465;-0 .5[...]

sresidualize 375 Frequency 0.1 1 10 0.01 100 Magnitude 0.0001 0.001 0.01 1e-05 0.1 original system residualized system truncated system residualization error truncation error[...]

376 Chapter 6. F unctio n Reference See also rifd , simtransform , sresidualize , orderstate modalstate .[...]

6.2. X µ SUBR OUTINES AND UTILITIES 377 6.2 X µ Subroutines and Utilities Several subroutines ma y also b e of interest to the user. These subroutines typically perform self cont ained part s of a calculat ion. They ma y be of interest t o those developing new robust con trol algorithms or wishing to study the calculation details of the algorit h[...]

blkbal 379 blkbal Syn t ax d = blkbal(M) Description Balances a square ma trix assuming only scalar blo cks. The Osb or ne metho d (growth rate: n 2 ) is used for large systems and the P erron method (gro wth rate: n 3 ) for smaller systems. The Perron method will exactl y calculate m u for p ositiv e matrices.[...]

hinfcalc 381 hinfcalc Syn t ax [X,Y,f,h,Ric fail,HX,HY,HXmin,HYmin ] = ... hinfcalc(p,nmeas,nco n,g,epr, { keyword s } )[...]

382 Chapter 6. F unctio n Reference Pa r a m e t e r Li s t Inputs: p Generalized interconnection structur e ( Dynamic Sys- tem ) nmeas measuremen t v ector dimension. ncon con trol v ector dimension. g H ∞ norm of sub o ptimal controller to b e ca lculated. Re- ferred to in the lit erature as gamma. epr T oler ance for deter mining when the Hami[...]

hinfcalc 383 Description F orm and solve th e Riccati equations for the H ∞ con trol problem. X and Y are the resulting Riccati solutio ns. THIS FUNCTION IS INTENDED ONL Y AS A SUBR OUTINE CALLED BY THE HINFSYN FUNCTION. *** NO ERROR C HECKIN G ***[...]

powermu 385 po w e r m u Syn t ax [lbnd,delta,errstat] = powermu(M,blk,rp,cp) Description Low er bound p ow er algorithm based on th e work of Andy Pac k ard. The v ector naming roughly corresponds to that in his t hesis. *** NO ERROR C HECKIN G ***[...]

riccati eig 387 riccati eig Syn t ax [x1,x2,stat,Heig min] = riccati eig(H,epp) Pa r a m e t e r Li s t Inputs: H Hamiltonian matrix. epp T olera nce for detecting proximit y of eigenv alues to the jω axis. Outputs: x1,x2 Basis v ector s for stable subspace. See description b elow. stat Statu s flag. 0 Stable subspace calculated. 1 F a ilure to d[...]

388 Chapter 6. F unctio n Reference v ariab le, X = x 2 x − 1 1 , is the stabilizing solution to the Ricca ti equation. If H has jω axis eigenv alues then no stabilizing solution exists and the function retur ns a failure status. If any eigenv alue of H is within epp o f the jω axis it is co nsidered to lie on the jω axis and no solution is fo[...]

riccati schur 389 riccati sc h ur Syn t ax [x1,x2,stat,Heig min] = riccati schur(H,epp) Pa r a m e t e r Li s t Inputs: H Hamiltonian matrix. epp T olera nce for detecting proximit y of eigenv alues to the jω axis. Outputs: x1,x2 Basis v ector s for stable subspace. See description b elow. stat Statu s flag. 0 Stable subspace calculated. 1 F a il[...]

390 Chapter 6. F unctio n Reference v ariab le, X = x 2 x − 1 1 , is the stabilizing solution to the Ricca ti equation. If H has jω axis eigenv alues then no stabilizing solution exists and the function retur ns a failure status. If any eigenv alue of H is within epp o f the jω axis it is co nsidered to lie on the jω axis and no solution is fo[...]

App endices A T ranslation Bet ween Ma tlab µ -T o ols and X µ This app endix outlines the functional equiv alences betw een the Ma tlab µ -T o ols and Xmath X µ . The ob ject ive is t o pro vide a smooth transiti on for users mo ving from µ -T o ols to X µ . W e will assume that the rea der is familiar with Ma tlab µ -T o ols and the genera[...]

392 APPEND ICES The functionally equiv alen t commands will b e listed, in each sub-section, for conv enience. A more detailed discussion is given to illustrate the mor e subtle differences in the mu and D - K iteration asp ects. Mor e imp ortantly , the Hima t demo is av ailable in bot h µ -T oo ls and X µ . F or a fast start on moving bet ween[...]

A. TRANSLA TION BETWE EN MA TLA B µ -TOOLS AND X µ 393 Description µ -T o ols F unction Xmath/X µ equiv alen t form system pck system decompose system unpck abcd form system nd2sys transfer f unction data ob jects form system zp2sys transfer f unction data ob jects decompose system sys2pss dynamic system data ob jects form system pss2sys dynami[...]

394 APPEND ICES Subblo cks: s electing input & outputs In µ -T o ols the function sel selects rows and columns fro m a va r y i n g matrix or inputs and outputs from a system matrix. In Xma th these can b e o btained by specifying row and column indexes. More flexibility of selecting par ts of a pdm can be obtained b y using the indexlist fun[...]

A. TRANSLA TION BETWE EN MA TLA B µ -TOOLS AND X µ 395 Note that the transpose and conjugate tran sp ose op erators are defined differen tly f or Ma tl ab and Xmat h. system / Dynamic System F unctions The following functions p erfo rm useful manipula tions to, or information a b out, the sta te of a system or Dynamic System . Description µ -T[...]

396 APPEND ICES Description µ -T o ols F unction Xmath/X µ equiv alen t peak norm pkvnorm norm , max absolute v alue vabs abs diagonal m atrix vdiag diagonal round down wards vfloor round round upw ar ds vceil round imagin ary part vimag imag real part vreal real complex c onjugate vconj conj norm vnorm norm determinant vdet det eigen v alues vei[...]

A. TRANSLA TION BETWE EN MA TLA B µ -TOOLS AND X µ 397 Miscellaneou s Utilities Several utilities a re pr ovided in µ -T o ols. These are subroutines used b y other µ -T ools funct ions wh ich ma y be of more general use. Description µ -T o ols function Xma th/X µ equiv alen t complex rando m num b er crand randpdm fit system to data fitsys [...]

398 APPEND ICES A.3 System R esp onse F unctions Creating Ti me Domain Signal s The Xmath pdm data ob ject allows the crea tion of time domain signals via standard and operators. Description µ -T o ols function Xma th/X µ equiv alen t cosine wa veform cos tr cos sine wa veform sin tr sin stair-step w av eform step tr gstep general w av eform sigg[...]

A. TRANSLA TION BETWE EN MA TLA B µ -TOOLS AND X µ 399 information. A.4 System In terconnection Simple interconnection has already been outlined in the a ugmentation section above. The mo re complicated in ter connection functions are a lmost identical in µ -T o o ls and X µ . The only r eal difference is the calling syntax of sysic . Descript[...]

400 APPEND ICES Description µ -T o ols function Xma th/X µ equiv alen t H 2 norm calcula tion h2norm h2norm H ∞ norm calcula tion hinfnorm hinfnorm H 2 controller syn thesis h2syn h2syn H ∞ controller syn thesis hinfsyn hinfsyn The ma jor syntactical diff erence is that the X µ fu nction s do not return the closed loop system. This is easil[...]

A. TRANSLA TION BETWE EN MA TLA B µ -TOOLS AND X µ 401 Description µ -T o ols function Xma th/X µ equiv alen t structured singula r v alue mu mu D scale deco ding unwrapd not required pe rturbatio n deco ding unwrapp not required blo ck norm calculations blknorm blknorm rational perturbat ion dypert mkpert random p er turbations randel randpert[...]

402 APPEND ICES g1g = freq(g1,omega) [mubnds,Dmagdata] = mu(g1g,blk) F rom this, frequency do main D -scales a re obtained, a nd a ra tional appr oximation is obtained via musynfit . [Dsys,Dinvsys] = musynfit(Dmagdata,blk,nmea s,nctrls,w eight,g1g) A difference b et ween the X µ and µ -T o ols implementations of musynfit is that, in the X µ cas[...]

A. TRANSLA TION BETWE EN MA TLA B µ -TOOLS AND X µ 403 The adv an tage of this is that in order to restart, or reproduce, an iterat ion, one need only sav e the previous con troller. The µ -T o ols approach requ ires saving the rational approximation to the previous D -scales. The controller is a more applicable data ob ject to sa ve an d the sa[...]

T echnical Support and Professional Ser vices Visit the following sections of the National Instruments Web site at ni.com for technical support an d professional services: • Support —Online technical support resources at ni.com/support include the following: – Self-Help Resources —For immediate answ ers and solu tions, visit the award-winni[...]