Control System Toolbox For Use with MATLAB ® Computation Visualization Programming User’s Guide Version 4.
How to Contact The MathWorks: ☎ 508-647-7000 Phone 508-647-7001 Fax The MathWorks, Inc. 24 Prime Park Way Natick, MA 01760-1500 Mail http://www.mathworks.com Web Anonymous FTP server Newsgroup PHONE FAX ✉ MAIL INTERNET ftp.mathworks.com comp.soft-sys.matlab @ support@mathworks.com suggest@mathworks.com bugs@mathworks.com doc@mathworks.com subscribe@mathworks.com service@mathworks.com info@mathworks.
Contents Preface Installation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Getting Started . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Typographic Conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Quick Start 1 LTI Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . MIMO Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
System Interconnections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-19 Control Design Tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-20 The Root Locus Design GUI . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-23 LTI Models 2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . LTI Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
LTI Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Generic Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Model-Specific Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Setting LTI Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Accessing Property Values Using get . . . . . . . . . . . . . . . . . . . . Direct Property Referencing . . . . . . . . . . . . . . . . . . . . . . . . . . .
Operations on LTI Models 3 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-2 Precedence and Property Inheritance . . . . . . . . . . . . . . . . . . 3-3 Extracting and Modifying Subsystems . . . . . . . . . . . . . . . . . . Referencing FRD Models Through Frequencies . . . . . . . . . . . . . Referencing Channels by Name . . . . . . . . . . . . . . . . . . . . . . . . . Resizing LTI Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Resampling of Discrete-Time Models . . . . . . . . . . . . . . . . . . . 3-27 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-28 Arrays of LTI Models 4 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . When to Collect a Set of Models in an LTI Array . . . . . . . . . . . Restrictions for LTI Models Collected in an Array . . . . . . . . . . Where to Find Information on LTI Arrays . . . . . . . . . . . . . . . .
Dimension Requirements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Special Cases for Operations on LTI Arrays . . . . . . . . . . . . . . Examples of Operations on LTI Arrays with Single LTI Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Examples: Arithmetic Operations on LTI Arrays and SISO Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Other Operations on LTI Arrays . . . . . . . . . . . . . . . . . . . . . . . .
Displaying Response Characteristics on a Plot . . . . . . . . . . . . . 6-9 Importing Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-11 Zooming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6-12 The LTI Viewer Menus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The File Menu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Importing a New Model into the LTI Viewer Workspace . .
Specifying the Simulink Model Portion for Analysis . . . . . . . . Adding Input Point or Output Point Blocks to the Diagram Removing Input Points and Output Points . . . . . . . . . . . . . Specifying Open- Versus Closed-Loop Analysis Models . . . Setting the Operating Conditions . . . . . . . . . . . . . . . . . . . . . . . Modifying the Block Parameters . . . . . . . . . . . . . . . . . . . . . . . . Performing Linear Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . .
The Root Locus Design GUI 8 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8-2 A Servomechanism Example . . . . . . . . . . . . . . . . . . . . . . . . . . . 8-4 Controller Design Using the Root Locus Design GUI . . . . . 8-6 Opening the Root Locus Design GUI . . . . . . . . . . . . . . . . . . . . . 8-6 Importing Models into the Root Locus Design GUI . . . . . . . . . . 8-7 Opening the Import LTI Design Model Window . . . . . . . . . .
Clearing Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8-46 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8-48 Design Case Studies 9 Yaw Damper for a 747 Jet Transport . . . . . . . . . . . . . . . . . . . . 9-3 Open-Loop Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9-6 Root Locus Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9-9 Washout Filter Design . .
Reliable Computations 10 Conditioning and Numerical Stability . . . . . . . . . . . . . . . . . 10-4 Conditioning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10-4 Numerical Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10-6 Choice of LTI Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10-8 State Space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10-8 Transfer Function . . . . .
Conversion to Transfer Function . . . . . . . . . . . . . . . . . . . Example 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Example 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Example 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Example 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Creation of Zero-Pole-Gain Models . . . . . . . . . . . . . . . . . .
Preface Installation . . . . . . . . . . . . . . . . . . . . . 3 Getting Started . . . . . . . . . . . . . . . . . . . 4 Typographic Conventions . . . . . . . . . . . . . . .
Preface MATLAB® has a rich collection of functions immediately useful to the control engineer or system theorist. Complex arithmetic, eigenvalues, root-finding, matrix inversion, and FFTs are just a few examples of MATLAB’s important numerical tools. More generally, MATLAB’s linear algebra, matrix computation, and numerical analysis capabilities provide a reliable foundation for control system engineering as well as many other disciplines.
Installation Installation Instructions for installing the Control System Toolbox can be found in the MATLAB Installation Guide for your platform. We recommend that you store the files from this toolbox in a directory named control off the main matlab directory. To determine if the Control System Toolbox is already installed on your system, check for a subdirectory named control within the main toolbox directory or folder. Five demonstration files are available. The demonstration M-file called ctrldemo.
Preface Getting Started If you are a new user, begin with Chapters 2 through 5 to learn: • How to specify and manipulate linear time-invariant models • How to analyze such models and plot their time and frequency response If you are an experienced toolbox user, see: • The New Features Guide for details on the latest release • Chapter 1 for an overview of some product features • Chapter 4 to learn about LTI arrays • Chapter 6 for an introduction to the LTI Viewer GUI • Chapter 8 for an introduction to the
Typographic Conventions Typographic Conventions To Indicate This Guide Uses Example Example code Monospace type To assign the value 5 to A, enter A = 5 Function names Monospace type The cos function finds the cosine of each array element. Function syntax Monospace type for text that must appear as shown. The magic function uses the syntax M = magic(n) Monospace italics for components you can replace with any variable. Keys Boldface Press the Return key.
Preface 6
1 Quick Start LTI Models . . . MIMO Systems . Model Conversion LTI Arrays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-3 1-4 1-5 1-6 LTI Properties . . . . . . . . . . . . . . . . . . . 1-7 Model Characteristics . . . . . . . . . . . . . . . 1-10 Operations on LTI Models . . . . . . . . . . . . . . 1-11 Continuous/Discrete Conversions . . . . . . . . . .
1 Quick Start This chapter provides a quick overview of some features of the Control Systems Toolbox.
LTI Models LTI Models You can specify linear time-invariant (LTI) systems as transfer function (TF) models, zero/pole/gain (ZPK) models, state-space (SS) models, or frequency response data (FRD) model. You can construct the corresponding models using the constructor functions.
1 Quick Start To create discrete-time models, append the sample time Ts to the previous calling sequences. sys sys sys sys = = = = tf(num,den,Ts) zpk(z,p,k,Ts) ss(a,b,c,d,Ts) frd(response,frequency,Ts) For more information, see “Discrete-Time Models” on page 2-20. For example, type sys = zpk(0.5,[–0.1 0.3],1,0.05) and MATLAB responds with Zero/pole/gain: (z–0.5) --------------(z+0.1) (z–0.3) Sampling time: 0.
LTI Models where each SISO entry is characterized by its numerator and denominator. Cell arrays provide an ideal means to specify the resulting arrays of numerators and denominators; see, “MIMO Transfer Function Models” on page 2-10 for more information. For example, num = {0.5,[1 1]} den = {[1 0],[1 2]} H = tf(num,den) % 1-by-2 cell array of numerators % 1-by-2 cell array of denominators creates the one-output/two-input transfer function H ( s ) = 0.
1 Quick Start MATLAB derives the zero/pole/gain representation of the transfer function h Zero/pole/gain: 1 ------(s+1)^2 LTI Arrays You can now create multidimensional arrays of LTI models and manipulate them as a single entity. LTI arrays are useful to perform batch analysis on an entire set of models. For more information, see Chapter 4, “Arrays of LTI Models.
LTI Properties LTI Properties In addition to the model data, the TF, ZPK, FRD, and SS objects can store extra information, such as the system sample time, time delays, and input or output names. The various pieces of information that can be attached to an LTI object are called the LTI properties. For information on LTI properties, type ltiprops See also “LTI Properties” on page 2-26 and “Time Delays” on page 2-45.
1 Quick Start Type: get(sys) a: -1 b: 1 c: 1 d: 0 e: [] StateName: {''} Ts: 0.5 InputDelay: 0 OutputDelay: 0 ioDelayMatrix: 0 InputName: {''} OutputName: {''} InputGroup: {0x2 cell} OutputGroup: {0x2 cell} Notes: {} UserData: [] You can also use set and get to access/modify LTI properties in a Handle Graphics® fashion; see “Setting LTI Properties” on page 2-30 for more information. For example, give names to the input and output of the SISO state-space model sys.
LTI Properties MATLAB returns ans = 3.
1 Quick Start Model Characteristics The Control System Toolbox contains commands to query such model characteristics as the I/O dimensions, poles, zeros, and DC gain. See “General Model Characteristics” on page 5-2 for more information. These commands apply to both continuous- and discrete-time model. Their LTI-based syntax is summarized below (with sys denoting an arbitrary LTI model).
Operations on LTI Models Operations on LTI Models You can perform simple matrix operations, such as addition, multiplication, or concatenation on LTI models. See Chapter 3, “Operations on LTI Models” for more information. Thanks to MATLAB object-oriented programming capabilities, these operations assume appropriate functionalities when applied to LTI models. For example, addition performs a parallel interconnection.
1 Quick Start extracts the subsystem mapping the first input to the third output. Note that row indices select the outputs while column indices select the inputs. Similarly, sys(3,1) = tf(1,[1 0]) redefines the relation between the first input and third output as an integrator.
Continuous/Discrete Conversions Continuous/Discrete Conversions The commands c2d, d2c, and d2d perform continuous to discrete, discrete to continuous, and discrete to discrete (resampling) conversions, respectively. sysd = c2d(sysc,Ts) % discretization w/ sample period Ts sysc = d2c(sysd) % equivalent continuous-time model sysd1= d2d(sysd,Ts) % resampling at the period Ts See “Continuous/Discrete Conversions of LTI Models” on page 3-20 for more information.
1 Quick Start Time and Frequency Response The following commands produce various time and frequency response plots for LTI models (see “Time and Frequency Response” on page 5-9 for more information).
Time and Frequency Response produces the Bode plot shown below. To superimpose and compare the responses of several LTI systems, use the syntax bode(sys1,sys2,sys3,...) You can also control the plot style by specifying a color/linestyle/marker for each system, much as with the plot command; see “Plotting and Comparing Multiple Systems” on page 5-13 for more information.
1 Quick Start These commands automatically determine an appropriate simulation horizon or frequency range based on the system dynamics. To override the default range, type step(sys,tfinal) bode(sys,{wmin,wmax}) 1-16 % final time = tfinal % freq.
The LTI Viewer The LTI Viewer You can also analyze time and frequency domain responses using the LTI Viewer; see “The LTI Viewer” on page 6-1 for more information.
1 Quick Start For more detail on the use of the LTI Viewer and how it can be integrated into a Simulink diagram, see Chapter 6, “The LTI Viewer.
System Interconnections System Interconnections You can derive LTI models for various system interconnections ranging from simple series connections to more complex block diagrams; see “Model Interconnection Functions” on page 3-16 for more information. Related commands include append(sys1,sys2,...
1 Quick Start Control Design Tools The Control System Toolbox supports three mainstream control design methodologies: gain selection from root locus, pole placement, and linear-quadratic-Gaussian (LQG) regulation. The first two techniques are covered by the rlocus and place commands. The LQG design tools include commands to compute the LQ-optimal state-feedback gain (lqr, dlqr, and lqry), to design the Kalman filter (kalman), and to form the resulting LQG regulator (lqgreg).
Control Design Tools Plant + u+d d + 100 -----------------------------s 2 + s + 100 y u + F(s) yn + n LQG regulator Figure 1-1: Simple Regulation Loop The following commands design the optimal LQG regulator F ( s ) for this problem.
1 Quick Start To validate the design, close the loop with feedback and compare the open- and closed-loop impulse responses with impulse. % Close loop clsys = feedback(sys,F,+1) % Note positive feedback % Open- vs. closed-loop impulse responses impulse(sys,'r--',clsys,'b-') Impulse Response 10 8 6 Amplitude 4 2 0 −2 −4 −6 −8 0 2 4 6 Time (sec.
The Root Locus Design GUI The Root Locus Design GUI You can also design a compensator using the Root Locus Design Graphical User Interface (GUI). See Chapter 8, “The Root Locus Design GUI” for more information.
1 Quick Start 1-24
2 LTI Models Introduction . . . . . . . . . . . . . . . . . . . . 2-2 Creating LTI Models . . . . . . . . . . . . . . . . 2-8 LTI Properties . . . . . . . . . . . . . . . . . . . 2-26 Model Conversion . . . . . . . . . . . . . . . . . 2-42 Time Delays . . . . . . . . . . . . . . . . . . . . 2-45 Simulink Block for LTI Systems . . . . . . . . . . . 2-57 References . . . . . . . . . . . . . . . . . . . . .
2 LTI Models Introduction The Control System Toolbox offers extensive tools to manipulate and analyze linear time-invariant (LTI) models. It supports both continuous- and discrete-time systems. Systems can be single-input/single-output (SISO) or multiple-input/multiple-output (MIMO). In addition, you can store several LTI models in an array under a single variable name. See Chapter 4, “Arrays of LTI Models” for information on LTI arrays.
Introduction Using LTI Models in the Control System Toolbox You can manipulate TF, SS, and ZPK models using the arithmetic and model interconnection operations described in Chapter 3, “Operations on LTI Models,” and analyze them using the model analysis functions, such as bode and step, described in Chapter 5, “Model Analysis Tools.” FRD models can be manipulated and analyzed in much the same way you analyze the other model types, but analysis is restricted to frequency-domain methods.
2 LTI Models Creating an LTI Object: An Example An LTI object of the type TF, ZPK, SS, or FRD is created whenever you invoke the corresponding constructor function, tf, zpk, ss, or frd. For example, P = tf([1 2],[1 1 10]) creates a TF object, P, that stores the numerator and denominator coefficients of the transfer function s+2 P ( s ) = -------------------------s 2 + s + 10 See “Creating LTI Models” on page 2-8 for methods for creating all of the LTI object types.
Introduction Precedence Rules Operations like addition and commands like feedback operate on more than one LTI model at a time. If these LTI models are represented as LTI objects of different types (for example, the first operand is TF and the second operand is SS), it is not obvious what type (for example, TF or SS) the resulting model should be. Such type conflicts are resolved by precedence rules. Specifically, TF, ZPK, SS, and FRD objects are ranked according to the precedence hierarchy.
2 LTI Models For example, the parallel connection of two LTI systems sys1 and sys2 can be expressed as sys = sys1 + sys2 because parallel connection amounts to adding the transfer matrices. Similarly, subsystems of a given LTI model sys can be extracted using matrix-like subscripting.
Introduction Table 2-1: Creating LTI Models or Getting Data From Them (Continued) Command Description set Set LTI model properties. size Get output/input/array dimensions or model order. ss Create a state-space model. ssdata, dssdata Retrieve state-space data (respectively, descriptor state-space data) or convert it to cell array format. tf Create a transfer function. tfdata Retrieve transfer function data. zpk Create a zero-pole-gain model. zpkdata Retrieve zero-pole-gain data.
2 LTI Models Creating LTI Models The functions tf, zpk, ss, and frd create transfer function models, zero-pole-gain models, state-space models, and frequency response data models, respectively. These functions take the model data as input and produce TF, ZPK, SS, or FRD objects that store this data in a single MATLAB variable. This section shows how to create continuous or discrete, SISO or MIMO LTI models with tf, zpk, ss, and frd.
Creating LTI Models where num and den are row vectors listing the coefficients of the polynomials n ( s ) and d ( s ), respectively, when these polynomials are ordered in descending powers of s. The resulting variable h is a TF object containing the numerator and denominator data.
2 LTI Models MIMO Transfer Function Models MIMO transfer functions are two-dimensional arrays of elementary SISO transfer functions. There are several ways to specify MIMO transfer function models, including: • Concatenation of SISO transfer function models • Using tf with cell array arguments Consider the rational transfer matrix s–1 -----------s+1 . H(s) = s+2 --------------------------2 s + 4s + 5 You can specify H ( s ) by concatenation of its SISO entries.
Creating LTI Models For example, for the rational transfer matrix H ( s ) , the two cell arrays N and D should contain the row-vector representations of the polynomial entries of N( s ) = s–1 s+2 D( s ) = s+1 s2 + 4s + 5 You can specify this MIMO transfer matrix H ( s ) by typing N = {[1 –1];[1 2]}; % cell array for N(s) D = {[1 1];[1 4 5]}; % cell array for D(s) H = tf(N,D) MATLAB responds with Transfer function from input to output...
2 LTI Models while E = tf creates an empty transfer function. Zero-Pole-Gain Models This section explains how to specify continuous-time SISO and MIMO zero-pole-gain models. The specification for discrete-time zero-pole-gain models is a simple extension of the continuous-time case. See “Discrete-Time Models” on page 2-20. SISO Zero-Pole-Gain Models Continuous-time SISO zero-pole-gain models are of the form ( s – z 1 ) ...
Creating LTI Models produces Zero/pole/gain: –2 s -------------------(s–2) (s^2 – 2s + 2) You can also specify zero-pole-gain models as rational expressions in the Laplace variable s by: 1 Defining the variable s as a ZPK model s = zpk('s') 2 Entering the transfer function as a rational expression in s. For example, once s is defined with zpk, H = –2s/((s – 2)*(s^2 + 2*s + 2)); returns the same ZPK model as h = zpk([0], [2 –1–i –1+i ], –2); Note: You need only define the ZPK variable s once.
2 LTI Models where • Z is the p-by-m cell array of zeros (Z{i,j} = zeros of H ij ( s ) ) • P is the p-by-m cell array of poles (P{i,j} = poles of H ij ( s ) ) • K is the p-by-m matrix of gains (K(i,j) = gain of H ij ( s ) ) For example, typing Z = {[],–5;[1–i 1+i] []}; P = {0,[–1 –1];[1 2 3],[]}; K = [–1 3;2 0]; H = zpk(Z,P,K) creates the two-input/two-output zero-pole-gain model H(s) = –1 -----s 3(s + 5) -------------------2 (s + 1) 2 ( s 2 – 2s + 2 ) ---------------------------------------------
Creating LTI Models Use the command ss to create state-space models sys = ss(A,B,C,D) For a model with Nx states, Ny outputs, and Nu inputs • A is an Nx-by-Nx real-valued matrix. • B is an Nx-by-Nu real-valued matrix. • C is an Ny-by-Nx real-valued matrix. • D is an Ny-by-Nu real-valued matrix. This produces an SS object sys that stores the state-space matrices A, B, C, and D. For models with a zero D matrix, you can use D = 0 (zero) as a shorthand for a zero matrix of the appropriate dimensions.
2 LTI Models to which MATLAB responds a = x1 x2 x1 0 –5.00000 x1 x2 u1 0 3.00000 y1 x1 0 y1 u1 0 x2 1.00000 –2.00000 b = c = x2 1.00000 d = In addition to the A, B, C, and D matrices, the display of state-space models includes state names, input names, and output names. Default names (here, x1, x2, u1, and y1) are displayed whenever you leave these unspecified. See “LTI Properties” on page 2-26 for more information on how to specify state, input, or output names.
Creating LTI Models it is often desirable to work with the descriptor form when the E matrix is poorly conditioned with respect to inversion. The function dss is the counterpart of ss for descriptor state-space models. Specifically, sys = dss(A,B,C,D,E) creates a continuous-time DSS model with matrix data A,B,C,D,E. For example, consider the dynamical model dω J -------- + Fω = T dt y=ω with vector ω of angular velocities.
2 LTI Models Here w i is the input frequency of each sinusoid, i = 1 ... n, and G(w) = G ( w ) exp ( j ∠G ( w ) ) . The steady state output response of this system satisfies y i ( t ) = G ( w i ) sin ( w i t + ∠G ( w i ) ) ; i = 1…n A frequency response data (FRD) object is a model form you can use to store frequency response data (complex frequency response, along with a corresponding vector of frequency points) that you obtain either through simulations or experimentally.
Creating LTI Models For example, the MAT-file LTIexamples.mat contains a frequency vector freq, and a corresponding complex frequency response data vector respG. To load this frequency-domain data and construct an FRD model, type load LTIexamples sys = frd(respG,freq) Continuous-time frequency response with 1 output and 1 input at 5 frequency points. From input 1 to: Frequency(rad/s) ---------------1 2 4 5 output 1 -------–0.812505 –0.000312i –0.092593 –0.462963i –0.075781 –0.001625i –0.043735 –0.
2 LTI Models Discrete-Time Models Creating discrete-time models is very much like creating continuous-time models, except that you must also specify a sampling period or sample time for discrete-time models. The sample time value should be scalar and expressed in seconds. You can also use the value –1 to leave the sample time unspecified. To specify discrete-time LTI models using tf, zpk, ss, or frd, simply append the desired sample time value Ts to the list of inputs.
Creating LTI Models produces Transfer function: z – 0.2 ------z + 0.3 Sampling time: unspecified Note: Do not simply omit Ts in this case. This would make h a continuous-time transfer function. If you forget to specify the sample time when creating your model, you can still set it to the correct value by reassigning the LTI property Ts. See “Sample Time” on page 2-34 for more information on setting this property.
2 LTI Models Similarly, z = zpk('z', 0.1); H = [z/(z+0.1)/(z+0.2) ; (z^2+0.2*z+0.1)/(z^2+0.2*z+0.01)] produces the single-input, two-output ZPK model Zero/pole/gain from input to output... z #1: --------------(z+0.1) (z+0.2) #2: (z^2 + 0.2z + 0.1) -----------------(z+0.1)^2 Sampling time: 0.1 Note that: • The syntax z = tf('z') is equivalent to z = tf('z',–1) and leaves the sample time unspecified. The same applies to z = zpk('z').
Creating LTI Models produces the transfer function z + 0.5 --------------------------2 z + 2z + 3 which differs from H ( z – 1 ) by a factor z . To avoid such convention clashes, the Control System Toolbox offers a separate function filt dedicated to the DSP-like specification of transfer functions. Its syntax is h = filt(num,den) for discrete transfer functions with unspecified sample time, and h = filt(num,den,Ts) to further specify the sample time Ts.
2 LTI Models Data Retrieval The functions tf, zpk, ss, and frd pack the model data and sample time in a single LTI object. Conversely, the following commands provide convenient data retrieval for any type of TF, SS, or ZPK model sys, or FRD model sysfr.
Creating LTI Models displays the coefficients of the numerator and denominator of the first input channel. ans = 0 ans = 1 1 –1 2 10 Note that the same result is obtained using H.num{1,1}, H.den{1,1} See “Direct Property Referencing” on page 2-33 for more information about this syntax.
2 LTI Models LTI Properties The previous section shows how to create LTI objects that encapsulate the model data and sample time. You also have the option to attribute additional information, such as the input names or notes on the model history, to LTI objects. This section gives a complete overview of the LTI properties, i.e., the various pieces of information that can be attached to the TF, ZPK, SS, and FRD objects. Type help ltiprops for online help on available LTI properties.
LTI Properties The sample time property Ts keeps track of the sample time (in seconds) of discrete-time systems. By convention, Ts is 0 (zero) for continuous-time systems, and Ts is –1 for discrete-time systems with unspecified sample time. Ts is always a scalar, even for MIMO systems. The InputDelay, OutputDelay, and ioDelayMatrix properties allow you to specify time delays in the input or output channels, or for each input/output pair. Their default value is zero (no delay).
2 LTI Models Model-Specific Properties The remaining LTI properties are specific to one of the four model types (TF, ZPK, SS, or FRD). For single LTI models, these are summarized in the following four tables. The property values differ for LTI arrays. See set on page 11-193 for more information on these values.
LTI Properties Table 2-7: SS-Specific Properties (Continued) Property Name Description Property Value e Descriptor E matrix 2-D real matrix StateName State names Cell vector of strings Table 2-8: FRD-Specific Properties Property Name Description Property Value Frequency Frequency data points Real-valued vector ResponseData Frequency response Complex-valued multidimensional array Units Units for frequency String 'rad/s' or 'Hz' Most of these properties are dedicated to storing the mod
2 LTI Models Setting LTI Properties There are three ways to specify LTI property values: • You can set properties when creating LTI models with tf, zpk, ss, or frd. • You can set or modify the properties of an existing LTI model with set. • You can also set property values using structure-like assignments. This section discusses the first two options. See “Direct Property Referencing” on page 2-33 for details on the third option.
LTI Properties For example, you can specify the delay directly when you create the model, and then use the set command to assign InputName, OutputName, and Notes to sys. sys = tf(1,[1 1],'Inputdelay',0.3); set(sys,'inputname','energy','outputname','temperature',... 'notes','A simple heater model') Finally, you can also use the set command to obtain a listing of all settable properties for a given LTI model type, along with valid values for these properties.
2 LTI Models where the string PropertyName is either the full property name, or any abbreviation with enough characters to identify the property uniquely. For example, typing h = tf(100,[1 5 100],'inputname','voltage',... 'outputn','current',... 'notes','A simple circuit') get(h,'notes') produces ans = 'A simple circuit' To display all of the properties of an LTI model sys (and their values), use the syntax get(sys).
LTI Properties Direct Property Referencing An alternative way to query/modify property values is by structure-like referencing. Recall that LTI objects are basic MATLAB structures except for the additional flag that marks them as TF, ZPK, SS, or FRD objects (see page 2-3). The field names for LTI objects are the property names, so you can retrieve or modify property values with the structure-like syntax. PropertyValue = sys.PropertyName% gets property value sys.
2 LTI Models produces ans = 'u' Any valid syntax for structures extends to LTI objects. For example, given the TF model h ( p ) = 1 ⁄ p h = tf(1,[1,0],'variable','p'); you can reset the numerator to p + 2 by typing h.num{1} = [1 2]; or equivalently, with h.num{1}(2) = 2; Additional Insight into LTI Properties By reading this section, you can learn more about using the Ts, InputName, OutputName, InputGroup, and OutputGroup LTI properties through a set of examples.
LTI Properties This sets the Ts property to the value 0.5, as is confirmed by h.Ts ans = 0.5000 For continuous-time models, the sample time property Ts is 0 by convention. For example, type h = tf(1,[1 0]); get(h,'Ts') ans = 0 To leave the sample time of a discrete-time LTI model unspecified, set Ts to – 1. For example, h = tf(1,[1 –1],–1) produces Transfer function: 1 ----z – 1 Sampling time: unspecified The same result is obtained by using the Variable property.
2 LTI Models Note that tf(0.1,[1 –1],0.1) + tf(1,[1 0.5],0.5) returns an error message. ??? Error using ==> lti/plus In SYS1+SYS2, both models must have the same sample time. Caution: Resetting the sample time of a continuous-time LTI model sys from zero to a nonzero value does not discretize the original model sys. The command set(sys,'Ts',0.1) only affects the Ts property and does not alter the remaining model data.
LTI Properties This produces Transfer function from input "thrust" to output "velocity": 1 -----p + 10 Note how the display reflects the input and output names and the variable selection. In the MIMO case, use cell vectors of strings to specify input or output channel names. For example, type num = {3 , [1 2]}; den = {[1 10] , [1 0]}; H = tf(num,den); % H(s) has one output and two inputs set(H,'inputname',{'temperature' ; 'pressure'}) The specified input names appear in the display of H.
2 LTI Models To see how input and output groups (I/O groups) work: 1 Create a random state-space model with one state, three inputs, and three outputs. 2 Assign the first two inputs to a group named controls, the first output to a group named temperature, and the last two outputs to a group named measurements.
LTI Properties To do this, type h = rss(1,3,3); set(h, 'InputGroup',{[1 2] 'controls'}) set(h, 'OutputGroup', {[1] 'temperature'; [2 3] 'measurements'}) h and MATLAB returns a state-space model of the following form. a = x1 x1 –0.64884 x1 u1 0.12533 y1 y2 y3 x1 1.1909 1.1892 0 y1 y2 y3 u1 0.32729 0 0 b = u2 0 u3 0 u2 0 0 2.1832 u3 –0.1364 0 0 c = d = I/O Groups: Group Name controls temperature measurements I/O I O O Channel(s) 1,2 1 2,3 Continuous-time model.
2 LTI Models Notice that the middle column of the I/O group listing indicates whether the group is an input group (I) or an output group (O). In general, to specify M input groups (or output groups), you need an M-by-2 cell array organized as follows.
LTI Properties You can use regular cell array syntax for accessing or modifying I/O group components. For example, to delete the first output group, temperature, type h.OutputGroup(1,:) = [] ans = [1x2 double] 'measurements' Similarly, you can add or delete channels from an existing input or output group. Recalling that input group channels are stored in the first column of the corresponding cell array, to add channel three to the input group controls, type h.inputgroup{1,1} = [h.
2 LTI Models Model Conversion There are four LTI model types you can use with the Control System Toolbox: TF, ZPK, SS, and FRD. This section shows how to convert models from one type to the other. Explicit Conversion Model conversions are performed by tf, ss, zpk, and frd.
Model Conversion for continuous-time models, and –1 H ( z ) = D + C ( zI – A ) B for discrete-time models. Automatic Conversion Some algorithms operate only on one type of LTI model. For example, the algorithm for zero-order-hold discretization with c2d can only be performed on state-space models. Similarly, commands like tfdata expect one particular type of LTI models (TF). For convenience, such commands automatically convert LTI models to the appropriate or required model type.
2 LTI Models balancing automatically when you convert any TF or ZPK model to state space using ss. • Conversions to the transfer function representation using tf may incur a loss of accuracy. As a result, the transfer function poles may noticeably differ from the poles of the original zero-pole-gain or state-space model. • Conversions to state space are not uniquely defined in the SISO case, nor are they guaranteed to produce a minimal realization in the MIMO case.
Time Delays Time Delays Using the ioDelayMatrix, InputDelay, and OutputDelay properties of LTI objects, you can specify delays in both continuous- and discrete-time LTI models. With these properties, you can, for example, represent: • LTI models with independent delays for each input/output pair. For example, the continuous-time model with transfer function H(s) = e –0.1s 10 2 --s e e – 0.3s – 0.
2 LTI Models • Interconnections of continuous-time delay systems as long as the resulting transfer function from input j to output i is of the form exp ( – sτ ij ) h ij ( s ) where h ij ( s ) is a rational function of s • Padé approximation of time delays (pade) Specifying Input/Output Delays Using the ioDelayMatrix property, you can specify frequency-domain models with independent delays in each entry of the transfer function.
Time Delays This creates the LTI model with the following transfer function. H(s) = exp ( – sτ ij ) r ij ( s ) Here r ij ( s ) is the ( i, j ) entry of –1 R ( s ) = D + C ( sI – A ) B Note: State-space models with I/O delays have only a frequency-domain interpretation. They cannot, in general, be described by state-space equations with delayed inputs and outputs. Distillation Column Example This example is adapted from [2] and illustrates the use of I/O delays in process modeling.
2 LTI Models Enriched vapor Cooling water Condensate Distillate Feed Reflux Reboiler Vapor Steam flow Bottom liquid Bottom products Figure 2-3: Distillation Column Schematically, the distillation process functions as follows: • Steam flows into the reboiler and vaporizes the bottom liquid. This vapor is reinjected into the column and mixes with the feed • Methanol, being more volatile than water, tends to concentrate in the vapor moving upward.
Time Delays The regulated output variables are: • Percentage X D of methanol in the distillate • Percentage X B of methanol in the bottom products. The goal is to maximize X D by adjusting the reflux flow rate R and the steam flow rate S in the reboiler. To obtain a linearized model around the steady-state operating conditions, the transient responses to pulses in steam and reflux flow are fitted by first-order plus delay models.
2 LTI Models The resulting TF model is displayed as Transfer function from input "R" to output... 12.8 Xd: exp(–1*s) * ---------16.7 s + 1 Xb: 6.6 exp(–7*s) * ---------10.9 s + 1 Transfer function from input "S" to output... –18.9 Xd: exp(–3*s) * -------21 s + 1 Xb: –19.4 exp(–3*s) * ---------14.
Time Delays Note that the 0.1 second delay is on the input in the first model, and on the output in the second model. InputDelay and OutputDelay Properties When the state trajectory is of interest, you should use the InputDelay and OutputDelay properties to distinguish between delays on the inputs and delays on the outputs in state-space models. For example, you can accurately specify the two models above by M1 = ss(–1,1,1,0,'inputdelay',0.1) M2 = ss(–1,1,1,0,'outputdelay',0.
2 LTI Models The resulting model is displayed as Transfer function from input to output... 1 #1: exp(–0.1*s) * s #2: 2 exp(–0.1*s) * ----s + 1 By comparison, to produce an equivalent transfer function using the ioDelayMatrix property, you would need to type H = [1/s ; 2/(s+1)]; H.iodelay = [0.1 ; 0.1]; Notice that the 0.1 second delay is repeated twice in the I/O delay matrix. More generally, for a TF, ZPK, or FRD model with input delays [ α 1, ..., α m ] and output delays [ β 1, ...
Time Delays produces the discrete-time transfer function Transfer function: 1 z^(–3) * ----------------z^2 + 0.5 z + 0.2 Sampling time: 0.1 Notice the z^(–3) factor reflecting the three-sampling-period delay on the input. Mapping Discrete-Time Delays to Poles at the Origin Since discrete-time delays are equivalent to additional poles at z = 0 , they can be easily absorbed into the transfer function denominator or the state-space equations.
2 LTI Models absorbs the input delay in H1 into the transfer function denominator to produce the third-order transfer function Transfer function: 1 --------z^3 – z^2 Sampling time: 1 Note that H2.inputdelay now returns 0 (zero). Retrieving Information About Delays There are several ways to retrieve time delay information from a given LTI model sys: • Use property display commands to inspect the values of the ioDelayMatrix, InputDelay, and OutputDelay properties. For example, sys.
Time Delays channel delays. The resulting model has a minimum number of delays. When this minimization takes place: • All or part of the I/O delay matrix is absorbed into the input and output delay vectors. This minimizes the total number of I/O delays. • The input delays are transferred to output delays (or vice-versa), so as to minimize the overall number of input and output channel delays.
2 LTI Models type pade(H,[],[],[1 1;2 1]) where H is the TF representation of H ( s ) defined in the “Distillation Column Example” on page 2-47. This command produces a rational transfer function. Transfer function from input "R" to output... –12.8 s + 25.6 Xd: --------------------16.7 s^2 + 34.4 s + 2 Xb: 6.6 s^2 – 5.657 s + 1.616 --------------------------------------10.9 s^3 + 10.34 s^2 + 3.527 s + 0.2449 Transfer function from input "S" to output... 18.9 s – 12.
Simulink Block for LTI Systems Simulink Block for LTI Systems You can incorporate LTI objects into Simulink diagrams using the LTI System block shown below. This mask is linked to an LTI block in a Simulink diagram. Double-click on the block in your Simulink diagram to display or modify model information. The LTI System block can be accessed either by typing ltiblock at the MATLAB prompt or by selecting Control System Toolbox from the Blocksets and Toolboxes section of the main Simulink library.
2 LTI Models transfer functions or zero-pole-gain models, as it depends on the choice of state coordinates used by the realization algorithm. As a result, you cannot enter nonzero initial states when you supply TF or ZPK models to LTI blocks in a Simulink diagram. Note: • For MIMO systems, the input delays stored in the LTI object must be either all positive or all zero. • LTI blocks in a Simulink diagram cannot be used for FRD models or LTI arrays.
References References [1] Dorf, R.C. and R.H. Bishop, Modern Control Systems, Addison-Wesley, Menlo Park, CA, 1998. [2] Wood, R.K. and M.W. Berry, “Terminal Composition Control of a Binary Distillation Column,” Chemical Engineering Science, 28 (1973), pp. 1707-1717.
2 LTI Models 2-60
3 Operations on LTI Models Introduction . . . . . . . . . . . . . . . . . . . . 3-2 Precedence and Property Inheritance . . . . . . . . 3-3 Extracting and Modifying Subsystems . . Referencing FRD Models Through Frequencies Referencing Channels by Name . . . . . . . Resizing LTI Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-5 3-7 3-8 3-9 Arithmetic Operations . . . Addition and Subtraction . . . Multiplication . . . . . . . .
3 Operations on LTI Models Introduction You can perform basic matrix operations such as addition, multiplication, or concatenation on LTI models. Such operations are “overloaded,” which means that they use the same syntax as they do for matrices, but are adapted so as to apply to the LTI model context. These overloaded operations and their interpretation in this context are discussed in this chapter. You can read about discretization methods in this chapter as well.
Precedence and Property Inheritance Precedence and Property Inheritance You can apply operations to LTI models of different types. The resulting type is then determined by the rules discussed in “Precedence Rules” on page 2-5. For example, if sys1 is a transfer function and sys2 is a state-space model, then the result of their addition sys = sys1 + sys2 is a state-space model, since state-space models have precedence over transfer function models.
3 Operations on LTI Models • In general, when two LTI models sys1 and sys2 are combined using operations such as +, *, [,], [;], append, and feedback, the resulting model inherits its I/O names and I/O groups from sys1 and sys2. However, conflicting I/O names or I/O groups are not inherited. For example, the InputName property for sys1 + sys2 is left unspecified if sys1 and sys2 have different InputName property values.
Extracting and Modifying Subsystems Extracting and Modifying Subsystems Subsystems relate subsets of the inputs and outputs of a system. The transfer matrix of a subsystem is a submatrix of the system transfer matrix. For example, if sys is a system with two inputs, three outputs, and I/O relation y = Hu then H ( 3, 1 ) gives the relation between the first input and third output. y 3 = H ( 3,1 ) u 1 Accordingly, we use matrix-like subindexing to extract this subsystem.
3 Operations on LTI Models Note: • sys, the LTI model that has had a portion reassigned, retains its original model type (TF, ZPK, SS, or FRD) regardless of the model type of NewSubSys. • If NewSubSys is an FRD model, then sys must also be an FRD model. Furthermore, their frequencies must match. • Subsystem assignment does not reassign any I/O names or I/O group names of NewSubSys that are already assigned to NewSubSys. • Reassigning parts of a MIMO state-space model generally increases its order.
Extracting and Modifying Subsystems To extract the transfer function T 11 ( s ) from the first input to the first output, type T(1,1) Transfer function: 1 ------s + 0.1 Next reassign T 11 ( s ) to 1 ⁄ ( s + 0.5 ) and modify the second input channel of T by typing T(1,1) = tf(1,[1 0.5]); T(:,2) = [ 1 ; tf(0.4,[1 0]) ] Transfer function from input 1 to output... 1 #1: ------s + 0.5 #2: s – 1 ------------s^2 + 2 s + 2 Transfer function from input 2 to output... #1: 1 #2: 0.
3 Operations on LTI Models I/O channel or group names) as a keyword. There are two ways you can specify FRD models using frequencies: • Using integers to index into the frequency vector of the FRD model • Using a Boolean (logical) expression to specify desired frequency points in an FRD model For example, if sys is an FRD model with five frequencies, (e.g., sys.Frequency=[1 1.1 1.2 1.3 1.4]), then you can create a new FRD model sys2 by indexing into the frequencies of sys as follows.
Extracting and Modifying Subsystems is equivalent to sys(1,[2 4 5] Similarly, if pressure is the name assigned to an output channel of the LTI model sys, then sys('pressure',1) = tf(1, [1 1]) reassigns the subsystem from the first input of sys to the output labeled pressure. You can reference a set of channels by input or output name by using a cell array of strings for the names.
3 Operations on LTI Models For state-space models, both concatenation and subsystem assignment increase the model order because they assume that sys and h have independent states. If you intend to keep the same state matrix and merely update the input-to-state or state-to-output relations, use set instead and modify the corresponding state-space data directly.
Arithmetic Operations Arithmetic Operations You can apply almost all arithmetic operations to LTI models, including those shown below. Operation Description + Addition – Subtraction * Multiplication / Right matrix divide \ Left matrix divide inv Matrix inversion ' Pertransposition .
3 Operations on LTI Models represents the parallel interconnection shown below. sys1 y1 + u y + sys2 y2 sys If sys1 and sys2 are two state-space models with data A 1, B 1, C 1, D 1 and A 2, B 2, C 2, D 2 , the state-space data associated with sys1 + sys2 is A1 0 0 A2 , B1 B2 , C1 C 2 , D1 + D2 Scalar addition is also supported and behaves as follows: if sys1 is MIMO and sys2 is SISO, sys1 + sys2 produces a system with the same dimensions as sys1 whose ijth entry is sys1(i,j) + sys2.
Arithmetic Operations Multiplication Multiplication of two LTI models connects them in series. Specifically, sys = sys1 * sys2 returns an LTI model sys for the series interconnection shown below. v u sys2 sys1 y Notice the reverse orders of sys1 and sys2 in the multiplication and block diagram.
3 Operations on LTI Models The resulting inverse model is of the same type as sys. Related operations include: • Left division sys1\sys2, which is equivalent to inv(sys1)*sys2 • Right division sys1/sys2, which is equivalent to sys1*inv(sys2) For a state-space model sys with data A, B, C, D , inv(sys) is defined only when D is a square invertible matrix, in which case its state-space data is –1 A – BD C , BD –1 , –1 –D C , D –1 Transposition You can transpose an LTI model sys using sys.
Arithmetic Operations You can use pertransposition to obtain the Hermitian (conjugate) transpose of the frequency response of a given system. The frequency response of the pertranspose of H ( s ), G ( s ) = [ H ( – s ) ] T , is the Hermitian transpose of the frequency response of H ( s ): G ( jw ) = H ( jw ) H .
3 Operations on LTI Models Model Interconnection Functions The Control System Toolbox provides a number of functions to help with the model building process. These include model interconnection functions to perform I/O concatenation ([,], [;], and append), general parallel and series connections (parallel and series), and feedback connections (feedback and lft). These functions are useful to model open- and closed-loop systems.
Model Interconnection Functions Concatenation of LTI Models LTI model concatenation is done in a manner similar to the way you concatenate matrices in MATLAB, using sys = [sys1 , sys2] % horizontal concatenation sys = [sys1 ; sys2] % vertical concatenation sys = append(sys1,sys2)% block diagonal appending In I/O terms, horizontal and vertical concatenation have the following block-diagram interpretations (with H 1 and H 2 denoting the transfer matrices of sys1 and sys2).
3 Operations on LTI Models Use append(sys1,sys2) to specify the block-decoupled LTI model interconnection. u1 sys1 y1 sys1 0 0 sys2 u2 sys2 y2 Appended Models Transfer Function See append on page 11-12 for more information on this command. Feedback and Other Interconnection Functions The following LTI model interconnection functions are useful for specifying closed- and open-loop model configurations: • feedback puts two LTI models with compatible dimensions in a feedback configuration.
Model Interconnection Functions For example, if sys1 has m inputs and p outputs, while sys2 has p inputs and m outputs, then the negative feedback configuration of these two LTI models + u sys1 y – sys2 is realized with feedback(sys1,sys2) This specifies the LTI model with m inputs and p outputs whose I/O map is ( I + sys1 ⋅ sys2 ) – 1 sys1 See Chapter 11, “Reference,” for more information on feedback, series, parallel, lft, and connect.
3 Operations on LTI Models Continuous/Discrete Conversions of LTI Models The function c2d discretizes continuous-time TF, SS, or ZPK models. Conversely, d2c converts discrete-time TF, SS, or ZPK models to continuous time. Several discretization/interpolation methods are supported, including zero-order hold (ZOH), first-order hold (FOH), Tustin approximation with or without frequency prewarping, and matched poles and zeros.
Continuous/Discrete Conversions of LTI Models u( t ) = u[ k ] , kT s ≤ t ≤ ( k + 1 )T s The signal u ( t ) is then fed to the continuous system H ( s ) , and the resulting output y ( t ) is sampled every T s seconds to produce y [ k ] . Conversely, given a discrete system H d ( z ) , the d2c conversion produces a continuous system H ( s ) whose ZOH discretization coincides with H d ( z ) .
3 Operations on LTI Models and you get back the original discrete-time system (up to canceling the pole/zero pair at z=–0.5): Zero/pole/gain: (z+0.5) --------(z+0.5)^2 Sampling time: 0.1 First-Order Hold First-order hold (FOH) differs from ZOH by the underlying hold mechanism. To turn the input samples u [ k ] into a continuous input u ( t ) , FOH uses linear interpolation between samples.
Continuous/Discrete Conversions of LTI Models H d ( z ) = H ( s' ) , where 2 z–1 s' = ------ -----------Ts z + 1 Similarly, the d2c conversion relies on the inverse correspondence H ( s ) = H d ( z' ), where 1 + sT s ⁄ 2 z' = -------------------------1 – sT s ⁄ 2 Tustin with Frequency Prewarping This variation of the Tustin approximation uses the correspondence ω z–1 s' = --------------------------------- -----------tan ( ωT s ⁄ 2 ) z + 1 H d ( z ) = H ( s' ) , This change of variable ensures the mat
3 Operations on LTI Models Discretization of Systems with Delays You can also use c2d to discretize SISO or MIMO continuous-time models with time delays. If Ts is the sampling period used for discretization: • A delay of tau seconds in the continuous-time model is mapped to a delay of k sampling periods in the discretized model, where k = fix(tau/Ts).
Continuous/Discrete Conversions of LTI Models The step responses of the continuous and discretized models are compared in the figure below. This plot was produced by the command step(h,'--',hd,'-') Note: The Tustin and matched pole/zero methods are accurate only for delays that are integer multiples of the sampling period. It is therefore preferable to use the zoh and foh discretization methods for models with delays.
3 Operations on LTI Models vectors when it is possible to reduce the total number of I/O delays. The resulting model has a minimum number of such delays.
Resampling of Discrete-Time Models Resampling of Discrete-Time Models You can resample a discrete-time TF, SS, or ZPK model sys1 by typing sys2 = d2d(sys1,Ts) The new sampling period Ts does not have to be an integer multiple of the original sampling period. For example, typing h1 = tf([1 0.4],[1 –0.7],0.1); h2 = d2d(h1,0.25); resamples h1 at the sampling period of 0.25 seconds, rather than 0.1 seconds.
3 Operations on LTI Models References [1] Åström, K.J. and B. Wittenmark, Computer-Controlled Systems: Theory and Design, Prentice-Hall, 1990, pp. 48–52. [2] Franklin, G.F., J.D. Powell, and M.L. Workman, Digital Control of Dynamic Systems, Second Edition, Addison-Wesley, 1990.
4 Arrays of LTI Models Introduction . . . . . . . . . . . . . . . When to Collect a Set of Models in an LTI Array . Restrictions for LTI Models Collected in an Array Where to Find Information on LTI Arrays . . . . . . . . . . . . . . . . . . . . . . . 4-2 4-2 4-2 4-3 The Concept of an LTI Array . . . . . . . . . . . . 4-4 Higher Dimensional Arrays of LTI Models . . . . . . . . 4-6 Dimensions, Size, and Shape of an LTI Array . . . . . 4-7 size and ndims . . . . . . . . . . . . . . . . . . . .
4 Arrays of LTI Models Introduction You can use LTI arrays to collect a set of LTI models into a single MATLAB variable. You then use this variable to manipulate or analyze the entire collection of models in a vectorized fashion. You access the individual models in the collection through indexing rather than by individual model names.
Introduction Where to Find Information on LTI Arrays The next two sections give examples that illustrate the concept of an LTI array, its dimensions, and size. To read about how to build an LTI array, go to “Building LTI Arrays” on page 4-12. The remainder of the chapter is devoted to indexing and operations on LTI Arrays. You can also apply the analysis functions in the Control System Toolbox to LTI arrays. See Chapter 5, “Model Analysis Tools,” for more information on these functions.
4 Arrays of LTI Models The Concept of an LTI Array To visualize the concept of an LTI array, consider the set of five transfer function models shown below. In this example, each model has two inputs and two outputs. They differ by parameter variations in the individual model components. 1.1 -----------s+1 0 1.3 ----------------s + 1.1 0 1.11 ----------------s + 1.2 0 1.15 ----------------s + 1.3 0 1.09 ----------------s + 1.4 0 0 1 -----------s+5 0 1 ----------------s + 5.
The Concept of an LTI Array Just as you might collect a set of two-by-two matrices in a multidimensional array, you can collect this set of five transfer function models as a list in an LTI array under one variable name, say, sys. Each element of the LTI array is an LTI model. Individual models in the LTI array sys are accessed via indexing. The general form for the syntax you use to access data in an LTI array is sys(Outputs,Inputs,Models) The first index selects the output channels.
4 Arrays of LTI Models Higher Dimensional Arrays of LTI Models You can also collect a set of models in a two-dimensional array. The following diagram illustrates a 2-by-3 array of six, two-output, one-input models called m2d. m2d(:,:,1,3) Each entry in this 2-by-3 array of models is a two-output, one-input transfer function. 3.42 -------------------s + 2.84 7.29 m2d(:,:,1,1) m2d(:,:,2,1) 3.36 ----------------s + 2.9 7.23 m2d(:,:,1,2) m2d(:,:,2,2 ) 3.4 -------------------s + 2.86 7.
Dimensions, Size, and Shape of an LTI Array Dimensions, Size, and Shape of an LTI Array The dimensions and size of a single LTI model are determined by the output and input channels. An array of LTI models has additional quantities that determine its dimensions, size, and shape.
4 Arrays of LTI Models . The next figure illustrates the concepts of dimension and size for the LTI array m2d, a 2-by-3 array of one-input, two-output transfer function models. t ar r ay di m ens ion is 2 The length of the second array dimension is 3. m2d(:,:,1,2) m2d(:,:,1,3) m2d(:,:,2,2 ) m2d(:,:,2,1) m2d(:,:,2,3) T he le n gt h of th e fi rs m2d(:,:,1,1) 3.4 -------------------s + 2.86 7.27 3.45 -------------------s + 2.81 7.32 m2d(:,:,2,3) Output dimension (length 2) 3.
Dimensions, Size, and Shape of an LTI Array Five related quantities are pertinent to understanding the array dimensions: • N, the number of models in the LTI array • K, the number of array dimensions • S 1 S 2 …S K, the list of lengths of the array dimensions - S i is the number of models along the i th dimension. • S 1 – by – S 2 – by – … – by – S K , the configuration of the models in the array - The configuration determines the shape of the array.
4 Arrays of LTI Models Note: • By convention, a single LTI model is treated as a 1-by-1 array of models. For single LTI models, size returns only the I/O dimensions [Ny Nu]. • For LTI arrays, size always returns at least two array dimensions. For example, the size of a 2-by-1 LTI array in [Ny Nu 2 1] • size ignores trailing singleton dimensions beyond the second array dimension. For example, size returns [Ny Nu 2 3] for a 2-by-3-by-1-by-1 LTI array of models with Ny outputs and Nu inputs.
Dimensions, Size, and Shape of an LTI Array Notice that size returns a vector whose entries correspond to the length of each of the four dimensions of m2d: two outputs and one input in a 2-by-3 array of models. Type ndims(m2d) ans = 4 to see that there are indeed four dimensions attributed to this LTI array. reshape Use reshape to reorganize the arrangement (array configuration) of the models of an existing LTI array.
4 Arrays of LTI Models Building LTI Arrays There are several ways to build LTI arrays: • Using a for loop to assign each model in the array • Using stack to concatenate LTI models into an LTI array • Using tf, zpk, ss, or frd In addition, you can use the command rss to generate LTI arrays of random state-space models. Generating LTI Arrays Using rss A convenient way to generate arrays of state-space models with the same number of states in each model is to use rss. The syntax is rss(N,P,M,sdim1,...
Building LTI Arrays Suppose, based on measured input and output data, you estimate confidence intervals [ω 1,ω 2] , and [ζ 1,ζ 2] for each of the parameters, ω and ζ. All of the possible combinations of the confidence limits for these model parameter values give rise to a set of four SISO models.
4 Arrays of LTI Models The first two colon indices ( : ) select all I/O channels from the I/O dimensions of H. The third index of H refers to the first array dimension ( ζ), while the fourth index is for the second array dimension (ω). Suppose the limits of the ranges of values for ζ and ω are [0.66,0.76] and [1.2,1.5], respectively. Enter these at the command line. zeta = [0.66,0.75]; w = [1.2,1.
Building LTI Arrays For the purposes of efficient computation, you can initialize an LTI array to zero, and then reassign the entire array to the values you want to specify. The general syntax for zero assignment of LTI arrays is Lengths of the output/input dimensions Lengths of the array dimensions sys = tf(zeros(Ny,Nu,S1,...,SK)) sys = zpk(zeros(Ny,Nu,S1,...,SK)) sys = ss(zeros(Ny,Nu,S1,...,SK,Nx)) sys = frd(zeros(Ny,Nu,Nf,S1,...
4 Arrays of LTI Models When you concatenate several models or LTI arrays along the jth array dimension, such as in stack(j,sys1,sys2,...,sysn) • The lengths of the I/O dimensions of sys1,...,sysn must all match. • The lengths of all but the jth array dimension of sys1,...,sysn must match. For example, if two TF models sys1 and sys2 have the same number of inputs and outputs, sys = stack(1,sys1,sys2) concatenates them into a 2-by-1 array of models.
Building LTI Arrays % Set up the LTI array using stack. COL1 = stack(1,H11,H21); % The first column of the 2-by-2 array COL2 = stack(1,H12,H22); % The second column of the 2-by-2 array H = stack(2, COL1, COL2); % Concatenate the two columns of models. Notice that this result is very different from the single MIMO LTI model returned by H = [H11,H12;H21,H22]; Building LTI Arrays Using tf, zpk, ss, and frd You can also build LTI arrays using the tf, zpk, ss, and frd constructors.
4 Arrays of LTI Models where • Both zeros and poles are multidimensional cell arrays whose cell entries contain the vectors of zeros and poles for each I/O pair of each model in the LTI array. • gains is a multidimensional array containing the scalar gains for each I/O pair of each model in the array. • The dimensions (and their lengths) of zeros, poles, and gains, determine those of the LTI array, sys.
Building LTI Arrays where • N s is the maximum of the number of states in each model in the array. • N u is the number of inputs in each model. • N y is the number of outputs in each model. • S 1, S 2, …, S K are the lengths of the array dimensions.
4 Arrays of LTI Models Indexing Into LTI Arrays You can index into LTI arrays in much the same way as you would for multidimensional arrays to: • Access models • Extract subsystems • Reassign parts of an LTI array • Delete parts of an LTI array When you index into an LTI array sys, the indices should be organized according to the following format sys(Outputs, Inputs, n 1, …, n K ) where • Outputs are indices that select output channels. • Inputs are indices that select input channels.
Indexing Into LTI Arrays For example, if sys is a 5-by-2 array of state-space models defined by sys = rss(4,3,2,5,2); you can access (and display) the model located in the (3,2) position of the array sys by typing sys(:,:,3,2) If sys is a 5-by-2 array of 3-output, 2-input FRD models, with frequency vector [1,2,3,4,5], then you can access the response data corresponding to the middle frequency (3 rad/s), of the model in the (3,1) position by typing sys(:,:,3,1,'frequency',3.
4 Arrays of LTI Models selects the first two input channels, and the first output channel in each model of the LTI array A, and returns the resulting 5-by-2 array of one-output, two-input subsystems. You can also combine model selection with I/O selection within an LTI array.
Indexing Into LTI Arrays Similarly, the commands sys(:,:,3,2) = sys(:,:,4,1); sys(1,2,3,2) = 0; reassign the entire model in the (3,2) position of the LTI array sys and the (1,2) subsystem of this model, respectively. LTI Arrays of SS Models with Differing Numbers of States You must use an entire LTI model for reassignment if you have an LTI array sys of state-space models for which: • The numbers of states in each model in sys is not constant.
4 Arrays of LTI Models Similarly, sys(:,:,[3 4],:) = [] deletes the third and fourth rows of this two-dimensional array of models.
Operations on LTI Arrays Operations on LTI Arrays Using LTI arrays, you can apply almost all of the basic model operations that work on single LTI models to entire sets of models at once. These basic operations, discussed in Chapter 3, “Operations on LTI Models,” include: • The arithmetic operations: +, –, *, /,\,',.
4 Arrays of LTI Models Example: Addition of Two LTI Arrays The following diagram illustrates the addition of two 3-by-1 LTI arrays sys1+sys2. sys2(:,:,3) sys1(:,:,3) 1 -----------s+2 + sys1(:,:,2) 1 ----------------s + 2.5 sys1 2s + 6.5 -------------------------------2 s + 6.5s + 9 = sys2(:,:,2) + sys1(:,:,1) 1 ----------------s + 2.9 1 ----------------s + 4.5 sys(:,:,3) 2.1 -----------s+4 sys(:,:,2) = 3.15 + 9.25 ----------------------------------2 s + 6.
Operations on LTI Arrays Note that: • Each model in sys1 and sys2 must have the same number of inputs and outputs. This is required for the addition of two LTI arrays. • The lengths of the array dimensions of sys1 and sys2 must match. Dimension Requirements In general, when you apply any of these basic operations to two or more LTI arrays: • The I/O dimensions of each of the LTI arrays must be compatible with the requirements of the operation. • The lengths of array dimensions must match.
4 Arrays of LTI Models where sys, the result of the operation, is an LTI array with the same array dimensions as sys1. You can use shortcuts for coding sys = op(sys1,sys2) in the following cases: • For operations that apply to LTI arrays, sys2 does not have to be an array. It can be a single LTI model (or a gain matrix) whose I/O dimensions satisfy the compatibility requirements for op (with those of each of the models in sys1).
Operations on LTI Arrays You can do this efficiently by first setting up an LTI array h containing the SISO models 1 ⁄ ( s + τ ) and then using concatenation to form the LTI array H of MIMO LTI models H τ ( s ), τ ∈ { 1.1, 1.2, 1.3 }. To do this, type tau = [1.1 1.2 1.3]; for i=1:3 % Form LTI array h of SISO models.
4 Arrays of LTI Models adds a single SISO transfer function model to each entry in each model of the LTI array of MIMO models [h,h]. Finally, G = rand(1,1,3,1); sys = G + [h,h] adds the array of scalars to each entry of each MIMO model in the LTI array [h,h] on a model-by-model basis. This last command is equivalent to the following for loop.
5 Model Analysis Tools General Model Characteristics . . . . . . . . . . . 5-2 Model Dynamics . . . . . . . . . . . . . . . . . . 5-4 State-Space Realizations . . . . . . . . . . . . . . 5-7 Time and Frequency Response . . . Time Responses . . . . . . . . . . . Frequency Response . . . . . . . . . Plotting and Comparing Multiple Systems Customizing the Plot Display . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5 Model Analysis Tools General Model Characteristics General model characteristics include the model type, I/O dimensions, and continuous or discrete nature. Related commands are listed in the table below. These commands operate on continuous- or discrete-time LTI models or arrays of LTI models of any type. General Model Characteristics Commands 5-2 class Display model type ('tf', 'zpk', 'ss', or 'frd'). hasdelay Test true if LTI model has any type of delay.
General Model Characteristics This example illustrates the use of some of these commands. See the related reference pages for more details. H = tf({1 [1 –1]},{[1 0.1] [1 2 10]}) Transfer function from input 1 to output: 1 ------s + 0.1 Transfer function from input 2 to output: s – 1 -------------s^2 + 2 s + 10 class(H) ans = tf size(H) Transfer function with 2 input(s) and 1 output(s).
5 Model Analysis Tools Model Dynamics The Control System Toolbox offers commands to determine the system poles, zeros, DC gain, norms, etc. You can apply these commands to single LTI models or LTI arrays. The following table gives an overview of these commands. Model Dynamics covar Covariance of response to white noise. damp Natural frequency and damping of system poles. dcgain Low-frequency (DC) gain. dsort Sort discrete-time poles by magnitude. esort Sort continuous-time poles by real part.
Model Dynamics Here is an example of model analysis using some of these commands. h = tf([4 8.4 30.8 60],[1 4.12 17.4 30.8 60]) Transfer function: 4 s^3 + 8.4 s^2 + 30.8 s + 60 --------------------------------------s^4 + 4.12 s^3 + 17.4 s^2 + 30.8 s + 60 pole(h) ans = –1.7971 –1.7971 –0.2629 –0.2629 + – + – 2.2137i 2.2137i 2.7039i 2.7039i zero(h) ans = –0.0500 + 2.7382i –0.0500 – 2.7382i –2.0000 dcgain(h) ans = 1 [ninf,fpeak] = norm(h,inf)% peak gain of freq. response ninf = 1.
5 Model Analysis Tools These functions also operate on LTI arrays and return arrays. For example, the poles of a three dimensional LTI array sysarray are obtained as follows. sysarray = tf(rss(2,1,1,3)) Model sysarray(:,:,1,1) ======================= Transfer function: -0.6201 s - 1.905 --------------------s^2 + 5.672 s + 7.405 Model sysarray(:,:,2,1) ======================= Transfer function: 0.4282 s^2 + 0.3706 s + 0.04264 ------------------------------s^2 + 1.056 s + 0.
State-Space Realizations State-Space Realizations The following functions are useful to analyze, perform state coordinate transformations on, and derive canonical state-space realizations for single state-space LTI models or LTI arrays of state-space models. State-Space Realizations canon Canonical state-space realizations. ctrb Controllability matrix. ctrbf Controllability staircase form. gram Controllability and observability gramians. obsv Observability matrix.
5 Model Analysis Tools state-space computations. Consequently, it is wise to use them only for analysis purposes and not in control design algorithms.
Time and Frequency Response Time and Frequency Response The Control System Toolbox contains a set of commands that provide the basic time and frequency domain analysis tools required for control system engineering. These commands apply to any kind of LTI model (TF, ZPK, or SS, continuous or discrete, SISO or MIMO). You can only apply the frequency domain analysis tools FRDs.
5 Model Analysis Tools where sys is any continuous or discrete LTI model or LTI array. For MIMO models, these commands produce an array of plots with one plot per I/O channel (or per output for initial and lsim). For example, h = [tf(10,[1 2 10]) , tf(1,[1 1])] step(h) produces the following plot. The simulation horizon is automatically determined based on the model dynamics.
Time and Frequency Response Note: When specifying a time vector t = [0:dt:tf], remember the following constraints on the spacing dt between time samples: • For discrete systems, dt should match the system sample time. • Continuous systems are first discretized using zero-order hold and dt as sampling period, and step simulates the resulting discrete system. As a result, you should pick dt small enough to capture the main features of the continuous transient response.
5 Model Analysis Tools In addition, the function margin determines the gain and phase margins for a given SISO open-loop model. These functions can be applied to single LTI models or LTI arrays. Table 5-1: Frequency Response Function Name Description bode Computes the Bode plot. evalfr Computes the frequency response at a single complex frequency (not for FRD models). freqresp Computes the frequency response for a set of frequencies. margin Computes gain and phase margins.
Time and Frequency Response For example, bode(sys,{0.1 , 100}) draws the Bode plot between 0.1 and 100 radians/second. You can also specify a particular vector of frequency points as in w = logspace(–1,2,100) bode(sys,w) The logspace command generates a vector w of logarithmically spaced –1 2 frequencies starting at 10 = 0.1 rad/s and ending at 10 = 100 rad/s. See the reference page for linspace for linearly spaced frequency vectors.
5 Model Analysis Tools several LTI models on a single plot. To do so, invoke the corresponding command line function using the list sys1,..., sysN of models as the inputs. step(sys1,sys2,...,sysN) impulse(sys1,sys2,...,sysN) ... bode(sys1,sys2,...,sysN) nichols(sys1,sys2,...,sysN) ... All models in the argument lists of any of the response plotting functions (except for sigma) must have the same number of inputs and outputs.
Time and Frequency Response The following example compares a continuous model with its zero-order-hold discretization. sysc = tf(1000,[1 10 1000]) sysd = c2d(sysc,0.2) % ZOH sampled at 0.
5 Model Analysis Tools bode(sysc,'--',sysd,'-') % compare Bode responses A comparison of the continuous and discretized responses reveals a drastic undersampling of the continuous system. Specifically, there are hidden oscillations in the discretized time response and aliasing conceals the continuous-time resonance near 300 rad/sec.
Time and Frequency Response Customizing the Plot Display You can plot data generated by several response analysis functions applied to one or several LTI models, as well as your own data.
5 Model Analysis Tools For example, the following sequence of commands displays the Bode plot, step response, pole/zero map, and some additional data in a single figure window. h = tf([4 8.4 30.8 60],[1 4.12 17.4 30.
Time and Frequency Response Another example is subplot(221) bode(h) subplot(222) step(h) subplot(223) pzmap(h) subplot(224) plot(rand(1, 100)) % any data can go here title('Some noise') Note: Each of the plots generated by response analysis functions in these figures (here, bode, step, and pzmap) has its own right-click menu (similar to those in the LTI viewer. For more information, see “The Right-Click Menus” on page 6-18.
5 Model Analysis Tools Model Order Reduction You can derive reduced-order models with the following commands. Model Order Reduction balreal Input/output balancing. minreal Minimal realization or pole/zero cancellation. modred State deletion in I/O balanced realization. sminreal Structurally minimal realization Use minreal to delete uncontrollable or unobservable state dynamics in state-space models, or cancel pole/zero pairs in transfer functions or zero-pole-gain models.
6 The LTI Viewer Introduction . . . . . . . . . . . . . . . . . . . . 6-2 Getting Started Using the LTI Viewer: An Example . . 6-4 The LTI Viewer Menus . . . . . . . . . . . . . . . 6-15 The Right-Click Menus . . . . . . . . . . . . . . . 6-18 The LTI Viewer Tools Menu . . . . . . . . . . . . . 6-39 Simulink LTI Viewer . . . . . . . . . . . . . . . .
6 The LTI Viewer Introduction The LTI Viewer is a graphical user interface for viewing and manipulating the response plots of LTI models.
Introduction • Zoom in on or out from the individual displayed plots • Toggle the grid on or off on a plot • Select which I/O channels the LTI Viewer displays for MIMO models in each plot • For a given plot type, select how the LTI Viewer displays the I/O channels for MIMO models • Select which models of an LTI array you want displayed in the LTI Viewer by indexing into dimensions or model characteristics • Control plot characteristics such as the ranges for time and frequency used in various types of plot
6 The LTI Viewer Getting Started Using the LTI Viewer: An Example This section contains a brief introduction to the LTI Viewer through an example that leads you through the following steps: 1 Load two LTI models into the LTI Viewer, initialized with the step responses and Bode plots of both models.
Getting Started Using the LTI Viewer: An Example Transfer function: 2 s^3 + 1.2 s^2 + 15.1 s + 7.5 ---------------------------------------s^4 + 2.12 s^3 + 10.2 s^2 + 15.1 s + 7.5 Transfer function: 1.2 s^3 + 1.12 s^2 + 9.1 s + 7.5 ---------------------------------------s^4 + 1.32 s^3 + 10.12 s^2 + 9.1 s + 7.5 Initializing the LTI Viewer with Multiple Plots For a given LTI model, you can use the LTI Viewer to simultaneously display multiple response plot types, such as the Bode plot and the step response.
6 The LTI Viewer Plot Type Description pzmap Plot of poles and zeros sigma Singular values of the frequency response step Step response Note: When you initialize the LTI Viewer with lsim or initial, these plot types require some extra arguments. For more information on the syntax for calling ltiview, see ltiview on page 11-133.
Getting Started Using the LTI Viewer: An Example Plot region for the step response The File menu has several items, including ones that allow you to import models and print plots. The Tools menu items allow you to reconfigure the plot arrangement and set general plot and linestyle preferences for this open LTI Viewer. Right-click on either plot region to see the response plot menus. Two response curves are plotted on each of the two plot regions.
6 The LTI Viewer 1 Right-click anywhere in the plot region of the step response plots. This opens the following menu list in the plot region. Figure 6-1: The Right-Click Menu for SISO Models 2 Place your mouse pointer on the Characteristics menu item, and select Settling Time with your left mouse button. Figure 6-2: The Step Response Characteristics Submenu 3 Right-click anywhere in the plot region of the Bode plots to open a right-click menu.
Getting Started Using the LTI Viewer: An Example Your LTI Viewer should now look like this. Notice that there is one settling time or peak magnitude marker for each LTI model displayed in the LTI Viewer. Displaying Response Characteristics on a Plot To display the values of any plot characteristic marked on a plot: 1 Click on the marker 2 Hold the left or right mouse button down to read the values off the plot.
6 The LTI Viewer Hold the mouse button down on the marker to display the values. Note that you can: • Use either the right or the left mouse button when you select a marker on a plot. • Left-click anywhere on a particular plot line to see the response values of that plot at that point. • Right-click anywhere on a plot line to see I/O and model information.
Getting Started Using the LTI Viewer: An Example Importing Models If the closed-loop models Gcl1 and Gcl2 do not meet your specifications, you may want to design another compensator at the command line, and import the resulting closed-loop model Gcl3 for comparison: 1 Select Import from the File menu. This opens a browser listing all of the LTI models currently available in the MATLAB workspace.
6 The LTI Viewer Note: A given LTI Viewer can only be used to analyze models with the same number of inputs and outputs. If you want to analyze models with different numbers of inputs or outputs, you must import these into separate LTI Viewers. See “Opening a New LTI Viewer” on page 6-16 for more information. Zooming With three models loaded into the LTI Viewer, you may want to zoom in on one region of a given plot, in order to inspect the response behavior in that region more closely.
Getting Started Using the LTI Viewer: An Example a Point your mouse to any corner of the rectangle of the region you want to zoom in on. b Left-Click there, and hold the mouse button down. c Drag the mouse pointer until the rectangle covers the region you want to zoom in on. d Release your mouse. For this example, zoom in around the region near 4.5 seconds on the step response plot. Your step response plot looks like this as you select the zoom region.
6 The LTI Viewer After releasing the mouse on the zoom region, the LTI Viewer looks like this. Notice that you’ve only zoomed on the step response plot; the Bode plot remains unchanged. Note: To zoom out, i.e., to revert back to the original coordinate limits that were in place before you zoomed, follow the steps for zooming again, only this time select Out from the Zoom menu.
The LTI Viewer Menus The LTI Viewer Menus The LTI Viewer has three main menus: • File • Tools • Help The File menu provides features pertinent to bringing data in and out of the LTI Viewer. The Help provides help on the LTI Viewer features. The File and Help menus are covered in this section. The Tools menu allows you to control certain features common to all of the plots. You can read about the Tools menu items in “The LTI Viewer Tools Menu” on page 6-39.
6 The LTI Viewer For directions for loading LTI models into the LTI Viewer workspace when you open it, see “Initializing the LTI Viewer with Multiple Plots” on page 6-5. For directions for importing models into the workspace of an open LTI Viewer, see “Importing Models” on page 6-11. Opening a New LTI Viewer The New Viewer option in the File menu enables you to initialize a new LTI Viewer. This is the same as typing ltiview at the MATLAB prompt.
The LTI Viewer Menus The first submenu, Overview, opens the help text describing how to use the LTI Viewer menus and right-click menus that control the LTI Viewer. The remaining help menu items pertain to LTI Viewer controls you access from the Tools menu: the Response Preferences and Linestyle Preferences windows. These windows provide additional tools for manipulating the system responses.
6 The LTI Viewer The Right-Click Menus You can access most of the controls for the individual response plots displayed by the LTI Viewer through the right-click menus located in any plot region. There is one right-click menu per plot region displayed on the LTI Viewer. These menus vary slightly, depending on the model dimensions and plot type: • The menu items that appear on the plot regions of the responses for SISO models are the basic right-click menu items.
The Right-Click Menus This is the right-click menu for SISO models. These menu items control the LTI Viewer plots for all models: • Plot Type—You can choose which plot type you want displayed from this list of nine plot types. The check mark next to the plot type indicates your selection for the type of plot displayed. You can select from any of the items on this list. Figure 6-5: Plot Type Submenu • Systems—The Systems submenu lists the models in the LTI Viewer workspace.
6 The LTI Viewer are in the LTI Viewer workspace, but their responses are not displayed on the plot associated with the open menu. - You can select any model in the list with your mouse to toggle on (or off) the visibility of its response curve in the selected plot region. • Characteristics—You can toggle on and off the option to display a marker for various response characteristics for each plot type. For more information, see “Displaying Response Characteristics on a Plot” on page 6-9.
The Right-Click Menus Note: To multiselect submenu items (such as in the Characteristics or the Systems menus), re-open the right-click menu for each submenu item selection. The Right-Click Menu for MIMO Models When you load a MIMO model into the LTI Viewer, the right-click menu has a few more options than the right-click menu for SISO models does.
6 The LTI Viewer Your LTI Viewer looks like this. This plot region is displayed as a grid of four distinct plots: one for each I/O response of ssF8. Figure 6-8: Nyquist Plots of the Four I/O Responses in ssF8 Notice that the I/O names for this model appear on the Nyquist plot. Each of the four plots displayed represents the I/O response from a single input to a single output.
The Right-Click Menus Right-click on any part of the plot region (anywhere on the grid of plots). This opens the following menu. The right-click menu for MIMO models has two extra menu items: Axes Grouping and Select I/Os. Figure 6-9: The Right-Click Menu for MIMO Models The Axes Grouping Submenu The Axes Grouping submenu is as follows. Default Axes Grouping setting: Each I/O response is displayed individually.
6 The LTI Viewer You can use the Axes Grouping submenu to reconfigure the grouping of these I/O response curves with the following submenu items. • Inputs: The response curves from all of the inputs to a given output are plotted in the same portion of the plot region. There are as many separate portions of the plot region displayed as there are outputs. The plot region is divided into two portions: one for each output.
The Right-Click Menus • Outputs: The response curves from a given input to all of the outputs are plotted in the same portion of the plot region. There are as many separate portions of the plot region displayed as there are inputs. The plot region is divided into two portions: one for each input.
6 The LTI Viewer • All: All of the I/O response curves are displayed (grouped) in a single plot region. The responses appear all on the same plot region. Figure 6-12: Axes Grouping: All The Select I/Os Menu Item The LTI Viewer initially displays all of the I/O response curves from each input channel to each output channel. You can select the Select I/Os menu item to customize the display with respect to the input and output channels.
The Right-Click Menus When you select Select I/Os from the right-click menu, the following window opens. Select one of these input channel names to display only I/O responses from the selected input. Select all to display the responses from all I/O channels. Select one of these output channel names to display only I/O responses to the selected output. Hold down the Shift key while selecting individual I/O response markers to customize the display.
6 The LTI Viewer With the Axes Grouping set to None, the display looks like this. Note: To reset the Axes Grouping to None, open the right-click menu on the plot region, and select None. The Right-Click Menu for LTI Arrays When you load an LTI array into the LTI Viewer, all models in the LTI array are initially displayed.
The Right-Click Menus For a given LTI array in the LTI Viewer workspace, you can use this interface to display the plots of a subset of models in the LTI array, using either or both of the following options: • Indexing into the array dimensions • Indexing into the array through design specification criteria In order to have access to right-click menu item for LTI arrays, you must have at least one LTI array loaded in the LTI Viewer workspace.
6 The LTI Viewer To display the responses of only some of the models in the LTI array, you must first complete the following two steps: 1 Right-click anywhere in the plot region to open the following right-click menu. Figure 6-13: Right-Click Menu for LTI Arrays 2 Select the Select from LTI Array menu item.
The Right-Click Menus This opens the Model Selector for LTI Arrays window in the (default) Index into Dimensions setup. This tab lists all LTI arrays in the LTI Viewer. You only apply selection criteria to one LTI array at a time. Each listbox corresponds to a dimension of the LTI array. The number of entries in each listbox is the same as the number of models along the corresponding dimension. Each numbered entry in a given listbox represents the indices for each dimension of the selected LTI array.
6 The LTI Viewer Once you have selected the name of an LTI array in the Model Selector for LTI Arrays window, you can select models in the LTI array whose response plots you want displayed using either or both of the following: • Index into the array dimensions of the LTI array (Index into Dimensions) • Index into LTI array using design specification criteria (Bound on Characteristics) Indexing into the Array Dimensions of an LTI Array To index into the array dimensions of an LTI Array: 1 Select Index in
The Right-Click Menus For example to display only the first row of models in the 2-by-3 LTI array m2d, either: • Select the first index in the first listbox (corresponding to the first dimension of the LTI array) with your mouse. • Type the vector [1] in the textbox below the first listbox. The following figure depict the Model Selector for LTI Arrays window for selecting to display the responses of the first row of models in the LTI array m2d.
6 The LTI Viewer The plots of only three models (as opposed to the six models in the LTI array) are shown here. Figure 6-16: Step Response of the First Row of Models in m2d There are a variety of ways you can index into the dimensions of an LTI array using the textboxes located below each listbox. You can type both logical expressions, or ones that define indices directly.
The Right-Click Menus Indexing into the LTI Array Using Design Specification Criteria You can also use several plot-specific design criteria to select those models in the LTI array whose responses you want displayed. You index into the LTI array through these design criteria (response plot characteristics) using Boolean expressions. To plot selected models by indexing into the LTI array in this manner: 1 Select show all in the Index into Dimensions set-up, and select Apply.
6 The LTI Viewer 4 Position your mouse pointer in the textbox next to the design specification characteristic. 5 Type a MATLAB relational expression in the textbox, using $ as a variable name in the expression. Note that for arrays of MIMO models, you can consider $ to be an Ny-by-Nu matrix, if each model in the LTI array has Ny outputs and Nu inputs. This allows you to specify different requirements on different I/O channels. 6 Press Apply to implement your indexing selection.
The Right-Click Menus The step response of the model in the (2,3) position of the LTI array is displayed. Figure 6-19: Step Response of the Model with the Maximum Rise Time You can also use any logical expression in variables defined in the MATLAB workspace to index into a specific design criterion. For example, typing $(2,1) < 7.25 & $(1,1) > 1.2 next to Steady State (after unchecking Rise Time), displays the responses of any models for which the steady-state response has a value less than 7.
6 The LTI Viewer at the MATLAB command line, and $ > any(any(Q)) in the Model Selector for LTI Arrays window. This displays only the plots of those models for which the required bound is not satisfied on any of the I/O channels. You must use the any command twice, once for each I/O dimension. Typing $>Q would only display the plots of those models for which is bound is not satisfied on all of the I/O channels.
The LTI Viewer Tools Menu The LTI Viewer Tools Menu Three preferences windows provide additional options for customizing the LTI Viewer display. You can access these windows from the Tools menu.
6 The LTI Viewer load this model). With the Available LTI Viewer Configurations window open: 1 Select the radio button for the two-plot configuration. 2 Use the pull-down tab next to 1. to set the first plot to nyquist. 3 Use the pull-down tab next to 2. to set the second plot to bode. 4 Select OK.
The LTI Viewer Tools Menu Select Define to define your own time or frequency vectors in the editable text boxes. Select Define to define your own vertical axis range. Edit these text boxes to customize the response characteristics for step response plots. Use these radio buttons to change units on Bode and sigma plots. Always press OK or Apply to execute window options.
6 The LTI Viewer value. You can use the Time vector (sec.) portion of the Time Domain field shown below to do this. The Time vector (sec.) field accepts one, two, or three arguments, separated by colons and surrounded by square brackets: • [Tf] specifies only the final time. • [Ti:Tf] specifies the initial and final time. • [Ti:dt:Tf] specifies the initial and final time and provides the incremental step dt to use when generating the time vector.
The LTI Viewer Tools Menu The Frequency vector (rad/sec.) field also provides you with the option to recalculate a new frequency vector for each frequency response type. When this checkbox is selected along with Generate automatically, a new frequency vector and response is calculated each time you toggle between different frequency responses, e.g.from Bode to Nyquist.
6 The LTI Viewer rise time. For example, you can change the value for the settling time to 5% as shown below. Figure 6-22: Changing the Settling Time Percentage Value to 5% Changing the Frequency Domain Plot Units In addition to providing options for specifying the frequencies used in the frequency responses, the Frequency Domain field allows you to choose the units used for Bode plots and singular value (sigma) plots.
The LTI Viewer Tools Menu After selecting Linestyle Preferences from the Tools menu, the following Linestyle Preferences window opens. Distinguish multiple response curves according to any of these four characteristics. If you select a line property, you can use the arrow buttons next to one of the three selected line property listboxes to rearrange the order of the entries in that listbox. Distinguish multiple response curves by either of these three line properties, or not at all.
6 The LTI Viewer You can designate that the chosen preference (color, marker, or linestyle) distinguish the plotted response curves by any (or all) of the following. • Systems: The line properties vary with the models. • Inputs: The line properties vary with the input channels. • Outputs: The line properties vary with the output channels. • Channels: The line properties vary with the I/O channels.
The LTI Viewer Tools Menu 2 Select either: - Apply to keep the Linestyle Preferences window open when you apply these changes - OK to apply the changes and close the Linestyle Preferences window The Order in which Line Properties are Assigned You can determine the order in which the line properties are applied to models (or inputs, outputs, or I/O channels) by referring to the order of the line properties in the listboxes.
6 The LTI Viewer Simulink LTI Viewer If you have Simulink, you can use the Simulink LTI Viewer, a version of the LTI Viewer that performs linear analysis on any portion of a Simulink model. The Simulink LTI Viewer features: • Drag-and-drop blocks that identify the location for the inputs and outputs of the portion of a Simulink model you want to analyze. • The ability to specify the operating conditions about which the Simulink model is linearized for analysis in the LTI Viewer.
Simulink LTI Viewer Notice that the title of this Simulink model is vdp, and that it contains static nonlinearities. A Sample Analysis Task Suppose you want to: • Analyze the Bode plot of the linear response between the input x2 to the Product block and the output x1 of the second Integrator block. • Determine the effect of changing the value of the Gain block labeled Mu on this response.
6 The LTI Viewer The basic procedure for carrying out this type of analysis is outlined below: 1 Open the Simulink LTI Viewer. 2 Specify your analysis model: a Specify the portion of the Simulink model you want to analyze. This involves using special Simulink blocks to locate the inputs and outputs of this analysis model on your Simulink diagram. b Set the operating conditions for linear analysis (optional).
Simulink LTI Viewer When you select Linear Analysis, two new windows open: an LTI Viewer window and a Simulink diagram called Model_Inputs_and_Outputs containing two blocks: Input Point and Output Point. The following figure depicts how to open the Simulink LTI Viewer.
6 The LTI Viewer Select Linear Analysis from the Simulink model Tools menu to open the Simulink LTI Viewer. This set of Simulink blocks opens when you select Linear Analysis. Use these blocks to specify the inputs and outputs of the portion of the Simulink model you want to analyze. Your Tools menu may differ from this one, depending on your options.
Simulink LTI Viewer • The title bar shows the name of the Simulink model to which it is linked. • The menu bar contains an additional menu called Simulink that contains the following items: - Get Linearized Model linearizes the Simulink model and imports the resulting linearized analysis model to the LTI Viewer. Each time you select this menu item, a new version of the linearized analysis model is added to the Simulink LTI Viewer workspace.
6 The LTI Viewer To set up the analysis model for the vdp Simulink model: 1 Insert an Input Point block on the line labeled x2 going into the Product block. 2 Insert an Output Point block on the line entering the Outport block, Out1. This results in the following diagram.
Simulink LTI Viewer Grab an Input Point block, drag it, and release it on the Simulink model here. .
6 The LTI Viewer Keep the following in mind when using the Input Point and Output Point blocks to specify analysis models: • You must place at least one Input Point block and one Output Point block somewhere in the diagram in order to specify an analysis model. • You can place the Input Point and Output Point blocks on any scalar or vector signal line in the Simulink model, with the exception of signal lines connected to any block in the Power System Blockset.
Simulink LTI Viewer Consider the following simple diagram. Delete this line to isolate the plant, P, for open loop analysis. Based on the location of the Input Point and Output Point blocks, you might think that the analysis model specified by these blocks is simply the plant model, P. However, due to the feedback loop, this analysis model is actually the closed-loop transfer function P ⁄ ( 1 + PK ) .
6 The LTI Viewer 1 Select Set Operating Point in the Simulink menu. This opens the Operating Point window. When you select this radio button, the Simulink LTI Viewer uses zero state values for each linearization. The default setting for the Operating Point window uses initial states from the Simulink diagram. Select either of the other two radio buttons to set all state values for linearization to zero, or choose some other values.
Simulink LTI Viewer 4 Select OK. A dialog box opens. a Selecting Yes closes the Operating Point window and linearizes the model automatically. The new operating conditions you selected remain in effect for any future linearizations. b Selecting No closes the Operating Point window without linearizing the model. The new operating conditions you selected remain in effect for any future linearizations.
6 The LTI Viewer For this example, we use the Zero initial states setting, shown in the figure below. Choose this menu on the Simulink LTI Viewer to set the operating conditions. Edit the operating point values here. Simulink path name for the state variables Selecting OK opens a dialog box, asking you to accept the changes you made to the operating conditions and linearize the model with these operating conditions.
Simulink LTI Viewer • While the Operating Point window is in the Initial state in Simulink diagram mode, the linearization values used by the Simulink LTI Viewer are updated as you change any state values in your Simulink diagram. • To use the MATLAB command line to change Simulink diagram initial state values - Select Parameters under the Simulation menu on your Simulink diagram. - Choose the Workspace I/O tab in the Simulation Parameters window.
6 The LTI Viewer Importing a Linearized Analysis Model to the LTI Viewer To linearize and load your first analysis model for the van der Pol example into the LTI Viewer, select Get Linearized Model on the Simulink menu on the LTI Viewer. This produces the following response plot: The LTI Viewer displays the selected response plot for all of the models in the LTI Viewer workspace. These are listed in the Systems menu available by right-clicking in the plot region. The default plot type is step response.
Simulink LTI Viewer You can view the linearized models in the LTI Viewer workspace or change the plot type using the right-click menus. See “The Right-Click Menus” on page 6-18 for more information. Analyzing the Bode Plot of the Linearized Analysis Model To analyze the Bode plot for this model, 1 Right-click anywhere in the plot region. 2 Select the Plot Type menu, and then the Bode submenu item.
6 The LTI Viewer As is shown below, the linearized model for the new value of Mu appears as the last item in the Systems submenu, and the Bode plots for both models are displayed. After reselecting Get Linearized Model, the Systems submenu contains two model names with different version numbers. The Bode plots for both models are displayed in different colors. As you might expect, changing the gain alters the Bode plot.
Simulink LTI Viewer Notice the following about the models listed in the Systems list on the right-click menu: • The names reflect the title of the Simulink model. • The version number is incremented every time Get Linearized Model is selected. • The LTI Viewer simultaneously displays the response plots of all of the checked models listed. Saving Analysis Models The analysis models obtained each time you select Get Linearized Model are stored only in the Simulink LTI Viewer workspace.
6 The LTI Viewer To export analysis models from the LTI Viewer workspace: 1 Highlight the models you want to save in the Export List. (You can multiselect models on the list by holding down the control key while selecting list items). 2 Save the selected models by clicking on the appropriate button on this GUI. You can export the linearized models to either: a A MAT-file b The MATLAB workspace Note: If you save models to a MAT-file, you are prompted to name the file.
7 Control Design Tools Root Locus Design . . . . . . . . . . . . . . . . . 7-3 Pole Placement . . . . . . State-Feedback Gain Selection State Estimator Design . . . Pole Placement Tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7-5 7-5 7-5 7-6 LQG Design . . . . . . . Optimal State-Feedback Gain Kalman State Estimator . . . LQG Regulator . . . . . . LQG Design Tools . . . . . . . . . . . . . . . . . . . . . . . . . . .
7 Control Design Tools We use the term control system design to refer to the process of selecting feedback gains in a closed-loop control system. Most design methods are iterative, combining parameter selection with analysis, simulation, and physical insight. The Control System Toolbox offers functions to: • Design the control system gains using either classical root locus techniques or modern pole placement and LQG techniques. • Close the loop for simulation and validation purposes.
Root Locus Design Root Locus Design The root locus method is used to describe the trajectories of the closed-loop poles of a feedback system as one parameter varies over a continuous range of values. Typically, the root locus method is used to tune a feedback gain so as to specify the closed-loop poles of a SISO control system. Consider, for example, the tracking loop r + P( s ) y – k H(s) where P ( s ) is the plant, H ( s ) is the sensor dynamics, and k is a scalar gain to be adjusted.
7 Control Design Tools If you are interested in just the root locus plot, use the command rlocus. This command takes one argument: a SISO model of the open loop system, created with ss, tf, or zpk. In the tracking loop depicted on the previous page, this model would represent P ( s )H ( s ) . You can also use the function rlocfind to select a point on the root locus plot and determine the corresponding gain k . The following table summarizes the commands for root locus design.
Pole Placement Pole Placement The closed-loop pole locations have a direct impact on time response characteristics such as rise time, settling time, and transient oscillations. This suggests the following method for tuning the closed-loop behavior: 1 Based on the time response specifications, select desirable locations for the closed-loop poles. 2 Compute feedback gains that achieve these locations. This design technique is known as pole placement.
7 Control Design Tools · ξ = Aξ + Bu + L ( y – Cξ – Du ) The estimator poles are the eigenvalues of A – LC , which can be arbitrarily assigned by proper selection of the estimator gain matrix L . As a rule of thumb, the estimator dynamics should be faster than the controller dynamics (eigenvalues of A – BK ).
Pole Placement Pole Placement acker SISO pole placement. estim Form state estimator given estimator gain. place MIMO pole placement. reg Form output-feedback compensator given state-feedback and estimator gains. The function acker is limited to SISO systems and should only be used for systems with a small number of states. The function place is a more general and numerically robust alternative to acker. Caution: Pole placement can be badly conditioned if you choose unrealistic pole locations.
7 Control Design Tools LQG Design Linear-Quadratic-Gaussian (LQG) control is a modern state-space technique for designing optimal dynamic regulators. It enables you to trade off regulation performance and control effort, and to take into account process and measurement noise. Like pole placement, LQG design requires a state-space model of the plant (use ss to convert other LTI models to state space).
LQG Design Optimal State-Feedback Gain In LQG control, the regulation performance is measured by a quadratic performance criterion of the form J( u ) = ∞ ∫0 { x T T T Qx + 2x Nu + u Ru } dt The weighting matrices Q, N, R are user specified and define the trade-off between regulation performance (how fast x ( t ) goes to zero) and control effort. The first design step seeks a state-feedback law u = – Kx that minimizes the cost function J ( u ) .
7 Control Design Tools LQG Regulator To form the LQG regulator, simply connect the Kalman filter and LQ-optimal gain K as shown below: w Plant y u u –K x̂ + Kalman filter v yv + LQG regulator This regulator has state-space equations · x̂ = A – LC – ( B – LD )K x̂ + Ly v u = – Kx̂ LQG Design Tools The Control System Toolbox contains functions to perform the three LQG design steps outlined above.
LQG Design discrete problems as well as the design of discrete LQG regulators for continuous plants. LQG Design care Solve continuous-time algebraic Riccati equations. dare Solve discrete-time algebraic Riccati equations. dlqr LQ-optimal gain for discrete systems. kalman Kalman estimator. kalmd Discrete Kalman estimator for continuous plant. lqgreg Form LQG regulator given LQ gain and Kalman filter. lqr LQ-optimal gain for continuous systems. lqrd Discrete LQ gain for continuous plant.
7 Control Design Tools 7-12
8 The Root Locus Design GUI Introduction . . . . . . . . . . . . . . . . . . . . 8-2 A Servomechanism Example . . . . . . . . . . . . 8-4 Controller Design Using the Root Locus Design GUI Opening the Root Locus Design GUI . . . . . . . . . Importing Models into the Root Locus Design GUI . . . Changing the Gain Set Point and Zooming . . . . . . Displaying System Responses . . . . . . . . . . . Designing a Compensator to Meet Specifications . . . . Saving the Compensator and Models . . . . . . . . . . . .
8 The Root Locus Design GUI Introduction The root locus method is used to describe the trajectories in the complex plane of the closed-loop poles of a SISO feedback system as one parameter (usually a gain) varies over a continuous range of values. Along with the MATLAB command line function rlocus, the Control System Toolbox provides the Root Locus Design graphical user interface (GUI) for implementing root locus methods on single-input/single-output (SISO) LTI models defined using zpk, tf, or ss.
Introduction - Editing the compensator pole and zero locations - Opening the LTI Viewer • Saving the compensator and models to the workspace or the disk Other important features listed below and not covered through this example are described in the section, “Additional Root Locus Design GUI Features” on page 8-38: • Specifying design models: general concepts • Getting help with the Root Locus Design GUI • Erasing compensator poles and zeros • Listing poles and zeros • Printing the root locus • Drawing a Si
8 The Root Locus Design GUI A Servomechanism Example A simple version of an electrohydraulic servomechanism model consists of: • A push-pull amplifier (a pair of electromagnets) • A sliding spool in a vessel of high pressure hydraulic fluid • Valve openings in the vessel to allow for fluid to flow • A central chamber with a piston-driven ram to deliver force to a load • A symmetrical fluid return vessel The force on the spool is proportional to the current in the electromagnet coil.
A Servomechanism Example A closed-loop model for the electrohydraulic valve position control can be set up as follows. Desired Ram Position + – K(s) G servo ( s ) Compensator Linearized Servomechanism Ram Position Figure 8-1: Feedback Control for an Electrohydraulic Servomechanism K(s) represents the compensator for you to design. This compensator can be either a gain or a more general LTI system.
8 The Root Locus Design GUI Controller Design Using the Root Locus Design GUI In this section, we use the servomechanism example to describe some of the main Root Locus Design GUI features and operations: • Opening the Root Locus Design GUI • Importing models into the Root Locus Design GUI: - Opening the Import LTI Design Model window - Choosing a feedback structure - Specifying the design model • Changing the gain set point and zooming: - Dragging closed-loop poles to change the gain set point - Zooming
Controller Design Using the Root Locus Design GUI This brings up the following GUI. These are the main menus for importing/ exporting of models, and editing them. You can also perform discrete/continuous conversion. Compensator description: The default compensator is K=1. Root locus toolbar buttons to drag, add, or erase compensator poles and zeros, or to drag the closed-loop poles Plot region to display the root locus Toggle between positive/negative feedback.
8 The Root Locus Design GUI There are four ways to import SISO LTI models into the Root Locus Design GUI: • Load a model from the MATLAB workspace. • Load a model from a MAT-file on your disk. • Load SISO LTI blocks from an open or saved Simulink diagram. • Create models using tf, ss, or zpk within the GUI. You can also use any combination of these methods to import models for root locus analysis and design.
Controller Design Using the Root Locus Design GUI Opening the Import LTI Design Model Window To import the linearized electrohydraulic servomechanism model into the Root Locus Design GUI, first open the Import LTI Design Model window. To do this, select the Import Model menu item in the File menu. The Import LTI Design Model window that appears on your screen is as follows.
8 The Root Locus Design GUI Diagram of feedback structure Model listbox: lists available LTI models or blocks. Edit the name of the overall design model. This name appears on the title bar of the Root Locus Design GUI. F,P, and H compose the design model. Toggle button for the feedback structure SISO LTI models can be imported from either the workspace or disk. SISO LTI blocks can be imported from saved or open Simulink models.
Controller Design Using the Root Locus Design GUI The Feedback Structure portion of the Import LTI Design Model window shows the current selection for the closed-loop structure. The Other button toggles the location of the compensator between the two configurations shown above. For this example you want the compensator in the forward path. Specifying the Design Model The SISO LTI models in either feedback configuration are coded as follows: • F represents a pre-filter. • P is the plant model.
8 The Root Locus Design GUI To specify the design model components for this servomechanism example: 1 Load the linearized servomechanism model Gservo into the plant P, by first selecting it and then selecting the arrow button. The Import LTI Design Model window looks like this after you specify both the plant, P, and the model name as Gservo. Use the arrow buttons to transfer selected models from the listbox to the design model components, in this case, to the plant, P. 2 Press the OK button.
Controller Design Using the Root Locus Design GUI The root locus of the design model is displayed in the plot region of the GUI. Your Root Locus Design GUI looks like this. The feedback structure selected in the Import LTI Design Model window The (red) squares are the closed-loop poles corresponding to the value of the compensator gain. . Notice that the design model name appears in the title bar.
8 The Root Locus Design GUI unstable if the gain set point is increased so as to place any of the closed-loop poles in the right-half of the complex plane. You can test the limits of stability by either: • Dragging the closed-loop poles around the locus to change the gain set point • Changing the gain set point manually In this example we do the former.
Controller Design Using the Root Locus Design GUI The following two figures capture this procedure. As you get close to a closedloop pole on the locus, the mouse pointer becomes a hand. The poles appear to be on the imaginary axis. The right-most pair of closed-loop poles seem to be on the imaginary axis. Actually, they are only close. Let’s use the zoom controls to improve this result.
8 The Root Locus Design GUI The Zoom controls are shown below. Zoom in the y-direction. Zoom in the x-direction. Zoom in both the x- and y- directions. Zoom out to full view. Once it is selected, you can operate any of the first three Zoom buttons in one of two ways. • Use your mouse to rubberband around the area on the plot region you want to focus on. Rubberbanding involves clicking and holding the mouse down, while dragging the mouse pointer around the region of interest on the screen.
Controller Design Using the Root Locus Design GUI The use of zooming to find the limits of stability on the root locus is depicted in the following three figures. After you drag that pole close to the imaginary axis, use the left-most Zoom button to zoom in. Move the closed-loop pole to the imaginary axis to fine-tune your gain adjustment.
8 The Root Locus Design GUI The limiting gain set point for stability Use the right-most Zoom button to zoom out. You can also store this zoom focus setting using the Axes settings described in the next section. The critical gain value for stability is approximately 43.9.
Controller Design Using the Root Locus Design GUI Storing and Retrieving Axes Limits You can store and change axes limits in two ways: • Using the Axes settings toolbar • Using the Set Axes Preferences menu item from the Tools menu For more information on square axes and/or equal axes, type help axis at the command line. To use the Axes settings toolbar to save the current axes limits, select the left-most button shown below. Save current axes limits. Retrieve a set of saved axes limits.
8 The Root Locus Design GUI You can also set or revise axes limits and other axes preferences in the Root Locus Axes Preferences window. To open this window, select Set Axes Preferences from the Tools menu. You can also use this window to change the colors of the root locus plot and the compensator poles and zeros, and the color and type of marker used for the closed-loop poles. Type of marker used to designate the closed-loop poles Select Apply to implement changes and keep this window open.
Controller Design Using the Root Locus Design GUI Change the gain to 20 by editing the text box next to Gain, and pressing the Enter key. Notice that the locations of the closed-loop poles on the root locus are recalculated for the new gain set point. Your Root Locus Design GUI looks like this now.
8 The Root Locus Design GUI Check the Step box, as shown below. Once you check the Step box, an LTI Viewer window opens. The plot type for this LTI Viewer is the step response and this cannot be changed. This closed-loop response does not meet the desired settling time requirement (.05 seconds or less). You can design a compensator so that you do meet the required specifications.
Controller Design Using the Root Locus Design GUI To increase the gain to about 33 using your mouse: 1 Hold your mouse button down as you click on any of the red squares. 2 Drag the square in the direction of increasing gain, without allowing any of the closed-loop poles to enter the right half of the complex plane. If, when you start to move the red square, you see the gain value in the Gain box decreasing, drag it in the opposite direction. 3 Release the mouse button when the gain is near 33.
8 The Root Locus Design GUI The step response plot on the dynamically linked LTI Viewer automatically updates when you release the mouse button. As you may have noticed, the response time decreases with increasing gain, while the overshoot increases. Here we no longer meet the overshoot requirement. Since this gain is already relatively large, it’s likely that we will not be able to meet both design requirements using only a gain for the compensator.
Controller Design Using the Root Locus Design GUI Let’s place approximate design region boundaries on our root locus plot based on our design specifications. To do so, select the Add Grid/Boundary menu item in the Tools menu. Our design specifications require that the (2 percent) settling time be less than .05 seconds, and the maximum overshoot less than 5 percent. For second-order systems, the overshoot requirement can be translated directly as a requirement on the damping ratio of about .7 (see [2]).
8 The Root Locus Design GUI After you press OK, the Root Locus Design GUI calculates and displays the specified boundaries. Damping ratio boundary Settling time boundary The design region lies within this wedge to the left of the settling time boundary. Not all four branches of the root locus are within the design region. Let’s try adding poles and zeros to our compensator to see if we can meet the design criteria.
Controller Design Using the Root Locus Design GUI Once you have the gain specified, you can add poles or zeros to the compensator. You can add poles and zeros to the compensator (or remove them) in two ways: • Use buttons on the root locus toolbar section of the GUI for pole/zero placement. • Use the Edit Compensator menu item on the Tools menu. The root locus toolbar is convenient for graphically placing compensator poles and zeros so that you meet the design specifications for a given gain set point.
8 The Root Locus Design GUI The default mode for the toolbar is the drag mode. In this mode, you can: • Click on a specific location on the root locus to place a closed-loop pole there (and consequently reassign the gain set point). • Drag any of the closed-loop poles along its branch of the root locus (also reassigning the gain set point). • Drag any of the compensator poles or zeros around the complex plane to change the compensator.
Controller Design Using the Root Locus Design GUI Try placing a pair of complex poles just above the right-most real closed-loop pole. The resulting root locus plot looks like this. Add compensator poles here. Similarly, to add a pair of complex zeros to the compensator: 1 Select the add zero button from the root locus toolbar (third button from the left). 2 Click on the plot region where you would like to add one of the complex zeros.
8 The Root Locus Design GUI Try adding a pair of complex zeros just to the left of and a little bit below the complex closed-loop poles closest to the imaginary axis. The resulting root locus and step response plots are as follows. Add compensator zeros here. If your step response is unstable, lower the gain. In this example, the resulting step response is stable, but it still doesn’t meet the design criteria. As you can see, the compensator design process can involve some trial and error.
Controller Design Using the Root Locus Design GUI Editing Compensator Pole and Zero Locations In this section we use the Edit Compensator window to design a compensator with the following characteristics: • Gain: 9.7 • Poles: – 110 ± 140i • Zeros: – 70 ± 270i You can access the Edit Compensator window in one of three ways: • Double-click on any of the Current Compensator text. • Click on the (blue) compensator block in the feedback structure on the GUI. • Select Edit Compensator from the Tools menu.
8 The Root Locus Design GUI You can use the Edit Compensator window to: • Edit the locations of compensator poles and zeros. • Add compensator poles and zeros. • Delete compensator poles and zeros. • Change the name of the compensator (This name is used when exporting the compensator). For this example, edit the poles and zeros to be at – 110 ± 140i, and – 70 ± 270i, respectively. Your GUI looks like this. If, in addition, you want to add poles to the compensator: 1 Click on Add Pole.
Controller Design Using the Root Locus Design GUI With the gain set point at 9.7, your root locus looks as follows. Note: Whenever the numerator or denominator are too long for the Current Compensator field, a default numK or denK is displayed. To display the actual numerator or denominator, resize the Root Locus Design GUI.
8 The Root Locus Design GUI To check where the closed-loop roots are with respect to the boundary of the design region, let’s zoom in a bit. We don’t have all of the closed-loop poles in the design region, but let’s see if we meet the step response specifications. The updated LTI Viewer plot of the step response looks like this.
Controller Design Using the Root Locus Design GUI The step response looks good. However, to be certain that you’ve met your design specifications, you can check the settling time and peak response characteristics from the right-click menu. To do this: 1 Right-click anywhere in the plot region. 2 Select Settling Time from the Characteristics menu. 3 Set your mouse pointer to the settling time marker and click to display the value. The LTI Viewer will look like this.
8 The Root Locus Design GUI Saving the Compensator and Models Now that you have successfully designed your compensator, you may want to save your design parameters for future implementation. You can do this by selecting Export from the File menu on the Root Locus Design GUI. The window shown below opens. Use your mouse to select or multiselect models to save. Save to the disk as a file with a .mat extension. Save to the workspace.
Controller Design Using the Root Locus Design GUI the compensator is now in the workspace, in the variable named comp. Then type comp to see that this variable is stored in zpk format.
8 The Root Locus Design GUI Additional Root Locus Design GUI Features This section describes several features of the Root Locus Design GUI not covered in the servomechanism example.
Additional Root Locus Design GUI Features Designating the Model Source The source of your design model data is indicated in the Import From field shown below. If you are loading your models from MAT-files or Simulink models saved on your disk, you are prompted to enter the name of the file in the editable text box below the Import From radio buttons. You can also select the Browse button to search for and select the file.
8 The Root Locus Design GUI Using the Help Menu Click on the Help menu and you find that it contains five menu items: • Main Help • Edit Compensator • Convert Model • Add Grid/Boundary • Set Axis Preferences The first menu item, Main Help, opens a help window that describes how to use the controls located on the Root Locus Design GUI. The remaining menu items provide additional information on the features you can access from the Tools menu.
Additional Root Locus Design GUI Features Erasing Compensator Poles and Zeros You can delete any compensator poles and zeros in one of two ways: • Check the associated Delete box on the Edit Compensator... window obtained from the Tools menu. • Click on the erase button (fourth from left on the root locus toolbar) and click on the pole or zero to delete.
8 The Root Locus Design GUI To view the current closed-loop poles: • Select List Closed-Loop Poles from the Tools menu. If you select this menu, a window listing the closed-loop poles associated with the current gain set point opens. • Click on the OK button to close the window. Note: You cannot edit the closed-loop pole locations listed in this window. To view the design model poles and zeros, you can either: • Select the List Model Poles/Zeros item from the Tools menu.
Additional Root Locus Design GUI Features The following window opens. For each component of the design model, this window tells you: • The component name • The type of LTI model that the component is (that is, ss, zpk, or tf) • The component’s pole and zero locations • The associated numerator and denominator of the model transfer function You can edit any of the model names in the textboxes provided. Selecting OK closes this window and implements the changes you made to design model names.
8 The Root Locus Design GUI Note: The names you enter in this window are only used when you generate a Simulink diagram of the closed-loop structure. To list the denominator, (respectively, the numerator, including the gain) associated with a given component of the design model, in the List Model Poles/Zeros window, click on Poles, (respectively, Zeros) of that component. When you do, a window providing the selected information opens.
Additional Root Locus Design GUI Features If you answered Yes, a Simulink diagram such as this appears. Pre-filter, F The root locus design compensator, K The plant, P (Gservo) Sensor dynamics, H The Simulink diagram is linked to the workspace, not to the Root Locus Design GUI. If you change your compensator design in the Root Locus Design GUI, you must export it to the workspace in order to reflect these changes in the Simulink diagram.
8 The Root Locus Design GUI The Convert Model/Compensator window, shown below, opens. Choose the discretization method here. Click on OK to convert the model to your specifications and close this window. Enter the sampling time. Critical frequency is only required for the Tustin method. Once you select a discretization method, sampling time, and critical frequency (if required), click on the OK button. The Convert Model/Compensator window closes and the discretization is performed.
Additional Root Locus Design GUI Features To test these features, select Clear Model. Notice that your compensator is not altered, and its poles and zeros remain plotted in the root locus plot region of the GUI. Now, select Clear Compensator. The Root Locus Design GUI returns to the state it was in when you first opened it and you are ready to start a new design.
8 The Root Locus Design GUI References [1] Clark, R.N., Control System Dynamics, Cambridge University Press, 1996. [2] Dorf, R.C. and R.H. Bishop, Modern Control Systems, 8th Edition, Addison Wesley, 1997.
9 Design Case Studies Yaw Damper for a 747 Jet Transport Open-Loop Analysis . . . . . . . . Root Locus Design . . . . . . . . . Washout Filter Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9-3 9-6 9-9 9-14 Hard-Disk Read/Write Head Controller . . . . . . . . 9-20 LQG Regulation . . . . . . Process and Disturbance Models LQG Design for the x-Axis . . . LQG Design for the y-Axis . . . Cross-Coupling Between Axes . MIMO LQG Design . . . . . . . . . . .
9 Design Case Studies This chapter contains four detailed case studies of control system design and analysis using the Control System Toolbox. Case Study 1 A yaw damper for a 747 jet transport aircraft that illustrates the classical design process. Case Study 2 A hard-disk read/write head controller that illustrates classical digital controller design. Case Study 3 LQG regulation of the beam thickness in a steel rolling mill.
Yaw Damper for a 747 Jet Transport Yaw Damper for a 747 Jet Transport This case study demonstrates the tools for classical control design by stepping through the design of a yaw damper for a 747 jet transport aircraft. The jet model during cruise flight at MACH = 0.8 and H = 40,000 ft. is A = [–0.0558 0.5980 –3.0500 0 B = [ 0.0729 –4.7500 1.5300 0 C = [0 0 1 0 D = [0 0 0 0]; –0.9968 –0.1150 0.3880 0.0805 0.0802 –0.0318 –0.4650 1.0000 0.0415 0 0 0]; 0.0000 0.00775 0.
9 Design Case Studies You can display the LTI model sys by typing sys. MATLAB responds with a = beta yaw roll phi beta -0.0558 0.598 -3.05 0 yaw -0.9968 -0.115 0.388 0.0805 beta yaw roll phi rudder 0.00729 -0.475 0.153 0 aileron 0 0.00775 0.143 0 yaw bank angle beta 0 0 yaw 1 0 rudder 0 0 aileron 0 0 roll 0.0802 -0.0318 -0.465 1 phi 0.0415 0 0 0 roll 0 0 phi 0 1 b = c = d = yaw bank angle Continuous-time model. The model has two inputs and two outputs.
Yaw Damper for a 747 Jet Transport Compute the open-loop eigenvalues and plot them in the s -plane. damp(sys) Eigenvalue -7.28e-003 -5.63e-001 -3.29e-002 + 9.47e-001i -3.29e-002 - 9.47e-001i Damping Freq. (rad/s) 1.00e+000 1.00e+000 3.48e-002 3.48e-002 7.28e-003 5.63e-001 9.47e-001 9.47e-001 pzmap(sys) This model has one pair of lightly damped poles. They correspond to what is called the “Dutch roll mode.
9 Design Case Studies Suppose you want to design a compensator that increases the damping of these poles, so that the resulting complex poles have a damping ratio ζ > 0.35 with natural frequency ω n < 1 rad/sec. You can do this using the Control System toolbox analysis tools. Open-Loop Analysis First perform some open-loop analysis to determine possible control strategies. Start with the time response (you could use step or impulse here).
Yaw Damper for a 747 Jet Transport concerned about the behavior during the first few seconds rather than the first few minutes. Next look at the response over a smaller time frame of 20 seconds. impulse(sys,20) Look at the plot from aileron (input 2) to bank angle (output 2). The aircraft is oscillating around a nonzero bank angle. Thus, the aircraft is turning in response to an aileron impulse. This behavior will prove important later in this case study.
9 Design Case Studies Typically, yaw dampers are designed using the yaw rate as sensed output and the rudder as control input. Look at the corresponding frequency response (input 1 to output 1). bode(sys(1,1)) From this Bode diagram, you can see that the rudder has significant effect around the lightly damped Dutch roll mode (that is, near ω = 1 rad/sec). To make the design easier, select the subsystem from rudder to yaw rate.
Yaw Damper for a 747 Jet Transport Root Locus Design Since the simplest compensator is a static gain, first try to determine appropriate gain values using the root locus technique.
9 Design Case Studies This is the root locus for negative feedback and shows that the system goes unstable almost immediately. If, instead, you use positive feedback, you may be able to keep the system stable. rlocus(–sys11) sgrid This looks better. Just using simple feedback you can achieve a damping ratio of ζ = 0.45 . You can graphically select some pole locations and determine the corresponding gain with rlocfind.
Yaw Damper for a 747 Jet Transport The '+' marks on the previous figure show a possible selection. The corresponding gain and closed-loop poles are k k = 2.7615 damp(poles) Eigenvalue -1.01e+000 -3.03e-001 + 6.18e-001i -3.03e-001 - 6.18e-001i -3.33e-001 Damping Freq. (rad/s) 1.00e+000 4.41e-001 4.41e-001 1.00e+000 1.01e+000 6.89e-001 6.89e-001 3.33e-001 Next, form the closed-loop system so that you can analyze this design.
9 Design Case Studies Plot the closed-loop impulse response for a duration of 20 seconds. impulse(cl11,20) The response settles quickly and does not oscillate much. Now close the loop on the original model and see how the response from the aileron looks.
Yaw Damper for a 747 Jet Transport feedback with index vectors selecting this input/output pair). At the MATLAB prompt, type cloop = feedback(sys,–k,1,1); damp(cloop) % closed-loop poles Eigenvalue -3.33e-001 -3.03e-001 + 6.18e-001i -3.03e-001 - 6.18e-001i -1.01e+000 Damping Freq. (rad/s) 1.00e+000 4.41e-001 4.41e-001 1.00e+000 3.33e-001 6.89e-001 6.89e-001 1.
9 Design Case Studies Plot the MIMO impulse response. impulse(cloop,20) Look at the plot from aileron (input 2) to bank angle (output 2). When you move the aileron, the system no longer continues to bank like a normal aircraft. You have over-stabilized the spiral mode. The spiral mode is typically a very slow mode and allows the aircraft to bank and turn without constant aileron input. Pilots are used to this behavior and will not like your design if it does not allow them to fly normally.
Yaw Damper for a 747 Jet Transport s H ( s ) = -----------s+a The washout filter places a zero at the origin, which constrains the spiral mode pole to remain near the origin. We choose a = 0.333 for a time constant of three seconds and use the root locus technique to select the filter gain k . First specify the fixed part s ⁄ ( s + a ) of the washout by H = zpk(0,–0.
9 Design Case Studies and draw another root locus for this open-loop model. rlocus(–oloop) sgrid Now the maximum damping is about ζ = 0.25 .
Yaw Damper for a 747 Jet Transport The selected pole locations are marked by '+' on the root locus above. The resulting gain value and closed-loop dynamics are found by typing k k = 2.5832 damp(poles) Eigenvalue -1.37e+000 -1.76e-001 + 6.75e-001i -1.76e-001 - 6.75e-001i -4.74e-001 -3.88e-003 Damping Freq. (rad/s) 1.00e+000 2.52e-001 2.52e-001 1.00e+000 1.00e+000 1.37e+000 6.98e-001 6.98e-001 4.74e-001 3.
9 Design Case Studies Look at the closed-loop response from rudder to yaw rate. cl11 = feedback(oloop,–k); impulse(cl11,20) The response settles nicely but has less damping than your previous design. Finally, you can verify that the washout filter has fixed the spiral mode problem. First form the complete washout filter kH ( s ) (washout + gain).
Yaw Damper for a 747 Jet Transport Then close the loop around the first I/O pair of the MIMO model sys and simulate the impulse response. cloop = feedback(sys,WOF,1,1); % Final closed-loop impulse response impulse(cloop,20) The bank angle response (output 2) due to an aileron impulse (input 2) now has the desired nearly constant behavior over this short time frame.
9 Design Case Studies Hard-Disk Read/Write Head Controller Hard Disk Drive Disk Platen Ω Disk Drive Motor Read/Write Head l θ Solenoid This case study demonstrates the ability to perform classical digital control design by going through the design of a computer hard-disk read/write head position controller.
Hard-Disk Read/Write Head Controller 2 Using the values J = 0.01 kg m , C = 0.004 Nm/(rad/sec), K = 10 Nm/rad, and K i = 0.05 Nm/rad, form the transfer function description of this system. At the MATLAB prompt, type J = num den H = .01; C = 0.004; K = 10; Ki = .05; = Ki; = [J C K]; tf(num,den) MATLAB responds with Transfer function: 0.05 ----------------------0.01 s^2 + 0.004 s + 10 The task here is to design a digital controller that provides accurate positioning of the read/write head.
9 Design Case Studies You can compare the Bode plots of the continuous and discretized models with bode(H,'-',Hd,'--') 9-22
Hard-Disk Read/Write Head Controller To analyze the discrete system, plot its step response, type step(Hd) The system oscillates quite a bit. This is probably due to very light damping. You can check this by computing the open-loop poles. Type % Open–loop poles of discrete model damp(Hd) Eigenvalue 9.87e–01 + 1.57e–01i 9.87e–01 – 1.57e–01i Magnitude 9.99e–01 9.99e–01 Equiv. Damping 6.32e–03 6.32e–03 Equiv. Freq. 3.16e+01 3.
9 Design Case Studies The simplest compensator is just a gain, so try the root locus technique to select an appropriate feedback gain. rlocus(Hd) As shown in the root locus, the poles quickly leave the unit circle and go unstable. You need to introduce some lead or a compensator with some zeros. Try the compensator z+a D ( z ) = -----------z+b with a = – 0.85 and b = 0 .
Hard-Disk Read/Write Head Controller The corresponding open-loop model u D(z) Hd( z ) Compensator Plant y is obtained by the series connection D = zpk(0.85,0,1,Ts) oloop = Hd * D Now see how this compensator modifies the open-loop frequency response.
9 Design Case Studies The plant response is the dashed line and the open-loop response with the compensator is the solid line. The plot above shows that the compensator has shifted up the phase plot (added lead) in the frequency range ω > 10 rad/sec.
Hard-Disk Read/Write Head Controller Now try the root locus again with the plant and compensator as open loop. rlocus(oloop) zgrid This time, the poles stay within the unit circle for some time (the lines drawn by zgrid show the damping ratios from ζ = 0 to 1 in steps of 0.1). This plot shows a set of poles '+' selected using rlocfind. At the MATLAB prompt, type [k,poles] = rlocfind(oloop) k k = 4.
9 Design Case Studies Type ddamp(poles,Ts) to see the equivalent damping and natural frequency for each of the eigenvalues. To analyze this design, form the closed-loop system and plot the closed-loop step response. cloop = feedback(oloop,k); step(cloop) This response depends on your closed loop set point. The one shown here is relatively fast and settles in about 0.07 seconds. Therefore, this closed loop disk drive system has a seek time of about 0.07 seconds.
Hard-Disk Read/Write Head Controller Now look at the robustness of your design. The most common classical robustness criteria are the gain and phase margins. Use the function margin to determine these margins. With output arguments, margin returns the gain and phase margins as well as the corresponding crossover frequencies. Without output argument, margin plots the Bode response and displays the margins graphically.
9 Design Case Studies The command produces the plot shown below. Bode Diagrams Gm=11.6 dB (Wcg=296.1); Pm=43.1 deg. (Wcp=106.5) 100 50 Phase (deg); Magnitude (dB) 0 −50 −100 100 0 −100 −200 −300 0 10 1 2 10 10 3 10 Frequency (rad/sec) This design is robust and can tolerate a 11 dB gain increase or a 40 degree phase lag in the open-loop system without going unstable.
LQG Regulation LQG Regulation This case study demonstrates the use of the LQG design tools in a process control application. The goal is to regulate the horizontal and vertical thickness of the beam produced by a hot steel rolling mill. This example is adapted from [1]. The full plant model is MIMO and the example shows the advantage of direct MIMO LQG design over separate SISO designs for each axis. Type milldemo at the command line to run this demonstration interactively.
9 Design Case Studies This shape is impressed by two pairs of rolling cylinders (one per axis) positioned by hydraulic actuators. The gap between the two cylinders is called the roll gap. rolling mill stand incoming beam shaped beam rolling cylinders x-axis The objective is to maintain the beam thickness along the x- and y-axes within the quality assurance tolerances.
LQG Regulation This leads to the following generic model for each axis of the rolling process.
9 Design Case Studies The model data for each axis is summarized below. Model Data for the x-Axis 8 2.4 × 10 H x ( s ) = -----------------------------------2 2 s + 72s + 90 4 10 F ix ( s ) = -------------------s + 0.05 4 3 × 10 s F ex ( s ) = ----------------------------------------2 2 s + 0.125s + 6 g x = 10 –6 Model Data for the y-Axis 8 7.8 × 10 H y ( s ) = -----------------------------------2 2 s + 71s + 88 4 2 × 10 F iy ( s ) = -------------------s + 0.
LQG Regulation Start with the x -axis. First specify the model components as transfer function objects. % Hydraulic actuator (with input "u-x") Hx = tf(2.4e8,[1 72 90^2],'inputname','u-x') % Input thickness/hardness disturbance model Fix = tf(1e4,[1 0.05],'inputn','w-ix') % Rolling eccentricity model Fex = tf([3e4 0],[1 0.125 6^2],'inputn','w-ex') % Gain from force to thickness gap gx = 1e–6; Next build the open-loop model shown in Figure 9-1 above.
9 Design Case Studies The variable Px now contains an open-loop state-space model complete with input and output names. Px.inputname ans = 'u-x' 'w-ex' 'w-ix' Px.outputname ans = 'x-gap' 'x-force' The second output 'x-force' is the rolling force measurement. The LQG regulator will use this measurement to drive the hydraulic actuator and reduce disturbance-induced thickness variations δ x .
LQG Regulation high-frequency content of δ x with the low-pass filter 30 ⁄ ( s + 30 ) and use the filtered value in the LQ performance criterion. lpf = tf(30,[1 30]) % Connect low-pass filter to first output of Px Pxdes = append(lpf,1) * Px set(Pxdes,'outputn',{'x-gap*' 'x-force'}) % Design the state-feedback gain using LQRY and q=1, r=1e–4 kx = lqry(Pxdes(1,1),1,1e–4) Note: lqry expects all inputs to be commands and all outputs to be measurements.
9 Design Case Studies Let’s look at the regulator Bode response between 0.1 and 1000 rad/sec. bode(Regx,{0.1 1000}) The phase response has an interesting physical interpretation. First, consider an increase in input thickness. This low-frequency disturbance boosts both output thickness and rolling force. Because the regulator phase is approximately 0o at low frequencies, the feedback loop then adequately reacts by increasing the hydraulic force to offset the thickness increase.
LQG Regulation Next, compare the open- and closed-loop responses from disturbance to thickness gap. Use feedback to close the loop. To help specify the feedback connection, look at the I/O names of the plant Px and regulator Regx. Px.inputname ans = 'u-x' 'w-ex' 'w-ix' Regx.outputname ans = 'u-x' Px.outputname ans = 'x-gap' 'x-force' Regx.inputname ans = 'x-force' This indicates that you must connect the first input and second output of Px to the regulator.
9 Design Case Studies You are now ready to compare the open- and closed-loop Bode responses from disturbance to thickness gap. bode(Px(1,2:3),'--',clx(1,2:3),'-',{0.1 100}) Bode Diagrams From: w−ex From: w−ix 0 −20 −40 Phase (deg); Magnitude (dB) −60 −80 −100 0 −50 To: x−gap −100 −150 −200 −250 −300 −1 10 0 10 1 10 2 −1 10 10 0 10 1 10 2 10 Frequency (rad/sec) The dashed lines show the open-loop response.
LQG Regulation simulation, and derive equivalent discrete white noise inputs for this sampling rate. dt = 0.01 t = 0:dt:50 % time samples % Generate unit-covariance driving noise wx = [w-ex;w-ix]. % Equivalent discrete covariance is 1/dt wx = sqrt(1/dt) * randn(2,length(t)) lsim(Px(1,2:3),':',clx(1,2:3),'-',wx,t) The dotted lines correspond to the open-loop response. In this simulation, the LQG regulation reduces the peak thickness variation by a factor 4.
9 Design Case Studies LQG Design for the y-Axis The LQG design for the y -axis (regulation of the y thickness) follows the exact same steps as for the x -axis. % Specify model components Hy = tf(7.8e8,[1 71 88^2],'inputn','u-y') Fiy = tf(2e4,[1 0.05],'inputn','w-iy') Fey = tf([1e5 0],[1 0.19 9.4^2],'inputn','w-ey') gy = 0.
LQG Regulation Compare the open- and closed-loop response to the white noise input disturbances. dt = 0.01 t = 0:dt:50 wy = sqrt(1/dt) * randn(2,length(t)) lsim(Py(1,2:3),':',cly(1,2:3),'-',wy,t) The dotted lines correspond to the open-loop response. The simulation results are comparable to those for the x -axis. Cross-Coupling Between Axes The x / y thickness regulation, is a MIMO problem. So far you have treated each axis separately and closed one SISO loop at a time.
9 Design Case Studies process exhibits some degree of cross-coupling between axes. Physically, an increase in hydraulic force along the x -axis compresses the material, which in turn boosts the repelling force on the y -axis cylinders. The result is an increase in y -thickness and an equivalent (relative) decrease in hydraulic force along the y-axis. The coupling between axes is as follows.
LQG Regulation Accordingly, the thickness gaps and rolling forces are related to the outputs δ x, f x, ... of the x- and y-axis models by δx δy fx 0 0 g yx g x δx 0 1 g xy g y 0 δy 0 0 1 – g yx fx 0 0 – g xy 1 fy fy = 1 cross-coupling matrix Let’s see how the previous “decoupled” LQG design fares when cross-coupling is taken into account. To build the two-axes model shown in Figure 9-2, append the models Px and Py for the x - and y-axes.
9 Design Case Studies You are now ready to simulate the open- and closed-loop responses to the driving white noises wx (for the x-axis) and wy (for the y-axis). wxy = [wx ; wy] lsim(Pc(1:2,3:6),':',cl(1:2,3:6),'-',wxy,t) The response reveals a severe deterioration in regulation performance along the x-axis (the peak thickness variation is about four times larger than in the simulation without cross-coupling).
LQG Regulation MIMO LQG Design Start with the complete two-axis state-space model Pc derived above. The model inputs and outputs are Pc.inputname ans = 'u-x' 'u-y' 'w-ex' 'w-ix' 'w_ey' 'w_iy' P.outputname ans = 'x-gap' 'y-gap' 'x-force' 'y-force' As earlier, add low-pass filters in series with the 'x-gap' and 'y-gap' outputs to penalize only low-frequency thickness variations. Pdes = append(lpf,lpf,eye(2)) * Pc Pdes.outputn = Pc.
9 Design Case Studies The resulting LQG regulator RegMIMO has two inputs and two outputs. RegMIMO.inputname ans = 'x-force' 'y-force' RegMIMO.outputname ans = 'u-x' 'u-y' Plot its singular value response (principal gains).
LQG Regulation Next, plot the open- and closed-loop time responses to the white noise inputs (using the MIMO LQG regulator for feedback). % Form the closed-loop model cl = feedback(Pc,RegMIMO,1:2,3:4,+1); % Simulate with LSIM using same noise inputs lsim(Pc(1:2,3:6),':',cl(1:2,3:6),'-',wxy,t) The MIMO design is a clear improvement over the separate SISO designs for each axis. In particular, the level of x / y thickness variation is now comparable to that obtained in the decoupled case.
9 Design Case Studies Kalman Filtering This final case study illustrates the use of the Control System Toolbox for Kalman filter design and simulation. Both steady-state and time-varying Kalman filters are considered. Consider the discrete plant x [ n + 1 ] = Ax [ n ] + B ( u [ n ] + w [ n ] ) y [ n ] = Cx [ n ] with additive Gaussian noise w [ n ] on the input u [ n ] and data A = [1.1269 1.0000 0 –0.4940 0 1.0000 0.1129 0 0]; B = [–0.3832 0.5919 0.
Kalman Filtering In these equations: • x̂ [ n n – 1 ] is the estimate of x [ n ] given past measurements up to y v [ n – 1 ] • x̂ [ n n ] is the updated estimate based on the last measurement y v [ n ] Given the current estimate x̂ [ n n ] , the time update predicts the state value at the next sample n + 1 (one-step-ahead predictor). The measurement update then adjusts this prediction based on the new measurement y v [ n + 1 ] . The correction term is a function of the innovation, that is, the discrepancy.
9 Design Case Studies This is done by % Note: set sample time to –1 to mark model as discrete Plant = ss(A,[B B],C,0,–1,'inputname',{'u' 'w'},... 'outputname','y'); Assuming that Q = R = 1 , you can now design the discrete Kalman filter by Q = 1; R = 1; [kalmf,L,P,M] = kalman(Plant,Q,R); This returns a state-space model kalmf of the filter as well as the innovation gain M M = 3.7980e–01 8.1732e–02 –2.
Kalman Filtering Because you are interested in the output estimate y e , keep only the first output of kalmf. Type kalmf = kalmf(1,:); kalmf a = x1_e x2_e x3_e x1_e 0.76830 0.62020 –0.08173 x2_e –0.49400 0 1.00000 x1_e x2_e x3_e u –0.38320 0.59190 0.51910 y 0.35860 0.37980 0.08173 y_e x1_e 0.62020 x2_e 0 y_e u 0 y 0.37980 x3_e 0.11290 0 0 b = c = x3_e 0 d = Sampling time: unspecified Discrete-time system.
9 Design Case Studies The block diagram below shows how to generate both true and filtered outputs. u Plant Process noise y Kalman filter yv ye Sensor noise y You can construct a state-space model of this block diagram with the functions parallel and feedback. First build a complete plant model with u, w, v as inputs and y and y v (measurements) as outputs. a b c d P = = = = = A; [B B 0*B]; [C;C]; [0 0 0;0 0 1]; ss(a,b,c,d,–1,'inputname',{'u' 'w' 'v'},...
Kalman Filtering Finally, close the sensor loop by connecting the plant output y v to the filter input y v with positive feedback. % Close loop around input #4 and output #2 SimModel = feedback(sys,1,4,2,1) % Delete yv from I/O list SimModel = SimModel([1 3],[1 2 3]) The resulting simulation model has w, v, u as inputs and y, y e as outputs. SimModel.input ans = 'w' 'v' 'u' SimModel.output ans = 'y' 'y_e' You are now ready to simulate the filter behavior.
9 Design Case Studies and compare the true and filtered responses graphically. subplot(211), plot(t,y,'--',t,ye,'-'), xlabel('No. of samples'), ylabel('Output') subplot(212), plot(t,y–yv,'-.',t,y–ye,'-'), xlabel('No. of samples'), ylabel('Error') The first plot shows the true response y (dashed line) and the filtered output y e (solid line). The second plot compares the measurement error (dash-dot) with the estimation error (solid). This plot shows that the noise level has been significantly reduced.
Kalman Filtering The error covariance before filtering (measurement error) is MeasErrCov MeasErrCov = 1.1138 while the error covariance after filtering (estimation error) is only EstErrCov EstErrCov = 0.2722 Time-Varying Kalman Filter The time-varying Kalman filter is a generalization of the steady-state filter for time-varying systems or LTI systems with nonstationary noise covariance.
9 Design Case Studies with x̂ [ n n – 1 ] and x̂ [ n n ] as defined on page 9-50, and in the following. T Q [ n ] = E ( w [ n ]w [ n ] ) T R [ n ] = E ( v [ n ]v [ n ] ) T P [ n n ] = E( { x [ n] – x [ n n ] }{ x [ n] – x [ n n ] } ) T P[n n – 1 ] = E( { x[n ] – x[n n – 1 ] }{ x[ n ] – x[ n n – 1 ]} ) For simplicity, we have dropped the subscripts indicating the time dependence of the state-space matrices.
Kalman Filtering you can implement the time-varying filter with the following for loop.
9 Design Case Studies You can now compare the true and estimated output graphically. subplot(211), plot(t,y,'--',t,ye,'-') subplot(212), plot(t,y-yv,'-.',t,y-ye,'-') The first plot shows the true response y (dashed line) and the filtered response y e (solid line). The second plot compares the measurement error (dash-dot) with the estimation error (solid).
Kalman Filtering The time-varying filter also estimates the covariance errcov of the estimation error y – y e at each sample. Plot it to see if your filter reached steady state (as you expect with stationary input noise). subplot(211) plot(t,errcov), ylabel('Error covar') From this covariance plot, you can see that the output covariance did indeed reach a steady state in about five samples. From then on, your time-varying filter has the same performance as the steady-state version.
9 Design Case Studies Compare with the estimation error covariance derived from the experimental data. Type EstErr = y–ye; EstErrCov = sum(EstErr.*EstErr)/length(EstErr) EstErrCov = 0.2718 This value is smaller than the theoretical value errcov and close to the value obtained for the steady-state design. Finally, note that the final value M [ n ] and the steady-state value M of the innovation gain matrix coincide. Mn, M Mn = 0.3798 0.0817 –0.2570 M = 0.3798 0.0817 –0.
Kalman Filtering References [1] Grimble, M.J., Robust Industrial Control: Optimal Design Approach for Polynomial Systems, Prentice Hall, 1994, p. 261 and pp. 443–456.
9 Design Case Studies 9-64
10 Reliable Computations Conditioning and Numerical Stability . . . . . . . . 10-4 Conditioning . . . . . . . . . . . . . . . . . . . . 10-4 Numerical Stability . . . . . . . . . . . . . . . . . . 10-6 Choice of LTI Model State Space . . . . . Transfer Function . . Zero-Pole-Gain Models Scaling Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10-8 . 10-8 . 10-8 10-14 . . . . . . . . . . . . . . . . . . . . .
10 Reliable Computations When working with low-order SISO models (less than five states), computers are usually quite forgiving and insensitive to numerical problems. You generally won't encounter any numerical difficulties and MATLAB will give you accurate answers regardless of the model or conversion method you choose. For high order SISO models and MIMO models, however, the finite-precision arithmetic of a computer is not so forgiving and you must exercise caution.
At the same time, it is important to appreciate the limitations of computer analyses. By following a few guidelines, you can avoid certain tools and models when they are likely to get you into trouble. The following sections try to illustrate, through examples, some of the numerical pitfalls to be avoided. We also encourage you to get the most out of the good algorithms by ensuring, if possible, that your models give rise to problems that are well-conditioned.
10 Reliable Computations Conditioning and Numerical Stability Two of the key concepts in numerical analysis are the conditioning of problems and the stability of algorithms. Conditioning Consider the linear system Ax = b given by A = 0.7800 0.9130 0.5630 0.6590 b = 0.2170 0.2540 The true solution is x = [1, –1]' and you can calculate it approximately using MATLAB. x = A\b x = 1.0000 –1.0000 format long, x x = 0.99999999991008 –0.
Conditioning and Numerical Stability Notice how much the small change in the data is magnified in the solution. One way to measure the magnification factor is by means of the quantity A A –1 called the condition number of A with respect to inversion. The condition number determines the loss in precision due to roundoff errors in Gaussian elimination and can be used to estimate the accuracy of results obtained from matrix inversion and linear equation solution.
10 Reliable Computations row of A. This perturbed matrix has n distinct eigenvalues λ 1, ..., λ n with λ k = 1 ⁄ 2 exp ( j2πk ⁄ n ) . Thus, you can see that this small perturbation in the n data has been magnified by a factor on the order of 2 to result in a rather large perturbation in the solution (the eigenvalues of A). Further details and related examples are to be found in [7].
Conditioning and Numerical Stability Using row 1 as the pivot row (i.e., subtracting 1000 times row 1 from row 2) you arrive at the equivalent triangular system. 0.001 1.000 x 1 0 – 1000 x 2 = 1.000 – 1000 Note that the coefficient multiplying x 2 in the second equation should be – 1001, but because of roundoff, becomes – 1000 . As a result, the second equation yields x 2 = 1.000 , a good approximation, but now back-substitution in the first equation 0.001x 1 = 1.000 – ( 1.000 ) ( 1.
10 Reliable Computations Choice of LTI Model Now turn to the implications of the results in the last section on the linear modeling techniques used for control engineering. The Control System Toolbox includes the following types of LTI models that are applicable to discussions of computational reliability: • State space • Transfer function, polynomial form • Transfer function, factored zero-pole-gain form The following subsections show that state space is most preferable for numerical computations.
Choice of LTI Model A major difficulty is the extreme sensitivity of the roots of a polynomial to its coefficients. This example is adapted from Wilkinson, [6] as an illustration. Consider the transfer function 1 1 H ( s ) = ------------------------------------------------------------- = ---------------------------------------------------------20 19 ( s + 1 ) ( s + 2 )... ( s + 20 ) s + 210s + ... + 20! The A matrix of the companion realization of H ( s ) is A = 0 0 : 0 – 20! 1 0 : 0 . 0 1 . ... ... .
10 Reliable Computations very little. This is true in general. Different roots have different sensitivities to different perturbations. Computed roots may then be quite meaningless for a polynomial, particularly high-order, with imprecisely known coefficients. Finding all the roots of a polynomial (equivalently, the poles of a transfer function or the eigenvalues of a matrix in controllable or observable canonical form) is often an intrinsically sensitive problem.
Choice of LTI Model and the conversion from zero-pole-gain to state space incurs no loss of accuracy in the poles. format long e [sort(eig(a1)) sort(eig(a2))] ans = 9.999999999998792e-01 2.000000000001984e+00 3.000000000475623e+00 3.999999981263996e+00 5.000000270433721e+00 5.999998194359617e+00 7.000004542844700e+00 8.000013753274901e+00 8.999848908317270e+00 1.000059459550623e+01 1.099854678336595e+01 1.200255822210095e+01 1.299647702454549e+01 1.400406940833612e+01 1.499604787386921e+01 1.
10 Reliable Computations Its eigenvectors and eigenvalues are given as follows. [v,d] = eig(A) v = 0.7071 –0.7071 0.0000 –0.0000 –0.0000 0.0000 0.7071 –0.7071 –0.3162 –0.3162 0.6325 0.6325 0.6325 0.6325 0.3162 0.3162 0 2.0000 0 0 0 0 5.0000 0 0 0 0 10.0000 d = 1.0000 0 0 0 The condition number (with respect to inversion) of the eigenvector matrix is cond(v) ans = 1.000 Now convert a state-space model with the above A matrix to transfer function form, and back again to state-space form.
Choice of LTI Model Note also that the eigenvectors have changed. [vc,dc] = eig(Ac) vc = –0.5017 –0.8026 –0.3211 –0.0321 0.2353 0.7531 0.6025 0.1205 0.0510 0.4077 0.8154 0.4077 0.0109 0.1741 0.6963 0.6963 dc = 10.0000 0 0 0 0 5.0000 0 0 0 0 2.0000 0 0 0 0 1.0000 The condition number of the new eigenvector matrix cond(vc) ans = 34.5825 is thirty times larger. The phenomenon illustrated above is not unusual.
10 Reliable Computations • The pole locations are very sensitive to the coefficients of the denominator polynomial. • The balanced companion form produced by ss, while better than the standard companion form, often results in ill-conditioned eigenproblems, especially with higher-order systems. The above statements hold even for systems with distinct poles, but are particularly relevant when poles are multiple.
Scaling Scaling State space is the preferred model for LTI systems, especially with higher order models. Even with state-space models, however, accurate results are not guaranteed, because of the finite-word-length arithmetic of the computer. A well-conditioned problem is usually a prerequisite for obtaining accurate results. You should generally normalize or scale the ( A, B, C, D ) matrices of a system to improve their conditioning.
10 Reliable Computations excellent discussion of scaling is given in the introduction to the LINPACK Users’ Guide, [1]. Choose scaling based upon physical insight to the problem at hand. If you choose not to scale, and for many small problems scaling is not necessary, be aware that this choice affects the accuracy of your answers. Finally, note that the function ssbal performs automatic scaling of the state vector.
Summary Summary This chapter has described numerous things that can go wrong when performing numerical computations. You won’t encounter most of these difficulties when you solve practical lower-order problems. The problems described here pertain to all computer analysis packages. MATLAB has some of the best algorithms available, and, where possible, notifies you when there are difficulties.
10 Reliable Computations References [1] Dongarra, J.J., J.R. Bunch, C.B. Moler, and G.W. Stewart, LINPACK Users Guide, SIAM Publications, Philadelphia, PA, 1978. [2] Franklin, G.F. and J.D. Powell, Digital Control of Dynamic Systems, Addison-Wesley, 1980. [3] Kailath, T., Linear Systems, Prentice-Hall, 1980. [4] Laub, A.J., “Numerical Linear Algebra Aspects of Control Design Computations,” IEEE Transactions on Automatic Control, Vol. AC-30, No. 2, February 1985, pp. 97-108. [5] Wilkinson, J.H.
11 Reference
11 Reference This chapter contains detailed descriptions of all Control System Toolbox functions. It begins with a list of functions grouped by subject area and continues with the reference entries in alphabetical order. Information is also available through the online Help facility.
Category Tables Category Tables Table 11-1: LTI Models Function Name Description drss Generate random discrete state-space model. dss Create descriptor state-space model. filt Create discrete filter with DSP convention. frd Create a frequency response data (FRD) model. frdata Retrieve data from an FRD model. get Query LTI model properties. ltimodels Information on LTI models ltiprops Information on all LTI properties. set Set LTI or response object properties.
11 Reference Table 11-2: Model Characteristics Function Name Description class Display model type ('tf', 'zpk', 'ss', or 'frd'). hasdelay Test true if LTI model has any type of delay. isa Test true if LTI model is of specified type. isct Test true for continuous-time models. isdt Test true for discrete-time models. isempty Test true for empty LTI models. isproper Test true for proper LTI models. issiso Test true for SISO models. ndims Get the number of model/array dimensions.
Category Tables Table 11-3: Model Conversion (Continued) Function Name Description residue Provide partial fraction expansion (see Using MATLAB). ss Convert to a state space model. tf Convert to a transfer function model. zpk Convert to a zero-pole-gain model. Table 11-4: Model Order Reduction Function Name Description balreal Calculate an I/O balanced realization. minreal Calculate minimal realization or pole/zero cancellation. modred Delete states in I/O balanced realization.
11 Reference Table 11-5: State-Space Realizations (Continued) Function Name Description ss2ss State coordinate transformation. ssbal Diagonal balancing of state-space realizations. Table 11-6: Model Dynamics Function Name Description damp Calculate natural frequency and damping. dcgain Calculate low-frequency (DC) gain. covar Calculate covariance of response to white noise. dsort Sort discrete-time poles by magnitude. esort Sort continuous-time poles by real part.
Category Tables Table 11-7: Model Building (Continued) Function Name Description connect Connect the subsystems of a block-diagonal model according to an interconnection scheme of your choice. conv Convolve two polynomials (see Using MATLAB). drmodel, drss Generate random discrete-time model. feedback Calculate the feedback connection of models. lft Calculate the star product (LFT interconnection). ord2 Generate second-order model. pade Compute the Padé approximation of time delays.
11 Reference Table 11-8: Time Response (Continued) Function Name Description ltiview Open the LTI Viewer for linear response analysis. step Calculate step response. Table 11-9: Frequency Response 11-8 Function Name Description bode Calculate bode plot. evalfr Evaluate response at single complex frequency. freqresp Evaluate frequency response for selected frequencies. linspace Create a vector of evenly spaced frequencies. logspace Create a vector of logarithmically spaced frequencies.
Category Tables Table 11-9: Frequency Response (Continued) Function Name Description sigma Calculate singular value plot. zgrid Superimpose z-plane grid on root locus or pole/zero map. Table 11-10: Pole Placement Function Name Description acker Calculate SISO pole placement design. place Calculate MIMO pole placement design. estim Form state estimator given estimator gain. reg Form output-feedback compensator given state-feedback and estimator gains.
11 Reference Table 11-11: LQG Design (Continued) Function Name Description kalmd Calculate the discrete Kalman estimator for continuous models. lqgreg Form LQG regulator given LQ gain and Kalman filter. Table 11-12: Equation Solvers Function Name Description care Solve continuous-time algebraic Riccati equations. dare Solve discrete-time algebraic Riccati equations. lyap Solve continuous-time Lyapunov equations. dlyap Solve discrete-time Lyapunov equations.
acker Purpose 11acker Pole placement design for single-input systems Syntax k = acker(A,b,p) Description Given the single-input system · x = Ax + bu and a vector p of desired closed-loop pole locations, acker (A,b,p)uses Ackermann’s formula [1] to calculate a gain vector k such that the state feedback u = – kx places the closed-loop poles at the locations p. In other words, the eigenvalues of A – bk match the entries of p (up to ordering).
append Purpose 11append Group LTI models by appending their inputs and outputs Syntax sys = append(sys1,sys2,...,sysN) Description append appends the inputs and outputs of the LTI models sys1,...,sysN to form the augmented model sys depicted below. u1 sys1 y1 u2 sys2 y2 : : uN yN sysN sys For systems with transfer functions H 1 ( s ) ,..., H N ( s ) , the resulting system sys has the block-diagonal transfer function H1 ( s ) 0 : 0 0 .. 0 H2 ( s ) . : . .. .
append · A1 0 x1 B 0 u1 x1 = + 1 · 0 A2 x2 0 B2 u2 x2 y1 y2 Arguments = C1 0 x1 0 C2 x2 + D1 0 u1 0 D2 u2 The input arguments sys1,..., sysN can be LTI models of any type. Regular matrices are also accepted as a representation of static gains, but there should be at least one LTI object in the input list. The LTI models should be either all continuous, or all discrete with the same sample time.
append c = y1 y2 y3 x1 1.00000 0 0 x2 0 0 3.00000 y1 y2 y3 u1 0 0 0 u2 0 10.00000 0 d = u3 0 0 4.00000 Continuous-time system.
augstate Purpose 11augstate Append the state vector to the output vector Syntax asys = augstate(sys) Description Given a state-space model sys with equations · x = Ax + Bu y = Cx + Du (or their discrete-time counterpart), augstate appends the states x to the outputs y to form the model · x = Ax + Bu y C D x+ = u x I 0 This command prepares the plant so that you can use the feedback command to close the loop on a full-state feedback u = – K x .
balreal Purpose 11balreal Input/output balancing of state-space realizations Syntax sysb = balreal(sys) [sysb,g,T,Ti] = balreal(sys) Description sysb = balreal(sys) produces a balanced realization sysb of the LTI model sys with equal and diagonal controllability and observability gramians (see gram for a definition of gramian). balreal handles both continuous and discrete systems. If sys is not a state-space model, it is first and automatically converted to state space using ss.
balreal which indicates that the last two states of sysb are weakly coupled to the input and output. You can then delete these states by sysr = modred(sysb,[2 3],'del') to obtain the following first-order approximation of the original system. zpk(sysr) Zero/pole/gain: 1.0001 -------(s+4.97) Compare the Bode responses of the original and reduced-order models.
balreal Algorithm Consider the model · x = Ax + Bu y = Cx + Du with controllability and observability gramians W c and W o . The state coordinate transformation x = Tx produces the equivalent model · –1 x = TAT x + TBu –1 y = CT x + Du and transforms the gramians to T W c = TW c T , Wo = T –T WoT –1 The function balreal computes a particular similarity transformation T such that W c = W o = diag ( g ) See [1,2] for details on the algorithm. Limitations The LTI model sys must be stable.
bode Purpose Syntax 11bode Compute the Bode frequency response of LTI models bode(sys) bode(sys,w) bode(sys1,sys2,...,sysN) bode(sys1,sys2,...,sysN,w) bode(sys1,'PlotStyle1',...,sysN,'PlotStyleN') [mag,phase,w] = bode(sys) Description bode computes the magnitude and phase of the frequency response of LTI models. When invoked without left-hand arguments, bode produces a Bode plot on the screen. The magnitude is plotted in decibels (dB), and the phase in degrees.
bode uses red dashed lines for the first system sys1 and green 'x' markers for the second system sys2. When invoked with left-hand arguments [mag,phase,w] = bode(sys) [mag,phase] = bode(sys,w) return the magnitude and phase (in degrees) of the frequency response at the frequencies w (in rad/sec). The outputs mag and phase are 3-D arrays with the frequency as the last dimension (see “Arguments” below for details).
bode 2 s + 0.1s + 7.5 H ( s ) = -------------------------------------------4 3 2 s + 0.12s + 9s by typing g = tf([1 0.1 7.5],[1 0.12 9 0 0]); bode(g) Bode Diagrams 40 20 Phase (deg); Magnitude (dB) 0 −20 −40 0 −50 −100 −150 −200 −1 10 0 10 1 10 Frequency (rad/sec) To plot the response on a wider frequency range, for example, from 0.1 to 100 rad/sec, type bode(g,{0.
bode You can also discretize this system using zero-order hold and the sample time T s = 0.5 second, and compare the continuous and discretized responses by typing gd = c2d(g,0.5) bode(g,'r',gd,'b--') Bode Diagrams 50 Phase (deg); Magnitude (dB) 0 −50 −100 0 −50 −100 −150 −200 −250 −300 −1 10 0 10 Frequency (rad/sec) Algorithm For continuous-time systems, bode computes the frequency response by evaluating the transfer function H ( s ) on the imaginary axis s = jω .
bode structure. The reduction to Hessenberg form provides a good compromise between efficiency and reliability. See [1] for more details on this technique. For discrete-time systems, the frequency response is obtained by evaluating the transfer function H ( z ) on the unit circle. To facilitate interpretation, the upper-half of the unit circle is parametrized as z=e jωT s , π 0 ≤ ω ≤ ω N = -----Ts where T s is the sample time. ω N is called the Nyquist frequency.
c2d Purpose 11c2d Discretize continuous-time systems Syntax sysd = c2d(sys,Ts) sysd = c2d(sys,Ts,method) Description sysd = c2d(sys,Ts) discretizes the continuous-time LTI model sys using zero-order hold on the inputs and a sample time of Ts seconds. sysd = c2d(sys,Ts,method) gives access to alternative discretization schemes. The string method selects the discretization method among the following: 'zoh' Zero-order hold. The control inputs are assumed piecewise constant over the sampling period Ts.
c2d with input delay T d = 0.35 second. To discretize this system using the triangle approximation with sample time T s = 0.1 second, type H = tf([1 –1],[1 4 5],'inputdelay',0.35) Transfer function: s - 1 exp(-0.35*s) * ------------s^2 + 4 s + 5 Hd = c2d(H,0.1,'foh') Transfer function: 0.0115 z^3 + 0.0456 z^2 – 0.0562 z – 0.009104 --------------------------------------------z^6 – 1.629 z^5 + 0.6703 z^4 Sampling time: 0.
c2d The next command compares the continuous and discretized step responses. step(H,'-',Hd,'--') See Also d2c d2d References [1] Franklin, G.F., J.D. Powell, and M.L. Workman, Digital Control of Dynamic Systems, Second Edition, Addison-Wesley, 1990.
canon Purpose 11canon Compute canonical state-space realizations Syntax csys = canon(sys,'type') [csys,T] = canon(sys,'type') Description canon computes a canonical state-space model for the continuous or discrete LTI system sys. Two types of canonical forms are supported.
canon For state-space models sys, [csys,T] = canon(a,b,c,d,'type') also returns the state coordinate transformation T relating the original state vector x and the canonical state vector x c. x c = Tx This syntax returns T=[] when sys is not a state-space model. Algorithm Transfer functions or zero-pole-gain models are first converted to state space using ss. The transformation to modal form uses the matrix P of eigenvectors of the A matrix.
care Purpose Syntax 11care Solve continuous-time algebraic Riccati equations (CARE) [X,L,G,rr] = care(A,B,Q) [X,L,G,rr] = care(A,B,Q,R,S,E) [X,L,G,report] = care(A,B,Q,...,'report') [X1,X2,L,report] = care(A,B,Q,...,'implicit') Description [X,L,G,rr] = care(A,B,Q) computes the unique solution X of the algebraic Riccati equation T T Ric ( X ) = A X + XA – XBB X + Q = 0 T such that A – BB X has all its eigenvalues in the open left-half plane.
care Alternatively, [X1,X2,L,report] = care(A,B,Q,...,'implicit') also turns off error messages but now returns X in implicit form. –1 X = X2 X 1 Note that this syntax returns report = 0 when successful. Examples Example 1 Given A = –3 2 1 1 B = 0 1 C = 1 –1 R = 3 you can solve the Riccati equation T –1 T T A X + XA – XBR B X + C C = 0 by a = [–3 2;1 1] b = [0 ; 1] c = [1 –1] r = 3 [x,l,g] = care(a,b,c'*c,r) This yields the solution x x = 0.5895 1.8216 1.8216 8.
care Finally, note that the variable l contains the closed-loop eigenvalues eig(a– b*g). l l = –3.5026 –1.
care Q S T >0 S R See Also dare lyap References [1] Arnold, W.F., III and A.J. Laub, “Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations,” Proc. IEEE, 72 (1984), pp. 1746–1754.
chgunits Purpose 11chgunits Convert the frequency units of an FRD model Syntax sys = chgunits(sys,units) Description sys = chgunits(sys,units) converts the units of the frequency points stored in an FRD model, sys to units, where units is either of the strings 'Hz' or 'rad/s'. This operation changes the assigned frequencies by applying the appropriate (2*pi) scaling factor, and the 'Units' property is updated. If the 'Units' field already matches units, no conversion is made.
connect Purpose 11connect Derive state-space model from block diagram description Syntax sysc = connect(sys,Q,inputs,outputs) Description Complex dynamical systems are often given in block diagram form. For systems of even moderate complexity, it can be quite difficult to find the state-space model required in order to bring certain analysis and design tools into use. Starting with a block diagram description, you can use append and connect to construct a state-space model of the system.
connect external outputs are outputs 2 and 7 of sys, then inputs and outputs should be set to inputs = [1 2 15]; outputs = [2 7]; The final model sysc has these particular inputs and outputs. Since it is easy to make a mistake entering all the data required for a large model, be sure to verify your model in as many ways as you can. Here are some suggestions: • Make sure the poles of the unconnected model sys match the poles of the various blocks in the diagram.
connect Given the matrices of the state-space model sys2 A = [ –9.0201 17.7791 –1.6943 3.2138 ]; B = [ –.5112 .5362 –.002 –1.8470]; C = [ –3.2897 2.4544 –13.5009 18.0745]; D = [–.5476 –.1410 –.6459 .2958 ]; Define the three blocks as individual LTI models. sys1 = tf(10,[1 5],'inputname','uc') sys2 = ss(A,B,C,D,'inputname',{'u1' 'u2'},... 'outputname',{'y1' 'y2'}) sys3 = zpk(–1,–2,2) Next append these blocks to form the unconnected model sys.
connect c = ? y1 y2 ? x1 2.5 0 0 0 x2 0 -3.2897 -13.501 0 x3 0 2.4544 18.075 0 x4 0 0 0 -1.4142 ? y1 y2 ? uc 0 0 0 0 u1 0 -0.5476 -0.6459 0 u2 0 -0.141 0.2958 0 ? 0 0 0 2 d = Continuous-time system. Note that the ordering of the inputs and outputs is the same as the block ordering you chose. Unnamed inputs or outputs are denoted by ?. To derive the overall block diagram model from sys, specify the interconnections and the external inputs and outputs.
connect You can obtain a state-space model for the overall interconnection by typing sysc = connect(sys,Q,inputs,outputs) a = x1 x2 x3 x4 x1 -5 0.84223 -2.9012 0.65708 x1 x2 x3 x4 uc 4 0 0 0 y1 y2 x1 -0.22148 0.46463 y1 y2 uc 0 0 x2 0 0.076636 -33.029 -11.996 x3 0 5.6007 45.164 16.06 x4 0 0.47644 -1.6411 -1.6283 x3 5.6568 11.356 x4 -0.12529 0.26283 b = u1 0 -0.076001 -1.5011 -0.57391 c = x2 -5.6818 -8.4826 d = u1 -0.66204 -0.40582 Continuous-time system.
connect References [1] Edwards, J.W., “A Fortran Program for the Analysis of Linear Continuous and Sampled-Data Systems,” NASA Report TM X56038, Dryden Research Center, 1976.
covar Purpose 11covar Output and state covariance of a system driven by white noise Syntax [P,Q] = covar(sys,W) Description covar calculates the stationary covariance of the output y of an LTI model sys driven by Gaussian white noise inputs w . This function handles both continuous- and discrete-time cases.
covar and MATLAB returns p = 30.3167 You can compare this output of covar to simulation results. randn('seed',0) w = sqrt(5)∗randn(1,1000); % 1000 samples % Simulate response to w with LSIM: y = lsim(sys,w); % Compute covariance of y values psim = sum(y .∗ y)/length(w); This yields psim = 32.6269 The two covariance values p and psim do not agree perfectly due to the finite simulation horizon. Algorithm Transfer functions and zero-pole-gain models are first converted to state space with ss.
covar T and P is given by P = CQC + DWD T Note that P is well defined for nonzero D in the discrete case. Limitations The state and output covariances are defined for stable systems only. For continuous systems, the output response covariance P is finite only when the D matrix is zero (strictly proper system). See Also dlyap lyap References [1] Bryson, A.E. and Y.C. Ho, Applied Optimal Control, Hemisphere Publishing, 1975, pp. 458-459.
ctrb Purpose 11ctrb Form the controllability matrix Syntax Co = ctrb(A,B) Co = ctrb(sys) Description ctrb computes the controllability matrix for state-space systems. For an n-by-n matrix A and an n-by-m matrix B, ctrb(A,B) returns the controllability matrix 2 Co = B AB A B … A n–1 B (11-1) where Co has n rows and nm columns. Co = ctrb(sys) calculates the controllability matrix of the state-space LTI object sys. This syntax is equivalent to executing Co = ctrb(sys.A,sys.
ctrb Limitations The calculation of Co may be ill-conditioned with respect to inversion. An indication of this can be seen from this simple example. A = 1 δ , 01 B = 1 δ This pair is controllable if δ ≠ 0 but if δ < eps , where eps is the relative machine precision. ctrb(A,B) returns B AB = 11 δ δ which is not full rank. For cases like these, it is better to determine the controllability of a system using ctrbf.
ctrbf Purpose Syntax Description 11ctrbf Compute the controllability staircase form [Abar,Bbar,Cbar,T,k] = ctrbf(A,B,C) [Abar,Bbar,Cbar,T,k] = ctrbf(A,B,C,tol) If the controllability matrix of ( A, B ) has rank r ≤ n , where n is the size of A , then there exists a similarity transformation such that T A = TAT , B = TB, C = CT T where T is unitary, and the transformed system has a staircase form, in which the uncontrollable modes, if there are any, are in the upper left corner.
ctrbf Example Compute the controllability staircase form for A = 1 4 1 –2 1 1 –1 –1 1 0 0 1 B = C = and locate the uncontrollable mode. [Abar,Bbar,Cbar,T,k]=ctrbf(A,B,C) Abar = –3.0000 –3.0000 0 2.0000 Bbar = 0.0000 1.4142 0.0000 –1.4142 Cbar = –0.7071 0.7071 0.7071 0.7071 T = –0.7071 0.7071 k = 1 0 0.7071 0.7071 The decomposed system Abar shows an uncontrollable mode located at –3 and a controllable mode located at 2.
ctrbf Algorithm ctrbf is an M-file that implements the Staircase Algorithm of [1]. See Also ctrb minreal References [1] Rosenbrock, M.M., State-Space and Multivariable Theory, John Wiley, 1970.
d2c Purpose 11d2c Convert discrete-time LTI models to continuous time Syntax sysc = d2c(sysd) sysc = d2c(sysd,method) Description d2c converts LTI models from discrete to continuous time using one of the following conversion methods: 'zoh' Zero-order hold on the inputs. The control inputs are assumed piecewise constant over the sampling period. 'tustin' Bilinear (Tustin) approximation to the derivative. 'prewarp' Tustin approximation with frequency prewarping.
d2c As with zero-order hold, the inverse discretization operation c2d(Hc,0.1,'tustin') gives back the original H ( z ) . Algorithm The 'zoh' conversion is performed in state space and relies on the matrix logarithm (see logm in Using MATLAB). Limitations The Tustin approximation is not defined for systems with poles at z = – 1 and is ill-conditioned for systems with poles near z = –1 . The zero-order hold method cannot handle systems with poles at z = 0 .
d2c Convert Hc back to discrete time by typing c2d(Hc,Ts) yielding Zero/pole/gain: (z+0.5) (z+0.2) ------------------------(z+0.5)^2 (z^2 + z + 0.4) Sampling time: 0.1 This discrete model coincides with H ( z ) after canceling the pole/zero pair at z = – 0.5 . See Also c2d d2d logm References [1] Franklin, G.F., J.D. Powell, and M.L. Workman, Digital Control of Dynamic Systems, Second Edition, Addison-Wesley, 1990.
d2d Purpose 11d2d Resample discrete-time LTI models or add input delays Syntax sys1 = d2d(sys,Ts) Description sys1 = d2d(sys,Ts) resamples the discrete-time LTI model sys to produce an equivalent discrete-time model sys1 with the new sample time Ts (in seconds). The resampling assumes zero-order hold on the inputs and is equivalent to consecutive d2c and c2d conversions. sys1 = c2d(d2c(sys),Ts) Example Consider the zero-pole-gain model z – 0.7 H ( z ) = ----------------z – 0.5 with sample time 0.
damp Purpose 11damp Compute damping factors and natural frequencies Syntax [Wn,Z] = damp(sys) [Wn,Z,P] = damp(sys) Description damp calculates the damping factor and natural frequencies of the poles of an LTI model sys. When invoked without lefthand arguments, a table of the eigenvalues in increasing frequency, along with their damping factors and natural frequencies, is displayed on the screen.
damp Type damp(H) and MATLAB returns Eigenvalue Damping –1.00e+000 + 1.41e+000i –1.00e+000 – 1.41e+000i See Also eig esort,dsort pole pzmap zero Freq. (rad/s) 5.77e–001 5.77e–001 1.73e+000 1.
dare Purpose Syntax 11dare Solve discrete-time algebraic Riccati equations (DARE) [X,L,G,rr] = dare(A,B,Q,R) [X,L,G,rr] = dare(A,B,Q,R,S,E) [X,L,G,report] = dare(A,B,Q,...,'report') [X1,X2,L,report] = dare(A,B,Q,...
dare and L = eig(A–B*G,E). Two additional syntaxes are provided to help develop applications such as H ∞ -optimal control design. [X,L,G,report] = dare(A,B,Q,...,'report') turns off the error messages when the solution X fails to exist and returns a failure report instead.
dare References 11-56 [1] Arnold, W.F., III and A.J. Laub, “Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations,” Proc. IEEE, 72 (1984), pp. 1746–1754.
dcgain Purpose 11dcgain Compute low frequency (DC) gain of LTI system Syntax k = dcgain(sys) Description k = dcgain(sys) computes the DC gain k of the LTI model sys. Continuous Time The continuous-time DC gain is the transfer function value at the frequency s = 0 . For state-space models with matrices ( A, B, C, D ) , this value is –1 K = D – CA B Discrete Time The discrete-time DC gain is the transfer function value at z = 1 .
delay2z Purpose Replace delays of discrete-time TF, SS, or ZPK models by poles at z=0, or replace delays of FRD models by a phase shift Syntax sys = delay2z(sys) Description sys = delay2z(sys) maps all time delays to poles at z=0 for discrete-time TF, ZPK, or SS models sys. Specifically, a delay of k sampling periods is replaced by (1/z)^k in the transfer function corresponding to the model.
dlqr Purpose Syntax Description 11dlqr Design linear-quadratic (LQ) state-feedback regulator for discrete-time plant [K,S,e] = dlqr(a,b,Q,R) [K,S,e] = dlqr(a,b,Q,R,N) [K,S,e] = dlqr(a,b,Q,R,N) calculates the optimal gain matrix K such that the state-feedback law u [ n ] = – Kx [ n ] minimizes the quadratic cost function ∞ ∑ J( u ) = T T T ( x [ n ] Qx [ n ] + u [ n ] Ru [ n ] + 2x [ n ] Nu [ n ] ) n=1 for the discrete-time state-space mode l x [ n + 1 ] = Ax [ n ] + Bu [ n ] The default value N=
dlqr lqr lqrd lqry 11-60 State-feedback LQ regulator for continuous plant Discrete LQ regulator for continuous plant State-feedback LQ regulator with output weighting
dlyap Purpose 11dlyap Solve discrete-time Lyapunov equations Syntax X = dlyap(A,Q) Description dlyap solves the discrete-time Lyapunov equation T A XA – X + Q = 0 where A and Q are n -by- n matrices. The solution X is symmetric when Q is symmetric, and positive definite when Q is positive definite and A has all its eigenvalues inside the unit disk. Diagnostics The discrete-time Lyapunov equation has a (unique) solution if the eigenvalues α 1, α 2, ...
drmodel, drss Purpose Syntax 11drmodel, drss Generate stable random discrete test models sys sys sys sys = = = = drss(n) drss(n,p) drss(n,p,m) drss(n,p,m,s1,...sn) [num,den] = drmodel(n) [A,B,C,D] = drmodel(n) [A,B,C,D] = drmodel(n,p,m) Description sys = drss(n) produces a random n-th order stable model with one input and one output, and returns the model in the state-space object sys. drss(n,p) produces a random n-th order stable model with one input and p outputs.
drmodel, drss Example Generate a random discrete LTI system with three states, two inputs, and two outputs. sys = drss(3,2,2) a = x1 x2 x3 x1 0.38630 –0.23390 –0.03412 x2 –0.21458 –0.15220 0.11394 x1 x2 x3 u1 0.98833 0 0.42350 u2 0.51551 0.33395 0.43291 y1 y2 x1 0.22595 0 x2 0.76037 0 y1 y2 u1 0 0.78333 u2 0.68085 0.46110 x3 –0.09914 –0.06572 –0.22618 b = c = x3 0 0 d = Sampling time: unspecified Discrete-time system.
dsort Purpose 11dsort Sort discrete-time poles by magnitude Syntax s = dsort(p) [s,ndx] = dsort(p) Description dsort sorts the discrete-time poles contained in the vector p in descending order by magnitude. Unstable poles appear first. When called with one lefthand argument, dsort returns the sorted poles in s. [s,ndx] = dsort(p) also returns the vector ndx containing the indices used in the sort. Example Sort the following discrete poles. p = –0.2410 + 0.5573i –0.2410 – 0.5573i 0.1503 –0.
dss Purpose Syntax 11dss Specify descriptor state-space models sys = dss(a,b,c,d,e) sys = dss(a,b,c,d,e,Ts) sys = dss(a,b,c,d,e,ltisys) sys = dss(a,b,c,d,e,'Property1',Value1,...,'PropertyN',ValueN) sys = dss(a,b,c,d,e,Ts,'Property1',Value1,...,'PropertyN',ValueN) Description sys = dss(a,b,c,d,e) creates the continuous-time descriptor state-space model · Ex = Ax + Bu y = Cx + Du The E matrix must be nonsingular. The output sys is an SS model storing the model data (see “LTI Objects” on page 2-3).
dss Example The command sys = dss(1,2,3,4,5,'td',0.1,'inputname','voltage',... 'notes','Just an example') creates the model · 5x = x + 2u y = 3x + 4u with a 0.1 second input delay. The input is labeled 'voltage', and a note is attached to tell you that this is just an example.
dssdata Purpose 11dssdata Quick access to descriptor state-space data Syntax [a,b,c,d,e] = dssdata(sys) [a,b,c,d,e,Ts] = dssdata(sys) Description [a,b,c,d,e] = dssdata(sys) extracts the descriptor matrix data ( A, B, C, D, E ) from the state-space model sys. If sys is a transfer function or zero-pole-gain model, it is first converted to state space. Note that dssdata is then equivalent to ssdata because it always returns E = I .
esort Purpose 11esort Sort continuous-time poles by real part Syntax s = esort(p) [s,ndx] = esort(p) Description esort sorts the continuous-time poles contained in the vector p by real part. Unstable eigenvalues appear first and the remaining poles are ordered by decreasing real parts. When called with one left-hand argument, s = esort(p) returns the sorted eigenvalues in s. [s,ndx] = esort(p) returns the additional argument ndx, a vector containing the indices used in the sort.
esort See Also dsort, sort eig pole pzmap zero Sort system poles Calculate eigenvalues and eigenvectors Compute system poles Pole-zero map Compute (transmission) zeros 11-69
estim Purpose 11estim Form state estimator given estimator gain Syntax est = estim(sys,L) est = estim(sys,L,sensors,known) Description est = estim(sys,L) produces a state/output estimator est given the plant state-space model sys and the estimator gain L. All inputs w of sys are assumed stochastic (process and/or measurement noise), and all outputs y are measured. The estimator est is returned in state-space form (SS object).
estim · x̂ = Ax̂ + B 2 u + L ( y – C 2 x̂ – D 22 u ) ŷ = x̂ C2 I x̂ + D 22 u 0 u (known) ŷ est y (sensors) x̂ estim handles both continuous- and discrete-time cases. You can use the functions place (pole placement) or kalman (Kalman filtering) to design an adequate estimator gain L . Note that the estimator poles (eigenvalues of A – LC ) should be faster than the plant dynamics (eigenvalues of A ) to ensure accurate estimation.
evalfr Purpose 11evalfr Evaluate frequency response at a single (complex) frequency Syntax frsp = evalfr(sys,f) Description frsp = evalfr(sys,f) evaluates the transfer function of the TF, SS, or ZPK model sys at the complex number f. For state-space models with data ( A, B, C, D ) , the result is –1 H ( f ) = D + C ( fI – A ) B evalfr is a simplified version of freqresp meant for quick evaluation of the response at a single point.
feedback Purpose Syntax Description 11feedback Feedback connection of two LTI models sys = feedback(sys1,sys2) sys = feedback(sys1,sys2,sign) sys = feedback(sys1,sys2,feedin,feedout,sign) sys = feedback(sys1,sys2) returns an LTI model sys for the negative feedback interconnection. + u sys1 y – sys2 The closed-loop model sys has u as input vector and y as output vector. The LTI models sys1 and sys2 must be both continuous or both discrete with identical sample times.
feedback computes a closed-loop model sys for the more general feedback loop. v z + sys1 u y – sys2 sys The vector feedin contains indices into the input vector of sys1 and specifies which inputs u are involved in the feedback loop. Similarly, feedout specifies which outputs y of sys1 are used for feedback. The resulting LTI model sys has the same inputs and outputs as sys1 (with their order preserved).
feedback Examples Example 1 + torque G velocity – H To connect the plant 2 2s + 5s + 1 G ( s ) = -----------------------------2 s + 2s + 3 with the controller 5(s + 2) H ( s ) = -------------------s + 10 using negative feedback, type G = tf([2 5 1],[1 2 3],'inputname','torque',... 'outputname','velocity'); H = zpk(–2,–10,5) Cloop = feedback(G,H) and MATLAB returns Zero/pole/gain from input "torque" to output "velocity": 0.18182 (s+10) (s+2.281) (s+0.2192) ----------------------------------(s+3.
feedback Example 2 Consider a state-space plant P with five inputs and four outputs and a state-space feedback controller K with three inputs and two outputs.
filt Purpose Syntax 11filt Specify discrete transfer functions in DSP format sys = filt(num,den) sys = filt(num,den,Ts) sys = filt(M) sys = filt(num,den,'Property1',Value1,...,'PropertyN',ValueN) sys = filt(num,den,Ts,'Property1',Value1,...
filt MIMO transfer functions are regarded as arrays of SISO transfer functions (one per I/O channel), each of which is characterized by its numerator and denominator. The input arguments num and den are then cell arrays of row vectors such that: • num and den have as many rows as outputs and as many columns as inputs. • Their ( i, j ) entries num{i,j} and den{i,j} specify the numerator and denominator of the transfer function from input j to output i.
frd Purpose Syntax 11frd Create a frequency response data (FRD) object or convert another model type to an FRD model sys sys sys sys = = = = frd(response,frequency) frd(response,frequency,Ts) frd frd(response,frequency,ltisys) sysfrd = frd(sys,frequency) sysfrd = frd(sys,frequency,'Units',units) Description sys = frd(response,frequency) creates an FRD model sys from the frequency response data stored in the multidimensional array response.
frd sysfrd = frd(sys,frequency,'Units',units)converts an FRD model from a TF, SS, or ZPK model while specifying the units for frequency to be units ('rad/s' or 'Hz'). Arguments When you specify a SISO or MIMO FRD model, or an array of FRD models, the input argument frequency is always a vector of length Nf, where Nf is the number of frequency data points in the FRD. The specification of the input argument response is summarized in the following table.
frd tf zpk Create transfer functions Create zero-pole-gain models 11-81
frdata Purpose 11frdata Quick access to data for a frequency response data object Syntax [response,freq] = frdata(sys) [response,freq,Ts] = frdata(sys) [response,freq] = frdata(sys,'v') Description [response,freq] = frdata(sys) returns the response data and frequency samples of the FRD model sys. For an FRD model with Ny outputs and Nu inputs at Nf frequencies: • response is an Ny-by-Nu-by-Nf multidimensional array where the (i,j) entry specifies the response from input j to output i.
frdata returns the FRD model data resp = 0.2040 + 0.4565i 2.4359 - 4.
freqresp Purpose 11freqresp Compute frequency response over grid of frequencies Syntax H = freqresp(sys,w) Description H = freqresp(sys,w) computes the frequency response of the LTI model sys at the real frequency points specified by the vector w. The frequencies must be in radians/sec. For single LTI Models, freqresp(sys,w) returns a 3-D array H with the frequency as the last dimension (see “Arguments” below). For LTI arrays of size [Ny Nu S1 ... Sn], freqresp(sys,w) returns a [Ny–by–Nu–by– S1–by–.
freqresp P( s ) = 0 1 -----------s+1 s–1 -----------s+2 1 at the frequencies ω = 1, 10, 100 . Type w = [1 10 100] H = freqresp(P,w) H(:,:,1) = 0 –0.2000+ 0.6000i 0.5000– 0.5000i 1.0000 H(:,:,2) = 0 0.9423+ 0.2885i 0.0099– 0.0990i 1.0000 H(:,:,3) = 0 0.9994+ 0.0300i 0.0001– 0.0100i 1.
freqresp Algorithm For transfer functions or zero-pole-gain models, freqresp evaluates the numerator(s) and denominator(s) at the specified frequency points. For continuous-time state-space models ( A, B, C, D ) , the frequency response is –1 D + C ( jω – A ) B , ω = ω 1, ..., ω N When numerically safe, A is diagonalized for fast evaluation of this expression at the frequencies ω 1, ..., ω N .
gensig Purpose 11gensig Generate test input signals for lsim Syntax [u,t] = gensig(type,tau) [u,t] = gensig(type,tau,Tf,Ts) Description [u,t] = gensig(type,tau) generates a scalar signal u of class type and with period tau (in seconds). The following types of signals are available. type = 'sin' Sine wave. type = 'square' Square wave. type = 'pulse' Periodic pulse. gensig returns a vector t of time samples and the vector u of signal values at these samples.
gensig Plot the resulting signal. plot(t,u) axis([0 30 –1 2]) 2 1.5 1 0.5 0 −0.
get Purpose 11get Access/query LTI property values Syntax Value = get(sys,'PropertyName') get(sys) Description Value = get(sys,'PropertyName') returns the current value of the property PropertyName of the LTI model sys. The string 'PropertyName' can be the full property name (for example, 'UserData') or any unambiguous case-insensitive abbreviation (for example, 'user').
get or query only about the numerator and sample time values by get(h,'num') ans = [1x2 double] and get(h,'ts') ans = 0.1000 Because the numerator data (num property) is always stored as a cell array, the first command evaluates to a cell array containing the row vector [0 1]. Remark An alternative to the syntax Value = get(sys,'PropertyName') is the structure-like referencing Value = sys.PropertyName For example, sys.Ts sys.a sys.
gram Purpose 11gram Compute controllability and observability gramians Syntax Wc = gram(sys,'c') Wo = gram(sys,'o') Description gram calculates controllability and observability gramians. You can use gramians to study the controllability and observability properties of state-space models and for model reduction [1,2]. They have better numerical properties than the controllability and observability matrices formed by ctrb and obsv.
gram to compute the gramians of a continuous or discrete system. The LTI model sys must be in state-space form. Algorithm The controllability gramian W c is obtained by solving the continuous-time Lyapunov equation T AW c + W c A + BB T = 0 or its discrete-time counterpart T AW c A – W c + BB T = 0 Similarly, the observability gramian W o solves the Lyapunov equation T T A Wo + W o A + C C = 0 in continuous time, and the Lyapunov equation T T A Wo A – Wo + C C = 0 in discrete time.
hasdelay Purpose 11hasdelay Test if an LTI model has time delays Syntax hasdelay(sys) Description hasdelay(sys) returns 1 (true) if the LTI model sys has input delays, output delays, or I/O delays, and 0 (false) otherwise.
impulse Purpose Syntax 11impulse Compute the impulse response of LTI models impulse(sys) impulse(sys,t) impulse(sys1,sys2,...,sysN) impulse(sys1,sys2,...,sysN,t) impulse(sys1,'PlotStyle1',...,sysN,'PlotStyleN') [y,t,x] = impulse(sys) Description impulse calculates the unit impulse response of a linear system. The impulse response is the response to a Dirac input δ ( t ) for continuous-time systems and to a unit pulse at t = 0 for discrete-time systems.
impulse As with bode or plot, you can specify a particular color, linestyle, and/or marker for each system, for example, impulse(sys1,'y:',sys2,'g--') See “Plotting and Comparing Multiple Systems” on page 5-13 and the bode entry in this chapter for more details.
impulse use the following commands. a = [–0.5572 –0.7814;0.7814 b = [1 –1;0 2]; c = [1.9691 6.4493]; sys = ss(a,b,c,0); impulse(sys) 0]; Impulse Response From: U1 From: U2 12 10 8 Amplitude 6 4 2 0 −2 −4 −6 0 5 10 15 200 5 10 15 Time (sec.) The left plot shows the impulse response of the first input channel, and the right plot shows the impulse response of the second input channel.
impulse Because this system has two inputs, y is a 3-D array with dimensions size(y) ans = 101 1 2 (the first dimension is the length of t). The impulse response of the first input channel is then accessed by y(:,:,1) Algorithm Continuous-time models are first converted to state space. The impulse response of a single-input state-space model · x = Ax + bu y = Cx is equivalent to the following unforced response with initial state b .
initial Purpose Syntax 11initial Compute the initial condition response of state-space models initial(sys,x0) initial(sys,x0,t) initial(sys1,sys2,...,sysN,x0) initial(sys1,sys2,...,sysN,x0,t) initial(sys1,'PlotStyle1',...,sysN,'PlotStyleN',x0) [y,t,x] = initial(sys,x0) Description initial calculates the unforced response of a state-space model with an initial condition on the states. · x = Ax , x ( 0 ) = x0 y = Cx This function is applicable to either continuous- or discrete-time models.
initial To plot the initial condition responses of several LTI models on a single figure, use initial(sys1,sys2,...,sysN,x0) initial(sys1,sys2,...,sysN,x0,t) (see impulse for details). When invoked with lefthand arguments, [y,t,x] = initial(sys,x0) [y,t,x] = initial(sys,x0,t) return the output response y, the time vector t used for simulation, and the state trajectories x. No plot is drawn on the screen. The array y has as many rows as time samples (length of t) and as many columns as outputs.
initial to the initial condition 1 0 x( 0) = a = [–0.5572 –0.7814;0.7814 c = [1.9691 6.4493]; x0 = [1 ; 0] 0]; sys = ss(a,[],c,[]); initial(sys,x0) Initial Condition Results 5 4 3 Amplitude 2 1 0 −1 −2 0 2 4 6 8 10 12 14 Time (sec.
inv Purpose 11inv Invert LTI systems Syntax isys = inv(sys) Description inv inverts the input/output relation y = G ( s )u –1 to produce the LTI system with the transfer matrix H ( s ) = G ( s ) . u = H ( s )y This operation is defined only for square systems (same number of inputs and outputs) with an invertible feedthrough matrix D . inv handles both continuous- and discrete-time systems.
inv You can verify that H * Hi is the identity transfer function (static gain I). Limitations Do not use inv to model feedback connections such as + G – H While it seems reasonable to evaluate the corresponding closed-loop transfer –1 function ( I + GH ) G as inv(1+g*h) * g this typically leads to nonminimal closed-loop models. For example, g = zpk([],1,1) h = tf([2 1],[1 0]) cloop = inv(1+g*h) * g yields a third-order closed-loop model with an unstable pole-zero cancellation at s = 1.
inv Use feedback or star to avoid such pitfalls.
isct, isdt Purpose 11isct, isdt Determine whether an LTI model is continuous or discrete Syntax boo = isct(sys) boo = isdt(sys) Description boo = isct(sys) returns 1 (true) if the LTI model sys is continuous and 0 (false) otherwise. sys is continuous if its sample time is zero, that is, sys.Ts=0. boo = isdt(sys) returns 1 (true) if sys is discrete and 0 (false) otherwise.
isempty Purpose 11isempty Test if an LTI model is empty Syntax boo = isempty(sys) Description isempty(sys) returns 1 (true) if the LTI model sys has no input or no output, and 0 (false) otherwise. Example Both commands isempty(tf) % tf by itself returns an empty transfer function isempty(ss(1,2,[],[])) return 1 (true) while isempty(ss(1,2,3,4)) returns 0 (false).
isproper Purpose Syntax Description 11isproper Test if an LTI model is proper boo = isproper(sys) isproper(sys) returns 1 (true) if the LTI model sys is proper and 0 (false) otherwise. State-space models are always proper. SISO transfer functions or zero-pole-gain models are proper if the degree of their numerator is less than or equal to the degree of their denominator. MIMO transfer functions are proper if all their SISO entries are proper.
issiso Purpose Syntax Description 11issiso Test if an LTI model is single-input/single-output (SISO) boo = issiso(sys) issiso(sys) returns 1 (true) if the LTI model sys is SISO and 0 (false) otherwise.
kalman Purpose 11kalman Design continuous- or discrete-time Kalman estimator Syntax [kest,L,P] = kalman(sys,Qn,Rn,Nn) [kest,L,P,M,Z] = kalman(sys,Qn,Rn,Nn) % discrete time only [kest,L,P] = kalman(sys,Qn,Rn,Nn,sensors,known) Description kalman designs a Kalman state estimator given a state-space model of the plant and the process and measurement noise covariance data. The Kalman estimator is the optimal solution to the following continuous or discrete estimation problems.
kalman the output and state estimates ŷ and x̂ .
kalman and generates optimal “current” output and state estimates ŷ [ n n ] and x̂ [ n n ] using all available measurements including y v [ n ] . The gain matrices L and M are derived by solving a discrete Riccati equation. The innovation gain M is used to update the prediction x̂ [ n n – 1 ] using the new measurement y v [ n ] .
kalman for more general plants sys where the known inputs u and stochastic inputs w are mixed together, and not all outputs are measured. The index vectors sensors and known then specify which outputs y of sys are measured and which inputs u are known. All other inputs are assumed stochastic. Example See examples on “Control Design Tools” on page 1-20, “LQG Design for the x-Axis” on page 9-34, and “Kalman Filtering” on page 9-50.
kalmd Purpose 11kalmd Design discrete Kalman estimator for continuous plant Syntax [kest,L,P,M,Z] = kalmd(sys,Qn,Rn,Ts) Description kalmd designs a discrete-time Kalman estimator that has response characteristics similar to a continuous-time estimator designed with kalman. This command is useful to derive a discrete estimator for digital implementation after a satisfactory continuous estimator has been designed.
kalmd lqgreg lqrd References Assemble LQG regulator Discrete LQ-optimal gain for continuous plant [1] Franklin, G.F., J.D. Powell, and M.L. Workman, Digital Control of Dynamic Systems, Second Edition, Addison-Wesley, 1990. [2] Van Loan, C.F., “Computing Integrals Involving the Matrix Exponential,” IEEE Trans. Automatic Control, AC-15, October 1970.
kalmd 11-114
lft Purpose 11lft Redheffer star product (linear fractional transformation) of two LTI models Syntax sys = lft(sys1,sys2) sys = lft(sys1,sys2,nu,ny) Description lft forms the star product or linear fractional transformation (LFT) of two LTI models or LTI arrays. Such interconnections are widely used in robust control techniques. sys = lft(sys1,sys2,nu,ny) forms the star product sys of the two LTI models (or LTI arrays) sys1 and sys2.
lft The abbreviated syntax sys = lft(sys1,sys2) produces: • The lower LFT of sys1 and sys2 if sys2 has fewer inputs and outputs than sys1. This amounts to deleting w 2 and z 2 in the above diagram. • The upper LFT of sys1 and sys2 if sys1 has fewer inputs and outputs than sys2. This amounts to deleting w 1 and z 1 in the above diagram.
lqgreg Purpose Syntax 11lqgreg Form LQG regulator given state-feedback gain and Kalman estimator rlqg = lqgreg(kest,k) rlqg = lqgreg(kest,k,'current') % discrete-time only rlqg = lqgreg(kest,k,controls) Description lqgreg forms the LQG regulator by connecting the Kalman estimator designed with kalman and the optimal state-feedback gain designed with lqr, dlqr, or lqry. The LQG regulator minimizes some quadratic cost function that trades off regulation performance and control effort.
lqgreg Process noise Plant y u u –K x̂ Kalman filter + yv + Measurement noise LQG regulator In discrete time, you can form the LQG regulator using either the prediction x̂ [ n n – 1 ] of x [ n ] based on measurements up to y v [ n – 1 ] , or the current state estimate x̂ [ n n ] based on all available measurements including y v [ n ] .
lqgreg • Continuous regulator for continuous plant: use lqr or lqry and kalman. • Discrete regulator for discrete plant: use dlqr or lqry and kalman. • Discrete regulator for continuous plant: use lqrd and kalmd. In discrete time, lqgreg produces the regulator u [ n ] = – Kx̂ [ n n – 1 ] by default (see “Description”). To form the “current” LQG regulator instead, use u [ n ] = – Kx̂ [ n n ] the syntax rlqg = lqgreg(kest,k,'current') This syntax is meaningful only for discrete-time problems.
lqgreg u ud yv Kalman –K estimator u x̂ LQG regulator Example See the examples “Control Design Tools” on page 1-20 and “LQG Regulation” on page 9-31.
lqr Purpose Syntax Description 11lqr Design linear-quadratic (LQ) state-feedback regulator for continuous plant [K,S,e] = lqr(A,B,Q,R) [K,S,e] = lqr(A,B,Q,R,N) [K,S,e] = lqr(A,B,Q,R,N) calculates the optimal gain matrix K such that the state-feedback law u = – Kx minimizes the quadratic cost function J( u ) = ∞ ∫0 ( x T T T Qx + u Ru + 2x Nu ) dt · for the continuous-time state-space model x = Ax + Bu The default value N=0 is assumed when N is omitted.
lqr See Also 11-122 care dlqr lqgreg lqrd lqry Solve continuous Riccati equations State-feedback LQ regulator for discrete plant Form LQG regulator Discrete LQ regulator for continuous plant State-feedback LQ regulator with output weighting
lqrd Purpose 11lqrd Design discrete LQ regulator for continuous plant Syntax [Kd,S,e] = lqrd(A,B,Q,R,Ts) [Kd,S,e] = lqrd(A,B,Q,R,N,Ts) Description lqrd designs a discrete full-state-feedback regulator that has response characteristics similar to a continuous state-feedback regulator designed using lqr. This command is useful to design a gain matrix for digital implementation after a satisfactory continuous state-feedback gain has been designed.
lqrd With the notation Φ( τ) = e Γ(τ) = Aτ τ ∫0 e Ad = Φ ( Ts ) , Aη B dη , Bd = Γ ( T s ) the discretized plant has equations x [ n + 1 ] = Ad x [ n ] + B d u [ n ] and the weighting matrices for the equivalent discrete cost function are Qd Nd T Nd Rd Ts = ∫0 T Q N Φ( τ) Γ(τ ) dτ T T 0 I Γ (τ) I N R Φ (τ) 0 The integrals are computed using matrix exponential formulas due to Van Loan (see [2]).
lqry Purpose 11lqry Linear-quadratic (LQ) state-feedback regulator with output weighting Syntax [K,S,e] = lqry(sys,Q,R) [K,S,e] = lqry(sys,Q,R,N) Description Given the plant · x = Ax + Bu y = Cx + Du or its discrete-time counterpart, lqry designs a state-feedback control u = – Kx that minimizes the quadratic cost function with output weighting J( u ) = ∞ ∫0 ( y T T T Qy + u Ru + 2y Nu ) dt (or its discrete-time counterpart).
lsim Purpose Syntax 11lsim Simulate LTI model response to arbitrary inputs lsim(sys,u,t) lsim(sys,u,t,x0) lsim(sys1,sys2,...,sysN,u,t) lsim(sys1,sys2,...,sysN,u,t,x0) lsim(sys1,'PlotStyle1',...,sysN,'PlotStyleN',u,t) [y,t,x] = lsim(sys,u,t,x0) Description lsim simulates the (time) response of continuous or discrete linear systems to arbitrary inputs. When invoked without left-hand arguments, lsim plots the response on the screen.
lsim simulates the responses of several LTI models to the same input history t,u and plots these responses on a single figure. As with bode or plot, you can specify a particular color, linestyle, and/or marker for each system, for example, lsim(sys1,'y:',sys2,'g--',u,t,x0) The multisystem behavior is similar to that of bode or step.
lsim Then simulate with lsim.
lsim Note that if you begin the simulation with a nonzero value for t(1), then the zero-initial condition response is shifted in time, as shown below. lsim(H,u,t+4) Algorithm Discrete-time systems are simulated with ltitr (state space) or filter (transfer function and zero-pole-gain). Continuous-time systems are discretized with c2d using either the 'zoh' or 'foh' method ('foh' is used for smooth input signals and 'zoh' for discontinuous signals such as pulses or square waves).
lsim To illustrate why resampling is sometimes necessary, consider the second-order model 2 ω - , H ( s ) = -----------------------------2 2 s + 2s + ω ω = 62.83 To simulate its response to a square wave with period 1 second, you can proceed as follows: w2 = 62.83^2 h = tf(w2,[1 2 w2]); t = 0:0.1:5; u = (rem(t,1)>=0.
lsim The response exhibits strong oscillations. Less obvious from this plot is the fact that lsim has resampled the input to reveal the oscillatory behavior. To see this, discretize H ( s ) using the sampling period 0.1 second (spacing in your t vector) and simulate the response of the discretized model: hd = c2d(h,0.
lsim The two responses look quite different. To clarify this discrepancy, superimpose the two plots by lsim(h,'b--',hd,'r-',u,t) The cause is now obvious: hd is undersampled and its response (solid line) masks the intersample oscillations of the continuous model H ( s ) . By comparing the suggested sampling dt=t(2)–t(1) against the system dynamics, lsim detects such undersampling and resamples the input to produce accurate continuous-time simulations.
ltiview Purpose Syntax 11ltiview Initialize an LTI Viewer for LTI system response analysis ltiview ltiview(plottype,sys) ltiview(plottype,sys,extras) ltiview(plottype,sys1,sys2,...sysN) ltiview(plottype,sys1,sys2,...sysN,extras) ltiview(plottype,sys1,PlotStyle1,sys2,PlotStyle2,...) Description ltiview when invoked without input arguments, initializes a new LTI Viewer for LTI system response analysis. Only frequency-domain analysis functions can be applied to FRDs.
ltiview extras is one or more input arguments as specified by the function named in plottype. These arguments may be required or optional, depending on the type of LTI response. For example, if plottype is 'step' then extras may be the desired final time, Tfinal, as shown below. ltiview('step',sys,Tfinal) However, if plottype is 'initial', the extras arguments must contain the initial conditions x0 and may contain other arguments, such as Tfinal.
lyap Purpose 11lyap Solve continuous-time Lyapunov equations Syntax X = lyap(A,Q) X = lyap(A,B,C) Description lyap solves the special and general forms of the Lyapunov matrix equation. Lyapunov equations arise in several areas of control, including stability theory and the study of the RMS behavior of systems. X = lyap(A,Q) solves the Lyapunov equation T AX + XA + Q = 0 where A and Q are square matrices of identical sizes. The solution X is a symmetric matrix if Q is.
lyap References [1] Bartels, R.H. and G.W. Stewart, “Solution of the Matrix Equation AX + XB = C,” Comm. of the ACM, Vol. 15, No. 9, 1972. [2] Bryson, A.E. and Y.C. Ho, Applied Optimal Control, Hemisphere Publishing, 1975. pp. 328–338.
margin Purpose 11margin Compute gain and phase margins and associated crossover frequencies Syntax [Gm,Pm,Wcg,Wcp] = margin(sys) [Gm,Pm,Wcg,Wcp] = margin(mag,phase,w) margin(sys) Description margin calculates the gain margin, phase margin, and associated crossover frequencies of SISO open-loop models. The gain and phase margins indicate the relative stability of the control system when the loop is closed.
margin Example You can compute the gain and phase margins of the open-loop discrete-time transfer function. Type hd = tf([0.04798 0.0464],[1 –1.81 0.9048],0.1) MATLAB responds with Transfer function: 0.04798 z + 0.0464 --------------------z^2 – 1.81 z + 0.9048 Sampling time: 0.1 Type [Gm,Pm,Wcg,Wcp] = margin(hd); [Gm,Pm,Wcg,Wcp] and MATLAB returns ans = 2.0517 11-138 13.5712 5.4374 4.
margin You can also display these margins graphically. margin(hd) Bode Diagrams Gm=6.3 dB (Wcg=5.4); Pm=13.6 deg. (Wcp=4.4) 50 Phase (deg); Magnitude (dB) 0 −50 −100 0 −50 −100 −150 −200 −250 −300 0 10 1 10 2 10 Frequency (rad/sec) Algorithm The phase margin is computed using H ∞ theory, and the gain margin by solving H ( jω ) = H ( jω ) for the frequency ω .
minreal Purpose 11minreal Minimal realization or pole-zero cancellation Syntax sysr = minreal(sys) sysr = minreal(sys,tol) Description sysr = minreal(sys) eliminates uncontrollable or unobservable state in state-space models, or cancels pole-zero pairs in transfer functions or zero-pole-gain models. The output sysr has minimal order and the same response characteristics as the original model sys. sysr = minreal(sys,tol) specifies the tolerance used for state elimination or pole-zero cancellation.
minreal See Also balreal modred sminreal Gramian-based input/output balancing Model order reduction Structured model reduction 11-141
modred Purpose 11modred Model order reduction Syntax rsys = modred(sys,elim) rsys = modred(sys,elim,'mdc') rsys = modred(sys,elim,'del') Description modred reduces the order of a continuous or discrete state-space model sys. This function is usually used in conjunction with balreal.
modred The last three diagonal entries of the balanced gramians are small, so eliminate the last three states with modred using both matched DC gain and direct deletion methods. hmdc = modred(hb,2:4,'mdc') hdel = modred(hb,2:4,'del') Both hmdc and hdel are first-order models. Compare their Bode responses against that of the original model h ( s ) .
modred The reduced-order model hdel is clearly a better frequency-domain approximation of h ( s ) . Now compare the step responses. step(h,'-',hmdc,'-.',hdel,'--') Step Response 0.3 0.25 Amplitude 0.2 0.15 0.1 0.05 0 −0.05 0 0.5 1 1.5 2 2.5 3 Time (sec.) While hdel accurately reflects the transient behavior, only hmdc gives the true steady-state response. Algorithm The algorithm for the matched DC gain method is as follows.
modred · A 11 A 12 x 1 B x1 = + 1 u · A 21 A 22 x 2 B2 x2 y = C 1 C 2 x + Du Next, the derivative of x 2 is set to zero and the resulting equation is solved for x 1 .
ndims Purpose 11ndims Provide the number of the dimensions of an LTI model or LTI array Syntax n = ndims(sys) Description n = ndims(sys) is the number of dimensions of an LTI model or an array of LTI models sys. A single LTI model has two dimensions (one for outputs, and one for inputs). An LTI array has 2+p dimensions, where p ≥ 2 is the number of array dimensions. For example, a 2-by-3-by-4 array of models has 2+3=5 dimensions.
ngrid Purpose 11ngrid Superimpose a Nichols chart on a Nichols plot Syntax ngrid Description ngrid superimposes Nichols chart grid lines over the Nichols frequency response of a SISO LTI system. The range of the Nichols grid lines is set to encompass the entire Nichols frequency response. The chart relates the complex number H ⁄ ( 1 + H ) to H , where H is any complex number.
ngrid Type nichols(H) ngrid Nichols Charts 40 0 dB 30 0.25 dB 0.
nichols Purpose Syntax 11nichols Compute Nichols frequency response of LTI models nichols(sys) nichols(sys,w) nichols(sys1,sys2,...,sysN) nichols(sys1,sys2,...,sysN,w) nichols(sys1,'PlotStyle1',...,sysN,'PlotStyleN') [mag,phase,w] = nichols(sys) [mag,phase] = nichols(sys,w) Description nichols computes the frequency response of an LTI model and plots it in the Nichols coordinates.
nichols When invoked with left-hand arguments, [mag,phase,w] = nichols(sys) [mag,phase] = nichols(sys,w) return the magnitude and phase (in degrees) of the frequency response at the frequencies w (in rad/sec). The outputs mag and phase are 3-D arrays similar to those produced by bode (see bode on page 11-19). They have dimensions ( number of outputs ) × ( number of inputs ) × ( length of w ) Remark If sys is an FRD model, nichols(sys,w), w can only include frequencies in sys.frequency.
nichols Nichols Charts 40 0 dB 30 0.25 dB 0.5 dB Open−Loop Gain (dB) 20 1 dB −1 dB 3 dB 10 −3 dB 6 dB −6 dB 0 −10 −12 dB −20 −20 dB −30 −40 −500 −450 −400 −350 −300 −250 −200 −150 −100 −50 −40 dB Open−Loop Phase (deg) Algorithm See bode.
norm Purpose Syntax 11norm Compute LTI model norms norm(sys) norm(sys,2) norm(sys,inf) norm(sys,inf,tol) [ninf,fpeak] = norm(sys) Description norm computes the H 2 or L ∞ norm of a continuous- or discrete-time LTI model.
norm The discrete-time counterpart is H(z ) Usage ∞ = jθ max σ max ( H ( e ) ) θ ∈ [ 0, π ] norm(sys) or norm(sys,2) both return the H 2 norm of the TF, SS, or ZPK model sys. This norm is infinite in the following cases: • sys is unstable. • sys is continuous and has a nonzero feedthrough (that is, nonzero gain at the frequency ω = ∞ ). Note that norm(sys) produces the same result as sqrt(trace(covar(sys,1))) norm(sys,inf) computes the infinity norm of any type of LTI model sys.
norm Compute its infinity norm by typing [ninf,fpeak] = norm(H,inf) ninf = 2.5488 fpeak = 3.0844 These values are confirmed by the Bode plot of H ( z ) .
norm MATLAB returns ans = 8.1268 Algorithm norm uses the same algorithm as covar for the H 2 norm, and the algorithm of [1] for the infinity norm. sys is first converted to state space. See Also bode freqresp sigma References Bode plot Frequency response computation Singular value plot [1] Bruisma, N.A. and M. Steinbuch, “A Fast Algorithm to Compute the H ∞ -Norm of a Transfer Function Matrix,” System Control Letters, 14 (1990), pp. 287–293.
nyquist Purpose Syntax 11nyquist Compute Nyquist frequency response of LTI models nyquist(sys) nyquist(sys,w) nyquist(sys1,sys2,...,sysN) nyquist(sys1,sys2,...,sysN,w) nyquist(sys1,'PlotStyle1',...,sysN,'PlotStyleN') [re,im,w] = nyquist(sys) [re,im] = nyquist(sys,w) Description nyquist calculates the Nyquist frequency response of LTI models. When invoked without left-hand arguments, nyquist produces a Nyquist plot on the screen.
nyquist When invoked with left-hand arguments [re,im,w] = nyquist(sys) [re,im] = nyquist(sys,w) return the real and imaginary parts of the frequency response at the frequencies w (in rad/sec). re and im are 3-D arrays with the frequency as last dimension (see “Arguments” below for details). Remark If sys is an FRD model, nyquist(sys,w), w can only include frequencies in sys.frequency.
nyquist 2 2s + 5s + 1 H ( s ) = -----------------------------2 s + 2s + 3 H = tf([2 5 1],[1 2 3]) nyquist(H) Nyquist Diagrams 2 1.5 1 Imaginary Axis 0.5 0 −0.5 −1 −1.5 −2 −1 −0.5 0 0.5 1 1.5 2 Real Axis See Also 11-158 bode evalfr freqresp ltiview nichols sigma Bode plot Response at single complex frequency Frequency response computation LTI system viewer Nichols plot Singular value plot 2.
obsv Purpose 11obsv Form the observability matrix Syntax Ob = obsv(A,B) Ob = obsv(sys) Description obsv computes the observability matrix for state-space systems. For an n-by-n matrix A and a p-by-n matrix C, obsv(A,C) returns the observability matrix C CA Ob = CA : CA 2 n–1 with n columns and np rows. Ob = obsv(sys) calculates the observability matrix of the state-space model sys. This syntax is equivalent to executing Ob = obsv(sys.A,sys.C) The model is observable if Ob has full rank n.
obsv MATLAB responds with unob = 0 See Also 11-160 obsvf Compute the observability staircase form
obsvf Purpose 11obsvf Compute the observability staircase form Syntax [Abar,Bbar,Cbar,T,k] = obsvf(A,B,C) [Abar,Bbar,Cbar,T,k] = obsvf(A,B,C,tol) Description If the observability matrix of (A,C) has rank r ≤ n , where n is the size of A, then there exists a similarity transformation such that T A = TAT , B = TB, C = CT T where T is unitary and the transformed system has a staircase form with the unobservable modes, if any, in the upper left corner.
obsvf Example Form the observability staircase form of A = 1 4 1 –2 1 1 –1 –1 1 0 0 1 B = C = by typing [Abar,Bbar,Cbar,T,k] = obsvf(A,B,C) Abar = 1 4 Bbar = 1 1 Cbar = 1 0 T = 1 0 k = 2 1 –2 1 –1 0 1 0 1 0 Algorithm obsvf is an M-file that implements the Staircase Algorithm of [1] by calling ctrbf and using duality. See Also ctrbf obsv References 11-162 Compute the controllability staircase form Calculate the observability matrix [1] Rosenbrock, M.M.
ord2 Purpose Syntax Description 11ord2 Generate continuous second-order systems [A,B,C,D] = ord2(wn,z) [num,den] = ord2(wn,z) [A,B,C,D] = ord2(wn,z) generates the state-space description (A,B,C,D) of the second-order system 1 h ( s ) = --------------------------------------------2 2 s + 2ζω n s + ω n given the natural frequency wn ( ω n ) and damping factor z ( ζ ). Use ss to turn this description into a state-space object.
pade Purpose Syntax 11pade Compute the Padé approximation of models with time delays [num,den] = pade(T,N) pade(T,N) sysx = pade(sys,N) sysx = pade(sys,NI,NO,Nio) Description pade approximates time delays by rational LTI models. Such approximations are useful to model time delay effects such as transport and computation delays within the context of continuous-time systems. The Laplace transform of an time delay of T seconds is exp ( – sT ) .
pade sysx = pade(sys,NI,NO,Nio) specifies independent approximation orders for each input, output, and I/O delay. These approximation orders are given by the arrays of integers NI, NO, and Nio, such that: • NI(j) is the approximation order for the j-th input channel. • NO(i) is the approximation order for the i-th output channel. • Nio(i,j) is the approximation order for the I/O delay from input j to output i.
pade Example Compute a third-order Padé approximation of a 0.1 second I/O delay and compare the time and frequency responses of the true delay and its approximation. To do this, type pade(0.1,3) Step response of 3rd−order Pade approximation 1.5 Amplitude 1 0.5 0 −0.5 −1 0 0.02 0.04 0.06 0.08 0.1 0.12 Time (secs) 0.14 0.16 0.18 0.2 Phase response 0 Phase (deg.
pade References [1] Golub, G. H. and C. F. Van Loan, Matrix Computations, Johns Hopkins University Press, Baltimore, 1989, pp. 557–558.
parallel Purpose 11parallel Parallel connection of two LTI models Syntax sys = parallel(sys1,sys2) sys = parallel(sys1,sys2,inp1,inp2,out1,out2) Description parallel connects two LTI models in parallel. This function accepts any type of LTI model. The two systems must be either both continuous or both discrete with identical sample time. Static gains are neutral and can be specified as regular matrices. sys = parallel(sys1,sys2) forms the basic parallel connection shown below.
parallel sys = parallel(sys1,sys2,inp1,inp2,out1,out2) forms the more general parallel connection. sys v1 y1 u1 + + u u2 v2 z1 sys1 sys2 y y2 z2 The index vectors inp1 and inp2 specify which inputs u 1 of sys1 and which inputs u 2 of sys2 are connected. Similarly, the index vectors out1 and out2 specify which outputs y 1 of sys1 and which outputs y 2 of sys2 are summed. The resulting model sys has [ v 1 ; u ; v 2 ] as inputs and [ z 1 ; y ; z 2 ] as outputs.
place Purpose 11place Pole placement design Syntax K = place(A,B,p) [K,prec,message] = place(A,B,p) Description Given the single- or multi-input system · x = Ax + Bu and a vector p of desired self-conjugate closed-loop pole locations, place computes a gain matrix K such that the state feedback u = – Kx places the closed-loop poles at the locations p. In other words, the eigenvalues of A – BK match the entries of p (up to the ordering).
place Algorithm place uses the algorithm of [1] which, for multi-input systems, optimizes the choice of eigenvectors for a robust solution. We recommend place rather than acker even for single-input systems. In high-order problems, some choices of pole locations result in very large gains. The sensitivity problems attached with large gains suggest caution in the use of pole placement techniques. See [2] for results from numerical testing.
pole Purpose 11pole Compute the poles of an LTI system Syntax p = pole(sys) Description pole computes the poles p of the SISO or MIMO LTI model sys. Algorithm For state-space models, the poles are the eigenvalues of the A matrix, or the generalized eigenvalues of A – λE in the descriptor case. For SISO transfer functions or zero-pole-gain models, the poles are simply the denominator roots (see roots).
pzmap Purpose 11pzmap Compute the pole-zero map of an LTI model Syntax pzmap(sys) [p,z] = pzmap(sys) Description pzmap(sys) plots the pole-zero map of the continuous- or discrete-time LTI model sys. For SISO systems, pzmap plots the transfer function poles and zeros. For MIMO systems, it plots the system poles and transmission zeros. The poles are plotted as x’s and the zeros are plotted as o’s.
pzmap Example Plot the poles and zeros of the continuous-time system. 2 2s + 5s + 1 H ( s ) = -----------------------------2 s + 2s + 3 H = tf([2 5 1],[1 2 3]) pzmap(H) Pole−zero map 1.5 1 Imag Axis 0.5 0 −0.5 −1 −1.5 −2.5 −2 −1.5 −1 −0.5 Real Axis Algorithm pzmap uses a combination of pole and zero.
pzmap sgrid, zgrid zero Plot lines of constant damping and natural frequency Compute system (transmission) zeros 11-175
reg Purpose 11reg Form regulator given state-feedback and estimator gains Syntax rsys = reg(sys,K,L) rsys = reg(sys,K,L,sensors,known,controls) Description rsys = reg(sys,K,L) forms a dynamic regulator or compensator rsys given a state-space model sys of the plant, a state-feedback gain matrix K, and an estimator gain matrix L. The gains K and L are typically designed using pole placement or LQG techniques. The function reg handles both continuous- and discrete-time cases.
reg This regulator should be connected to the plant using positive feedback. Plant y u u State Estimator –K x̂ y Regulator rsys = reg(sys,K,L,sensors,known,controls) handles more general regulation problems where: • The plant inputs consist of controls u , known inputs u d , and stochastic inputs w . • Only a subset y of the plant outputs is measured. The index vectors sensors, known, and controls specify y , u d , and u as subsets of the outputs and inputs of sys.
reg Example Given a continuous-time state-space model sys = ss(A,B,C,D) with seven outputs and four inputs, suppose you have designed: • A state-feedback controller gain K using inputs 1, 2, and 4 of the plant as control inputs • A state estimator with gain L using outputs 4, 7, and 1 of the plant as sensors, and input 3 of the plant as an additional known input You can then connect the controller and estimator and form the complete regulation system by controls = [1,2,4]; sensors = [4,7,1]; known = [3];
reshape Purpose 11reshape Change the shape of an LTI array Syntax sys = reshape(sys,s1,s2,...,sk) sys = reshape(sys,[s1 s2 ... sk]) Description sys = reshape(sys,s1,s2,...,sk) (or, equivalently, sys = reshape(sys,[s1 s2 ... sk])) reshapes the LTI array sys into an s1-by-s2-by...-sk array of LTI models. Equivalently, sys = reshape(sys,[s1 s2 ... sk]) reshapes the LTI array sys into an s1-by-s2-by...-sk array of LTI models. With either syntax, there must be s1*s2*...*sk models in sys to begin with.
rlocfind Purpose 11rlocfind Select feedback gain from root locus plot Syntax [k,poles] = rlocfind(sys) [k,poles] = rlocfind(sys,p) Description rlocfind returns the feedback gain associated with a particular set of poles on the root locus. rlocfind works with both continuous- and discrete-time SISO systems. [k,poles] = rlocfind(sys) is used for interactive gain selection from the root locus plot of the SISO system sys generated by rlocus.
rlocfind where p is the complex point you supply as an input argument to rlocfind (or where you point to with the mouse), and n and d are respectively the numerator and denominator polynomials of the transfer function associated with sys. The poles associated with the gain k are then computed as poles = roots(k*n(s) + d(s)) Limitations rlocfind assumes that a root locus is in the current figure window when this function is called without second input argument.
rlocus Purpose Syntax 11rlocus Evans root locus rlocus(sys) rlocus(sys,k) [r,k] = rlocus(sys) r = rlocus(sys,k) Description rlocus computes the Evans root locus of a SISO open-loop model. The root locus gives the closed-loop pole trajectories as a function of the feedback gain k (assuming negative feedback). Root loci are used to study the effects of varying feedback gains on closed-loop pole locations. In turn, these locations provide indirect information on the time and frequency responses.
rlocus the closed-loop poles are the roots of d( s ) + k n( s ) = 0 rlocus adaptively selects a set of positive gains k to produce a smooth plot. Alternatively, rlocus(sys,k) uses the user-specified vector k of gains to plot the root locus. When invoked with output arguments, [r,k] = rlocus(sys) r = rlocus(sys,k) return the vector k of selected gains and the complex root locations r for these gains. The matrix r has length(k) columns and its jth column lists the closed-loop roots for the gain k(j).
rlocus 2 2s + 5s + 1 h ( s ) = -----------------------------2 s + 2s + 3 h = tf([2 5 1],[1 2 3]); rlocus(h) 1.5 1 Imag Axis 0.5 0 −0.5 −1 −1.5 −3 −2.5 −2 −1.5 −1 Real Axis −0.5 0 0.5 1 For examples, see “Root Locus Design” on page 9-9 and “Hard-Disk Read/Write Head Controller” on page 9-20.
rltool Purpose 11rltool Initialize the Root Locus Design GUI Syntax rltool rltool(sys) rltool(sys,comp) rltool(sys,comp,LocationFlag,FeedbackSign) Description When invoked without input arguments, rltool initializes a new Root Locus Design GUI for interactive compensator design. This GUI allows you to design a single-input/single-output (SISO) compensator using root locus techniques.
rltool The Root Locus Design GUI looks like this. These are the main menus for importing/exporting of models, and editing them. You can also perform discrete/ continuous conversion. Compensator description: The default compensator is K=1.
rltool This tool can be applied to SISO LTI systems whose feedback structure is in one of the following two configurations. In either configuration, F is a pre-filter, P is the plant model, H is the sensor dynamics, and K is the compensator to be designed. In terms of the GUI design procedure, once you specify them, F, P, and H are fixed in the feedback structure. This triple, along with the feedback structure, is called the design model.
rltool model (sys) in the forward loop of a negative unity feedback system, as shown in the diagram below. u comp sys y +_ In this case, F and H are taken to be 1, while P is sys. If you want to include F and H in the design model after loading rltool(sys) or rltool(sys,comp), select the Import Model menu item from the File menu of the Root Locus Design GUI to load F and H. rltool(sys,comp,LocationFlag,FeedbackSign) allows you to override the default compensator location and feedback sign.
rmodel, rss Purpose Syntax 11rmodel, rss Generate stable random continuous test models sys = rss(n) sys = rss(n,p) sys = rss(n,p,m) sys = rss(n,p,m,s1,...,sn) [num,den] = rmodel(n) [A,B,C,D] = rmodel(n) [A,B,C,D] = rmodel(n,p,m) Description rss(n) produces a stable random n-th order model with one input and one output and returns the model in the state-space object sys.
rmodel, rss Example Obtain a stable random continuous LTI model with three states, two inputs, and two outputs by typing sys = rss(3,2,2) a = x1 x2 x3 x1 –0.54175 0.09729 0.08304 x2 0.09729 –0.89491 0.58707 x1 x2 x3 u1 –0.88844 0 –0.07162 u2 –2.41459 –0.69435 –1.39139 y1 y2 x1 0.32965 0.59854 x2 0.14718 –0.10144 y1 y2 u1 –0.87631 0 u2 –0.32758 0 x3 0.08304 0.58707 –1.95271 b = c = x3 0 0.02805 d = Continuous-time system.
series Purpose 11series Series connection of two LTI models Syntax sys = series(sys1,sys2) sys = series(sys1,sys2,outputs1,inputs2) Description series connects two LTI models in series. This function accepts any type of LTI model. The two systems must be either both continuous or both discrete with identical sample time. Static gains are neutral and can be specified as regular matrices. sys = series(sys1,sys2) forms the basic series connection shown below.
series The index vectors outputs1 and inputs2 indicate which outputs y 1 of sys1 and which inputs u 2 of sys2 should be connected. The resulting model sys has u as input and y as output. Example Consider a state-space system sys1 with five inputs and four outputs and another system sys2 with two inputs and three outputs. Connect the two systems in series by connecting outputs 2 and 4 of sys1 with inputs 1 and 2 of sys2.
set Purpose Syntax 11set Set or modify LTI model properties set(sys,'Property',Value) set(sys,'Property1',Value1,'Property2',Value2,...) set(sys,'Property') set(sys) Description set is used to set or modify the properties of an LTI model (see “LTI Properties” on page 2-26 for background on LTI properties). Like its Handle Graphics counterpart, set uses property name/property value pairs to update property values.
set Note that set does not require any output argument. Check the result with get by typing get(sys) a: 1 b: 2 c: 3 d: 0 e: [] StateName: {''} Ts: 0 InputDelay: 0.1 OutputDelay: 0 ioDelayMatrix: 0 InputName: {'torque'} OutputName: {''} InputGroup: {0x2 cell} OutputGroup: {0x2 cell} Notes: {} UserData: -2 Property Values 11-194 The following table lists the admissible values for each LTI property. N u and N y denotes the number of inputs and outputs of the underlying LTI model.
set Table 11-15: LTI Properties Property Name Admissible Property Values Ts • 0 (zero) for continuous-time systems • Sample time in seconds for discrete-time systems • –1 or [] for discrete systems with unspecified sample time Note: Resetting the sample time property does not alter the model data. Use c2d, d2c, or d2d for discrete/continuous and discrete/discrete conversions.
set Table 11-15: LTI Properties (Continued) Property Name Admissible Property Values OutputDelay Output delays specified with • Nonnegative real numbers for continuous-time models (seconds) • Integers for discrete-time models (number of sample periods) • Scalar when N y = 1 or system has uniform output delay • Vector of length N y to specify independent delay times for each output channel • Array of size N y -by- N u -by- S 1 -by-. . .
set Table 11-17: TF Model Properties Property Name Admissible Property Values num, den • Real-valued row vectors for the coefficients of the numerator or denominator polynomials in the SISO case. List the coefficients in descending powers of the variable s or z by default, and in ascending –1 powers of q = z when the Variable property is set to 'q' or 'z^–1' (see note below).
set Table 11-19: FRD Model Properties Property Name Admissible Property Values Frequency Real-valued vector of length N f -by-1, where N f is the number of frequencies Response • N y -by- N u -by- N f -dimensional array of complex data for single LTI models • N y -by- N u -by- N f -by- S 1 -by- …-by- S K -dimensional array for LTI arrays Units Remark String 'rad/s' (default), or 'Hz' For discrete-time transfer functions, the convention used to represent the numerator and denominator depends on the
set –1 1 + 2z –1 h ( z ) = --------------------------------------- = zh ( z ) –1 –2 1 + 3z + 4z Because the resulting transfer functions are different, make sure to use the convention consistent with your choice of variable.
sgrid Purpose 11sgrid Generate an s-plane grid of constant damping factors and natural frequencies Syntax sgrid sgrid(z,wn) Description sgrid generates a grid of constant damping factors from zero to one in steps of 0.1 and natural frequencies from zero to 10 rad/sec in steps of one rad/sec, and plots the grid over the current axis. If the current axis contains a continuous s-plane root locus diagram or pole-zero map, sgrid draws the grid over the plot.
sgrid 1.5 1 Imag Axis 0.5 0 −0.5 −1 −1.5 −3 Limitations −2.5 −2 −1.5 −1 Real Axis −0.5 0 0.5 1 sgrid plots the grid over the current axis regardless of whether the axis contains a root locus diagram or pole-zero map. Therefore, if the current axes contains, for example, a step response, you may superimpose a meaningless s-plane grid over the plot. In addition, sgrid does not rescale the axis limits of the current axis.
sigma Purpose Syntax 11sigma Singular values of the frequency response of LTI models sigma(sys) sigma(sys,w) sigma(sys,w,type) sigma(sys1,sys2,...,sysN) sigma(sys1,sys2,...,sysN,w) sigma(sys1,sys2,...,sysN,w,type) sigma(sys1,'PlotStyle1',...,sysN,'PlotStyleN') [sv,w] = sigma(sys) sv = sigma(sys,w) Description sigma calculates the singular values of the frequency response of an LTI model. For an FRD model, sys, sigma computes the singular values of sys.Response at the frequencies, sys.frequency.
sigma w = {wmin,wmax}. To use particular frequency points, set w to the corresponding vector of frequencies. Use logspace to generate logarithmically spaced frequency vectors. The frequencies must be specified in rad/sec. sigma(sys,[],type) or sigma(sys,w,type) plots the following modified singular value responses: –1 type = 1 Singular values of the frequency response H the frequency response of sys. type = 2 Singular values of the frequency response I + H .
sigma Remark If sys is an FRD model, sigma(sys,w), w can only include frequencies in sys.frequency. Example Plot the singular value responses of 0 H(s) = s+1 -----------s+5 and I + H ( s ) .
sigma You can do this by typing H = [0 tf([3 0],[1 1 10]) ; tf([1 1],[1 5]) tf(2,[1 6])] subplot(211) sigma(H) subplot(212) sigma(H,[],2) Singular Values Singular Values (dB) 20 0 −20 −40 −60 −2 10 −1 10 0 10 Frequency (rad/sec) 1 10 2 10 Singular Values Singular Values (dB) 20 10 0 −10 −20 −30 −40 −2 10 Algorithm −1 10 0 10 Frequency (rad/sec) 1 10 2 10 sigma uses the svd function in MATLAB to compute the singular values of a complex matrix.
sigma nichols nyquist 11-206 Nichols plot Nyquist plot
size Purpose 11size Provide the output/input/array dimensions of LTI models, the model order of TF, SS, and ZPK models, and the number of frequencies of FRD models Syntax size(sys) d = size(sys) Ny = size(sys,1) Nu = size(sys,2) Sk = size(sys,2+k) Ns = size(sys,'order') Nf = size(sys,'frequency') Description When invoked without output arguments, size(sys) returns a vector of the number of outputs and inputs for a single LTI model.
size Example Consider the random LTI array of state-space models sys = rss(5,3,2,3); Its dimensions are obtained by typing size(sys) 3x1 array of state-space models Each model has 3 outputs, 2 inputs, and 5 states.
sminreal Purpose 11sminreal Perform model reduction based on structure Syntax msys = sminreal(sys) Description msys = sminreal(sys) eliminates the states of the state-space model sys that don’t affect the input/output response. All of the states of the resulting state-space model msys are also states of sys and the input/output response of msys is equivalent to that of sys. sminreal eliminates only structurally non minimal states, i.e.
sminreal all of the states of sys, including those of sys2 are retained.
ss Purpose Syntax 11ss Specify state-space models or convert an LTI model to state space sys sys sys sys = = = = ss(a,b,c,d) ss(a,b,c,d,Ts) ss(d) ss(a,b,c,d,ltisys) sys = ss(a,b,c,d,'Property1',Value1,...,'PropertyN',ValueN) sys = ss(a,b,c,d,Ts,'Property1',Value1,...,'PropertyN',ValueN) sys_ss = ss(sys) sys_ss = ss(sys,'minimal') Description ss is used to create real-valued state-space models (SS objects) or to convert transfer function or zero-pole-gain models to state space.
ss with sample time Ts (in seconds). Set Ts = –1 or Ts = [] to leave the sample time unspecified. sys = ss(d) specifies a static gain matrix D and is equivalent to sys = ss([],[],[],d) sys = ss(a,b,c,d,ltisys) creates a state-space model with generic LTI properties inherited from the LTI model ltisys (including the sample time). See “Generic Properties” on page 2-26 for an overview of generic LTI properties.
ss Examples Example 1 The command sys = ss(A,B,C,D,0.05,'statename',{'position' 'velocity'},... 'inputname','force',... 'notes','Created 10/15/96') creates a discrete-time model with matrices A, B, C, D and sample time 0.05 second. This model has two states labeled position and velocity, and one input labeled force (the dimensions of A, B, C, D should be consistent with these numbers of states and inputs). Finally, a note is attached with the date of creation of the model.
ss The resulting state-space model order has order three, the minimum number of states needed to represent H(s). This can be seen directly by factoring H(s) as the product of a first order system with a second order one. 1 ------------ 0 H(s) = s + 2 0 1 See Also 11-214 dss frd get set ssdata tf zpk s+1 ----------------------2 s +s+1 2 s +3 ----------------------2 s +s+1 Specify descriptor state-space models. Specify FRD models or convert to an FRD. Get properties of LTI models.
ss2ss Purpose 11ss2ss State coordinate transformation for state-space models Syntax sysT = ss2ss(sys,T) Description Given a state-space model sys with equations · x = Ax + Bu y = Cx + Du (or their discrete-time counterpart), ss2ss performs the similarity transformation x = Tx on the state vector x and produces the equivalent state-space model sysT with equations.
ssbal Purpose 11ssbal Balance state-space models using a diagonal similarity transformation Syntax [sysb,T] = ssbal(sys) [sysb,T] = ssbal(sys,condT) Description Given a state-space model sys with matrices ( A, B, C, D ) , [sysb,T] = ssbal(sys) computes a diagonal similarity transformation T and a scalar α such that TAT αCT –1 –1 TB ⁄ α 0 has approximately equal row and column norms.
ssbal Balance this model with ssbal by typing ssbal(sys) a = x1 x2 x3 x1 1 0 2560 x1 x2 x3 u1 0.125 0.5 32 y1 x1 0.8 y1 u1 0 x2 2500 100 64 x3 0.39063 1562.5 0 x2 20 x3 3.125 b = c = d = Continuous-time system. Direct inspection shows that the range of numerical values has been compressed by a factor 100 and that the B and C matrices now have nearly equal norms. Algorithm ssbal uses the MATLAB function balance to compute T and α .
ssdata Purpose 11ssdata Quick access to state-space model data Syntax [a,b,c,d] = ssdata(sys) [a,b,c,d,Ts] = ssdata(sys) Description [a,b,c,d] = ssdata(sys) extracts the matrix (or multidimensional array) data ( A, B, C, D ) from the state-space model (LTI array) sys. If sys is a transfer function or zero-pole-gain model (LTI array), it is first converted to state space. See Table 11-16, “State-Space Model Properties,” on page 11-196 for more information on the format of state-space model data.
stack Purpose 11stack Build an LTI array by stacking LTI models or LTI arrays along array dimensions of an LTI array Syntax sys = stack(arraydim,sys1,sys2,...) Description sys = stack(arraydim,sys1,sys2,...) produces an array of LTI models sys by stacking (concatenating) the LTI models (or LTI arrays) sys1,sys2,... along the array dimension arraydim. All models must have the same number of inputs and outputs (the same I/O dimensions). The I/O dimensions are not counted in the array dimensions.
step Purpose Syntax 11step Step response of LTI systems step(sys) step(sys,t) step(sys1,sys2,...,sysN) step(sys1,sys2,...,sysN,t) step(sys1,'PlotStyle1',...,sysN,'PlotStyleN') [y,t,x] = step(sys) Description step calculates the unit step response of a linear system. Zero initial state is assumed in the state-space case. When invoked with no output arguments, this function plots the step response on the screen. step(sys) plots the step response of an arbitrary LTI model sys.
step All systems must have the same number of inputs and outputs but may otherwise be a mix of continuous- and discrete-time systems. This syntax is useful to compare the step responses of multiple systems. You can also specify a distinctive color, linestyle, and/or marker for each system. For example, step(sys1,'y:',sys2,'g--') plots the step response of sys1 with a dotted yellow line and the step response of sys2 with a green dashed line.
step · x1 – 0.5572 – 0.7814 x 1 1 = + · x2 0.7814 0 0 x2 y = 1.9691 6.4493 –1 u1 2 u2 x1 x2 a = [–0.5572 –0.7814;0.7814 b = [1 –1;0 2]; c = [1.9691 6.4493]; sys = ss(a,b,c,0); step(sys) 0]; Step Response From: U1 From: U2 12 10 8 Amplitude 6 4 2 0 −2 −4 0 5 10 15 200 5 10 15 20 Time (sec.) The left plot shows the step response of the first input channel, and the right plot shows the step response of the second input channel.
step Algorithm Continuous-time models are converted to state space and discretized using zero-order hold on the inputs. The sampling period is chosen automatically based on the system dynamics, except when a time vector t = 0:dt:Tf is supplied (dt is then used as sampling period).
tf Purpose Syntax 11tf Specify transfer functions or convert LTI model to transfer function form sys sys sys sys = = = = tf(num,den) tf(num,den,Ts) tf(M) tf(num,den,ltisys) sys = tf(num,den,'Property1',Value1,...,'PropertyN',ValueN) sys = tf(num,den,Ts,'Property1',Value1,...
tf If all SISO entries of a MIMO transfer function have the same denominator, you can set den to the row vector representation of this common denominator. See “Examples” for more details. sys = tf(num,den,Ts) creates a discrete-time transfer function with sample time Ts (in seconds). Set Ts = –1 or Ts = [] to leave the sample time unspecified. The input arguments num and den are as in the continuous-time case and must list the numerator and denominator coefficients in descending powers of z .
tf • s = tf('s') to specify a TF model using a rational function in the Laplace variable, s. • z = tf('z',Ts) to specify a TF model with sample time Ts using a rational function in the discrete-time variable, z. Once you specify either of these variables, you can specify TF models directly as rational expressions in the variable s or z by entering your transfer function as a rational expression in either s or z.
tf To do this, type num = {[1 1] ; 1} den = {[1 2 2] ; [1 0]} H = tf(num,den,'inputn','current',... 'outputn',{'torque' 'ang. velocity'},... 'variable','p') Transfer function from input "current" to output... p + 1 torque: ------------p^2 + 2 p + 2 ang. velocity: 1 p Note how setting the 'variable' property to 'p' causes the result to be displayed as a transfer function of the variable p .
tf with common denominator d ( z ) = z + 0.3 and sample time of 0.2 seconds. nums = {1 [1 0];[–1 2] 3} Ts = 0.2 H = tf(nums,[1 0.3],Ts) % Note: row vector for common den. d(z) Example 4 Compute the transfer function of the state-space model with the following data.
tf engineers use the z variable and order the numerator and denominator terms in descending powers of z , for example, 2 z h ( z ) = --------------------------z 2 + 2z + 3 2 The polynomials z and z 2 + 2z + 3 are then specified by the row vectors [1 0 0] and [1 2 3], respectively. By contrast, DSP engineers prefer to write this transfer function as 1 –1 h ( z ) = --------------------------------------–1 –2 1 + 2z + 3z and specify its numerator as 1 (instead of [1 0 0]) and its denominator as [1 2 3].
tf uses the DSP convention and creates –1 1+z –1 h ( z ) = --------------------------------------- = zg ( z ) –1 –2 1 + 2z + 3z See also filt for direct specification of discrete transfer functions using the DSP convention. Note that tf stores data so that the numerator and denominator lengths are made equal.
tfdata Purpose 11tfdata Quick access to transfer function data Syntax [num,den] = tfdata(sys) [num,den] = tfdata(sys,'v') [num,den,Ts] = tfdata(sys) Description [num,den] = tfdata(sys) returns the numerator(s) and denominator(s) of the transfer function for the TF, SS or ZPK model (or LTI array of TF, SS or ZPK models) sys.
tfdata you can extract the numerator and denominator coefficients by typing [num,den] = tfdata(h,'v') num = 0 1 1 den = 1 2 5 This syntax returns two row vectors. If you turn h into a MIMO transfer function by typing H = [h ; tf(1,[1 1])] the command [num,den] = tfdata(H) now returns two cell arrays with the numerator/denominator data for each SISO entry. Use celldisp to visualize this data. Type celldisp(num) and MATLAB returns the numerator vectors of the entries of H.
tfdata tf zpkdata Specify transfer functions Quick access to zero-pole-gain data 11-233
totaldelay Purpose 11totaldelay Return the total combined I/O delays for an LTI model Syntax td = totaldelay(sys) Description td = totaldelay(sys) returns the total combined I/O delays for an LTI model sys. The matrix td combines contributions from the InputDelay, OutputDelay, and ioDelayMatrix properties, (see set on page 11-193 or type ltiprops for details on these properties).
zero Purpose 11zero Transmission zeros of LTI models Syntax z = zero(sys) [z,gain] = zero(sys) Description zero computes the zeros of SISO systems and the transmission zeros of MIMO systems. For a MIMO system with matrices ( A, B, C, D ) , the transmission zeros are the complex values λ for which the normal rank of A – λI B C D drops. z = zero(sys) returns the (transmission) zeros of the LTI model sys as a column vector.
zgrid Purpose 11zgrid Generate a z-plane grid of constant damping factors and natural frequencies Syntax zgrid zgrid(z,wn) Description zgrid generates a grid of constant damping factors from zero to one in steps of 0.1 and natural frequencies from zero to π in steps of π ⁄ 10 , and plots the grid over the current axis. If the current axis contains a discrete z-plane root locus diagram or pole-zero map, zgrid draws the grid over the plot without altering the current axis limits.
zgrid To see the z-plane grid on the root locus plot, type rlocus(H) zgrid axis('square') 1 0.8 0.6 0.4 Imag Axis 0.2 0 −0.2 −0.4 −0.6 −0.8 −1 −1 Limitations −0.8 −0.6 −0.4 −0.2 0 0.2 Real Axis 0.4 0.6 0.8 1 zgrid plots the grid over the current axis regardless of whether the axis contains a root locus diagram or pole-zero map. Therefore, if the current axis contains, for example, a step response, you may superimpose a meaningless z-plane grid over the plot.
zpk Purpose Syntax 11zpk Specify zero-pole-gain models or convert LTI model to zero-pole-gain form sys sys sys sys = = = = zpk(z,p,k) zpk(z,p,k,Ts) zpk(M) zpk(z,p,k,ltisys) sys = zpk(z,p,k,'Property1',Value1,...,'PropertyN',ValueN) sys = zpk(z,p,k,Ts,'Property1',Value1,...
zpk • s = zpk('s') to specify a ZPK model from a rational transfer function of the Laplace variable, s. • z = zpk('z',Ts) to specify a ZPK model with sample time Ts from a rational transfer function of the discrete-time variable, z. Once you specify either of these variables, you can specify ZPK models directly as real-valued rational expressions in the variable s or z. To create a MIMO zero-pole-gain model, specify the zeros, poles, and gain of each SISO entry of this model.
zpk Each pair specifies a particular LTI property of the model, for example, the input names or the input delay time. See set entry and the example below for details. Note that sys = zpk(z,p,k,'Property1',Value1,...,'PropertyN',ValueN) is a shortcut for the following sequence of commands. sys = zpk(z,p,k) set(sys,'Property1',Value1,...,'PropertyN',ValueN) Zero-Pole-Gain Models as Rational Expressions in s or z You can also use rational expressions to create a ZPK model.
zpk Example Example 1 Specify the following zero-pole-gain model. 1 ----------------z – 0.3 H(z ) = 2 ( z + 0.5 ) -----------------------------------------------------------( z – 0.1 + j ) ( z – 0.1 – j ) To do this, type z p k H = = = = {[] ; –0.5} {0.3 ; [0.1+i 0.
zpk Example 3 Create a discrete-time ZPK model from a rational expression in the variable z, by typing z = zpk('z',0.1); H = (z+.1)*(z+.2)/(z^2+.6*z+.09) Zero/pole/gain: (z+0.1) (z+0.2) --------------(z+0.3)^2 Sampling time: 0.1 Algorithm zpk uses the MATLAB function roots to convert transfer functions and the functions zero and pole to convert state-space models.
zpkdata Purpose 11zpkdata Quick access to zero-pole-gain data Syntax [z,p,k] = zpkdata(sys) [z,p,k] = zpkdata(sys,'v') [z,p,k,Ts,Td] = zpkdata(sys) Description [z,p,k] = zpkdata(sys) returns the zeros z, poles p, and gain(s) k of the zeropole-gain model sys. The outputs z and p are cell arrays with the following characteristics: • z and p have as many rows as outputs and as many columns as inputs.
zpkdata Example Given a zero-pole-gain model with two outputs and one input H = zpk({[0];[–0.5]},{[0.3];[0.1+i 0.1–i]},[1;2],–1) Zero/pole/gain from input to output... 1 #1: ------(z–0.3) #2: 2 (z+0.5) ------------------(z^2 – 0.2z + 1.01) Sampling time: unspecified you can extract the zero/pole/gain data embedded in H with [z,p,k] = zpkdata(H) z = [ 0] [-0.5000] p = [ 0.
zpkdata See Also get ssdata tfdata zpk Get properties of LTI models Quick access to state-space data Quick access to transfer function data Specify zero-pole-gain models 11-245
zpkdata 11-246
Index A acker 11-11 addition of LTI models 3-11 scalar 3-12 adjoint. See pertransposition algebraic loop 11-76 aliasing 5-13 analysis models specifying 6-56 See also Simulink LTI Viewer append 3-16, 3-17, 4-29, 11-12 array dimensions 4-7 arrays. See LTI arrays augstate 11-15 axes grouping menu 6-23 B balancing realizations 5-7, 11-16, 11-216 balreal 11-16 block diagram.
Index SS model, to 2-42 state-space, to 2-44, 11-212 TF model, to 2-42 ZPK model, to 2-42 covar 11-40 covariance error 9-56, 9-62, 11-110 noise 7-9, 11-110 output 11-40, 11-110 state 11-40, 11-110 ctrb 11-43 ctrbf 11-45 D d2c 11-48 d2d 3-27, 11-51 damp 11-52 damping 11-52, 11-200, 11-236 dare 11-54 dcgain 11-57, 11-58 delay2z 11-58 delays arithmetic operations 3-15 c2d and d2c conversions 3-25 combining 2-54, 11-234 conversion 11-58 conversion to SS 2-54 delay2z 11-58 discrete-time models 2-52 discretiza
Index control design 9-20 equivalent continuous poles 11-52 frequency 5-13, 11-23 Kalman estimator 9-50, 11-108 random 11-62 resampling 3-27 See also LTI models discretization 2-36, 3-20, 9-21, 11-24 available methods 11-24 delay systems 3-24 first-order hold 3-22 intersample behavior 11-129 matched poles/zeros 3-23 Tustin method 3-22 zero-order hold 3-20 dlqr 11-59 dlyap 11-61 drmodel 11-62 drss 11-62 dsort 11-64 DSP convention 11-77, 11-198, 11-228 dss 11-65 dssdata 11-67 dual.
Index Nyquist 5-13, 11-23 range 5-12 frequency response 2-18, 5-11 at single frequency (evalfr) 11-72 Bode plot 11-19 customized plots 5-17 discrete-time frequency 5-13, 11-23 freqresp 11-84 magnitude 11-19 MIMO 5-12, 11-19, 11-149, 11-156 Nichols chart (ngrid) 11-147 Nichols plot 11-149 Nyquist plot 11-156 phase 11-19 plotting 5-13, 11-19, 11-134 singular value plot 11-202 viewing the gain and phase margins 11-137 G gain 2-11 estimator gain 7-6, 7-9 feedthrough 2-28 low frequency (DC) 11-57 property LTI
Index square wave 11-87 input point block 6-56 See also Simulink LTI Viewer InputDelay. See delays InputGroup 2-26, 2-27 conflicts, naming 3-4 See also I/O groups InputName 2-34, 2-36 conflicts, naming 3-4 See also I/O names interconnection. See model building intersample behavior 11-129 inv 11-101 inversion 11-101 limitations 11-102 model 3-13 ioDelayMatrix. See delay isct 11-104 isdt 11-104 isempty 11-105 isproper 11-106 issiso 11-107 K Kalman filter.
Index extracting subsystems 4-21 indexing into 4-20 interconnection functions 4-25 LTI Viewer, model selector 6-28, 6-31 model dimensions 4-7 operations on 4-25 dimension requirements 4-27 special cases 4-27 reassigning parts of 4-22 SS models 4-23 shape, changing 11-179 size 4-7 stack 4-15, 11-219 LTI models addition 3-11 scalar 3-12 building 3-16 characteristics 5-2 comparing multiple models 5-13, 11-19, 11-94, 11-127, 11-134 concatenation effects on model order 3-10 horizontal 3-17 vertical 3-17 continu
Index model-specific properties 2-28 online help (ltiprops) 2-26 property names 2-26, 2-30, 11-89, 11-193 property values 2-26, 2-31, 11-89, 11-193 setting 2-30, 11-193, 11-212, 11-225, 11-240 sample time 3-3 variable property 3-4 LTI Viewer 11-133 axes grouping menu 6-23 command line initializing 6-5 file menu 6-15 frequency domain plot units 6-44 getting help 6-16 importing models 6-11, 6-15 multiselection and deselection 6-11 line properties, order 6-47 linestyle properties 6-46 LTI array selector 6-28
Index multiplication 3-13 scalar 3-13 multiselection of items in a window 6-11 N natural frequency 11-52 ndims 11-146 ngrid 11-147 Nichols chart 11-147 plot (nichols) 11-149 nichols 11-149 noise covariance 7-9, 11-110 measurement 7-8, 11-70 process 7-8, 11-70 white 5-9, 7-8, 11-40 norm 11-152 norms of LTI systems (norm) 11-152 Notes 2-27 numerator property 2-28, 11-197 specification 2-8, 2-10, 2-11, 2-22, 11-77 value 2-25, 11-90 numerical stability 10-6 Nyquist frequency 5-13, 11-23 plot (nyquist) 11-156
Index See also I/O, groups OutputName 2-34 conflicts, naming 3-4 See also I/O, names overshoot 5-9 P pade 11-164 Padé approximation (pade) 2-55, 11-164 parallel 11-168 parallel connection 3-12, 9-54, 11-168 pertransposition 3-14 phase margins 9-29, 11-19, 11-137 place 11-170 plot configuration, LTI Viewer 6-39 plotting customized plots 5-17 frequency response.
Index See also LTI Viewer response preferences, setting 6-40 response, I/O 3-5 Riccati equation 7-9 continuous (care) 11-29 discrete (dare) 11-54 for LQG design 11-110, 11-121 Η∞-like 11-31 stabilizing solution 11-29, 11-54 right-click menus LTI arrays 6-28 right-click menus, LTI Viewer 6-7 See also LTI Viewer, right-click menus rise time 5-9 rlocfind 11-180 rlocus 11-182 rltool 11-179, 11-185 rmodel 11-189 robustness 9-29 root locus design 7-3, 9-9, 9-24 plot (rlocus) 11-182 select gain from 11-180 See al
Index analysis models 6-50 clearing 6-53 open and closed loop 6-56 saving 6-65 specifying 6-53, 6-56 input point blocks 6-53 linearizing models 6-53, 6-63 opening 6-50 operating conditions, changing 6-61 operating conditions, setting 6-58 operating points 6-53 output point blocks 6-53 Simulink menu 6-53 specifying models for 6-51 sine wave 11-87 singular value plot (bode) 11-202 SISO 2-2, 5-2, 8-2, 11-107 size 11-207 sminreal 11-209 square wave 11-87 ss 2-15, 11-211 SS objects.
Index tfdata output, form of 2-24 time delays.
Index pole-zero map 11-173 property 2-28 transmission 11-235 zooming LTI Viewer 6-12 Root Locus Design GUI 8-15 zpk 2-12, 11-238 ZPK objects.
